diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 76 |
1 files changed, 54 insertions, 22 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 90c20b01d..77aa7342f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -66,9 +66,38 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) return k_comp; } +SpectralShiftFactor +SpectralKSpace::getSpectralShiftFactor( + const DistributionMapping& dm, const int i_dim, const int shift_type ) const +{ + // Initialize an empty ManagedVector in each box + SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); + // Loop over boxes + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + const ManagedVector<Real>& k = k_vec[i_dim][mfi]; + ManagedVector<Complex>& shift = shift_factor[mfi]; + + // Allocate shift coefficients + shift.resize( k.size() ); + + // Fill the shift coefficients + Real sign = 0; + switch (shift_type){ + case ShiftType::CenteredToNodal: sign = -1.; break; + case ShiftType::NodalToCentered: sign = 1.; + } + constexpr Complex I{0,1}; + for (int i=0; i<k.size(); i++ ){ + shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] ); + } + } + return shift_factor; +} + KVectorComponent SpectralKSpace::getModifiedKComponent( - const DistributionMapping& dm, const int i_dim, const int order ) const + const DistributionMapping& dm, const int i_dim, + const int n_order, const bool nodal ) const { // Initialize an empty ManagedVector in each box KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); @@ -90,30 +119,33 @@ SpectralKSpace::getModifiedKComponent( return modified_k_comp; } -SpectralShiftFactor -SpectralKSpace::getSpectralShiftFactor( - const DistributionMapping& dm, const int i_dim, const int shift_type ) const +/* TODO: Documentation: point to Fonberg paper ; explain recurrence relation + */ +Array<Real> +getFonbergStencilCoefficients( const int n_order, const bool nodal ) { - // Initialize an empty ManagedVector in each box - SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); - // Loop over boxes - for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ - const ManagedVector<Real>& k = k_vec[i_dim][mfi]; - ManagedVector<Complex>& shift = shift_factor[mfi]; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, "n_order should be even."); + const int m = n_order/2; + Array<Real> stencil_coef; + stencil_coef.resize( m+1 ); - // Allocate shift coefficients - shift.resize( k.size() ); - - // Fill the shift coefficients - Real sign = 0; - switch (shift_type){ - case ShiftType::CenteredToNodal: sign = -1.; break; - case ShiftType::NodalToCentered: sign = 1.; + // Coefficients for nodal (a.k.a. centered) finite-difference + if (nodal == true) { + stencil_coef[0] = -2.; // First coefficient + for (int n=1; n<m+1; n++){ // Get the other coefficients by recurrence + stencil_coef[n] = - (m+1-n)*1./(m+n)*stencil_coef[n-1]; } - constexpr Complex I{0,1}; - for (int i=0; i<k.size(); i++ ){ - shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] ); + } + // Coefficients for staggered finite-difference + else { + Real prod = 1.; + for (int k=1; k<m+1; k++){ + prod *= (m+k)*1./(4*k); + } + stencil_coef[0] = 4*m*prod*prod; // First coefficient + for (int n=1; n<m+1; n++){ // Get the other coefficients by recurrence + stencil_coef[n] = - ((2*n-3)*(m+1-n))*1./((2*n-1)*(m-1+n))*stencil_coef[n-1]; } } - return shift_factor; + return stencil_coef; } |