From 438590e9a37556ec5f0a029f0f88469e7ba3565a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 05:39:49 -0700 Subject: Change name of PsatdSolver to PsatdAlgorithm --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 29 ++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H new file mode 100644 index 000000000..a235315ae --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -0,0 +1,29 @@ +#ifndef WARPX_PSATD_ALGORITHM_H_ +#define WARPX_PSATD_ALGORITHM_H_ + +#include +#include + +using namespace amrex; +using namespace Gpu; + +/* TODO: Write documentation +*/ +class PsatdAlgorithm +{ + using SpectralCoefficients = FabArray>; + + public: + PsatdAlgorithm(const SpectralKSpace& spectral_kspace, + const DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const Real dt); + void pushSpectralFields( SpectralData& f ) const; + + private: + // Modified finite-order vectors + SpectralKVector modified_kx_vec, modified_ky_vec, modified_kz_vec; + SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; +}; + +#endif // WARPX_PSATD_ALGORITHM_H_ -- cgit v1.2.3 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 --- Source/FieldSolver/SpectralSolver/Make.package | 5 +- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 4 +- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 2 +- Source/FieldSolver/SpectralSolver/SpectralData.H | 45 ----- Source/FieldSolver/SpectralSolver/SpectralData.cpp | 184 -------------------- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 46 +++++ .../SpectralSolver/SpectralFieldData.cpp | 186 +++++++++++++++++++++ 7 files changed, 238 insertions(+), 234 deletions(-) delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralData.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralData.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralFieldData.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H') diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index 52328b259..50127914d 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,6 +1,7 @@ CEXE_headers += WarpX_Complex.H -CEXE_headers += SpectralData.H -CEXE_sources += SpectralData.cpp +CEXE_headers += SpectralSolver.H +CEXE_headers += SpectralFieldData.H +CEXE_sources += SpectralFieldData.cpp CEXE_headers += PsatdAlgorithm.H CEXE_sources += PsatdAlgorithm.cpp CEXE_headers += SpectralKSpace.H diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index a235315ae..7e94bdec6 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -2,7 +2,7 @@ #define WARPX_PSATD_ALGORITHM_H_ #include -#include +#include using namespace amrex; using namespace Gpu; @@ -18,7 +18,7 @@ class PsatdAlgorithm const DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const Real dt); - void pushSpectralFields( SpectralData& f ) const; + void pushSpectralFields( SpectralFieldData& f ) const; private: // Modified finite-order vectors diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 78a4ba191..acc74148e 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -84,7 +84,7 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, }; void -PsatdAlgorithm::pushSpectralFields( SpectralData& f ) const{ +PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ // Loop over boxes for ( MFIter mfi(f.Ex); mfi.isValid(); ++mfi ){ diff --git a/Source/FieldSolver/SpectralSolver/SpectralData.H b/Source/FieldSolver/SpectralSolver/SpectralData.H deleted file mode 100644 index 49764b96b..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralData.H +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef WARPX_PSATD_DATA_H_ -#define WARPX_PSATD_DATA_H_ - -#include -#include - -using namespace amrex; - -// Declare type for spectral fields -using SpectralField = FabArray>; - -/* 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 - */ -class SpectralData -{ - friend class PsatdAlgorithm; - -#ifdef AMREX_USE_GPU -// Add cuFFT-specific code -#else - using FFTplans = LayoutData; -#endif - - public: - SpectralData( const BoxArray& realspace_ba, - const BoxArray& spectralspace_ba, - const DistributionMapping& dm ); - ~SpectralData(); - void ForwardTransform( const MultiFab& mf, const int field_index ); - void BackwardTransform( MultiFab& mf, const int field_index ); - - private: - 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; - SpectralField& getSpectralField( const int field_index ); -}; - -#endif // WARPX_PSATD_DATA_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralData.cpp b/Source/FieldSolver/SpectralSolver/SpectralData.cpp deleted file mode 100644 index b3c902b20..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralData.cpp +++ /dev/null @@ -1,184 +0,0 @@ -#include - -using namespace amrex; - -SpectralData::SpectralData( const BoxArray& realspace_ba, - const BoxArray& spectralspace_ba, - const DistributionMapping& dm ) -{ - // 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 - } -} - - -SpectralData::~SpectralData() -{ - 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 -SpectralData::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 -SpectralData::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& -SpectralData::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 - } -} diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H new file mode 100644 index 000000000..2f6274e40 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -0,0 +1,46 @@ +#ifndef WARPX_SPECTRAL_FIELD_DATA_H_ +#define WARPX_SPECTRAL_FIELD_DATA_H_ + +#include +#include +#include + +using namespace amrex; + +// Declare type for spectral fields +using SpectralField = FabArray>; + +/* 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 + */ +class SpectralFieldData +{ + friend class PsatdAlgorithm; + +#ifdef AMREX_USE_GPU +// Add cuFFT-specific code +#else + using FFTplans = LayoutData; +#endif + + public: + SpectralFieldData( const BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const DistributionMapping& dm ); + ~SpectralFieldData(); + void ForwardTransform( const MultiFab& mf, const int field_index ); + void BackwardTransform( MultiFab& mf, const int field_index ); + + private: + 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; + SpectralField& getSpectralField( const int field_index ); +}; + +#endif // WARPX_SPECTRAL_FIELD_DATA_H_ 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 f9e7f5168ceca79eb1c68905787d5985b5bbe427 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 06:16:03 -0700 Subject: Do not import namespace in .H file --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 11 ++++------- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 14 ++++++-------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 22 +++++++++++----------- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 3 +++ 4 files changed, 24 insertions(+), 26 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 7e94bdec6..9fbdc7073 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -4,20 +4,17 @@ #include #include -using namespace amrex; -using namespace Gpu; - /* TODO: Write documentation */ class PsatdAlgorithm { - using SpectralCoefficients = FabArray>; + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; public: PsatdAlgorithm(const SpectralKSpace& spectral_kspace, - const DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const Real dt); + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const amrex::Real dt); void pushSpectralFields( SpectralFieldData& f ) const; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 2f6274e40..9f5b85a2b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -5,10 +5,8 @@ #include #include -using namespace amrex; - // Declare type for spectral fields -using SpectralField = FabArray>; +using SpectralField = amrex::FabArray< amrex::BaseFab >; /* Fields that will be stored in spectral space */ struct SpectralFieldIndex{ @@ -25,16 +23,16 @@ class SpectralFieldData #ifdef AMREX_USE_GPU // Add cuFFT-specific code #else - using FFTplans = LayoutData; + using FFTplans = amrex::LayoutData; #endif public: - SpectralFieldData( const BoxArray& realspace_ba, + SpectralFieldData( const amrex::BoxArray& realspace_ba, const SpectralKSpace& k_space, - const DistributionMapping& dm ); + const amrex::DistributionMapping& dm ); ~SpectralFieldData(); - void ForwardTransform( const MultiFab& mf, const int field_index ); - void BackwardTransform( MultiFab& mf, const int field_index ); + void ForwardTransform( const amrex::MultiFab& mf, const int field_index ); + void BackwardTransform( amrex::MultiFab& mf, const int field_index ); private: SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 04264e629..f61cffe14 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -4,10 +4,7 @@ #include #include -using namespace amrex; -using namespace Gpu; - -using SpectralKVector = LayoutData>; +using SpectralKVector = amrex::LayoutData< amrex::Gpu::ManagedVector >; /* TODO Documentation: class represent k space, and how it is distribution * for local FFT on each MPI: k spaces are not connected. @@ -15,18 +12,21 @@ using SpectralKVector = LayoutData>; class SpectralKSpace { public: - SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, const Real* dx ); - BoxArray spectralspace_ba; + SpectralKSpace( const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, const amrex::Real* dx ); + amrex::BoxArray spectralspace_ba; SpectralKVector kx_vec, ky_vec, kz_vec; - const Real* dx; + const amrex::Real* dx; }; void -AllocateAndFillKvector( ManagedVector& k, const Box& bx, const Real* dx, const int i_dim ); +AllocateAndFillKvector( amrex::Gpu::ManagedVector& k, + const amrex::Box& bx, const amrex::Real* dx, const int i_dim ); void -ComputeModifiedKVector( ManagedVector& modified_k, - const ManagedVector& k, - const Box& bx, const Real dx, const int norder ); +ComputeModifiedKVector( amrex::Gpu::ManagedVector& modified_k, + const amrex::Gpu::ManagedVector& k, + const amrex::Box& bx, const amrex::Real dx, + const int norder ); #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 03d62892f..0bd57c7ea 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -1,6 +1,9 @@ #include #include +using namespace amrex; +using namespace Gpu; + SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, const Real* realspace_dx ) -- cgit v1.2.3 From 551e934fdee50f2321076b0dd1882a74cc92fb30 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 10:43:18 -0700 Subject: Link spectral solver to the rest of the code --- Docs/source/running_cpp/parameters.rst | 10 ++++-- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 10 +++--- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 2 ++ Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 7 ++-- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 5 +-- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 40 +++++++++++----------- Source/FieldSolver/WarpXFFT.cpp | 9 ++++- Source/WarpX.H | 7 +++- Source/WarpX.cpp | 28 ++++++++------- 9 files changed, 72 insertions(+), 46 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 4b1496ed4..f8f3dcde2 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -108,7 +108,7 @@ user-defined constant and ``x`` and ``y`` are variables. The names are case sens define functions by intervals. User-defined constants can be used in parsed functions only (i.e., ``density_function(x,y,z)`` and ``field_function(X,Y,t)``, see below). User-defined constants can contain only letter, numbers and character _. -The name of each constant has to begin with a letter. The following names are used +The name of each constant has to begin with a letter. The following names are used by WarpX, and cannot be used as user-defined constants: `x`, `y`, `z`, `X`, `Y`, `t`. For example, parameters ``a0`` and ``z_plateau`` can be specified with: @@ -428,8 +428,14 @@ Numerics and algorithms * ``psatd.nox``, ``psatd.noy``, ``pstad.noz`` (`integer`) optional (default `16` for all) The order of accuracy of the spatial derivatives, when using the code compiled with a PSATD solver. +* ``psatd.hybrid_mpi_decomposition`` (`0` or `1`; default: 0) + Whether to use a different MPI decomposition for the particle-grid operations + (deposition and gather) and for the PSATD solver. If `1`, the FFT will + be performed over MPI groups. + * ``psatd.ngroups_fft`` (`integer`) - The number of MPI groups that are created for the FFT, when using the code compiled with a PSATD solver. + The number of MPI groups that are created for the FFT, when using the code compiled with a PSATD solver + (and only if `hybrid_mpi_decomposition` is `1`). The FFTs are global within one MPI group and use guard cell exchanges in between MPI groups. (If ``ngroups_fft`` is larger than the number of MPI ranks used, than the actual number of MPI ranks is used instead.) diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 9fbdc7073..375fa48bb 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -11,10 +11,12 @@ class PsatdAlgorithm using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; public: - PsatdAlgorithm(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const amrex::Real dt); + PsatdAlgorithm( const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const amrex::Real dt); + PsatdAlgorithm() = default; // Default constructor + PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; // Default move assignment void pushSpectralFields( SpectralFieldData& f ) const; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 9f5b85a2b..a74bfe22b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -30,6 +30,8 @@ class SpectralFieldData SpectralFieldData( const amrex::BoxArray& realspace_ba, const SpectralKSpace& k_space, const amrex::DistributionMapping& dm ); + 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 ); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index f61cffe14..a790be94f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -13,15 +13,16 @@ class SpectralKSpace { public: SpectralKSpace( const amrex::BoxArray& realspace_ba, - const amrex::DistributionMapping& dm, const amrex::Real* dx ); + const amrex::DistributionMapping& dm, + const amrex::Array dx ); amrex::BoxArray spectralspace_ba; SpectralKVector kx_vec, ky_vec, kz_vec; - const amrex::Real* dx; + amrex::Array dx; }; void AllocateAndFillKvector( amrex::Gpu::ManagedVector& k, - const amrex::Box& bx, const amrex::Real* dx, const int i_dim ); + const amrex::Box& bx, const amrex::Array 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 0bd57c7ea..2b1e7ee33 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 Real* realspace_dx ) + const Array realspace_dx ) { // Create the box array that corresponds to spectral space BoxList spectral_bl; // Create empty box list @@ -37,7 +37,8 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } void -AllocateAndFillKvector( ManagedVector& k, const Box& bx, const Real* dx, const int i_dim ) +AllocateAndFillKvector( ManagedVector& k, const Box& bx, + const Array 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 cde1fccd2..521e558ba 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -1,24 +1,39 @@ #ifndef WARPX_SPECTRAL_SOLVER_H_ #define WARPX_SPECTRAL_SOLVER_H_ +#include +#include +#include + /* \brief * TODO */ class SpectralSolver { private: - SpectralKSpace k_space; // Contains size of each box in spectral space, - // and corresponding values of the k vectors SpectralFieldData field_data; // Store field in spectral space // and perform the Fourier transforms PsatdAlgorithm algorithm; // Contains Psatd coefficients // and field update equation public: - SpectralSolver(); - void ForwardTransform( const MultiFab& mf, const int field_index ){ + 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 ) + { + // Initialize all structures using the same distribution mapping dm + // - Initialize k space (Contains size of each box in spectral space, + // and corresponding values of the k vectors) + const SpectralKSpace k_space = SpectralKSpace( realspace_ba, dm, dx ); + // - Initialize algorithm (coefficients) over this space + algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, dt ); + // - 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 BackwardTransform( MultiFab& mf, const int field_index ){ + void BackwardTransform( amrex::MultiFab& mf, const int field_index ){ field_data.BackwardTransform( mf, field_index ); }; void pushSpectralFields(){ @@ -26,19 +41,4 @@ class SpectralSolver }; }; -SpectralSolver::SpectralSolver( const SpectralKSpace& realspace_ba, - const DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const Real* dx, - const Real dt ) -{ - // Initialize all structures using the same distribution mapping dm - // - Initialize k space values - k_space = SpectralKSpace( realspace_ba, dm, dx ); - // - Initialize algorithm (coefficients) over this space - algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z ); - // - Initialize arrays for fields in Fourier space + FFT plans - field_data = SpectralFieldData( realspace_ba, k_space, dm ); -}; - #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index f09290f29..1b1da4d6f 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -126,6 +126,13 @@ WarpX::AllocLevelDataFFT (int lev) FFTDomainDecomposition(lev, ba_fp_fft, dm_fp_fft, ba_valid_fp_fft[lev], domain_fp_fft[lev], geom[lev].Domain()); + if (fft_hybrid_mpi_decomposition == false){ + // Allocate and initialize objects for the spectral solver + // (all use the same distribution mapping) + spectral_solver_fp[lev].reset( new SpectralSolver( ba_fp_fft, dm_fp_fft, + nox_fft, noy_fft, noz_fft, CellSize(lev), dt[lev] ) ); + } + // rho2 has one extra ghost cell, so that it's safe to deposit charge density after // pushing particle. @@ -426,7 +433,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) WARPX_TO_FORTRAN_ANYD(fab), WARPX_TO_FORTRAN_ANYD(fab), WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab)); + WARPX_TO_FORTRAN_ANYD(fab)); } else // Multiple FFT patches on this MPI rank diff --git a/Source/WarpX.H b/Source/WarpX.H index cc01a042e..d21e36b40 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -26,6 +26,7 @@ #ifdef WARPX_USE_PSATD #include +#include #endif #if defined(BL_USE_SENSEI_INSITU) @@ -538,7 +539,7 @@ private: amrex::Vector particle_plot_vars; amrex::Vector particle_plot_flags; - + #ifdef WARPX_USE_PSATD // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) // This includes data for the FFT guard cells (between FFT groups) @@ -586,6 +587,9 @@ private: amrex::Vector > > dataptr_fp_fft; amrex::Vector > > dataptr_cp_fft; + amrex::Vector> spectral_solver_fp; + amrex::Vector> spectral_solver_cp; + // Domain decomposition containing the valid part of the dual grids (i.e. without FFT guard cells) amrex::Vector ba_valid_fp_fft; amrex::Vector ba_valid_cp_fft; @@ -596,6 +600,7 @@ private: amrex::Vector comm_fft; amrex::Vector color_fft; + bool fft_hybrid_mpi_decomposition = false; int ngroups_fft = 4; int fftw_plan_measure = 1; int nox_fft = 16; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 51d9133be..cb8803c1b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -261,12 +261,12 @@ WarpX::ReadParameters () pp.query("verbose", verbose); pp.query("regrid_int", regrid_int); pp.query("do_subcycling", do_subcycling); - + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1, "Subcycling method 1 only works for 2 levels."); - + ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); - + pp.queryarr("B_external", B_external); pp.query("do_moving_window", do_moving_window); @@ -313,7 +313,7 @@ WarpX::ReadParameters () pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); if (do_boosted_frame_diagnostic) { - + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boosted frame diagnostic."); @@ -328,8 +328,8 @@ WarpX::ReadParameters () pp.query("do_boosted_frame_fields", do_boosted_frame_fields); pp.query("do_boosted_frame_particles", do_boosted_frame_particles); - - + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); @@ -352,7 +352,7 @@ WarpX::ReadParameters () filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; #if (AMREX_SPACEDIM == 3) filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; -#endif +#endif pp.query("serialize_ics", serialize_ics); pp.query("refine_plasma", refine_plasma); @@ -389,7 +389,7 @@ WarpX::ReadParameters () } } } - + if (plot_F){ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_dive_cleaning, "plot_F only works if warpx.do_dive_cleaning = 1"); @@ -432,7 +432,7 @@ WarpX::ReadParameters () // select which particle comps to write { pp.queryarr("particle_plot_vars", particle_plot_vars); - + if (particle_plot_vars.size() == 0) { particle_plot_flags.resize(PIdx::nattribs, 1); @@ -440,18 +440,18 @@ WarpX::ReadParameters () else { particle_plot_flags.resize(PIdx::nattribs, 0); - + for (const auto& var : particle_plot_vars) { particle_plot_flags[ParticleStringNames::to_index.at(var)] = 1; } } } - + pp.query("load_balance_int", load_balance_int); pp.query("load_balance_with_sfc", load_balance_with_sfc); pp.query("load_balance_knapsack_factor", load_balance_knapsack_factor); - + pp.query("do_dynamic_scheduling", do_dynamic_scheduling); pp.query("do_nodal", do_nodal); @@ -509,11 +509,14 @@ WarpX::ReadParameters () #ifdef WARPX_USE_PSATD { ParmParse pp("psatd"); + pp.query("hybrid_mpi_decomposition", fft_hybrid_mpi_decomposition); pp.query("ngroups_fft", ngroups_fft); pp.query("fftw_plan_measure", fftw_plan_measure); pp.query("nox", nox_fft); pp.query("noy", noy_fft); pp.query("noz", noz_fft); + // Override value + if (fft_hybrid_mpi_decomposition==false) ngroups_fft=ParallelDescriptor::NProcs(); } #endif @@ -1075,4 +1078,3 @@ WarpX::PicsarVersion () return std::string("Unknown"); #endif } - -- 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/PsatdAlgorithm.H') 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 7b0571362a6eb23427ea56b1d780378c6f7730d8 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 20 Apr 2019 21:36:19 -0700 Subject: Change type name --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 6 +++--- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 8 ++++---- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 8 ++++---- 3 files changed, 11 insertions(+), 11 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index b54f5c143..2c3d8abfd 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -21,10 +21,10 @@ class PsatdAlgorithm private: // Modified finite-order vectors - SpectralKVector modified_kx_vec, modified_kz_vec; + KVectorComponent modified_kx_vec, modified_kz_vec; #if (AMREX_SPACEDIM==3) - SpectralKVector modified_ky_vec; -#endif + KVectorComponent modified_ky_vec; +#endif SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 6589c6346..366bc2a79 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -4,7 +4,7 @@ #include #include -using SpectralKVector = amrex::LayoutData< amrex::Gpu::ManagedVector >; +using KVectorComponent = amrex::LayoutData< amrex::Gpu::ManagedVector >; /* TODO Documentation: class represent k space, and how it is distribution * for local FFT on each MPI: k spaces are not connected. @@ -15,14 +15,14 @@ class SpectralKSpace SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const amrex::RealVect dx ); - SpectralKVector& AllocateAndFillKVector( + KVectorComponent& AllocateAndFillKVector( const amrex::DistributionMapping& dm, const int i_dim ) const; - SpectralKVector& AllocateAndFillModifiedKVector( + KVectorComponent& AllocateAndFillModifiedKVector( const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; amrex::BoxArray spectralspace_ba; private: - amrex::Array k_vec; + amrex::Array k_vec; // 3D: array of 3 components, corresponding to kx, ky, kz // 2D: array of 2 components, corresponding to kx, kz amrex::RealVect dx; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 3c600221e..6e6115cec 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -29,11 +29,11 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } } -SpectralKVector& +KVectorComponent& SpectralKSpace::AllocateAndFillKVector( const DistributionMapping& dm, const int i_dim ) const { // Initialize an empty vector in each box - SpectralKVector k_vec = SpectralKVector(spectralspace_ba, dm); + KVectorComponent k_vec = KVectorComponent(spectralspace_ba, dm); // Loop over boxes for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; @@ -65,12 +65,12 @@ SpectralKSpace::AllocateAndFillKVector( const DistributionMapping& dm, const int return k_vec; } -SpectralKVector& +KVectorComponent& SpectralKSpace::AllocateAndFillModifiedKVector( const DistributionMapping& dm, const int i_dim, const int order ) const { // Initialize an empty vector in each box - SpectralKVector modified_k_vec = SpectralKVector( spectralspace_ba, dm ); + KVectorComponent modified_k_vec = KVectorComponent( spectralspace_ba, dm ); // Loop over boxes for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ const ManagedVector& k = k_vec[i_dim][mfi]; -- cgit v1.2.3 From eb5ee68612afe016afd47a621937ee57006c135c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 22 Apr 2019 11:28:08 -0700 Subject: Use Fonberg coefficients to calculate the modified k --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 6 +- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 64 +++++++++++----------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 6 +- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 21 +++++-- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 4 +- Source/FieldSolver/WarpXFFT.cpp | 2 +- 6 files changed, 59 insertions(+), 44 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 2c3d8abfd..16d27cab2 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -11,13 +11,13 @@ class PsatdAlgorithm using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; public: - PsatdAlgorithm( const SpectralKSpace& spectral_kspace, + PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, - const int norder_z, const amrex::Real dt); + const int norder_z, const bool nodal, const amrex::Real dt); PsatdAlgorithm() = default; // Default constructor PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; // Default move assignment - void pushSpectralFields( SpectralFieldData& f ) const; + void pushSpectralFields(SpectralFieldData& f) const; private: // Modified finite-order vectors diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 37d3b33ee..5ebb7144d 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -7,29 +7,29 @@ using namespace amrex; PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, - const int norder_z, const Real dt) + const int norder_z, const bool nodal, const Real dt) { const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.getModifiedKComponent( dm, 0, norder_x ); + modified_kx_vec = spectral_kspace.getModifiedKComponent(dm, 0, norder_x, nodal); #if (AMREX_SPACEDIM==3) - modified_ky_vec = spectral_kspace.getModifiedKComponent( dm, 1, norder_y ); - modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 2, norder_z ); + modified_ky_vec = spectral_kspace.getModifiedKComponent(dm, 1, norder_y, nodal); + modified_kz_vec = spectral_kspace.getModifiedKComponent(dm, 2, norder_z, nodal); #else - modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 1, norder_z ); + modified_kz_vec = spectral_kspace.getModifiedKComponent(dm, 1, norder_z, nodal); #endif // Allocate the arrays of coefficients - C_coef = SpectralCoefficients( ba, dm, 1, 0 ); - S_ck_coef = SpectralCoefficients( ba, dm, 1, 0 ); - X1_coef = SpectralCoefficients( ba, dm, 1, 0 ); - X2_coef = SpectralCoefficients( ba, dm, 1, 0 ); - X3_coef = SpectralCoefficients( ba, dm, 1, 0 ); + C_coef = SpectralCoefficients(ba, dm, 1, 0); + S_ck_coef = SpectralCoefficients(ba, dm, 1, 0); + X1_coef = SpectralCoefficients(ba, dm, 1, 0); + X2_coef = SpectralCoefficients(ba, dm, 1, 0); + X3_coef = SpectralCoefficients(ba, dm, 1, 0); // Fill them with the right values: // Loop over boxes - for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ const Box& bx = ba[mfi]; @@ -47,26 +47,26 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, Array4 X3 = X3_coef[mfi].array(); // Loop over indices within one box - ParallelFor( bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Calculate norm of vector const Real k_norm = std::sqrt( - std::pow( modified_kx[i], 2 ) + + std::pow(modified_kx[i], 2) + #if (AMREX_SPACEDIM==3) - std::pow( modified_ky[j], 2 ) + + std::pow(modified_ky[j], 2) + #endif - std::pow( modified_kz[k], 2 ) ); + std::pow(modified_kz[k], 2)); // Calculate coefficients constexpr Real c = PhysConst::c; constexpr Real ep0 = PhysConst::ep0; - if ( k_norm != 0 ){ - C(i,j,k) = std::cos( c*k_norm*dt ); - S_ck(i,j,k) = std::sin( c*k_norm*dt )/( c*k_norm ); + if (k_norm != 0){ + C(i,j,k) = std::cos(c*k_norm*dt); + S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); X1(i,j,k) = (1. - C(i,j,k))/(ep0 * c*c * k_norm*k_norm); - X2(i,j,k) = (1. - S_ck(i,j,k)/dt )/(ep0 * k_norm*k_norm); - X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt )/(ep0 * k_norm*k_norm); + X2(i,j,k) = (1. - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); + X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); } else { // Handle k_norm = 0, by using the analytical limit C(i,j,k) = 1.; S_ck(i,j,k) = dt; @@ -79,10 +79,10 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, }; void -PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ +PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ // Loop over boxes - for ( MFIter mfi(f.Ex); mfi.isValid(); ++mfi ){ + for (MFIter mfi(f.Ex); mfi.isValid(); ++mfi){ const Box& bx = f.Ex[mfi].box(); @@ -113,7 +113,7 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); // Loop over indices within one box - ParallelFor( bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Record old values of the fields to be updated @@ -148,24 +148,24 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ // Update E (see WarpX online documentation: theory section) Ex_arr(i,j,k) = C*Ex_old - + S_ck*( c2*I*(ky*Bz_old - kz*By_old) - inv_ep0*Jx ) - - I*( X2*rho_new - X3*rho_old )*kx; + + S_ck*(c2*I*(ky*Bz_old - kz*By_old) - inv_ep0*Jx) + - I*(X2*rho_new - X3*rho_old)*kx; Ey_arr(i,j,k) = C*Ey_old - + S_ck*( c2*I*(kz*Bx_old - kx*Bz_old) - inv_ep0*Jy ) - - I*( X2*rho_new - X3*rho_old )*ky; + + S_ck*(c2*I*(kz*Bx_old - kx*Bz_old) - inv_ep0*Jy) + - I*(X2*rho_new - X3*rho_old)*ky; Ez_arr(i,j,k) = C*Ez_old - + S_ck*( c2*I*(kx*By_old - ky*Bx_old) - inv_ep0*Jz ) - - I*( X2*rho_new - X3*rho_old )*kz; + + S_ck*(c2*I*(kx*By_old - ky*Bx_old) - inv_ep0*Jz) + - I*(X2*rho_new - X3*rho_old)*kz; // Update B (see WarpX online documentation: theory section) Bx_arr(i,j,k) = C*Bx_old - S_ck*I*(ky*Ez_old - kz*Ey_old) - + X1*I*(ky*Jz - kz*Jy ); + + X1*I*(ky*Jz - kz*Jy); By_arr(i,j,k) = C*By_old - S_ck*I*(kz*Ex_old - kx*Ez_old) - + X1*I*(kz*Jx - kx*Jz ); + + X1*I*(kz*Jx - kx*Jz); Bz_arr(i,j,k) = C*Bz_old - S_ck*I*(kx*Ey_old - ky*Ex_old) - + X1*I*(kx*Jy - ky*Jx ); + + X1*I*(kx*Jy - ky*Jx); }); } }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 2f0681c3b..fd77bb7b8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -25,7 +25,8 @@ class SpectralKSpace KVectorComponent getKComponent( const amrex::DistributionMapping& dm, const int i_dim ) const; KVectorComponent getModifiedKComponent( - const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; + 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; @@ -36,4 +37,7 @@ class SpectralKSpace amrex::RealVect dx; }; +amrex::Vector +getFonbergStencilCoefficients( const int n_order, const bool nodal ); + #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 77aa7342f..111643197 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -101,8 +101,13 @@ SpectralKSpace::getModifiedKComponent( { // Initialize an empty ManagedVector in each box KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); + + // Compute real-space stencil coefficients + Vector stencil_coef = getFonbergStencilCoefficients(n_order, nodal); + // Loop over boxes for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + Real delta_x = dx[i_dim]; const ManagedVector& k = k_vec[i_dim][mfi]; ManagedVector& modified_k = modified_k_comp[mfi]; @@ -111,9 +116,15 @@ SpectralKSpace::getModifiedKComponent( // Fill the modified k vector for (int i=0; i +Vector getFonbergStencilCoefficients( const int n_order, const bool nodal ) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, "n_order should be even."); const int m = n_order/2; - Array stencil_coef; + Vector stencil_coef; stencil_coef.resize( m+1 ); // Coefficients for nodal (a.k.a. centered) finite-difference diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 06d53a2a9..363ac7bd8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -19,14 +19,14 @@ 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::RealVect dx, const amrex::Real dt ) + const bool nodal, 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, // and corresponding values of the k vectors) const SpectralKSpace k_space = SpectralKSpace( realspace_ba, dm, dx ); // - Initialize algorithm (coefficients) over this space - algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, dt ); + algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, nodal, dt ); // - Initialize arrays for fields in Fourier space + FFT plans field_data = SpectralFieldData( realspace_ba, k_space, dm ); }; diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index e655bd11f..0997643ab 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -132,7 +132,7 @@ WarpX::AllocLevelDataFFT (int lev) 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, dx_vect, dt[lev] ) ); + nox_fft, noy_fft, noz_fft, do_nodal, 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 4a3265daf691f0a8fa7461322c7299016fd706bf Mon Sep 17 00:00:00 2001 From: Remi Lehe 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/PsatdAlgorithm.H') 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