diff options
Diffstat (limited to 'Source')
21 files changed, 389 insertions, 191 deletions
diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index ccbdd280c..186a8d45e 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -631,15 +631,16 @@ WarpX::WritePlotFile () const particle_varnames.push_back("theta"); #endif -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - particle_varnames.push_back("xold"); - particle_varnames.push_back("yold"); - particle_varnames.push_back("zold"); - - particle_varnames.push_back("uxold"); - particle_varnames.push_back("uyold"); - particle_varnames.push_back("uzold"); -#endif + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + particle_varnames.push_back("xold"); + particle_varnames.push_back("yold"); + particle_varnames.push_back("zold"); + + particle_varnames.push_back("uxold"); + particle_varnames.push_back("uyold"); + particle_varnames.push_back("uzold"); + } mypc->WritePlotFile(plotfilename, particle_plot_flags, particle_varnames); diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index 50127914d..b0ee54095 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,11 +1,12 @@ CEXE_headers += WarpX_Complex.H CEXE_headers += SpectralSolver.H +CEXE_sources += SpectralSolver.cpp CEXE_headers += SpectralFieldData.H CEXE_sources += SpectralFieldData.cpp -CEXE_headers += PsatdAlgorithm.H -CEXE_sources += PsatdAlgorithm.cpp CEXE_headers += SpectralKSpace.H CEXE_sources += SpectralKSpace.cpp +include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H deleted file mode 100644 index acefcc466..000000000 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef WARPX_PSATD_ALGORITHM_H_ -#define WARPX_PSATD_ALGORITHM_H_ - -#include <SpectralKSpace.H> -#include <SpectralFieldData.H> - -/* \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 <amrex::Real> >; - - public: - PsatdAlgorithm(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - 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; - void pushSpectralFields(SpectralFieldData& f) const; - - private: - // Modified finite-order vectors - KVectorComponent modified_kx_vec, modified_kz_vec; -#if (AMREX_SPACEDIM==3) - KVectorComponent modified_ky_vec; -#endif - SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; -}; - -#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package new file mode 100644 index 000000000..c62c21f44 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += SpectralBaseAlgorithm.H +CEXE_headers += PsatdAlgorithm.H +CEXE_sources += PsatdAlgorithm.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H new file mode 100644 index 000000000..0487e5226 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -0,0 +1,24 @@ +#ifndef WARPX_PSATD_ALGORITHM_H_ +#define WARPX_PSATD_ALGORITHM_H_ + +#include <SpectralBaseAlgorithm.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithm : public SpectralBaseAlgorithm +{ + public: + PsatdAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::Real dt); + // Redefine update equation from base class + virtual void pushSpectralFields(SpectralFieldData& f) const override final; + + private: + SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; +}; + +#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index ada7506c3..37892d35a 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -9,14 +9,9 @@ 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) -// 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)) -#else - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) -#endif + // Initialize members of base class + : SpectralBaseAlgorithm( spectral_kspace, dm, + norder_x, norder_y, norder_z, nodal ) { const BoxArray& ba = spectral_kspace.spectralspace_ba; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H new file mode 100644 index 000000000..602eb2473 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -0,0 +1,51 @@ +#ifndef WARPX_SPECTRAL_BASE_ALGORITHM_H_ +#define WARPX_SPECTRAL_BASE_ALGORITHM_H_ + +#include <SpectralKSpace.H> +#include <SpectralFieldData.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + * + * `SpectralBaseAlgorithm` is only a base class and cannot be used directly. + * Instead use its subclasses, which implement the specific field update + * equations for a given spectral algorithm. + */ +class SpectralBaseAlgorithm +{ + public: + // Member function that updates the fields in spectral space ; + // meant to be overridden in subclasses + virtual void pushSpectralFields(SpectralFieldData& f) const = 0; + // The destructor should also be a virtual function, so that + // a pointer to subclass of `SpectraBaseAlgorithm` actually + // calls the subclass's destructor. + virtual ~SpectralBaseAlgorithm() {}; + + protected: // Meant to be used in the subclasses + + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; + + // Constructor + SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool 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)) +#else + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) +#endif + {}; + + // Modified finite-order vectors + KVectorComponent modified_kx_vec, modified_kz_vec; +#if (AMREX_SPACEDIM==3) + KVectorComponent modified_ky_vec; +#endif +}; + +#endif // WARPX_SPECTRAL_BASE_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 9e31ac5b8..8e58aa1d8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -47,7 +47,8 @@ class SpectralFieldData SpectralField fields; // tmpRealField and tmpSpectralField store fields // right before/after the Fourier transform - SpectralField tmpRealField, tmpSpectralField; + SpectralField tmpSpectralField; // contains Complexs + amrex::MultiFab tmpRealField; // contains Reals 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 6e6cc124f..02fa2015f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -16,7 +16,7 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // 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); + tmpRealField = MultiFab(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); // By default, we assume the FFT is done from/to a nodal grid in real space @@ -48,7 +48,10 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // 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]; + // Note: the size of the real-space box and spectral-space box + // differ when using real-to-complex FFT. When initializing + // the FFT plan, the valid dimensions are those of the real-space box. + IntVect fft_size = realspace_ba[mfi].length(); #ifdef AMREX_USE_GPU // Add cuFFT-specific code #else @@ -56,23 +59,23 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, forward_plan[mfi] = // 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), + fftw_plan_dft_r2c_3d( fft_size[2], fft_size[1], fft_size[0], #else - fftw_plan_dft_2d( bx.length(1), bx.length(0), + fftw_plan_dft_r2c_2d( fft_size[1], fft_size[0], #endif - reinterpret_cast<fftw_complex*>( tmpRealField[mfi].dataPtr() ), + tmpRealField[mfi].dataPtr(), reinterpret_cast<fftw_complex*>( tmpSpectralField[mfi].dataPtr() ), - FFTW_FORWARD, FFTW_ESTIMATE ); + FFTW_ESTIMATE ); backward_plan[mfi] = // 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), + fftw_plan_dft_c2r_3d( fft_size[2], fft_size[1], fft_size[0], #else - fftw_plan_dft_2d( bx.length(1), bx.length(0), + fftw_plan_dft_c2r_2d( fft_size[1], fft_size[0], #endif reinterpret_cast<fftw_complex*>( tmpSpectralField[mfi].dataPtr() ), - reinterpret_cast<fftw_complex*>( tmpRealField[mfi].dataPtr() ), - FFTW_BACKWARD, FFTW_ESTIMATE ); + tmpRealField[mfi].dataPtr(), + FFTW_ESTIMATE ); #endif } } @@ -123,7 +126,7 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, realspace_bx.enclosedCells(); // Discard last point in nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4<const Real> mf_arr = mf[mfi].array(); - Array4<Complex> tmp_arr = tmpRealField[mfi].array(); + Array4<Real> tmp_arr = tmpRealField[mfi].array(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); @@ -194,7 +197,6 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // 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 { Array4<const Complex> field_arr = SpectralFieldData::fields[mfi].array(); Array4<Complex> tmp_arr = tmpSpectralField[mfi].array(); @@ -205,8 +207,6 @@ 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,field_index); @@ -218,8 +218,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, #elif (AMREX_SPACEDIM == 2) if (is_nodal_z==false) spectral_field_value *= zshift_arr[j]; #endif - // Copy field into temporary array (after normalization) - tmp_arr(i,j,k) = inv_N*spectral_field_value; + // Copy field into temporary array + tmp_arr(i,j,k) = spectral_field_value; }); } @@ -232,13 +232,18 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, #endif // Copy the temporary field `tmpRealField` to the real-space field `mf` + + // Normalize (divide by 1/N) since the FFT+IFFT results in a factor N { const Box realspace_bx = tmpRealField[mfi].box(); Array4<Real> mf_arr = mf[mfi].array(); - Array4<const Complex> tmp_arr = tmpRealField[mfi].array(); + Array4<const Real> tmp_arr = tmpRealField[mfi].array(); + // Normalization: divide by the number of points in realspace + const Real inv_N = 1./realspace_bx.numPts(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); + // Copy and normalize field + mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k); }); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index ad17e6423..ae7124b2e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -33,7 +33,9 @@ class SpectralKSpace const amrex::DistributionMapping& dm, const amrex::RealVect realspace_dx ); KVectorComponent getKComponent( - const amrex::DistributionMapping& dm, const int i_dim ) const; + const amrex::DistributionMapping& dm, + const amrex::BoxArray& realspace_ba, + const int i_dim, const bool only_positive_k ) const; KVectorComponent getModifiedKComponent( const amrex::DistributionMapping& dm, const int i_dim, const int n_order, const bool nodal ) const; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index ddb2020d8..2fe78cedd 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -28,19 +28,33 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // 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]; - 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 ); + IntVect fft_size = realspace_bx.length(); + // Because the spectral solver uses real-to-complex FFTs, we only + // need the positive k values along the fastest axis + // (first axis for AMReX Fortran-order arrays) in spectral space. + // This effectively reduces the size of the spectral space by half + // see e.g. the FFTW documentation for real-to-complex FFTs + IntVect spectral_bx_size = fft_size; + spectral_bx_size[0] = fft_size[0]/2 + 1; + // Define the corresponding box + Box spectral_bx = Box( IntVect::TheZeroVector(), + spectral_bx_size - IntVect::TheUnitVector() ); + spectral_bl.push_back( spectral_bx ); } spectralspace_ba.define( spectral_bl ); // Allocate the components of the k vector: kx, ky (only in 3D), kz + bool only_positive_k; for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++) { - k_vec[i_dim] = getKComponent( dm, i_dim ); + if (i_dim==0) { + // Real-to-complex FFTs: first axis contains only the positive k + only_positive_k = true; + } else { + only_positive_k = false; + } + k_vec[i_dim] = getKComponent(dm, realspace_ba, i_dim, only_positive_k); } } @@ -50,7 +64,9 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, */ KVectorComponent SpectralKSpace::getKComponent( const DistributionMapping& dm, - const int i_dim ) const + const BoxArray& realspace_ba, + const int i_dim, + const bool only_positive_k ) const { // Initialize an empty ManagedVector in each box KVectorComponent k_comp(spectralspace_ba, dm); @@ -65,21 +81,31 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, k.resize( N ); // Fill the k vector - const Real dk = 2*MathConst::pi/(N*dx[i_dim]); + IntVect fft_size = realspace_ba[mfi].length(); + const Real dk = 2*MathConst::pi/(fft_size[i_dim]*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."); - const int mid_point = (N+1)/2; - // Fill positive values of k (FFT conventions: first half is positive) - for (int i=0; i<mid_point; i++ ){ - k[i] = i*dk; - } - // Fill negative values of k (FFT conventions: second half is negative) - for (int i=mid_point; i<N; i++){ - k[i] = (i-N)*dk; + if (only_positive_k){ + // Fill the full axis with positive k values + // (typically: first axis, in a real-to-complex FFT) + for (int i=0; i<N; i++ ){ + k[i] = i*dk; + } + } else { + const int mid_point = (N+1)/2; + // Fill positive values of k + // (FFT conventions: first half is positive) + for (int i=0; i<mid_point; i++ ){ + k[i] = i*dk; + } + // Fill negative values of k + // (FFT conventions: second half is negative) + for (int i=mid_point; i<N; i++){ + k[i] = (i-N)*dk; + } } - // TODO: this will be different for the real-to-complex transform // TODO: this will be different for the hybrid FFT scheme } return k_comp; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 7444452af..d4019a9a3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -1,8 +1,7 @@ #ifndef WARPX_SPECTRAL_SOLVER_H_ #define WARPX_SPECTRAL_SOLVER_H_ -#include <SpectralKSpace.H> -#include <PsatdAlgorithm.H> +#include <SpectralBaseAlgorithm.H> #include <SpectralFieldData.H> /* \brief Top-level class for the electromagnetic spectral solver @@ -14,28 +13,17 @@ class SpectralSolver { public: - // Inline definition of the member functions of `SpectralSolver` + // Inline definition of the member functions of `SpectralSolver`, + // except the constructor (see `SpectralSolver.cpp`) // The body of these functions is short, since the work is done in the // underlying classes `SpectralFieldData` and `PsatdAlgorithm` - /* \brief Initialize the spectral solver */ + // Constructor SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { - // Initialize all structures using the same distribution mapping dm - - // - Initialize k space object (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) - const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); - // - Initialize the algorithm (coefficients) over this space - 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 ); - }; + const amrex::RealVect dx, const amrex::Real dt ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally @@ -59,14 +47,20 @@ class SpectralSolver /* \brief Update the fields in spectral space, over one timestep */ void pushSpectralFields(){ BL_PROFILE("SpectralSolver::pushSpectralFields"); - algorithm.pushSpectralFields( field_data ); + // Virtual function: the actual function used here depends + // on the sub-class of `SpectralBaseAlgorithm` that was + // initialized in the constructor of `SpectralSolver` + 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 + std::unique_ptr<SpectralBaseAlgorithm> algorithm; + // Defines field update equation in spectral space, + // and the associated coefficients. + // SpectralBaseAlgorithm is a base class ; this pointer is meant + // to point an instance of a *sub-class* defining a specific algorithm }; #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp new file mode 100644 index 000000000..c21c3cfb1 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -0,0 +1,35 @@ +#include <SpectralKSpace.H> +#include <SpectralSolver.H> +#include <PsatdAlgorithm.H> + +/* \brief Initialize the spectral Maxwell solver + * + * This function selects the spectral algorithm to be used, allocates the + * corresponding coefficients for the discretized field update equation, + * and prepares the structures that store the fields in spectral space. + */ +SpectralSolver::SpectralSolver( + const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::RealVect dx, const amrex::Real dt ) { + + // Initialize all structures using the same distribution mapping dm + + // - Initialize k space object (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) + const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); + + // - Select the algorithm depending on the input parameters + // Initialize the corresponding coefficients over k space + // TODO: Add more algorithms + selection depending on input parameters + // For the moment, this only uses the standard PsatdAlgorithm + algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + + // - Initialize arrays for fields in spectral space + FFT plans + field_data = SpectralFieldData( realspace_ba, k_space, dm ); + +}; diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 473ca645a..f246f9f54 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -81,17 +81,16 @@ extern "C" { #endif -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS void warpx_copy_attribs(const long* np, const amrex_real* xp, const amrex_real* yp, const amrex_real* zp, const amrex_real* uxp, const amrex_real* uyp, const amrex_real* uzp, amrex_real* xpold, amrex_real* ypold, amrex_real* zpold, amrex_real* uxpold, amrex_real* uypold, amrex_real* uzpold); -#endif - void WRPX_COPY_SLICE(const int* lo, const int* hi, - const amrex_real* tmp, const int* tlo, const int* thi, - amrex_real* buf, const int* blo, const int* bhi, - const int* ncomp, const int* i_boost, const int* i_lab); + + void WRPX_COPY_SLICE(const int* lo, const int* hi, + const amrex_real* tmp, const int* tlo, const int* thi, + amrex_real* buf, const int* blo, const int* bhi, + const int* ncomp, const int* i_boost, const int* i_lab); // Charge deposition void warpx_charge_deposition(amrex::Real* rho, diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 926523d82..95e2b4ec0 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -123,10 +123,6 @@ ifeq ($(USE_RZ),TRUE) DEFINES += -DWARPX_RZ endif -ifeq ($(STORE_OLD_PARTICLE_ATTRIBS),TRUE) - DEFINES += -DWARPX_STORE_OLD_PARTICLE_ATTRIBS -endif - ifeq ($(DO_ELECTROSTATIC),TRUE) include $(AMREX_HOME)/Src/LinearSolvers/C_to_F_MG/Make.package include $(AMREX_HOME)/Src/LinearSolvers/F_MG/FParallelMG.mak diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 0485f7b8f..a4df1f83a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -29,7 +29,26 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[i].reset(new LaserParticleContainer(amr_core,i, lasers_names[i-nspecies])); } - pc_tmp.reset(new PhysicalParticleContainer(amr_core)); + pc_tmp.reset(new PhysicalParticleContainer(amr_core)); + + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + for (int i = 0; i < nspecies + nlasers; ++i) + { + allcontainers[i]->AddRealComp("xold"); + allcontainers[i]->AddRealComp("yold"); + allcontainers[i]->AddRealComp("zold"); + allcontainers[i]->AddRealComp("uxold"); + allcontainers[i]->AddRealComp("uyold"); + allcontainers[i]->AddRealComp("uzold"); + } + pc_tmp->AddRealComp("xold"); + pc_tmp->AddRealComp("yold"); + pc_tmp->AddRealComp("zold"); + pc_tmp->AddRealComp("uxold"); + pc_tmp->AddRealComp("uyold"); + pc_tmp->AddRealComp("uzold"); + } } void diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index e31d43204..17e6d98d9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -184,6 +184,18 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } + AddOneParticle(0, 0, 0, x, y, z, attribs); } } @@ -455,16 +467,18 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::ux] = u[0]; attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - attribs[PIdx::xold] = x; - attribs[PIdx::yold] = y; - attribs[PIdx::zold] = z; - - attribs[PIdx::uxold] = u[0]; - attribs[PIdx::uyold] = u[1]; - attribs[PIdx::uzold] = u[2]; -#endif + + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); } @@ -695,15 +709,18 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - attribs[PIdx::xold] = x; - attribs[PIdx::yold] = y; - attribs[PIdx::zold] = z; - - attribs[PIdx::uxold] = u[0]; - attribs[PIdx::uyold] = u[1]; - attribs[PIdx::uzold] = u[2]; -#endif + // note - this will be slow on the GPU, need to revisit + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } ParticleType p; p.id() = ParticleType::NextID(); @@ -1666,20 +1683,20 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - auto& xpold = attribs[PIdx::xold]; - auto& ypold = attribs[PIdx::yold]; - auto& zpold = attribs[PIdx::zold]; - auto& uxpold = attribs[PIdx::uxold]; - auto& uypold = attribs[PIdx::uyold]; - auto& uzpold = attribs[PIdx::uzold]; - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - -#endif + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& xpold = pti.GetAttribs(particle_comps["xold"]); + auto& ypold = pti.GetAttribs(particle_comps["yold"]); + auto& zpold = pti.GetAttribs(particle_comps["zold"]); + auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); + auto& uypold = pti.GetAttribs(particle_comps["uyold"]); + auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); + + warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), + uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); + } warpx_particle_pusher(&np, xp.dataPtr(), @@ -1801,7 +1818,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real { BL_PROFILE("PhysicalParticleContainer::GetParticleSlice"); -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS // Assume that the boost in the positive z direction. #if (AMREX_SPACEDIM == 2) AMREX_ALWAYS_ASSERT(direction == 1); @@ -1864,12 +1880,12 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = attribs[PIdx::xold ]; - auto& yp_old = attribs[PIdx::yold ]; - auto& zp_old = attribs[PIdx::zold ]; - auto& uxp_old = attribs[PIdx::uxold]; - auto& uyp_old = attribs[PIdx::uyold]; - auto& uzp_old = attribs[PIdx::uzold]; + auto& xp_old = pti.GetAttribs(particle_comps["xold"]); + auto& yp_old = pti.GetAttribs(particle_comps["yold"]); + auto& zp_old = pti.GetAttribs(particle_comps["zold"]); + auto& uxp_old = pti.GetAttribs(particle_comps["uxold"]); + auto& uyp_old = pti.GetAttribs(particle_comps["uyold"]); + auto& uzp_old = pti.GetAttribs(particle_comps["uzold"]); const long np = pti.numParticles(); @@ -1919,10 +1935,6 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real } } } -#else - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false , -"ERROR: WarpX must be compiled with STORE_OLD_PARTICLE_ATTRIBS=TRUE to use the back-transformed diagnostics"); -#endif } int PhysicalParticleContainer::GetRefineFac(const Real x, const Real y, const Real z) diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 3ee4d87e5..a5acca281 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -225,20 +225,20 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - auto& xpold = attribs[PIdx::xold]; - auto& ypold = attribs[PIdx::yold]; - auto& zpold = attribs[PIdx::zold]; - auto& uxpold = attribs[PIdx::uxold]; - auto& uypold = attribs[PIdx::uyold]; - auto& uzpold = attribs[PIdx::uzold]; - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - -#endif + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& xpold = pti.GetAttribs(particle_comps["xold"]); + auto& ypold = pti.GetAttribs(particle_comps["yold"]); + auto& zpold = pti.GetAttribs(particle_comps["zold"]); + auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); + auto& uypold = pti.GetAttribs(particle_comps["uyold"]); + auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); + + warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), + uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); + } // Save the position and momenta, making copies Cuda::ManagedDeviceVector<Real> xp_save, yp_save, zp_save; diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 050060b47..275554cd8 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -14,9 +14,6 @@ struct PIdx #ifdef WARPX_RZ theta, // RZ needs all three position components #endif -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - xold, yold, zold, uxold, uyold, uzold, -#endif nattribs }; }; @@ -42,15 +39,7 @@ namespace ParticleStringNames {"Ez", PIdx::Ez }, {"Bx", PIdx::Bx }, {"By", PIdx::By }, - {"Bz", PIdx::Bz }, -#ifdef WARPX_STORE_OLD_PARTICLE_ATTRIBS - {"xold", PIdx::xold }, - {"yold", PIdx::yold }, - {"zold", PIdx::zold }, - {"uxold", PIdx::uxold}, - {"uyold", PIdx::uyold}, - {"uzold", PIdx::uzold}, -#endif + {"Bz", PIdx::Bz } }; } @@ -231,8 +220,25 @@ public: // split along axes (0) or diagonals (1) int split_type = 0; + using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; + using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddIntComp; + + void AddRealComp (const std::string& name, bool comm=true) + { + particle_comps[name] = NumRealComps(); + AddRealComp(comm); + } + + void AddIntComp (const std::string& name, bool comm=true) + { + particle_comps[name] = NumIntComps(); + AddIntComp(comm); + } + protected: + std::map<std::string, int> particle_comps; + int species_id; amrex::Real charge; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c52e0a6d0..2edd3c636 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -63,6 +63,32 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) SetParticleSize(); ReadParameters(); + // build up the map of string names to particle component numbers + particle_comps["w"] = PIdx::w; + particle_comps["ux"] = PIdx::ux; + particle_comps["uy"] = PIdx::uy; + particle_comps["uz"] = PIdx::uz; + particle_comps["Ex"] = PIdx::Ex; + particle_comps["Ey"] = PIdx::Ey; + particle_comps["Ez"] = PIdx::Ez; + particle_comps["Bx"] = PIdx::Bx; + particle_comps["By"] = PIdx::By; + particle_comps["Bz"] = PIdx::Bz; +#ifdef WARPX_RZ + particle_comps["theta"] = PIdx::theta; +#endif + + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + particle_comps["xold"] = PIdx::nattribs; + particle_comps["yold"] = PIdx::nattribs+1; + particle_comps["zold"] = PIdx::nattribs+2; + particle_comps["uxold"] = PIdx::nattribs+3; + particle_comps["uyold"] = PIdx::nattribs+4; + particle_comps["uzold"] = PIdx::nattribs+5; + + } + // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -204,6 +230,15 @@ WarpXParticleContainer::AddNParticles (int lev, #endif p.pos(1) = z[i]; #endif + + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x[i]); + particle_tile.push_back_real(particle_comps["yold"], y[i]); + particle_tile.push_back_real(particle_comps["zold"], z[i]); + } + particle_tile.push_back(p); } @@ -214,6 +249,14 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + } + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_RZ diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d4ddca90e..793b96db7 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -448,11 +448,25 @@ WarpX::ReadParameters () if (particle_plot_vars.size() == 0) { - particle_plot_flags.resize(PIdx::nattribs, 1); + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + particle_plot_flags.resize(PIdx::nattribs + 6, 1); + } + else + { + particle_plot_flags.resize(PIdx::nattribs, 1); + } } else { - particle_plot_flags.resize(PIdx::nattribs, 0); + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + { + particle_plot_flags.resize(PIdx::nattribs + 6, 0); + } + else + { + particle_plot_flags.resize(PIdx::nattribs, 0); + } for (const auto& var : particle_plot_vars) { |