From 949f25f10118d8b66fe827706904c49ab42178c8 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 18 Apr 2019 15:37:11 -0700 Subject: Add Spectral K space --- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 62 ++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp new file mode 100644 index 000000000..ab684444d --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -0,0 +1,62 @@ +#include +#include + +SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, + const DistributionMapping& dm, + const Real* realspace_dx ) +{ + // Create the box array that corresponds to spectral space + 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) + Box realspace_bx = realspace_ba[i]; + Box bx = Box( IntVect::TheZeroVector(), realspace_bx.bigEnd() - realspace_bx.smallEnd() ); + spectral_bl.push_back( bx ); + } + spectralspace_ba.define( spectral_bl ); + + // Allocate the 1D vectors + kx_vec = SpectralKVector( spectralspace_ba, dm ); + ky_vec = SpectralKVector( spectralspace_ba, dm ); + 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 ); + AllocateAndFillKvector( ky_vec[mfi], bx, dx, 1 ); + AllocateAndFillKvector( kz_vec[mfi], bx, dx, 2 ); + } + + // Store the cell size + dx = realspace_dx; +} + +void +AllocateAndFillKvector( ManagedVector& k, const Box& bx, const Real* dx, const int i_dim ) +{ + // Alllocate k to the right size + int N = bx.length( i_dim ); + k.resize( N ); + + // Fill the k vector + const Real dk = 2*MathConst::pi/(N*dx[i_dim]); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.smallEnd(i_dim) == 0, + "Expected box to start at 0, in spectral space."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.bigEnd(i_dim) == N-1, + "Expected different box end index in spectral space."); + // Fill positive values of k (FFT conventions: first half is positive) + for (int i=0; i<(N+1)/2; i++ ){ + k[i] = i*dk; + } + // Fill negative values of k (FFT conventions: second half is negative) + for (int i=(N+1)/2; i Date: Thu, 18 Apr 2019 16:31:30 -0700 Subject: Added modified k vector --- Source/FieldSolver/SpectralSolver/PsatdSolver.H | 6 ++-- Source/FieldSolver/SpectralSolver/PsatdSolver.cpp | 35 +++++++++++++++------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 8 +++-- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 16 ++++++++++ 4 files changed, 49 insertions(+), 16 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdSolver.H b/Source/FieldSolver/SpectralSolver/PsatdSolver.H index bb46da670..06a355021 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdSolver.H +++ b/Source/FieldSolver/SpectralSolver/PsatdSolver.H @@ -14,8 +14,10 @@ class PsatdSolver using SpectralCoefficients = FabArray>; public: - PsatdSolver( const SpectralKSpace& spectral_kspace, - const DistributionMapping& dm, Real dt ); + PsatdSolver(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: diff --git a/Source/FieldSolver/SpectralSolver/PsatdSolver.cpp b/Source/FieldSolver/SpectralSolver/PsatdSolver.cpp index 66409ca33..f8a3f7c6e 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdSolver.cpp @@ -4,20 +4,33 @@ using namespace amrex; -/* - * ba: BoxArray for spectral space - * dm: DistributionMapping for spectral space - */ + PsatdSolver::PsatdSolver(const SpectralKSpace& spectral_kspace, - const DistributionMapping& dm, const Real dt) + const DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const Real dt) { + const BoxArray& ba = spectral_kspace.spectralspace_ba; + // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.getModifiedKVector( 0 ); - modified_ky_vec = spectral_kspace.getModifiedKVector( 1 ); - modified_kz_vec = spectral_kspace.getModifiedKVector( 2 ); + modified_kx_vec = SpectralKVector( ba, dm ); + modified_ky_vec = SpectralKVector( ba, dm ); + modified_kz_vec = SpectralKVector( ba, dm ); + // Allocate and fill them by computing the modified vector + for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ + Box bx = ba[mfi]; + ComputeModifiedKVector( + modified_kx_vec[mfi], spectral_kspace.kx_vec[mfi], + bx, spectral_kspace.dx[0], norder_x ); + 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 ); + } // Allocate the arrays of coefficients - const BoxArray& ba = spectral_kspace.spectralspace_ba; C_coef = SpectralCoefficients( ba, dm, 1, 0 ); S_ck_coef = SpectralCoefficients( ba, dm, 1, 0 ); X1_coef = SpectralCoefficients( ba, dm, 1, 0 ); @@ -69,7 +82,7 @@ PsatdSolver::PsatdSolver(const SpectralKSpace& spectral_kspace, } }); } -} +}; void PsatdSolver::pushSpectralFields( SpectralData& f ) const{ @@ -155,4 +168,4 @@ PsatdSolver::pushSpectralFields( SpectralData& f ) const{ + X1*I*(kx*Jy - ky*Jx ); }); } -} +}; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index d1aea836e..04264e629 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -16,10 +16,7 @@ class SpectralKSpace { public: SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, const Real* dx ); - SpectralKVector& getModifiedKVector( const int i_dim ) const; BoxArray spectralspace_ba; - - private: SpectralKVector kx_vec, ky_vec, kz_vec; const Real* dx; }; @@ -27,4 +24,9 @@ class SpectralKSpace void AllocateAndFillKvector( ManagedVector& k, const Box& bx, const Real* dx, const int i_dim ); +void +ComputeModifiedKVector( ManagedVector& modified_k, + const ManagedVector& k, + const Box& bx, const Real dx, const int norder ); + #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index ab684444d..03d62892f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -60,3 +60,19 @@ AllocateAndFillKvector( ManagedVector& k, const Box& bx, const Real* dx, c // TODO: For real-to-complex, } + +void +ComputeModifiedKVector( ManagedVector& modified_k, + const ManagedVector& k, + const Box& bx, const Real dx, const int norder ) +{ + // Allocate modified_k to the right size + int N = k.size(); + modified_k.resize( N ); + + // For now, this simply copies the infinite order k + for (int i=0; i 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/SpectralKSpace.cpp') 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/SpectralKSpace.cpp') 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 e54e25a2cf23007bc8132b87079664fcf891fdfb Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 20:50:38 -0700 Subject: Fix NaN issue --- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 2b1e7ee33..d05748192 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -8,7 +8,10 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, const Array realspace_dx ) { - // Create the box array that corresponds to spectral space + // Store the cell size + dx = realspace_dx; + + // Create the box array that corresponds to spectral space BoxList spectral_bl; // Create empty box list // Loop over boxes and fill the box list for (int i=0; i < realspace_ba.size(); i++ ) { @@ -31,9 +34,6 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, AllocateAndFillKvector( ky_vec[mfi], bx, dx, 1 ); AllocateAndFillKvector( kz_vec[mfi], bx, dx, 2 ); } - - // Store the cell size - dx = realspace_dx; } void -- 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/SpectralKSpace.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 7104eee32f6f909d8a0b3abd11528b9e059a36d4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 20 Apr 2019 21:29:07 -0700 Subject: Change function interface --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 25 +---- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 23 ++--- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 106 +++++++++++---------- 3 files changed, 67 insertions(+), 87 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 75df01343..0b1dbc4da 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -12,30 +12,13 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = SpectralKVector( ba, dm ); + modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 0, norder_x ); #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 ){ - Box bx = ba[mfi]; - 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 ); + modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_y ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 2, norder_z ); #else - ComputeModifiedKVector( - modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], - bx, spectral_kspace.dx[1], norder_z ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_z ); #endif - } // Allocate the arrays of coefficients C_coef = SpectralCoefficients( ba, dm, 1, 0 ); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 16ec95cf5..6589c6346 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -15,22 +15,17 @@ class SpectralKSpace SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const amrex::RealVect dx ); + SpectralKVector& AllocateAndFillKVector( + const amrex::DistributionMapping& dm, const int i_dim ) const; + SpectralKVector& AllocateAndFillModifiedKVector( + const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; amrex::BoxArray spectralspace_ba; - SpectralKVector kx_vec, kz_vec; -#if (AMREX_SPACEDIM==3) - SpectralKVector ky_vec; -#endif + + private: + 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; }; -void -AllocateAndFillKvector( amrex::Gpu::ManagedVector& k, - const amrex::Box& bx, const amrex::RealVect dx, const int i_dim ); - -void -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 76e8aef15..3c600221e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -11,7 +11,7 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Store the cell size dx = realspace_dx; - // Create the box array that corresponds to spectral space + // Create the box array that corresponds to spectral space BoxList spectral_bl; // Create empty box list // Loop over boxes and fill the box list for (int i=0; i < realspace_ba.size(); i++ ) { @@ -23,66 +23,68 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } spectralspace_ba.define( spectral_bl ); - // 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 + // Allocate the components of the k vector: kx, ky (only in 3D), kz + for (int i_dim=0; i_dim& k, const Box& bx, - const RealVect dx, const int i_dim ) +SpectralKVector& +SpectralKSpace::AllocateAndFillKVector( const DistributionMapping& dm, const int i_dim ) const { - // Alllocate k to the right size - int N = bx.length( i_dim ); - k.resize( N ); + // Initialize an empty vector in each box + SpectralKVector k_vec = SpectralKVector(spectralspace_ba, dm); + // Loop over boxes + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + Box bx = spectralspace_ba[mfi]; + ManagedVector& k = k_vec[mfi]; - // Fill the k vector - const Real dk = 2*MathConst::pi/(N*dx[i_dim]); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.smallEnd(i_dim) == 0, - "Expected box to start at 0, in spectral space."); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.bigEnd(i_dim) == N-1, - "Expected different box end index in spectral space."); - // Fill positive values of k (FFT conventions: first half is positive) - for (int i=0; i<(N+1)/2; i++ ){ - k[i] = i*dk; - } - // Fill negative values of k (FFT conventions: second half is negative) - for (int i=(N+1)/2; i& modified_k, - const ManagedVector& k, - const Box& bx, const Real dx, const int norder ) +SpectralKVector& +SpectralKSpace::AllocateAndFillModifiedKVector( + const DistributionMapping& dm, const int i_dim, const int order ) const { - // Allocate modified_k to the right size - int N = k.size(); - modified_k.resize( N ); + // Initialize an empty vector in each box + SpectralKVector modified_k_vec = SpectralKVector( 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]; - // For now, this simply copies the infinite order k - for (int i=0; i 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/SpectralKSpace.cpp') 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 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/SpectralKSpace.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/SpectralKSpace.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/SpectralKSpace.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: Mon, 22 Apr 2019 09:56:04 -0700 Subject: Add calculation of Fonberg coefficients --- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 76 +++++++++++++++------- 1 file changed, 54 insertions(+), 22 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 90c20b01d..77aa7342f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -66,9 +66,38 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) return k_comp; } +SpectralShiftFactor +SpectralKSpace::getSpectralShiftFactor( + const DistributionMapping& dm, const int i_dim, const int shift_type ) const +{ + // Initialize an empty ManagedVector in each box + SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); + // Loop over boxes + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + const ManagedVector& 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.; break; + case ShiftType::NodalToCentered: sign = 1.; + } + constexpr Complex I{0,1}; + for (int i=0; i +getFonbergStencilCoefficients( const int n_order, const bool nodal ) { - // Initialize an empty ManagedVector in each box - SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); - // Loop over boxes - for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ - const ManagedVector& k = k_vec[i_dim][mfi]; - ManagedVector& shift = shift_factor[mfi]; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, "n_order should be even."); + const int m = n_order/2; + Array stencil_coef; + stencil_coef.resize( m+1 ); - // Allocate shift coefficients - shift.resize( k.size() ); - - // Fill the shift coefficients - Real sign = 0; - switch (shift_type){ - case ShiftType::CenteredToNodal: sign = -1.; break; - case ShiftType::NodalToCentered: sign = 1.; + // Coefficients for nodal (a.k.a. centered) finite-difference + if (nodal == true) { + stencil_coef[0] = -2.; // First coefficient + for (int n=1; n 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/SpectralKSpace.cpp') 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 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/SpectralKSpace.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 14:25:34 -0700 Subject: Add comments, switch value of k --- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 22 ++++++++++---- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 34 +++++++++++++++++----- 2 files changed, 43 insertions(+), 13 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index c71cf2643..eff42e32c 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -5,15 +5,25 @@ #include #include -using KVectorComponent = amrex::LayoutData>; -using SpectralShiftFactor = amrex::LayoutData>; +// `KVectorComponent` and `SpectralShiftFactor` hold one 1D array +// ("ManagedVector") for each box ("LayoutData"). The arrays are +// only allocated if the corresponding box is owned by the local MPI rank. +using KVectorComponent = amrex::LayoutData< + amrex::Gpu::ManagedVector >; +using SpectralShiftFactor = amrex::LayoutData< + amrex::Gpu::ManagedVector >; +// Indicate the type of correction "shift" factor to apply +// when the FFT is performed from/to a cell-centered grid in real space. struct ShiftType { enum{ TransformFromCellCentered=0, TransformToCellCentered=1 }; }; -/* TODO Documentation: class represent k space, and how it is distribution - * for local FFT on each MPI: k spaces are not connected. +/* \brief Class that represents the spectral space. + * + * (Contains info about the size of the spectral space corresponding + * to each box in `realspace_ba`, as well as the value of the + * corresponding k coordinates) */ class SpectralKSpace { @@ -33,8 +43,8 @@ class SpectralKSpace private: amrex::Array k_vec; - // 3D: array of 3 components, corresponding to kx, ky, kz - // 2D: array of 2 components, corresponding to kx, kz + // 3D: k_vec is an Array of 3 components, corresponding to kx, ky, kz + // 2D: k_vec is an 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 4852c7d2b..90866cd3b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -5,10 +5,21 @@ using namespace amrex; using namespace Gpu; +/* \brief Initialize k space object. + * + * \param realspace_ba Box array that corresponds to the decomposition + * of the fields in real space (cell-centered ; includes guard cells) + * \param dm Indicates which MPI proc owns which box, in realspace_ba. + * \param realspace_dx Cell size of the grid in real space + */ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, const RealVect realspace_dx ) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + realspace_ba.ixType()==IndexType::TheCellType(), + "SpectralKSpace expects a cell-centered box."); + // Store the cell size dx = realspace_dx; @@ -16,8 +27,10 @@ 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, boxes in spectral space start at 0 in each direction - // and have the same number of points as the real space box + // For local FFTs, boxes in spectral space start at 0 in + // each direction and have the same number of points as the + // (cell-centered) real space box + // TODO: this will be different for the real-to-complex FFT // TODO: this will be different for the hybrid FFT scheme Box realspace_bx = realspace_ba[i]; Box bx = Box( IntVect::TheZeroVector(), @@ -32,13 +45,19 @@ 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` + */ KVectorComponent SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) const { // Initialize an empty ManagedVector in each box KVectorComponent k_comp = KVectorComponent(spectralspace_ba, dm); - // Loop over boxes + // Loop over boxes and allocate the corresponding ManagedVector + // for each box owned by the local MPI proc ("mfi.isValid") for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; ManagedVector& k = k_comp[mfi]; @@ -53,15 +72,16 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, "Expected box to start at 0, in spectral space."); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.bigEnd(i_dim) == N-1, "Expected different box end index in spectral space."); + const int mid_point = (N+1)/2; // Fill positive values of k (FFT conventions: first half is positive) - for (int i=0; i<(N+1)/2; i++ ){ + for (int i=0; i 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/SpectralKSpace.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 119d84674eaa0ca410a66bbc2a47a05408361cf6 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 25 Apr 2019 14:58:56 -0700 Subject: Correct shift --- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index c82e45c5e..7bed96fc4 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -114,8 +114,8 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, // Fill the shift coefficients Real sign = 0; switch (shift_type){ - case ShiftType::TransformFromCellCentered: sign = 1.; break; - case ShiftType::TransformToCellCentered: sign = -1.; + case ShiftType::TransformFromCellCentered: sign = -1.; break; + case ShiftType::TransformToCellCentered: sign = 1.; } constexpr Complex I{0,1}; for (int i=0; i Date: Fri, 26 Apr 2019 19:42:55 -0700 Subject: Correct naming of `dx` ; use initializer list --- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 2 +- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index eff42e32c..ad17e6423 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -31,7 +31,7 @@ class SpectralKSpace amrex::BoxArray spectralspace_ba; SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, - const amrex::RealVect dx ); + const amrex::RealVect realspace_dx ); KVectorComponent getKComponent( const amrex::DistributionMapping& dm, const int i_dim ) const; KVectorComponent getModifiedKComponent( diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 7bed96fc4..6df3ea8fa 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -15,14 +15,12 @@ using namespace Gpu; SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, const RealVect realspace_dx ) + : dx(realspace_dx) // Store the cell size as member `dx` { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( realspace_ba.ixType()==IndexType::TheCellType(), "SpectralKSpace expects a cell-centered box."); - // Store the cell size - dx = realspace_dx; - // Create the box array that corresponds to spectral space BoxList spectral_bl; // Create empty box list // Loop over boxes and fill the box list @@ -33,6 +31,7 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // TODO: this will be different for the real-to-complex FFT // TODO: this will be different for the hybrid FFT scheme Box realspace_bx = realspace_ba[i]; + Print() << realspace_bx.smallEnd() << " " << realspace_bx.bigEnd() << std::endl; Box bx = Box( IntVect::TheZeroVector(), realspace_bx.bigEnd() - realspace_bx.smallEnd() ); spectral_bl.push_back( bx ); -- cgit v1.2.3 From 759f7e1bae3db1c2f436800189a3780a8447841e Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Apr 2019 19:48:42 -0700 Subject: Cosmetic changes --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 15 ++++++------ .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 6 ++--- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 11 +++++---- Source/FieldSolver/SpectralSolver/WarpX_Complex.H | 27 ---------------------- 4 files changed, 16 insertions(+), 43 deletions(-) delete mode 100644 Source/FieldSolver/SpectralSolver/WarpX_Complex.H (limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 60e9d58c0..900b542be 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -9,17 +9,16 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, 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, nodal); +// Compute and assign the modified k vectors +: 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, nodal); - modified_kz_vec = spectral_kspace.getModifiedKComponent(dm, 2, norder_z, nodal); + 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, nodal); + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, norder_z, nodal)) #endif +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the arrays of coefficients C_coef = SpectralCoefficients(ba, dm, 1, 0); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 6df3ea8fa..d91891a30 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -53,7 +53,7 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) const { // Initialize an empty ManagedVector in each box - KVectorComponent k_comp = KVectorComponent(spectralspace_ba, dm); + KVectorComponent k_comp(spectralspace_ba, dm); // 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 ){ @@ -100,7 +100,7 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, const int shift_type ) const { // Initialize an empty ManagedVector in each box - SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); + SpectralShiftFactor shift_factor( spectralspace_ba, dm ); // 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 ){ @@ -144,7 +144,7 @@ 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(spectralspace_ba, dm); // Compute real-space stencil coefficients Vector stencil_coef = getFonbergStencilCoefficients(n_order, nodal); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 4ed44cb9f..7444452af 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -13,11 +13,6 @@ */ class SpectralSolver { - private: - SpectralFieldData field_data; // Store field in spectral space - // and perform the Fourier transforms - PsatdAlgorithm algorithm; // Contains Psatd coefficients - // and field update equation public: // Inline definition of the member functions of `SpectralSolver` // The body of these functions is short, since the work is done in the @@ -66,6 +61,12 @@ class SpectralSolver BL_PROFILE("SpectralSolver::pushSpectralFields"); algorithm.pushSpectralFields( field_data ); }; + + private: + SpectralFieldData field_data; // Store field in spectral space + // and perform the Fourier transforms + PsatdAlgorithm algorithm; // Contains Psatd coefficients + // and field update equation }; #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/SpectralSolver/WarpX_Complex.H b/Source/FieldSolver/SpectralSolver/WarpX_Complex.H deleted file mode 100644 index c898c5baa..000000000 --- a/Source/FieldSolver/SpectralSolver/WarpX_Complex.H +++ /dev/null @@ -1,27 +0,0 @@ -#ifndef WARPX_COMPLEX_H_ -#define WARPX_COMPLEX_H_ - -#include - -// Define complex type on GPU/CPU -#ifdef AMREX_USE_GPU - -#include -#include -using Complex = thrust::complex; -static_assert( sizeof(Complex) == sizeof(cuDoubleComplex), - "The complex types in WarpX and cuFFT do not match."); - -#else - -#include -#include -using Complex = std::complex; -static_assert( sizeof(Complex) == sizeof(fftw_complex), - "The complex types in WarpX and FFTW do not match."); - -#endif -static_assert(sizeof(Complex) == sizeof(amrex::Real[2]), - "Unexpected complex type."); - -#endif //WARPX_COMPLEX_H_ -- cgit v1.2.3