diff options
Diffstat (limited to 'Source')
24 files changed, 162 insertions, 163 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 4ef68cfc9..b7c344954 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -385,13 +385,6 @@ WarpX::OneStep_nosub (Real cur_time) SyncCurrent(); SyncRho(); - // Apply current correction in Fourier space: for periodic single-box global FFTs - // without guard cells, apply this after calling SyncCurrent - if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { - if (fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay)) - VayDeposition(); - } - // At this point, J is up-to-date inside the domain, and E and B are // up-to-date including enough guard cells for first step of the field // solve. @@ -933,25 +926,3 @@ WarpX::applyMirrors(Real time){ } } } - -// Compute current from Vay deposition in Fourier space -void -WarpX::VayDeposition () -{ -#ifdef WARPX_USE_PSATD - if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) - { - for (int lev = 0; lev <= finest_level; ++lev) - { - spectral_solver_fp[lev]->VayDeposition(lev, current_fp[lev]); - if (spectral_solver_cp[lev]) spectral_solver_cp[lev]->VayDeposition(lev, current_cp[lev]); - } - } else { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false, - "WarpX::VayDeposition: only implemented for spectral solver."); - } -#else - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false, - "WarpX::CurrentCorrection: requires WarpX build with spectral solver support."); -#endif -} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index fde6bcb02..9139fa433 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -112,15 +112,9 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm * base class \c SpectralBaseAlgorithm and cannot be overridden by further * derived classes. * - * \param[in] lev The mesh-refinement level * \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 ( - const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldData& field_data) override final; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 52813cb13..35429109a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -832,25 +832,13 @@ void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data) } void -PsatdAlgorithm::VayDeposition ( - const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data) { // Profiling BL_PROFILE("PsatdAlgorithm::VayDeposition()"); const SpectralFieldIndex& Idx = m_spectral_index; - // 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(lev, *current[0], Idx.Jx, 0, IntVect(1)); - field_data.ForwardTransform(lev, *current[1], Idx.Jy, 0, IntVect(1)); - field_data.ForwardTransform(lev, *current[2], Idx.Jz, 0, IntVect(1)); - - const amrex::IntVect& fill_guards = m_fill_guards; - // Loop over boxes for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi) { @@ -903,11 +891,6 @@ PsatdAlgorithm::VayDeposition ( else fields(i,j,k,Idx.Jz) = 0._rt; }); } - - // Backward Fourier transform of J - field_data.BackwardTransform(lev, *current[0], Idx.Jx, 0, fill_guards); - field_data.BackwardTransform(lev, *current[1], Idx.Jy, 0, fill_guards); - field_data.BackwardTransform(lev, *current[2], Idx.Jz, 0, fill_guards); } #endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H index 13d2f1d5d..6f8497f8f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.H @@ -65,14 +65,9 @@ class PsatdAlgorithmComoving : public SpectralBaseAlgorithm * base class \c SpectralBaseAlgorithm and cannot be overridden by further * derived classes. * - * \param[in] lev The mesh-refinement level * \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 (const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldData& field_data) override final; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp index 77d03b500..9ed33c801 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp @@ -498,9 +498,7 @@ void PsatdAlgorithmComoving::CurrentCorrection (SpectralFieldData& field_data) } void -PsatdAlgorithmComoving::VayDeposition (const int /*lev*/, - SpectralFieldData& /*field_data*/, - std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/) +PsatdAlgorithmComoving::VayDeposition (SpectralFieldData& /*field_data*/) { amrex::Abort("Vay deposition not implemented for comoving PSATD"); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H index 11a1ea33a..d3c519c94 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.H @@ -45,13 +45,9 @@ class PsatdAlgorithmGalileanRZ : public SpectralBaseAlgorithmRZ * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further * derived classes. * - * \param[in] lev mesh-refinement level * \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 (const int lev, SpectralFieldDataRZ& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldDataRZ& field_data) override final; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp index 831d1e89f..81674f396 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp @@ -360,9 +360,7 @@ PsatdAlgorithmGalileanRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) } void -PsatdAlgorithmGalileanRZ::VayDeposition (const int /*lev*/, - SpectralFieldDataRZ& /*field_data*/, - std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/) +PsatdAlgorithmGalileanRZ::VayDeposition (SpectralFieldDataRZ& /*field_data*/) { amrex::Abort("Vay deposition not implemented in RZ geometry"); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H index 414394ba9..0689a7e49 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H @@ -111,15 +111,9 @@ class PsatdAlgorithmJLinearInTime : public SpectralBaseAlgorithm * base class \c SpectralBaseAlgorithm and cannot be overridden by further * derived classes. * - * \param[in] lev The mesh-refinement level * \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 ( - const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldData& field_data) override final; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp index d06087484..e18055020 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp @@ -433,15 +433,12 @@ void PsatdAlgorithmJLinearInTime::CurrentCorrection (SpectralFieldData& field_da } void -PsatdAlgorithmJLinearInTime::VayDeposition ( - const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +PsatdAlgorithmJLinearInTime::VayDeposition (SpectralFieldData& field_data) { // Profiling BL_PROFILE("PsatdAlgorithmJLinearInTime::VayDeposition()"); - amrex::ignore_unused(lev, field_data, current); + amrex::ignore_unused(field_data); amrex::Abort("Vay deposition not implemented for multi-J PSATD algorithm"); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H index 243760857..630c26260 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.H @@ -63,14 +63,9 @@ class PsatdAlgorithmPml : public SpectralBaseAlgorithm * base class \c SpectralBaseAlgorithm and cannot be overridden by further * derived classes. * - * \param[in] lev The mesh-refinement level * \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 (const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldData& field_data) override final; private: SpectralFieldIndex m_spectral_index; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp index 7b8deec30..246609905 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp @@ -408,9 +408,7 @@ PsatdAlgorithmPml::CurrentCorrection (SpectralFieldData& /*field_data*/) } void -PsatdAlgorithmPml::VayDeposition (const int /*lev*/, - SpectralFieldData& /*field_data*/, - std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/) +PsatdAlgorithmPml::VayDeposition (SpectralFieldData& /*field_data*/) { amrex::Abort("Vay deposition not implemented for PML PSATD"); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H index a9892cb40..6ea63c5a2 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.H @@ -46,12 +46,8 @@ class PsatdAlgorithmPmlRZ : public SpectralBaseAlgorithmRZ * 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 (const int lev, - SpectralFieldDataRZ& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldDataRZ& field_data) override final; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp index 522fe1881..3b8609c03 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp @@ -166,9 +166,7 @@ PsatdAlgorithmPmlRZ::CurrentCorrection (SpectralFieldDataRZ& /* field_data */) } void -PsatdAlgorithmPmlRZ::VayDeposition (const int /* lev */, - SpectralFieldDataRZ& /*field_data*/, - std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/) +PsatdAlgorithmPmlRZ::VayDeposition (SpectralFieldDataRZ& /*field_data*/) { amrex::Abort("Vay deposition not implemented in RZ geometry PML"); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index 9bbb695d2..608da5fd5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -49,14 +49,9 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further * derived classes. * - * \param[in] lev The mesh-refinement level * \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 (const int lev, - SpectralFieldDataRZ& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; + virtual void VayDeposition (SpectralFieldDataRZ& field_data) override final; private: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index 5bb409269..d2faecaa0 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -479,9 +479,7 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) } void -PsatdAlgorithmRZ::VayDeposition (const int /* lev */, - SpectralFieldDataRZ& /*field_data*/, - std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/) +PsatdAlgorithmRZ::VayDeposition (SpectralFieldDataRZ& /*field_data*/) { 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 b4d63ec81..b4f70d5f5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -57,14 +57,9 @@ class SpectralBaseAlgorithm * (<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] lev The mesh-refinement level * \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 (const int lev, - SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0; + virtual void VayDeposition (SpectralFieldData& field_data) = 0; /** * \brief Compute spectral divergence of E diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H index 9e7719a49..a2e293dc2 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H @@ -50,14 +50,9 @@ class SpectralBaseAlgorithmRZ * (<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] lev The mesh-refinement level * \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 (const int lev, - SpectralFieldDataRZ& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0; + virtual void VayDeposition (SpectralFieldDataRZ& field_data) = 0; protected: // Meant to be used in the subclasses diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 0953996fa..47ad47150 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -86,13 +86,55 @@ class SpectralSolver const bool divb_cleaning); /** - * \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 ForwardTransform( const int lev, + * \brief Transform the component i_comp of the MultiFab mf to Fourier space, + * and store the result internally (in the spectral field specified by field_index) + * + * \param[in] lev mesh refinement level + * \param[in] mf MultiFab that is transformed to Fourier space (component i_comp) + * \param[in] field_index index of the spectral field that stores the FFT result + * \param[in] i_comp component of the MultiFab mf that is transformed to Fourier space + * \param[in] stag index type that needs to be used to perform the FFT operation + */ + void ForwardTransform (const int lev, + const amrex::MultiFab& mf, + const int field_index, + const int i_comp, + const amrex::IntVect& stag); + + /** + * \brief Overload of ForwardTransform used if stag is not passed + */ + AMREX_FORCE_INLINE + void ForwardTransform (const int lev, const amrex::MultiFab& mf, const int field_index, - const int i_comp=0 ); + const int i_comp) + { + ForwardTransform(lev, mf, field_index, i_comp, mf.ixType().toIntVect()); + } + + /** + * \brief Overload of ForwardTransform used if i_comp is not passed + */ + AMREX_FORCE_INLINE + void ForwardTransform (const int lev, + const amrex::MultiFab& mf, + const int field_index, + const amrex::IntVect& stag) + { + ForwardTransform(lev, mf, field_index, 0, stag); + } + + /** + * Overload of ForwardTransform used if both i_comp and stag are not passed + */ + AMREX_FORCE_INLINE + void ForwardTransform (const int lev, + const amrex::MultiFab& mf, + const int field_index) + { + ForwardTransform(lev, mf, field_index, 0, mf.ixType().toIntVect()); + } /** * \brief Transform spectral field specified by `field_index` back to @@ -139,9 +181,9 @@ class SpectralSolver * \param[in,out] current Array of unique pointers to \c MultiFab storing * the three components of the current density */ - void VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current) + void VayDeposition () { - algorithm->VayDeposition(lev, field_data, current); + algorithm->VayDeposition(field_data); } /** diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 3a9e23127..0f7804002 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -92,10 +92,11 @@ void SpectralSolver::ForwardTransform( const int lev, const amrex::MultiFab& mf, const int field_index, - const int i_comp ) + const int i_comp, + const amrex::IntVect& stag ) { WARPX_PROFILE("SpectralSolver::ForwardTransform"); - field_data.ForwardTransform( lev, mf, field_index, i_comp ); + field_data.ForwardTransform( lev, mf, field_index, i_comp, stag ); } void diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 4868bd443..30c922512 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -106,12 +106,8 @@ class SpectralSolverRZ * 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] lev The mesh refinement level - * \param[in,out] current Array of unique pointers to \c MultiFab storing - * the three components of the current density */ - void VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current); + void VayDeposition (); /** * \brief Copy spectral data from component \c src_comp to component \c dest_comp diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index 0d515100a..3f4bddd8b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -157,7 +157,7 @@ SpectralSolverRZ::CurrentCorrection () } void -SpectralSolverRZ::VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current) +SpectralSolverRZ::VayDeposition () { - algorithm->VayDeposition(lev, field_data, current); + algorithm->VayDeposition(field_data); } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 82ca52c4c..a9267e4cd 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -55,8 +55,7 @@ using namespace amrex; #ifdef WARPX_USE_PSATD namespace { - void - ForwardTransformVect ( + void ForwardTransformVect ( const int lev, #ifdef WARPX_DIM_RZ SpectralSolverRZ& solver, @@ -64,15 +63,35 @@ namespace { SpectralSolver& solver, #endif std::array<std::unique_ptr<amrex::MultiFab>,3>& vector_field, - const int compx, const int compy, const int compz) + const int compx, const int compy, const int compz, + const amrex::IntVect& stag_x, const amrex::IntVect& stag_y, const amrex::IntVect& stag_z) { #ifdef WARPX_DIM_RZ + amrex::ignore_unused(stag_x, stag_y, stag_z); solver.ForwardTransform(lev, *vector_field[0], compx, *vector_field[1], compy); + solver.ForwardTransform(lev, *vector_field[2], compz); #else - solver.ForwardTransform(lev, *vector_field[0], compx); - solver.ForwardTransform(lev, *vector_field[1], compy); + solver.ForwardTransform(lev, *vector_field[0], compx, stag_x); + solver.ForwardTransform(lev, *vector_field[1], compy, stag_y); + solver.ForwardTransform(lev, *vector_field[2], compz, stag_z); #endif - solver.ForwardTransform(lev, *vector_field[2], compz); + } + + AMREX_FORCE_INLINE + void ForwardTransformVect ( + const int lev, +#ifdef WARPX_DIM_RZ + SpectralSolverRZ& solver, +#else + SpectralSolver& solver, +#endif + std::array<std::unique_ptr<amrex::MultiFab>,3>& vector_field, + const int compx, const int compy, const int compz) + { + ForwardTransformVect(lev, solver, vector_field, compx, compy, compz, + (*vector_field[0]).ixType().toIntVect(), + (*vector_field[1]).ixType().toIntVect(), + (*vector_field[2]).ixType().toIntVect()); } void @@ -243,11 +262,41 @@ WarpX::PSATDForwardTransformJ () for (int lev = 0; lev <= finest_level; ++lev) { - ForwardTransformVect(lev, *spectral_solver_fp[lev], current_fp[lev], idx_jx, idx_jy, idx_jz); + // With Vay's deposition, J stores a temporary current (D) that is later modified + // in Fourier space: its staggering matches that of rho and not J + amrex::IntVect jx_stag = + (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? + WarpX::m_rho_nodal_flag : current_fp[lev][0]->ixType().toIntVect(); + + amrex::IntVect jy_stag = + (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? + WarpX::m_rho_nodal_flag : current_fp[lev][1]->ixType().toIntVect(); + + amrex::IntVect jz_stag = + (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? + WarpX::m_rho_nodal_flag : current_fp[lev][2]->ixType().toIntVect(); + + ForwardTransformVect(lev, *spectral_solver_fp[lev], current_fp[lev], + idx_jx, idx_jy, idx_jz, jx_stag, jy_stag, jz_stag); if (spectral_solver_cp[lev]) { - ForwardTransformVect(lev, *spectral_solver_cp[lev], current_cp[lev], idx_jx, idx_jy, idx_jz); + // With Vay's deposition, J stores a temporary current (D) that is later modified + // in Fourier space: its staggering matches that of rho and not J + jx_stag = + (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? + WarpX::m_rho_nodal_flag : current_cp[lev][0]->ixType().toIntVect(); + + jy_stag = + (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? + WarpX::m_rho_nodal_flag : current_cp[lev][1]->ixType().toIntVect(); + + jz_stag = + (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? + WarpX::m_rho_nodal_flag : current_cp[lev][2]->ixType().toIntVect(); + + ForwardTransformVect(lev, *spectral_solver_cp[lev], current_cp[lev], + idx_jx, idx_jy, idx_jz, jx_stag, jy_stag, jz_stag); } } @@ -348,6 +397,19 @@ void WarpX::PSATDCurrentCorrection () } } +void WarpX::PSATDVayDeposition () +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + spectral_solver_fp[lev]->VayDeposition(); + + if (spectral_solver_cp[lev]) + { + spectral_solver_cp[lev]->VayDeposition(); + } + } +} + void WarpX::PSATDPushSpectralFields () { @@ -476,6 +538,14 @@ WarpX::PushPSATD () PSATDBackwardTransformJ(); } + // Compute the current in Fourier space according to the Vay deposition scheme, and + // transform back to real space so that the Vay deposition is reflected in the diagnostics + if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) + { + PSATDVayDeposition(); + PSATDBackwardTransformJ(); + } + #ifdef WARPX_DIM_RZ if (pml_rz[0]) pml_rz[0]->PushPSATD(0); #endif diff --git a/Source/WarpX.H b/Source/WarpX.H index 460c017d3..2fc5b89b7 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -336,6 +336,9 @@ public: //! If true, all fields are evaluated on a nodal grid and all MultiFabs have a nodal index type static int do_nodal; + // Global rho nodal flag to know about rho index type when rho MultiFab is not allocated + amrex::IntVect m_rho_nodal_flag; + std::array<const amrex::MultiFab* const, 3> get_array_Bfield_aux (const int lev) const { return { @@ -944,24 +947,6 @@ private: void AddRhoFromFineLevelandSumBoundary (int lev, int icomp, int ncomp); void NodalSyncRho (int lev, PatchType patch_type, int icomp, int ncomp); - /** - * \brief Private function for current correction in Fourier space - * (equation (19) 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 CurrentCorrection of the spectral - * algorithm in use, via the public interface defined in the class SpectralSolver). - */ - 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 @@ -1463,6 +1448,11 @@ private: void PSATDCurrentCorrection (); /** + * \brief Vay deposition in Fourier space (https://doi.org/10.1016/j.jcp.2013.03.010) + */ + void PSATDVayDeposition (); + + /** * \brief Update all necessary fields in spectral space */ void PSATDPushSpectralFields (); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 971b13548..eb1ee803e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -192,6 +192,7 @@ int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; int WarpX::do_nodal = false; +amrex::IntVect m_rho_nodal_flag; int WarpX::do_similar_dm_pml = 1; @@ -1650,6 +1651,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm ncomps = n_rz_azimuthal_modes*2 - 1; #endif + // Set global rho nodal flag to know about rho index type when rho MultiFab is not allocated + m_rho_nodal_flag = rho_nodal_flag; + // set human-readable tag for each MultiFab auto const tag = [lev]( std::string tagname ) { tagname.append("[l=").append(std::to_string(lev)).append("]"); |