diff options
Diffstat (limited to 'Source')
23 files changed, 603 insertions, 19 deletions
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index 7340d568e..b019ad1f5 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -351,6 +351,7 @@ WarpXOpenPMDPlot::DumpToFile (WarpXParticleContainer* pc, currSpecies.setAttribute( "currentDeposition", [](){ switch( WarpX::current_deposition_algo ) { case CurrentDepositionAlgo::Esirkepov : return "Esirkepov"; + case CurrentDepositionAlgo::Vay : return "Vay"; default: return "directMorseNielson"; } }() ); diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 674680773..fda7100d4 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -296,6 +296,9 @@ WarpX::OneStep_nosub (Real cur_time) if ( !fft_periodic_single_box && current_correction ) amrex::Abort("\nCurrent correction does not guarantee charge conservation with local FFTs over guard cells:\n" "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation"); + if ( !fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ) + amrex::Abort("\nVay current deposition does not guarantee charge conservation with local FFTs over guard cells:\n" + "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation"); #endif #ifdef WARPX_QED @@ -310,6 +313,8 @@ WarpX::OneStep_nosub (Real cur_time) // without guard cells, apply this after calling SyncCurrent #ifdef WARPX_USE_PSATD if ( fft_periodic_single_box && current_correction ) CurrentCorrection(); + if ( fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ) + VayDeposition(); #endif @@ -792,19 +797,28 @@ WarpX::applyMirrors(Real time){ } } +// Apply current correction in Fourier space #ifdef WARPX_USE_PSATD void WarpX::CurrentCorrection () { for ( int lev = 0; lev <= finest_level; ++lev ) { - // Apply correction on fine patch spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] ); - if ( spectral_solver_cp[lev] ) - { - // Apply correction on coarse patch - spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] ); - } + if ( spectral_solver_cp[lev] ) spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] ); + } +} +#endif + +// Compute current from Vay deposition in Fourier space +#ifdef WARPX_USE_PSATD +void +WarpX::VayDeposition () +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + spectral_solver_fp[lev]->VayDeposition(current_fp[lev]); + if (spectral_solver_cp[lev]) spectral_solver_cp[lev]->VayDeposition(current_cp[lev]); } } #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H index f2ae61f33..45a35098b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H @@ -26,6 +26,20 @@ class AvgGalileanAlgorithm : public SpectralBaseAlgorithm const amrex::Array<amrex::Real, 3>& v_galilean, const amrex::Real dt); + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + private: SpectralRealCoefficients C_coef, S_ck_coef, C1_coef, C3_coef, S1_coef,S3_coef; SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef, Psi1_coef, Psi2_coef, Psi3_coef, Psi4_coef, A1_coef, A2_coef, Rhoold_coef, Rhonew_coef, Jcoef_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp index d9db745f1..0ba603031 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp @@ -369,3 +369,10 @@ AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ }); } }; + +void +AvgGalileanAlgorithm::VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +{ + amrex::Abort("Vay deposition not implemented for averaged Galilean PSATD"); +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H index 192c453ad..799eca8b1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H @@ -26,6 +26,20 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm const amrex::Array<amrex::Real, 3>& v_galilean, const amrex::Real dt); + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + private: SpectralRealCoefficients C_coef, S_ck_coef; SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp index 49e318e3c..e606b3232 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp @@ -243,4 +243,12 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe }); } } + +void +GalileanAlgorithm::VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +{ + amrex::Abort("Vay deposition not implemented for Galilean PSATD"); +} + #endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H index ef12ff12e..347012d5a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H @@ -34,6 +34,20 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm return SpectralPMLIndex::n_fields; } + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + private: SpectralRealCoefficients C_coef, S_ck_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp index b3f2aa93d..d2f087706 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp @@ -154,4 +154,12 @@ void PMLPsatdAlgorithm::InitializeSpectralCoefficients ( }); } }; + +void +PMLPsatdAlgorithm::VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +{ + amrex::Abort("Vay deposition not implemented for PML PSATD"); +} + #endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index b492d9a18..bfa9283e6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -53,6 +53,20 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm std::array<std::unique_ptr<amrex::MultiFab>,3>& current, const std::unique_ptr<amrex::MultiFab>& rho ) override final; + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + private: SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; amrex::Real m_dt; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index b8b64b80c..7f9fd3edb 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -256,7 +256,6 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, // Loop over indices within one box ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Record old values of the fields to be updated using Idx = SpectralFieldIndex; // Shortcuts for the values of J and rho @@ -297,4 +296,80 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, field_data.BackwardTransform( *current[1], Idx::Jy, 0 ); field_data.BackwardTransform( *current[2], Idx::Jz, 0 ); } + +void +PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) { + // Profiling + WARPX_PROFILE("PsatdAlgorithm::VayDeposition"); + + using Idx = SpectralFieldIndex; + + // Forward Fourier transform of D (temporarily stored in current): + // D is nodal and does not match the staggering of J, therefore we pass the + // actual staggering of D (IntVect(1)) to the ForwardTransform function + field_data.ForwardTransform(*current[0], Idx::Jx, 0, IntVect(1)); + field_data.ForwardTransform(*current[1], Idx::Jy, 0, IntVect(1)); + field_data.ForwardTransform(*current[2], Idx::Jz, 0, IntVect(1)); + + // Loop over boxes + for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi) { + + const amrex::Box& bx = field_data.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4<Complex> fields = field_data.fields[mfi].array(); + + // Extract pointers for the modified k vectors + const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + + // Loop over indices within one box + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + using Idx = SpectralFieldIndex; + + // Shortcuts for the values of D + const Complex Dx = fields(i,j,k,Idx::Jx); + const Complex Dy = fields(i,j,k,Idx::Jy); + const Complex Dz = fields(i,j,k,Idx::Jz); + + // Imaginary unit + constexpr Complex I = Complex{0._rt, 1._rt}; + + // Modified k vector values + const amrex::Real kx_mod = modified_kx_arr[i]; +#if (AMREX_SPACEDIM==3) + const amrex::Real ky_mod = modified_ky_arr[j]; + const amrex::Real kz_mod = modified_kz_arr[k]; +#else + constexpr amrex::Real ky_mod = 0._rt; + const amrex::Real kz_mod = modified_kz_arr[j]; +#endif + + // Compute Jx + if (kx_mod != 0._rt) fields(i,j,k,Idx::Jx) = I * Dx / kx_mod; + else fields(i,j,k,Idx::Jx) = 0._rt; + +#if (AMREX_SPACEDIM==3) + // Compute Jy + if (ky_mod != 0._rt) fields(i,j,k,Idx::Jy) = I * Dy / ky_mod; + else fields(i,j,k,Idx::Jy) = 0._rt; +#endif + + // Compute Jz + if (kz_mod != 0._rt) fields(i,j,k,Idx::Jz) = I * Dz / kz_mod; + else fields(i,j,k,Idx::Jz) = 0._rt; + + }); + } + + // Backward Fourier transform of J + field_data.BackwardTransform(*current[0], Idx::Jx, 0); + field_data.BackwardTransform(*current[1], Idx::Jy, 0); + field_data.BackwardTransform(*current[2], Idx::Jz, 0); +} #endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index 6a604155b..fc03f95f2 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -44,6 +44,20 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ std::array<std::unique_ptr<amrex::MultiFab>,3>& current, const std::unique_ptr<amrex::MultiFab>& rho ) override final; + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldDataRZ& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + private: bool coefficients_initialized; // Note that dt is saved to use in InitializeSpectralCoefficients diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index b094843ac..35a478f2c 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -204,3 +204,10 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data, { amrex::Abort("Current correction not implemented in RZ"); } + +void +PsatdAlgorithmRZ::VayDeposition (SpectralFieldDataRZ& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +{ + amrex::Abort("Vay deposition not implemented in RZ geometry"); +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index 7143b7aa9..2a34d21ab 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -48,6 +48,18 @@ class SpectralBaseAlgorithm const std::unique_ptr<amrex::MultiFab>& rho ) {}; /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This virtual function is pure and must be defined in derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0; + + /** * \brief Compute spectral divergence of E */ void ComputeSpectralDivE ( SpectralFieldData& field_data, diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H index ea50e979a..833a61aec 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H @@ -51,6 +51,18 @@ class SpectralBaseAlgorithmRZ const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This virtual function is pure and must be defined in derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + virtual void VayDeposition (SpectralFieldDataRZ& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0; + protected: // Meant to be used in the subclasses using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index b445054cc..f48272744 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -55,10 +55,17 @@ class SpectralFieldData SpectralFieldData() = default; // Default constructor SpectralFieldData& operator=(SpectralFieldData&& field_data) = default; ~SpectralFieldData(); - void ForwardTransform( const amrex::MultiFab& mf, - const int field_index, const int i_comp); - void BackwardTransform( amrex::MultiFab& mf, - const int field_index, const int i_comp); + + void ForwardTransform (const amrex::MultiFab& mf, const int field_index, + const int i_comp, const amrex::IntVect& stag); + AMREX_FORCE_INLINE + void ForwardTransform (const amrex::MultiFab& mf, const int field_index, const int i_comp) + { + ForwardTransform(mf, field_index, i_comp, mf.ixType().toIntVect()); + } + + void BackwardTransform (amrex::MultiFab& mf, const int field_index, const int i_comp); + // `fields` stores fields in spectral space, as multicomponent FabArray SpectralField fields; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 0f49e695b..e99e7d57b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -94,17 +94,16 @@ SpectralFieldData::~SpectralFieldData() * 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 ) +SpectralFieldData::ForwardTransform (const MultiFab& mf, const int field_index, + const int i_comp, const IntVect& stag) { // Check field index type, in order to apply proper shift in spectral space - const bool is_nodal_x = mf.is_nodal(0); + const bool is_nodal_x = (stag[0] == amrex::IndexType::NODE) ? true : false; #if (AMREX_SPACEDIM == 3) - const bool is_nodal_y = mf.is_nodal(1); - const bool is_nodal_z = mf.is_nodal(2); + const bool is_nodal_y = (stag[1] == amrex::IndexType::NODE) ? true : false; + const bool is_nodal_z = (stag[2] == amrex::IndexType::NODE) ? true : false; #else - const bool is_nodal_z = mf.is_nodal(1); + const bool is_nodal_z = (stag[1] == amrex::IndexType::NODE) ? true : false; #endif // Loop over boxes diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 91cfc84f1..ba805c790 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -84,6 +84,20 @@ class SpectralSolver algorithm->CurrentCorrection( field_data, current, rho ); }; + /** + * \brief Public interface to call the virtual function \c VayDeposition, + * declared in the base class SpectralBaseAlgorithm and defined in its + * derived classes, from objects of class SpectralSolver through the private + * unique pointer \c algorithm. + * + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + void VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current) + { + algorithm->VayDeposition(field_data, current); + } + private: void ReadParameters (); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 5762e1436..e1179d9cf 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -104,6 +104,20 @@ class SpectralSolverRZ algorithm->CurrentCorrection( field_data, current, rho ); }; + /** + * \brief Public interface to call the virtual function \c VayDeposition, + * declared in the base class SpectralBaseAlgorithmRZ and defined in its + * derived classes, from objects of class SpectralSolverRZ through the private + * unique pointer \c algorithm. + * + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + */ + void VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current) + { + algorithm->VayDeposition(field_data, current); + } + private: SpectralFieldDataRZ field_data; // Store field in spectral space diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index f434edff8..4e547e52c 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -15,6 +15,8 @@ #include <AMReX_Array4.H> #include <AMReX_REAL.H> +using namespace amrex::literals; + /** * \brief Current Deposition for thread thread_num * /param GetPosition : A functor for returning the particle position. @@ -541,4 +543,307 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, ); } +/** + * \brief Vay current deposition + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>) + * for thread \c thread_num: deposit \c D in real space and store the result in + * \c jx_fab, \c jy_fab, \c jz_fab + * + * \param[in] GetPosition Functor that returns the particle position + * \param[in] wp Pointer to array of particle weights + * \param[in] uxp Pointer to arrays of particle momentum along \c x + * \param[in] uyp Pointer to arrays of particle momentum along \c y + * \param[in] uzp Pointer to arrays of particle momentum along \c z + * \param[in] ion_lev Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle since \c q + is a scalar. For non-ionizable species, \c ion_lev is \c null + * \param[in,out] jx_fab FArrayBox of current density, either full array or tile + * \param[in,out] jy_fab FArrayBox of current density, either full array or tile + * \param[in,out] jz_fab FArrayBox of current density, either full array or tile + * \param[in] np_to_depose Number of particles for which current is deposited + * \param[in] dt Time step for particle level + * \param[in] dx 3D cell size + * \param[in] xyzmin 3D lower bounds of physical domain + * \param[in] lo Dimension-agnostic lower bounds of index domain + * \param[in] q Species charge + * \param[in] n_rz_azimuthal_modes Number of azimuthal modes in RZ geometry + */ +template <int depos_order> +void doVayDepositionShapeN (const GetParticlePosition& GetPosition, + const amrex::ParticleReal* const wp, + const amrex::ParticleReal* const uxp, + const amrex::ParticleReal* const uyp, + const amrex::ParticleReal* const uzp, + const int* const ion_lev, + amrex::FArrayBox& jx_fab, + amrex::FArrayBox& jy_fab, + amrex::FArrayBox& jz_fab, + const long np_to_depose, + const amrex::Real dt, + const std::array<amrex::Real,3>& dx, + const std::array<amrex::Real,3>& xyzmin, + const amrex::Dim3 lo, + const amrex::Real q, + const long n_rz_azimuthal_modes) +{ +#if (defined WARPX_DIM_RZ) + amrex::Abort("Vay deposition not implemented in RZ geometry"); +#endif + +#if !(defined WARPX_DIM_RZ) + + // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1 + const bool do_ionization = ion_lev; + + // Inverse cell volume in each direction + const amrex::Real dxi = 1._rt / dx[0]; + const amrex::Real dzi = 1._rt / dx[2]; +#if (defined WARPX_DIM_3D) + const amrex::Real dyi = 1._rt / dx[1]; +#endif + + // Inverse of time step + const amrex::Real invdt = 1._rt / dt; + + // Total inverse cell volume +#if (defined WARPX_DIM_XZ) + const amrex::Real invvol = dxi * dzi; +#elif (defined WARPX_DIM_3D) + const amrex::Real invvol = dxi * dyi * dzi; +#endif + + // Lower bound of physical domain in each direction + const amrex::Real xmin = xyzmin[0]; + const amrex::Real zmin = xyzmin[2]; +#if (defined WARPX_DIM_3D) + const amrex::Real ymin = xyzmin[1]; +#endif + + // Auxiliary constants +#if (defined WARPX_DIM_3D) + const amrex::Real onethird = 1._rt / 3._rt; + const amrex::Real onesixth = 1._rt / 6._rt; +#endif + + // Inverse of light speed squared + const amrex::Real invcsq = 1._rt / (PhysConst::c * PhysConst::c); + + // Arrays where D will be stored + amrex::Array4<amrex::Real> const& jx_arr = jx_fab.array(); + amrex::Array4<amrex::Real> const& jy_arr = jy_fab.array(); + amrex::Array4<amrex::Real> const& jz_arr = jz_fab.array(); + + // D must have nodal index type in each direction + amrex::IntVect const jx_type = amrex::IntVect(1); + amrex::IntVect const jy_type = amrex::IntVect(1); + amrex::IntVect const jz_type = amrex::IntVect(1); + + // Loop over particles and deposit (Dx,Dy,Dz) into jx_fab, jy_fab and jz_fab + amrex::ParallelFor(np_to_depose, [=] AMREX_GPU_DEVICE (long ip) + { + // Inverse of Lorentz factor gamma + const amrex::Real invgam = 1._rt / std::sqrt(1._rt + uxp[ip] * uxp[ip] * invcsq + + uyp[ip] * uyp[ip] * invcsq + + uzp[ip] * uzp[ip] * invcsq); + // Product of particle charges and weights + amrex::Real wq = q * wp[ip]; + if (do_ionization) wq *= ion_lev[ip]; + + // Current particle positions (in physical units) + amrex::ParticleReal xp, yp, zp; + GetPosition(ip, xp, yp, zp); + + // Particle velocities + const amrex::Real vx = uxp[ip] * invgam; + const amrex::Real vy = uyp[ip] * invgam; + const amrex::Real vz = uzp[ip] * invgam; + + // Particle current densities +#if (defined WARPX_DIM_XZ) + const amrex::Real wqy = wq * vy * invvol; +#endif + + // Current and old particle positions in grid units + amrex::Real const x_new = (xp - xmin) * dxi; + amrex::Real const x_old = x_new - vx * dt * dxi; +#if (defined WARPX_DIM_3D) + amrex::Real const y_new = (yp - ymin) * dyi; + amrex::Real const y_old = y_new - vy * dt * dyi; +#endif + amrex::Real const z_new = (zp - zmin) * dzi; + amrex::Real const z_old = z_new - vz * dt * dzi; + + // Shape factor arrays for current and old positions (nodal) + amrex::Real sx_new[depos_order+1] = {0._rt}; + amrex::Real sx_old[depos_order+1] = {0._rt}; +#if (defined WARPX_DIM_3D) + amrex::Real sy_new[depos_order+1] = {0._rt}; + amrex::Real sy_old[depos_order+1] = {0._rt}; +#endif + amrex::Real sz_new[depos_order+1] = {0._rt}; + amrex::Real sz_old[depos_order+1] = {0._rt}; + + // Compute shape factors for current positions + + // i_new leftmost grid point in x that the particle touches + // sx_new shape factor along x for the centering of each current + const int i_new = compute_shape_factor<depos_order>(sx_new, x_new); +#if (defined WARPX_DIM_3D) + // j_new leftmost grid point in y that the particle touches + // sy_new shape factor along y for the centering of each current + const int j_new = compute_shape_factor<depos_order>(sy_new, y_new); +#endif + // k_new leftmost grid point in z that the particle touches + // sz_new shape factor along z for the centering of each current + const int k_new = compute_shape_factor<depos_order>(sz_new, z_new); + + // Compute shape factors for old positions + + // i_old leftmost grid point in x that the particle touches + // sx_old shape factor along x for the centering of each current + const int i_old = compute_shape_factor<depos_order>(sx_old, x_old); +#if (defined WARPX_DIM_3D) + // j_old leftmost grid point in y that the particle touches + // sy_old shape factor along y for the centering of each current + const int j_old = compute_shape_factor<depos_order>(sy_old, y_old); +#endif + // k_old leftmost grid point in z that the particle touches + // sz_old shape factor along z for the centering of each current + const int k_old = compute_shape_factor<depos_order>(sz_old, z_old); + + // Deposit current into jx_arr, jy_arr and jz_arr + +#if (defined WARPX_DIM_XZ) + + for (int k=0; k<=depos_order; k++) { + for (int i=0; i<=depos_order; i++) { + + // Jx + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0), + wq * invvol * invdt * 0.5_rt * sx_new[i] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0), + - wq * invvol * invdt * 0.5_rt * sx_old[i] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0), + wq * invvol * invdt * 0.5_rt * sx_new[i] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0), + - wq * invvol * invdt * 0.5_rt * sx_old[i] * sz_old[k]); + + // Jy + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0), + wqy * 0.25_rt * sx_new[i] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + k_old + k, 0, 0), + wqy * 0.25_rt * sx_new[i] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + k_new + k, 0, 0), + wqy * 0.25_rt * sx_old[i] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0), + wqy * 0.25_rt * sx_old[i] * sz_old[k]); + + // Jz + amrex::Gpu::Atomic::Add(&jz_arr(lo.x + i_new + i, lo.y + k_new + k, 0, 0), + wq * invvol * invdt * 0.5_rt * sx_new[i] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jz_arr(lo.x+i_new+i,lo.y+k_old+k,0,0), + - wq * invvol * invdt * 0.5_rt * sx_new[i] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jz_arr(lo.x+i_old+i,lo.y+k_new+k,0,0), + wq * invvol * invdt * 0.5_rt * sx_old[i] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jz_arr(lo.x + i_old + i, lo.y + k_old + k, 0, 0), + - wq * invvol * invdt * 0.5_rt * sx_old[i] * sz_old[k]); + } + } + +#elif (defined WARPX_DIM_3D) + + for (int k=0; k<=depos_order; k++) { + for (int j=0; j<=depos_order; j++) { + for (int i=0; i<=depos_order; i++) { + + // Jx + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k), + wq * invvol * invdt * onethird * sx_new[i] * sy_new[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k), + - wq * invvol * invdt * onethird * sx_old[i] * sy_new[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k), + wq * invvol * invdt * onesixth * sx_new[i] * sy_old[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_old + j,lo.z + k_new + k), + - wq * invvol * invdt * onesixth * sx_old[i] * sy_old[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k), + wq * invvol * invdt * onesixth * sx_new[i] * sy_new[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k), + - wq * invvol * invdt * onesixth * sx_old[i] * sy_new[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k), + wq * invvol * invdt * onethird * sx_new[i] * sy_old[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jx_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k), + - wq * invvol * invdt * onethird * sx_old[i] * sy_old[j] * sz_old[k]); + + // Jy + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k), + wq * invvol * invdt * onethird * sx_new[i] * sy_new[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k), + - wq * invvol * invdt * onethird * sx_new[i] * sy_old[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k), + wq * invvol * invdt * onesixth * sx_old[i] * sy_new[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k), + - wq * invvol * invdt * onesixth * sx_old[i] * sy_old[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k), + wq * invvol * invdt * onesixth * sx_new[i] * sy_new[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k), + - wq * invvol * invdt * onesixth * sx_new[i] * sy_old[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k), + wq * invvol * invdt * onethird * sx_old[i] * sy_new[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jy_arr(lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k), + - wq * invvol * invdt * onethird * sx_old[i] * sy_old[j] * sz_old[k]); + + // Jz + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_new + j, lo.z + k_new + k), + wq * invvol * invdt * onethird * sx_new[i] * sy_new[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_new + j, lo.z + k_old + k), + - wq * invvol * invdt * onethird * sx_new[i] * sy_new[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_new + j, lo.z + k_new + k), + wq * invvol * invdt * onesixth * sx_old[i] * sy_new[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_new + j, lo.z + k_old + k), + - wq * invvol * invdt * onesixth * sx_old[i] * sy_new[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_old + j, lo.z + k_new + k), + wq * invvol * invdt * onesixth * sx_new[i] * sy_old[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_new + i, lo.y + j_old + j, lo.z + k_old + k), + - wq * invvol * invdt * onesixth * sx_new[i] * sy_old[j] * sz_old[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_old + j, lo.z + k_new + k), + wq * invvol * invdt * onethird * sx_old[i] * sy_old[j] * sz_new[k]); + + amrex::Gpu::Atomic::Add(&jz_arr( lo.x + i_old + i, lo.y + j_old + j, lo.z + k_old + k), + - wq * invvol * invdt * onethird * sx_old[i] * sy_old[j] * sz_old[k]); + } + } + } +#endif + } ); +#endif // #if !(defined WARPX_DIM_RZ) +} #endif // CURRENTDEPOSITION_H_ diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 90772129b..a708ca8e1 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -336,6 +336,26 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } + } else if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) { + if (WarpX::nox == 1){ + doVayDepositionShapeN<1>( + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes ); + } else if (WarpX::nox == 2){ + doVayDepositionShapeN<2>( + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes ); + } else if (WarpX::nox == 3){ + doVayDepositionShapeN<3>( + GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes ); + } } else { if (WarpX::nox == 1){ doDepositionShapeN<1>( diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 87c444f79..06440729c 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -53,7 +53,8 @@ struct ParticlePusherAlgo { struct CurrentDepositionAlgo { enum { Esirkepov = 0, - Direct = 1 + Direct = 1, + Vay = 2 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 9649a20c9..7d8af6023 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -34,6 +34,7 @@ const std::map<std::string, int> particle_pusher_algo_to_int = { const std::map<std::string, int> current_deposition_algo_to_int = { {"esirkepov", CurrentDepositionAlgo::Esirkepov }, {"direct", CurrentDepositionAlgo::Direct }, + {"vay", CurrentDepositionAlgo::Vay }, #ifdef WARPX_USE_PSATD {"default", CurrentDepositionAlgo::Direct } #else diff --git a/Source/WarpX.H b/Source/WarpX.H index b8405db73..5337d898f 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -611,6 +611,15 @@ private: */ void CurrentCorrection (); + /** + * \brief Private function for Vay deposition in Fourier space + * (equations (20)-(24) of https://doi.org/10.1016/j.jcp.2013.03.010): + * loops over the MR levels and applies the correction on the fine and coarse + * patches (calls the virtual method \c VayDeposition of the spectral + * algorithm in use, via the public interface defined in the class SpectralSolver). + */ + void VayDeposition (); + void ReadParameters (); /** This function queries deprecated input parameters and abort |