From 0e9c0209744cf032da4713a1041dd0ccbaa0b1c4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 06:04:40 -0700 Subject: Define SpectralSolver class --- .../SpectralSolver/SpectralFieldData.cpp | 186 +++++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100644 Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp new file mode 100644 index 000000000..25c8f438b --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -0,0 +1,186 @@ +#include + +using namespace amrex; + +SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const DistributionMapping& dm ) +{ + const BoxArray& spectralspace_ba = k_space.spectralspace_ba; + + // Allocate the arrays that contain the fields in spectral space + Ex = SpectralField(spectralspace_ba, dm, 1, 0); + Ey = SpectralField(spectralspace_ba, dm, 1, 0); + Ez = SpectralField(spectralspace_ba, dm, 1, 0); + Bx = SpectralField(spectralspace_ba, dm, 1, 0); + By = SpectralField(spectralspace_ba, dm, 1, 0); + Bz = SpectralField(spectralspace_ba, dm, 1, 0); + Jx = SpectralField(spectralspace_ba, dm, 1, 0); + Jy = SpectralField(spectralspace_ba, dm, 1, 0); + Jz = SpectralField(spectralspace_ba, dm, 1, 0); + rho_old = SpectralField(spectralspace_ba, dm, 1, 0); + rho_new = SpectralField(spectralspace_ba, dm, 1, 0); + + // Allocate temporary arrays - over different boxarrays + tmpRealField = SpectralField(realspace_ba, dm, 1, 0); + tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); + + // Allocate and initialize the FFT plans + forward_plan = FFTplans(spectralspace_ba, dm); + backward_plan = FFTplans(spectralspace_ba, dm); + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + Box bx = spectralspace_ba[mfi]; +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code +#else + // Create FFTW plans + forward_plan[mfi] = fftw_plan_dft_3d( + // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + bx.length(2), bx.length(1), bx.length(0), + reinterpret_cast( tmpRealField[mfi].dataPtr() ), + reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), + FFTW_FORWARD, FFTW_ESTIMATE ); + backward_plan[mfi] = fftw_plan_dft_3d( + // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + bx.length(2), bx.length(1), bx.length(0), + reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), + reinterpret_cast( tmpRealField[mfi].dataPtr() ), + FFTW_BACKWARD, FFTW_ESTIMATE ); + // TODO: Add 2D code + // TODO: Do real-to-complex transform instead of complex-to-complex + // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE +#endif + } +} + + +SpectralFieldData::~SpectralFieldData() +{ + for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code +#else + // Destroy FFTW plans + fftw_destroy_plan( forward_plan[mfi] ); + fftw_destroy_plan( backward_plan[mfi] ); +#endif + } +} + +/* TODO: Documentation + * Example: ForwardTransform( Efield_cp[0], SpectralFieldIndex::Ex ) + */ +void +SpectralFieldData::ForwardTransform( const MultiFab& mf, const int field_index ) +{ + // Loop over boxes + for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ + + // Copy the real-space field `mf` to the temporary field `tmpRealField` + // This ensures that all fields have the same number of points + // before the Fourier transform. + // As a consequence, the copy discards the *last* point of `mf` + // in any direction that has *nodal* index type. + { + Box bx = mf[mfi].box(); + const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); + Array4 mf_arr = mf[mfi].array(); + Array4 tmp_arr = tmpRealField[mfi].array(); + ParallelFor( realspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + tmp_arr(i,j,k) = mf_arr(i,j,k); + }); + } + + // Perform Fourier transform from `tmpRealField` to `tmpSpectralField` +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code ; make sure that this is done on the same + // GPU stream as the above copy +#else + fftw_execute( forward_plan[mfi] ); +#endif + + // Copy the spectral-space field `tmpSpectralField` to the appropriate field + // (specified by the input argument field_index ) + { + SpectralField& field = getSpectralField( field_index ); + Array4 field_arr = field[mfi].array(); + Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + field_arr(i,j,k) = tmp_arr(i,j,k); + }); + } + } +} + + +/* TODO: Documentation + */ +void +SpectralFieldData::BackwardTransform( MultiFab& mf, const int field_index ) +{ + // Loop over boxes + for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ + + // Copy the appropriate field (specified by the input argument field_index) + // to the spectral-space field `tmpSpectralField` + { + SpectralField& field = getSpectralField( field_index ); + Array4 field_arr = field[mfi].array(); + Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + tmp_arr(i,j,k) = field_arr(i,j,k); + }); + } + + // Perform Fourier transform from `tmpSpectralField` to `tmpRealField` +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code ; make sure that this is done on the same + // GPU stream as the above copy +#else + fftw_execute( backward_plan[mfi] ); +#endif + + // Copy the temporary field `tmpRealField` to the real-space field `mf` + // The copy does *not* fill the *last* point of `mf` + // in any direction that has *nodal* index type (but this point is + // in the guard cells and will be filled by guard cell exchange) + { + Box bx = mf[mfi].box(); + const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); + Array4 mf_arr = mf[mfi].array(); + Array4 tmp_arr = tmpRealField[mfi].array(); + ParallelFor( realspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + mf_arr(i,j,k) = tmp_arr(i,j,k).real(); + }); + } + } +} + + +SpectralField& +SpectralFieldData::getSpectralField( const int field_index ) +{ + switch(field_index) + { + case SpectralFieldIndex::Ex : return Ex; + case SpectralFieldIndex::Ey : return Ey; + case SpectralFieldIndex::Ez : return Ez; + case SpectralFieldIndex::Bx : return Bx; + case SpectralFieldIndex::By : return By; + case SpectralFieldIndex::Bz : return Bz; + case SpectralFieldIndex::Jx : return Jx; + case SpectralFieldIndex::Jy : return Jy; + case SpectralFieldIndex::Jz : return Jz; + case SpectralFieldIndex::rho_old : return rho_old; + case SpectralFieldIndex::rho_new : return rho_new; + default : return tmpSpectralField; // For synthax; should not occur in practice + } +} -- cgit v1.2.3 From eec0e25fc2f048a9278a7aa9911bd2e0abd4604b Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 11:09:14 -0700 Subject: Fix runtime errors --- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 12 +++++++----- Source/WarpX.cpp | 3 +++ 2 files changed, 10 insertions(+), 5 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 25c8f438b..66e1e3470 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -56,14 +56,16 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, SpectralFieldData::~SpectralFieldData() { - for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ + if (tmpRealField.size() > 0){ + for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + // Add cuFFT-specific code #else - // Destroy FFTW plans - fftw_destroy_plan( forward_plan[mfi] ); - fftw_destroy_plan( backward_plan[mfi] ); + // Destroy FFTW plans + fftw_destroy_plan( forward_plan[mfi] ); + fftw_destroy_plan( backward_plan[mfi] ); #endif + } } } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index cb8803c1b..7f831c542 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -206,6 +206,9 @@ WarpX::WarpX () dataptr_fp_fft.resize(nlevs_max); dataptr_cp_fft.resize(nlevs_max); + spectral_solver_fp.resize(nlevs_max); + spectral_solver_cp.resize(nlevs_max); + ba_valid_fp_fft.resize(nlevs_max); ba_valid_cp_fft.resize(nlevs_max); -- cgit v1.2.3 From 8bffe04d8eb683230a6b9919438adf9239c87fec Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 20:03:54 -0700 Subject: Transform rho_old and rho_new --- Source/FieldSolver/SpectralSolver/SpectralFieldData.H | 6 ++++-- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 10 ++++++---- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 10 ++++++---- Source/FieldSolver/WarpXFFT.cpp | 3 +++ Source/WarpX.cpp | 2 +- 5 files changed, 20 insertions(+), 11 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index a74bfe22b..2abe81889 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -33,8 +33,10 @@ class SpectralFieldData SpectralFieldData() = default; // Default constructor SpectralFieldData& operator=(SpectralFieldData&& field_data) = default; // Default move assignment ~SpectralFieldData(); - void ForwardTransform( const amrex::MultiFab& mf, const int field_index ); - void BackwardTransform( amrex::MultiFab& mf, const int field_index ); + void ForwardTransform( const amrex::MultiFab& mf, + const int field_index, const int i_comp ); + void BackwardTransform( amrex::MultiFab& mf, + const int field_index, const int i_comp ); private: SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 66e1e3470..ee6dbff6a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -73,7 +73,8 @@ SpectralFieldData::~SpectralFieldData() * Example: ForwardTransform( Efield_cp[0], SpectralFieldIndex::Ex ) */ void -SpectralFieldData::ForwardTransform( const MultiFab& mf, const int field_index ) +SpectralFieldData::ForwardTransform( const MultiFab& mf, + const int field_index, const int i_comp ) { // Loop over boxes for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ @@ -91,7 +92,7 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, const int field_index ) Array4 tmp_arr = tmpRealField[mfi].array(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - tmp_arr(i,j,k) = mf_arr(i,j,k); + tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); }); } @@ -122,7 +123,8 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, const int field_index ) /* TODO: Documentation */ void -SpectralFieldData::BackwardTransform( MultiFab& mf, const int field_index ) +SpectralFieldData::BackwardTransform( MultiFab& mf, + const int field_index, const int i_comp ) { // Loop over boxes for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ @@ -160,7 +162,7 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, const int field_index ) Array4 tmp_arr = tmpRealField[mfi].array(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k) = tmp_arr(i,j,k).real(); + mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); }); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 521e558ba..a471666b9 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -30,11 +30,13 @@ class SpectralSolver // - Initialize arrays for fields in Fourier space + FFT plans field_data = SpectralFieldData( realspace_ba, k_space, dm ); }; - void ForwardTransform( const amrex::MultiFab& mf, const int field_index ){ - field_data.ForwardTransform( mf, field_index ); + void ForwardTransform( const amrex::MultiFab& mf, + const int field_index, const int i_comp=0 ){ + field_data.ForwardTransform( mf, field_index, i_comp ); }; - void BackwardTransform( amrex::MultiFab& mf, const int field_index ){ - field_data.BackwardTransform( mf, field_index ); + void BackwardTransform( amrex::MultiFab& mf, + const int field_index, const int i_comp=0 ){ + field_data.BackwardTransform( mf, field_index, i_comp ); }; void pushSpectralFields(){ algorithm.pushSpectralFields( field_data ); diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 2e146b2a1..1f7fe6cd0 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -455,6 +455,9 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) solver.ForwardTransform(*current_fp_fft[lev][0], SpectralFieldIndex::Jx); solver.ForwardTransform(*current_fp_fft[lev][1], SpectralFieldIndex::Jy); solver.ForwardTransform(*current_fp_fft[lev][2], SpectralFieldIndex::Jz); + solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_old, 0); + solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_new, 1); + // TODO: Transform rho // Advance fields in spectral space diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 7f831c542..1bee1263f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -208,7 +208,7 @@ WarpX::WarpX () spectral_solver_fp.resize(nlevs_max); spectral_solver_cp.resize(nlevs_max); - + ba_valid_fp_fft.resize(nlevs_max); ba_valid_cp_fft.resize(nlevs_max); -- cgit v1.2.3 From 0b77b789043461c77780ab56ab46a1e564b98d6e Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 20:23:10 -0700 Subject: Normalize inverse FFT --- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index ee6dbff6a..595cbac16 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -154,15 +154,18 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // The copy does *not* fill the *last* point of `mf` // in any direction that has *nodal* index type (but this point is // in the guard cells and will be filled by guard cell exchange) + // Normalize (divide by 1/N) since the FFT result in a factor N { Box bx = mf[mfi].box(); const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); + // For normalization: divide by the number of points in the box + const Real inv_N = 1./(bx.length(0)*bx.length(1)*bx.length(2)); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); + mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k).real(); }); } } -- cgit v1.2.3 From d96e23d67f6b181657169453daef7c37cb62e635 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 21:50:52 -0700 Subject: Implement spectral solver in 2D --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 5 ++++- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 18 ++++++++++++++++++ .../FieldSolver/SpectralSolver/SpectralFieldData.cpp | 19 ++++++++++++------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 11 +++++++---- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 10 ++++++++-- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 2 +- Source/FieldSolver/WarpXFFT.cpp | 4 +++- 7 files changed, 53 insertions(+), 16 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 375fa48bb..b54f5c143 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -21,7 +21,10 @@ class PsatdAlgorithm private: // Modified finite-order vectors - SpectralKVector modified_kx_vec, modified_ky_vec, modified_kz_vec; + SpectralKVector modified_kx_vec, modified_kz_vec; +#if (AMREX_SPACEDIM==3) + SpectralKVector modified_ky_vec; +#endif SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; }; diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index acc74148e..75df01343 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -13,7 +13,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // Allocate the 1D vectors modified_kx_vec = SpectralKVector( ba, dm ); +#if (AMREX_SPACEDIM==3) modified_ky_vec = SpectralKVector( ba, dm ); +#endif modified_kz_vec = SpectralKVector( ba, dm ); // Allocate and fill them by computing the modified vector for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ @@ -21,12 +23,18 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, ComputeModifiedKVector( modified_kx_vec[mfi], spectral_kspace.kx_vec[mfi], bx, spectral_kspace.dx[0], norder_x ); +#if (AMREX_SPACEDIM==3) ComputeModifiedKVector( modified_ky_vec[mfi], spectral_kspace.ky_vec[mfi], bx, spectral_kspace.dx[1], norder_y ); ComputeModifiedKVector( modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], bx, spectral_kspace.dx[2], norder_z ); +#else + ComputeModifiedKVector( + modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], + bx, spectral_kspace.dx[1], norder_z ); +#endif } // Allocate the arrays of coefficients @@ -44,7 +52,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // Extract pointers for the k vectors const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); +#endif const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); // Extract arrays for the coefficients Array4 C = C_coef[mfi].array(); @@ -60,7 +70,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // Calculate norm of vector const Real k_norm = std::sqrt( std::pow( modified_kx[i], 2 ) + +#if (AMREX_SPACEDIM==3) std::pow( modified_ky[j], 2 ) + +#endif std::pow( modified_kz[k], 2 ) ); // Calculate coefficients @@ -112,7 +124,9 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ Array4 X3_arr = X3_coef[mfi].array(); // Extract pointers for the k vectors const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); // Loop over indices within one box @@ -134,7 +148,11 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ const Complex rho_new = rho_new_arr(i,j,k); // k vector values, and coefficients const Real kx = modified_kx_arr[i]; +#if (AMREX_SPACEDIM==3) const Real ky = modified_ky_arr[j]; +#else + constexpr Real ky = 0; +#endif const Real kz = modified_kz_arr[k]; constexpr Real c2 = PhysConst::c*PhysConst::c; constexpr Real inv_ep0 = 1./PhysConst::ep0; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 595cbac16..25c08a367 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -34,19 +34,24 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // Add cuFFT-specific code #else // Create FFTW plans - forward_plan[mfi] = fftw_plan_dft_3d( - // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order - bx.length(2), bx.length(1), bx.length(0), + forward_plan[mfi] = +#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), +#else + fftw_plan_dft_2d( bx.length(1), bx.length(0), +#endif reinterpret_cast( tmpRealField[mfi].dataPtr() ), reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), FFTW_FORWARD, FFTW_ESTIMATE ); - backward_plan[mfi] = fftw_plan_dft_3d( - // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order - bx.length(2), bx.length(1), bx.length(0), + backward_plan[mfi] = +#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), +#else + fftw_plan_dft_2d( bx.length(1), bx.length(0), +#endif reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), reinterpret_cast( tmpRealField[mfi].dataPtr() ), FFTW_BACKWARD, FFTW_ESTIMATE ); - // TODO: Add 2D code // TODO: Do real-to-complex transform instead of complex-to-complex // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index a790be94f..16ec95cf5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -14,15 +14,18 @@ class SpectralKSpace public: SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, - const amrex::Array dx ); + const amrex::RealVect dx ); amrex::BoxArray spectralspace_ba; - SpectralKVector kx_vec, ky_vec, kz_vec; - amrex::Array dx; + SpectralKVector kx_vec, kz_vec; +#if (AMREX_SPACEDIM==3) + SpectralKVector ky_vec; +#endif + amrex::RealVect dx; }; void AllocateAndFillKvector( amrex::Gpu::ManagedVector& k, - const amrex::Box& bx, const amrex::Array dx, const int i_dim ); + const amrex::Box& bx, const amrex::RealVect dx, const int i_dim ); void ComputeModifiedKVector( amrex::Gpu::ManagedVector& modified_k, diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index d05748192..76e8aef15 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -6,7 +6,7 @@ using namespace Gpu; SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, - const Array realspace_dx ) + const RealVect realspace_dx ) { // Store the cell size dx = realspace_dx; @@ -25,20 +25,26 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Allocate the 1D vectors kx_vec = SpectralKVector( spectralspace_ba, dm ); +#if (AMREX_SPACEDIM == 3) ky_vec = SpectralKVector( spectralspace_ba, dm ); +#endif kz_vec = SpectralKVector( spectralspace_ba, dm ); // Initialize the values on each box for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; AllocateAndFillKvector( kx_vec[mfi], bx, dx, 0 ); +#if (AMREX_SPACEDIM == 3) AllocateAndFillKvector( ky_vec[mfi], bx, dx, 1 ); AllocateAndFillKvector( kz_vec[mfi], bx, dx, 2 ); +#else + AllocateAndFillKvector( kz_vec[mfi], bx, dx, 1 ); +#endif } } void AllocateAndFillKvector( ManagedVector& k, const Box& bx, - const Array dx, const int i_dim ) + const RealVect dx, const int i_dim ) { // Alllocate k to the right size int N = bx.length( i_dim ); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index a471666b9..06d53a2a9 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -19,7 +19,7 @@ class SpectralSolver SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, - const amrex::Array dx, const amrex::Real dt ) + const amrex::RealVect dx, const amrex::Real dt ) { // Initialize all structures using the same distribution mapping dm // - Initialize k space (Contains size of each box in spectral space, diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 1f7fe6cd0..e655bd11f 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -129,8 +129,10 @@ WarpX::AllocLevelDataFFT (int lev) if (fft_hybrid_mpi_decomposition == false){ // Allocate and initialize objects for the spectral solver // (all use the same distribution mapping) + std::array dx = CellSize(lev); + RealVect dx_vect = RealVect( AMREX_D_DECL(dx[0], dx[1], dx[2]) ); spectral_solver_fp[lev].reset( new SpectralSolver( ba_fp_fft, dm_fp_fft, - nox_fft, noy_fft, noz_fft, CellSize(lev), dt[lev] ) ); + nox_fft, noy_fft, noz_fft, dx_vect, dt[lev] ) ); } // rho2 has one extra ghost cell, so that it's safe to deposit charge density after -- cgit v1.2.3 From d9415cd662ec6930569a35c3fc4c9040aae0514a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 22:14:48 -0700 Subject: Fix error with 2D spectral solver --- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 25c08a367..6b39e9fa4 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -167,7 +167,7 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); // For normalization: divide by the number of points in the box - const Real inv_N = 1./(bx.length(0)*bx.length(1)*bx.length(2)); + const Real inv_N = 1./bx.numPts(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k).real(); -- cgit v1.2.3 From 80787133443c9adebdf6a9c6cdc6538bb2bcd2df Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 20 Apr 2019 22:53:00 -0700 Subject: Added shift function --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 8 ++-- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 3 ++ .../SpectralSolver/SpectralFieldData.cpp | 8 ++++ Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 15 ++++-- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 55 +++++++++++++++++----- 5 files changed, 69 insertions(+), 20 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 0b1dbc4da..95c637907 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -12,12 +12,12 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 0, norder_x ); + modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 0, norder_x ); #if (AMREX_SPACEDIM==3) - modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_y ); - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 2, norder_z ); + modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_y ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 2, norder_z ); #else - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_z ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_z ); #endif // Allocate the arrays of coefficients diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 2abe81889..052fc045e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -42,6 +42,9 @@ class SpectralFieldData SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform FFTplans forward_plan, backward_plan; + // Factors that shift the fields between nodal and cell-centered, in spectral space + // (for each dimension of space) + amrex::Array shift_N2C, shift_C2N; SpectralField& getSpectralField( const int field_index ); }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 6b39e9fa4..78b7d7c67 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -25,6 +25,14 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, tmpRealField = SpectralField(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); + // Allocate the vectors that allow to shift between nodal and cell-centered + for (int i_dim=0; i_dim #include #include using KVectorComponent = amrex::LayoutData< amrex::Gpu::ManagedVector >; +using SpectralShiftFactor = amrex::LayoutData< amrex::Gpu::ManagedVector >; + +struct ShiftType { + enum{ NodalToCentered=0, CenteredToNodal=1 }; +}; /* TODO Documentation: class represent k space, and how it is distribution * for local FFT on each MPI: k spaces are not connected. @@ -12,14 +18,17 @@ using KVectorComponent = amrex::LayoutData< amrex::Gpu::ManagedVector k_vec; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 6e6115cec..94b1486af 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -1,5 +1,6 @@ #include #include +#include using namespace amrex; using namespace Gpu; @@ -25,19 +26,19 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Allocate the components of the k vector: kx, ky (only in 3D), kz for (int i_dim=0; i_dim& k = k_vec[mfi]; + ManagedVector& k = k_comp[mfi]; // Allocate k to the right size int N = bx.length( i_dim ); @@ -62,19 +63,19 @@ SpectralKSpace::AllocateAndFillKVector( const DistributionMapping& dm, const int // and distinguish the size of the local box and that of the global FFT // This will also be different for the real-to-complex transform } - return k_vec; + return k_comp; } -KVectorComponent& -SpectralKSpace::AllocateAndFillModifiedKVector( +KVectorComponent +SpectralKSpace::AllocateAndFillModifiedKComponent( const DistributionMapping& dm, const int i_dim, const int order ) const { - // Initialize an empty vector in each box - KVectorComponent modified_k_vec = KVectorComponent( spectralspace_ba, dm ); + // Initialize an empty ManagedVector in each box + KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); // Loop over boxes for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ const ManagedVector& k = k_vec[i_dim][mfi]; - ManagedVector& modified_k = modified_k_vec[mfi]; + ManagedVector& modified_k = modified_k_comp[mfi]; // Allocate modified_k to the same size as k modified_k.resize( k.size() ); @@ -86,5 +87,33 @@ SpectralKSpace::AllocateAndFillModifiedKVector( modified_k[i] = k[i]; } } - return modified_k_vec; + return modified_k_comp; +} + +SpectralShiftFactor +SpectralKSpace::AllocateAndFillSpectralShiftFactor( + 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& k = k_vec[i_dim][mfi]; + ManagedVector& 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.; + case ShiftType::NodalToCentered: sign = 1.; + } + constexpr Complex I{0,1}; + for (int i=0; i Date: Sun, 21 Apr 2019 06:16:35 -0700 Subject: Add break statements --- .../SpectralSolver/SpectralFieldData.cpp | 24 +++++++++++----------- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 2 +- 2 files changed, 13 insertions(+), 13 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 78b7d7c67..3fd7177e5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -30,7 +30,7 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, shift_C2N[i_dim] = k_space.AllocateAndFillSpectralShiftFactor( dm, i_dim, ShiftType::CenteredToNodal ); shift_N2C[i_dim] = k_space.AllocateAndFillSpectralShiftFactor( - dm, i_dim, ShiftType::NodalToCentered ); + dm, i_dim, ShiftType::NodalToCentered ); } // Allocate and initialize the FFT plans @@ -190,17 +190,17 @@ SpectralFieldData::getSpectralField( const int field_index ) { switch(field_index) { - case SpectralFieldIndex::Ex : return Ex; - case SpectralFieldIndex::Ey : return Ey; - case SpectralFieldIndex::Ez : return Ez; - case SpectralFieldIndex::Bx : return Bx; - case SpectralFieldIndex::By : return By; - case SpectralFieldIndex::Bz : return Bz; - case SpectralFieldIndex::Jx : return Jx; - case SpectralFieldIndex::Jy : return Jy; - case SpectralFieldIndex::Jz : return Jz; - case SpectralFieldIndex::rho_old : return rho_old; - case SpectralFieldIndex::rho_new : return rho_new; + case SpectralFieldIndex::Ex : return Ex; break; + case SpectralFieldIndex::Ey : return Ey; break; + case SpectralFieldIndex::Ez : return Ez; break; + case SpectralFieldIndex::Bx : return Bx; break; + case SpectralFieldIndex::By : return By; break; + case SpectralFieldIndex::Bz : return Bz; break; + case SpectralFieldIndex::Jx : return Jx; break; + case SpectralFieldIndex::Jy : return Jy; break; + case SpectralFieldIndex::Jz : return Jz; break; + case SpectralFieldIndex::rho_old : return rho_old; break; + case SpectralFieldIndex::rho_new : return rho_new; break; default : return tmpSpectralField; // For synthax; should not occur in practice } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 94b1486af..a7f726c8a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -107,7 +107,7 @@ SpectralKSpace::AllocateAndFillSpectralShiftFactor( // Fill the shift coefficients Real sign = 0; switch (shift_type){ - case ShiftType::CenteredToNodal: sign = -1.; + case ShiftType::CenteredToNodal: sign = -1.; break; case ShiftType::NodalToCentered: sign = 1.; } constexpr Complex I{0,1}; -- cgit v1.2.3 From 409775bcbe46b293a886808c24d404e3f37f547b Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 21 Apr 2019 07:00:59 -0700 Subject: Apply shift factors --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 8 ++-- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 6 ++- .../SpectralSolver/SpectralFieldData.cpp | 45 +++++++++++++++++----- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 9 ++--- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 8 ++-- 5 files changed, 52 insertions(+), 24 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 95c637907..37d3b33ee 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -12,12 +12,12 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 0, norder_x ); + modified_kx_vec = spectral_kspace.getModifiedKComponent( dm, 0, norder_x ); #if (AMREX_SPACEDIM==3) - modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_y ); - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 2, norder_z ); + modified_ky_vec = spectral_kspace.getModifiedKComponent( dm, 1, norder_y ); + modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 2, norder_z ); #else - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_z ); + modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 1, norder_z ); #endif // Allocate the arrays of coefficients diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 052fc045e..63a7c7520 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -43,8 +43,10 @@ class SpectralFieldData SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform FFTplans forward_plan, backward_plan; // Factors that shift the fields between nodal and cell-centered, in spectral space - // (for each dimension of space) - amrex::Array shift_N2C, shift_C2N; + SpectralShiftFactor xshift_N2C, xshift_C2N, zshift_N2C, zshift_C2N; +#if (AMREX_SPACEDIM==3) + SpectralShiftFactor yshift_N2C, yshift_C2N; +#endif SpectralField& getSpectralField( const int field_index ); }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 3fd7177e5..7d712ba03 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -25,13 +25,18 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, tmpRealField = SpectralField(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); - // Allocate the vectors that allow to shift between nodal and cell-centered - for (int i_dim=0; i_dim field_arr = field[mfi].array(); Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Complex* xshift_C2N_arr = xshift_C2N[mfi].dataPtr(); + const Complex* yshift_C2N_arr = yshift_C2N[mfi].dataPtr(); + const Complex* zshift_C2N_arr = zshift_C2N[mfi].dataPtr(); + // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - tmp_arr(i,j,k) = field_arr(i,j,k); + Complex spectral_field_value = field_arr(i,j,k); + // Apply proper shift in each dimension + if (is_nodal_x==false) spectral_field_value *= xshift_C2N_arr[i]; +#if (AMREX_SPACEDIM == 3) + if (is_nodal_y==false) spectral_field_value *= yshift_C2N_arr[j]; +#endif + if (is_nodal_z==false) spectral_field_value *= zshift_C2N_arr[k]; + // Copy field into temporary array + tmp_arr(i,j,k) = spectral_field_value; }); } @@ -170,7 +196,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // Normalize (divide by 1/N) since the FFT result in a factor N { Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + const Box realspace_bx = bx.enclosedCells(); + // `enclosedells` discards last point in each nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 90755bc0d..2f0681c3b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -22,13 +22,12 @@ class SpectralKSpace SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const amrex::RealVect dx ); - KVectorComponent AllocateAndFillKComponent( + KVectorComponent getKComponent( const amrex::DistributionMapping& dm, const int i_dim ) const; - KVectorComponent AllocateAndFillModifiedKComponent( + KVectorComponent getModifiedKComponent( const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; - SpectralShiftFactor AllocateAndFillSpectralShiftFactor( - const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; - + SpectralShiftFactor getSpectralShiftFactor( + const amrex::DistributionMapping& dm, const int i_dim, const int shift_type ) const; private: amrex::Array k_vec; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index a7f726c8a..90c20b01d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -26,12 +26,12 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Allocate the components of the k vector: kx, ky (only in 3D), kz for (int i_dim=0; i_dim Date: Sun, 21 Apr 2019 07:22:17 -0700 Subject: Fix a few bug --- .../SpectralSolver/SpectralFieldData.cpp | 39 ++++++++++++++++++---- 1 file changed, 32 insertions(+), 7 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 7d712ba03..68c9b7078 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -94,6 +94,15 @@ void SpectralFieldData::ForwardTransform( const MultiFab& mf, const int field_index, const int i_comp ) { + // Check field index type, in order to apply proper shift in spectral space + const bool is_nodal_x = mf.is_nodal(0); +#if (AMREX_SPACEDIM == 3) + const bool is_nodal_y = mf.is_nodal(1); + const bool is_nodal_z = mf.is_nodal(2); +#else + const bool is_nodal_z = mf.is_nodal(1); +#endif + // Loop over boxes for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ @@ -128,10 +137,24 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, SpectralField& field = getSpectralField( field_index ); Array4 field_arr = field[mfi].array(); Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Complex* xshift_C2N_arr = xshift_C2N[mfi].dataPtr(); +#if (AMREX_SPACEDIM == 3) + const Complex* yshift_C2N_arr = yshift_C2N[mfi].dataPtr(); +#endif + const Complex* zshift_C2N_arr = zshift_C2N[mfi].dataPtr(); + // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - field_arr(i,j,k) = tmp_arr(i,j,k); + Complex spectral_field_value = field_arr(i,j,k); + // Apply proper shift in each dimension + if (is_nodal_x==false) spectral_field_value *= xshift_C2N_arr[i]; +#if (AMREX_SPACEDIM == 3) + if (is_nodal_y==false) spectral_field_value *= yshift_C2N_arr[j]; +#endif + if (is_nodal_z==false) spectral_field_value *= zshift_C2N_arr[k]; + // Copy field into temporary array + field_arr(i,j,k) = spectral_field_value; }); } } @@ -162,20 +185,22 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, SpectralField& field = getSpectralField( field_index ); Array4 field_arr = field[mfi].array(); Array4 tmp_arr = tmpSpectralField[mfi].array(); - const Complex* xshift_C2N_arr = xshift_C2N[mfi].dataPtr(); - const Complex* yshift_C2N_arr = yshift_C2N[mfi].dataPtr(); - const Complex* zshift_C2N_arr = zshift_C2N[mfi].dataPtr(); + const Complex* xshift_N2C_arr = xshift_N2C[mfi].dataPtr(); +#if (AMREX_SPACEDIM == 3) + const Complex* yshift_N2C_arr = yshift_N2C[mfi].dataPtr(); +#endif + const Complex* zshift_N2C_arr = zshift_N2C[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = field_arr(i,j,k); // Apply proper shift in each dimension - if (is_nodal_x==false) spectral_field_value *= xshift_C2N_arr[i]; + if (is_nodal_x==false) spectral_field_value *= xshift_N2C_arr[i]; #if (AMREX_SPACEDIM == 3) - if (is_nodal_y==false) spectral_field_value *= yshift_C2N_arr[j]; + if (is_nodal_y==false) spectral_field_value *= yshift_N2C_arr[j]; #endif - if (is_nodal_z==false) spectral_field_value *= zshift_C2N_arr[k]; + if (is_nodal_z==false) spectral_field_value *= zshift_N2C_arr[k]; // Copy field into temporary array tmp_arr(i,j,k) = spectral_field_value; }); -- cgit v1.2.3 From bd11c215d0a2fab14e2a5c6a82a7a8a1b811e1c8 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 21 Apr 2019 08:15:57 -0700 Subject: Remove bug in Fourier transform --- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 68c9b7078..31390127e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -146,7 +146,7 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, const Box spectralspace_bx = tmpSpectralField[mfi].box(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Complex spectral_field_value = field_arr(i,j,k); + Complex spectral_field_value = tmp_arr(i,j,k); // Apply proper shift in each dimension if (is_nodal_x==false) spectral_field_value *= xshift_C2N_arr[i]; #if (AMREX_SPACEDIM == 3) -- cgit v1.2.3 From 37f932323e7e3381cf7eee72a9e45e0304754e10 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 23 Apr 2019 13:19:08 -0700 Subject: Switch and rename the shift factors --- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 10 ++-- .../SpectralSolver/SpectralFieldData.cpp | 67 +++++++++++++--------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 9 +-- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 58 ++++++++++--------- 4 files changed, 83 insertions(+), 61 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 63a7c7520..be1765ca0 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -31,7 +31,7 @@ class SpectralFieldData const SpectralKSpace& k_space, const amrex::DistributionMapping& dm ); SpectralFieldData() = default; // Default constructor - SpectralFieldData& operator=(SpectralFieldData&& field_data) = default; // Default move assignment + SpectralFieldData& operator=(SpectralFieldData&& field_data) = default; ~SpectralFieldData(); void ForwardTransform( const amrex::MultiFab& mf, const int field_index, const int i_comp ); @@ -42,10 +42,12 @@ class SpectralFieldData SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform FFTplans forward_plan, backward_plan; - // Factors that shift the fields between nodal and cell-centered, in spectral space - SpectralShiftFactor xshift_N2C, xshift_C2N, zshift_N2C, zshift_C2N; + // Correcting "shift" factors when performing FFT from/to + // a cell-centered grid in real space, instead of a nodal grid + SpectralShiftFactor xshift_FFTfromCell, xshift_FFTtoCell, + zshift_FFTfromCell, zshift_FFTtoCell; #if (AMREX_SPACEDIM==3) - SpectralShiftFactor yshift_N2C, yshift_C2N; + SpectralShiftFactor yshift_FFTfromCell, yshift_FFTtoCell; #endif SpectralField& getSpectralField( const int field_index ); }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 31390127e..ca5e338e5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -21,21 +21,32 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, rho_old = SpectralField(spectralspace_ba, dm, 1, 0); rho_new = SpectralField(spectralspace_ba, dm, 1, 0); - // Allocate temporary arrays - over different boxarrays + // Allocate temporary arrays - in real space and spectral space + // These arrays will store the data just before/after the FFT tmpRealField = SpectralField(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); - // Allocate the coefficients that allow to shift between nodal and cell-centered - xshift_C2N = k_space.getSpectralShiftFactor(dm, 0, ShiftType::CenteredToNodal); - xshift_N2C = k_space.getSpectralShiftFactor(dm, 0, ShiftType::NodalToCentered); + // By default, we assume the FFT is done from/to a nodal grid in real space + // It the FFT is performed from/to a cell-centered grid in real space, + // a correcting "shift" factor must be applied in spectral space. + xshift_FFTfromCell = k_space.getSpectralShiftFactor(dm, 0, + ShiftType::TransformFromCellCentered); + xshift_FFTtoCell = k_space.getSpectralShiftFactor(dm, 0, + ShiftType::TransformToCellCentered); #if (AMREX_SPACEDIM == 3) - yshift_C2N = k_space.getSpectralShiftFactor(dm, 1, ShiftType::CenteredToNodal); - yshift_N2C = k_space.getSpectralShiftFactor(dm, 1, ShiftType::NodalToCentered); - zshift_C2N = k_space.getSpectralShiftFactor(dm, 2, ShiftType::CenteredToNodal); - zshift_N2C = k_space.getSpectralShiftFactor(dm, 2, ShiftType::NodalToCentered); + yshift_FFTfromCell = k_space.getSpectralShiftFactor(dm, 1, + ShiftType::TransformFromCellCentered); + yshift_FFTtoCell = k_space.getSpectralShiftFactor(dm, 1, + ShiftType::TransformToCellCentered); + zshift_FFTfromCell = k_space.getSpectralShiftFactor(dm, 2, + ShiftType::TransformFromCellCentered); + zshift_FFTtoCell = k_space.getSpectralShiftFactor(dm, 2, + ShiftType::TransformToCellCentered); #else - zshift_C2N = k_space.getSpectralShiftFactor(dm, 1, ShiftType::CenteredToNodal); - zshift_N2C = k_space.getSpectralShiftFactor(dm, 1, ShiftType::NodalToCentered); + zshift_FFTfromCell = k_space.getSpectralShiftFactor(dm, 1, + ShiftType::TransformFromCellCentered); + zshift_FFTtoCell = k_space.getSpectralShiftFactor(dm, 1, + ShiftType::TransformToCellCentered); #endif // Allocate and initialize the FFT plans @@ -131,28 +142,30 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, fftw_execute( forward_plan[mfi] ); #endif - // Copy the spectral-space field `tmpSpectralField` to the appropriate field - // (specified by the input argument field_index ) + // Copy the spectral-space field `tmpSpectralField` to the appropriate + // field (specified by the input argument field_index ) + // and apply correcting shift factor if the real space data comes + // from a cell-centered grid in real space instead of a nodal grid. { SpectralField& field = getSpectralField( field_index ); Array4 field_arr = field[mfi].array(); Array4 tmp_arr = tmpSpectralField[mfi].array(); - const Complex* xshift_C2N_arr = xshift_C2N[mfi].dataPtr(); + const Complex* xshift_arr = xshift_FFTfromCell[mfi].dataPtr(); #if (AMREX_SPACEDIM == 3) - const Complex* yshift_C2N_arr = yshift_C2N[mfi].dataPtr(); + const Complex* yshift_arr = yshift_FFTfromCell[mfi].dataPtr(); #endif - const Complex* zshift_C2N_arr = zshift_C2N[mfi].dataPtr(); + const Complex* zshift_arr = zshift_FFTfromCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = tmp_arr(i,j,k); // Apply proper shift in each dimension - if (is_nodal_x==false) spectral_field_value *= xshift_C2N_arr[i]; + if (is_nodal_x==false) spectral_field_value *= xshift_arr[i]; #if (AMREX_SPACEDIM == 3) - if (is_nodal_y==false) spectral_field_value *= yshift_C2N_arr[j]; + if (is_nodal_y==false) spectral_field_value *= yshift_arr[j]; #endif - if (is_nodal_z==false) spectral_field_value *= zshift_C2N_arr[k]; + if (is_nodal_z==false) spectral_field_value *= zshift_arr[k]; // Copy field into temporary array field_arr(i,j,k) = spectral_field_value; }); @@ -179,28 +192,30 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // Loop over boxes for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ - // Copy the appropriate field (specified by the input argument field_index) - // to the spectral-space field `tmpSpectralField` + // Copy the spectral-space field `tmpSpectralField` to the appropriate + // field (specified by the input argument field_index ) + // and apply correcting shift factor if the field is to be transformed + // to a cell-centered grid in real space instead of a nodal grid. { SpectralField& field = getSpectralField( field_index ); Array4 field_arr = field[mfi].array(); Array4 tmp_arr = tmpSpectralField[mfi].array(); - const Complex* xshift_N2C_arr = xshift_N2C[mfi].dataPtr(); + const Complex* xshift_arr = xshift_FFTtoCell[mfi].dataPtr(); #if (AMREX_SPACEDIM == 3) - const Complex* yshift_N2C_arr = yshift_N2C[mfi].dataPtr(); + const Complex* yshift_arr = yshift_FFTtoCell[mfi].dataPtr(); #endif - const Complex* zshift_N2C_arr = zshift_N2C[mfi].dataPtr(); + const Complex* zshift_arr = zshift_FFTtoCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = field_arr(i,j,k); // Apply proper shift in each dimension - if (is_nodal_x==false) spectral_field_value *= xshift_N2C_arr[i]; + if (is_nodal_x==false) spectral_field_value *= xshift_arr[i]; #if (AMREX_SPACEDIM == 3) - if (is_nodal_y==false) spectral_field_value *= yshift_N2C_arr[j]; + if (is_nodal_y==false) spectral_field_value *= yshift_arr[j]; #endif - if (is_nodal_z==false) spectral_field_value *= zshift_N2C_arr[k]; + if (is_nodal_z==false) spectral_field_value *= zshift_arr[k]; // Copy field into temporary array tmp_arr(i,j,k) = spectral_field_value; }); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index fd77bb7b8..c71cf2643 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -5,11 +5,11 @@ #include #include -using KVectorComponent = amrex::LayoutData< amrex::Gpu::ManagedVector >; -using SpectralShiftFactor = amrex::LayoutData< amrex::Gpu::ManagedVector >; +using KVectorComponent = amrex::LayoutData>; +using SpectralShiftFactor = amrex::LayoutData>; struct ShiftType { - enum{ NodalToCentered=0, CenteredToNodal=1 }; + enum{ TransformFromCellCentered=0, TransformToCellCentered=1 }; }; /* TODO Documentation: class represent k space, and how it is distribution @@ -28,7 +28,8 @@ class SpectralKSpace const amrex::DistributionMapping& dm, const int i_dim, const int n_order, const bool nodal ) const; SpectralShiftFactor getSpectralShiftFactor( - const amrex::DistributionMapping& dm, const int i_dim, const int shift_type ) const; + const amrex::DistributionMapping& dm, const int i_dim, + const int shift_type ) const; private: amrex::Array k_vec; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 111643197..4852c7d2b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -16,10 +16,12 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, BoxList spectral_bl; // Create empty box list // Loop over boxes and fill the box list for (int i=0; i < realspace_ba.size(); i++ ) { - // For local FFTs, each box in spectral space starts at 0 in each direction - // and has the same number of points as the real space box (including guard cells) + // For local FFTs, boxes in spectral space start at 0 in each direction + // and have the same number of points as the real space box + // TODO: this will be different for the hybrid FFT scheme Box realspace_bx = realspace_ba[i]; - Box bx = Box( IntVect::TheZeroVector(), realspace_bx.bigEnd() - realspace_bx.smallEnd() ); + Box bx = Box( IntVect::TheZeroVector(), + realspace_bx.bigEnd() - realspace_bx.smallEnd() ); spectral_bl.push_back( bx ); } spectralspace_ba.define( spectral_bl ); @@ -31,7 +33,8 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } KVectorComponent -SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) const +SpectralKSpace::getKComponent( const DistributionMapping& dm, + const int i_dim ) const { // Initialize an empty ManagedVector in each box KVectorComponent k_comp = KVectorComponent(spectralspace_ba, dm); @@ -58,17 +61,16 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) for (int i=(N+1)/2; i getFonbergStencilCoefficients( const int n_order, const bool nodal ) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, "n_order should be even."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, + "n_order should be even."); const int m = n_order/2; - Vector stencil_coef; - stencil_coef.resize( m+1 ); + Vector coefs; + coefs.resize( m+1 ); // Coefficients for nodal (a.k.a. centered) finite-difference if (nodal == true) { - stencil_coef[0] = -2.; // First coefficient + coefs[0] = -2.; // First coefficient for (int n=1; n Date: Tue, 23 Apr 2019 16:12:11 -0700 Subject: Add comments ; copy only to valid box --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 7 +-- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 6 ++- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 15 ++++-- .../SpectralSolver/SpectralFieldData.cpp | 54 +++++++++++----------- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 51 ++++++++++++++++---- 5 files changed, 88 insertions(+), 45 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 16d27cab2..acefcc466 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -4,8 +4,9 @@ #include #include -/* TODO: Write documentation -*/ +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ class PsatdAlgorithm { using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; @@ -16,7 +17,7 @@ class PsatdAlgorithm const int norder_x, const int norder_y, const int norder_z, const bool nodal, const amrex::Real dt); PsatdAlgorithm() = default; // Default constructor - PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; // Default move assignment + PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; void pushSpectralFields(SpectralFieldData& f) const; private: diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 5ebb7144d..60e9d58c0 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -4,6 +4,7 @@ using namespace amrex; +/* \brief Initialize coefficients for the update equation */ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, @@ -28,7 +29,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, X3_coef = SpectralCoefficients(ba, dm, 1, 0); // Fill them with the right values: - // Loop over boxes + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ const Box& bx = ba[mfi]; @@ -78,6 +80,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, } }; +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ void PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index be1765ca0..c62513de6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -8,20 +8,23 @@ // Declare type for spectral fields using SpectralField = amrex::FabArray< amrex::BaseFab >; -/* Fields that will be stored in spectral space */ +/* Index for the fields that will be stored in spectral space */ struct SpectralFieldIndex{ enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new }; }; -/* \brief Class that stores the fields in spectral space - * and performs the spectral transforms to/from real space +/* \brief Class that stores the fields in spectral space, and performs the + * Fourier transforms between real space and spectral space */ class SpectralFieldData { friend class PsatdAlgorithm; + // Define the FFTplans type, which holds one fft plan per box + // (plans are only initialized for the boxes that are owned by + // the local MPI rank) #ifdef AMREX_USE_GPU -// Add cuFFT-specific code + // Add cuFFT-specific code #else using FFTplans = amrex::LayoutData; #endif @@ -40,7 +43,9 @@ class SpectralFieldData private: SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; - SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform + // tmpRealField and tmpSpectralField store fields + // right before/after the Fourier transform + SpectralField tmpRealField, tmpSpectralField; FFTplans forward_plan, backward_plan; // Correcting "shift" factors when performing FFT from/to // a cell-centered grid in real space, instead of a nodal grid diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index ca5e338e5..15652addc 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -2,6 +2,7 @@ using namespace amrex; +/* \brief Initialize fields in spectral space, and FFT plans */ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, const SpectralKSpace& k_space, const DistributionMapping& dm ) @@ -52,6 +53,8 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // Allocate and initialize the FFT plans forward_plan = FFTplans(spectralspace_ba, dm); backward_plan = FFTplans(spectralspace_ba, dm); + // Loop over boxes and allocate the corresponding plan + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; #ifdef AMREX_USE_GPU @@ -59,7 +62,8 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, #else // Create FFTW plans forward_plan[mfi] = -#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + // Swap dimensions: AMReX FAB are Fortran-order but FFTW is C-order +#if (AMREX_SPACEDIM == 3) fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), #else fftw_plan_dft_2d( bx.length(1), bx.length(0), @@ -68,7 +72,8 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), FFTW_FORWARD, FFTW_ESTIMATE ); backward_plan[mfi] = -#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + // Swap dimensions: AMReX FAB are Fortran-order but FFTW is C-order +#if (AMREX_SPACEDIM == 3) fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), #else fftw_plan_dft_2d( bx.length(1), bx.length(0), @@ -76,8 +81,6 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), reinterpret_cast( tmpRealField[mfi].dataPtr() ), FFTW_BACKWARD, FFTW_ESTIMATE ); - // TODO: Do real-to-complex transform instead of complex-to-complex - // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE #endif } } @@ -98,12 +101,13 @@ SpectralFieldData::~SpectralFieldData() } } -/* TODO: Documentation - * Example: ForwardTransform( Efield_cp[0], SpectralFieldIndex::Ex ) - */ +/* \brief Transform the component `i_comp` of MultiFab `mf` + * to spectral space, and store the corresponding result internally + * (in the spectral field specified by `field_index`) */ void SpectralFieldData::ForwardTransform( const MultiFab& mf, - const int field_index, const int i_comp ) + const int field_index, + const int i_comp ) { // Check field index type, in order to apply proper shift in spectral space const bool is_nodal_x = mf.is_nodal(0); @@ -123,8 +127,8 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, // As a consequence, the copy discards the *last* point of `mf` // in any direction that has *nodal* index type. { - Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + Box realspace_bx = mf[mfi].box(); // Copy the box + realspace_bx.enclosedCells(); // Discard last point in nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); @@ -174,8 +178,8 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, } -/* TODO: Documentation - */ +/* \brief Transform spectral field specified by `field_index` back to + * real space, and store it in the component `i_comp` of `mf` */ void SpectralFieldData::BackwardTransform( MultiFab& mf, const int field_index, const int i_comp ) @@ -193,9 +197,10 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ // Copy the spectral-space field `tmpSpectralField` to the appropriate - // field (specified by the input argument field_index ) + // field (specified by the input argument field_index) // and apply correcting shift factor if the field is to be transformed // to a cell-centered grid in real space instead of a nodal grid. + // Normalize (divide by 1/N) since the FFT+IFFT results in a factor N { SpectralField& field = getSpectralField( field_index ); Array4 field_arr = field[mfi].array(); @@ -207,6 +212,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, const Complex* zshift_arr = zshift_FFTtoCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); + // For normalization: divide by the number of points in the box + const Real inv_N = 1./spectralspace_bx.numPts(); ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = field_arr(i,j,k); @@ -216,8 +223,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, if (is_nodal_y==false) spectral_field_value *= yshift_arr[j]; #endif if (is_nodal_z==false) spectral_field_value *= zshift_arr[k]; - // Copy field into temporary array - tmp_arr(i,j,k) = spectral_field_value; + // Copy field into temporary array (after normalization) + tmp_arr(i,j,k) = inv_N*spectral_field_value; }); } @@ -230,22 +237,15 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, #endif // Copy the temporary field `tmpRealField` to the real-space field `mf` - // The copy does *not* fill the *last* point of `mf` - // in any direction that has *nodal* index type (but this point is - // in the guard cells and will be filled by guard cell exchange) - // Normalize (divide by 1/N) since the FFT result in a factor N + // in the *valid* part of the domain (guard cells are filled later, + // through guard cell exchange) { - Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); - // `enclosedells` discards last point in each nodal direction - AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); + const Box realspace_valid_bx = mfi.validbox(); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); - // For normalization: divide by the number of points in the box - const Real inv_N = 1./bx.numPts(); - ParallelFor( realspace_bx, + ParallelFor( realspace_valid_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k).real(); + mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); }); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 90866cd3b..c82e45c5e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -45,10 +45,9 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } } -/* For each box, in `spectralspace_ba`, which is owned - * by the local MPI proc (as indicated by the argument `dm`), - * compute the values of the corresponding k coordinate - * along the dimension specified by `i_dim` +/* For each box, in `spectralspace_ba`, which is owned by the local MPI rank + * (as indicated by the argument `dm`), compute the values of the + * corresponding k coordinate along the dimension specified by `i_dim` */ KVectorComponent SpectralKSpace::getKComponent( const DistributionMapping& dm, @@ -57,7 +56,7 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, // Initialize an empty ManagedVector in each box KVectorComponent k_comp = KVectorComponent(spectralspace_ba, dm); // Loop over boxes and allocate the corresponding ManagedVector - // for each box owned by the local MPI proc ("mfi.isValid") + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; ManagedVector& k = k_comp[mfi]; @@ -87,6 +86,15 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, return k_comp; } +/* For each box, in `spectralspace_ba`, which is owned by the local MPI rank + * (as indicated by the argument `dm`), compute the values of the + * corresponding correcting "shift" factor, along the dimension + * specified by `i_dim`. + * + * (By default, we assume the FFT is done from/to a nodal grid in real space + * It the FFT is performed from/to a cell-centered grid in real space, + * a correcting "shift" factor must be applied in spectral space.) + */ SpectralShiftFactor SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, const int i_dim, @@ -94,7 +102,8 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, { // Initialize an empty ManagedVector in each box SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); - // Loop over boxes + // Loop over boxes and allocate the corresponding ManagedVector + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ const ManagedVector& k = k_vec[i_dim][mfi]; ManagedVector& shift = shift_factor[mfi]; @@ -116,6 +125,19 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, return shift_factor; } +/* \brief For each box, in `spectralspace_ba`, which is owned by the local MPI + * rank (as indicated by the argument `dm`), compute the values of the + * corresponding finite-order modified k vector, along the + * dimension specified by `i_dim` + * + * The finite-order modified k vector is the spectral-space representation + * of a finite-order stencil in real space. + * + * \param n_order Order of accuracy of the stencil, in discretizing + * a spatial derivative + * \param nodal Whether the stencil is to be applied to a nodal or + staggered set of fields + */ KVectorComponent SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, const int i_dim, @@ -123,12 +145,13 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, const bool nodal ) const { // Initialize an empty ManagedVector in each box - KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); + KVectorComponent modified_k_comp = KVectorComponent(spectralspace_ba, dm); // Compute real-space stencil coefficients Vector stencil_coef = getFonbergStencilCoefficients(n_order, nodal); - // Loop over boxes + // Loop over boxes and allocate the corresponding ManagedVector + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Real delta_x = dx[i_dim]; const ManagedVector& k = k_vec[i_dim][mfi]; @@ -153,7 +176,12 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, return modified_k_comp; } -/* TODO: Documentation: point to Fonberg paper ; explain recurrence relation +/* Returns an array of coefficients, corresponding to the weight + * of each point in a finite-difference approximation (to order `n_order`) + * of a derivative. + * + * `nodal` indicates whether this finite-difference approximation is + * taken on a nodal grid or a staggered grid. */ Vector getFonbergStencilCoefficients( const int n_order, const bool nodal ) @@ -164,6 +192,11 @@ getFonbergStencilCoefficients( const int n_order, const bool nodal ) Vector coefs; coefs.resize( m+1 ); + // Note: there are closed-form formula for these coefficients, + // but they result in an overflow when evaluated numerically. + // One way to avoid the overflow is to calculate the coefficients + // by recurrence. + // Coefficients for nodal (a.k.a. centered) finite-difference if (nodal == true) { coefs[0] = -2.; // First coefficient -- cgit v1.2.3 From 8e5df6f8622a90c29ecfd4dad27971be76971cd9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 24 Apr 2019 10:35:05 -0700 Subject: Correct out-of-bounds with real-space box --- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 15652addc..291fe945e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -237,13 +237,11 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, #endif // Copy the temporary field `tmpRealField` to the real-space field `mf` - // in the *valid* part of the domain (guard cells are filled later, - // through guard cell exchange) { - const Box realspace_valid_bx = mfi.validbox(); + const Box realspace_bx = tmpRealField[mfi].box(); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); - ParallelFor( realspace_valid_bx, + ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); }); -- cgit v1.2.3