diff options
23 files changed, 198 insertions, 97 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index bf6772ee2..f943e5d26 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1829,6 +1829,12 @@ Numerics and algorithms Note that the update with and without rho is also supported in RZ geometry. +* ``psatd.J_in_time`` (``constant`` or ``linear``; default ``constant``) + This determines whether the current density is assumed to be constant or linear in time, within the time step over which the electromagnetic fields are evolved. + +* ``psatd.rho_in_time`` (``linear``; default ``linear``) + This determines whether the charge density is assumed to be linear in time, within the time step over which the electromagnetic fields are evolved. + * ``psatd.v_galilean`` (`3 floats`, in units of the speed of light; default ``0. 0. 0.``) Defines the Galilean velocity. A non-zero velocity activates the Galilean algorithm, which suppresses numerical Cherenkov instabilities (NCI) in boosted-frame simulations (see the section :ref:`Numerical Stability and alternate formulation in a Galilean frame <theory-boostedframe-galilean>` for more information). diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index e274b9058..9d1ae29a3 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -414,7 +414,7 @@ analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_multiJ] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -433,7 +433,7 @@ analysisOutputImage = Langmuir_multi_psatd_multiJ.png [Langmuir_multi_psatd_multiJ_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -699,7 +699,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_multiJ] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON @@ -718,7 +718,7 @@ analysisOutputImage = Langmuir_multi_2d_psatd_multiJ.png [Langmuir_multi_2d_psatd_multiJ_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 dim = 2 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON @@ -2719,7 +2719,7 @@ analysisRoutine = Examples/Tests/galilean/analysis.py [multi_J_rz_psatd] buildDir = . inputFile = Examples/Tests/multi_J/inputs_rz -runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 warpx.abort_on_warning_threshold=medium +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 warpx.abort_on_warning_threshold=medium psatd.J_in_time=linear dim = 2 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 840f9d825..4688c19e5 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -131,7 +131,7 @@ public: int ncell, int delta, amrex::IntVect ref_ratio, amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, int do_moving_window, int pml_has_particles, int do_pml_in_domain, - const bool do_multi_J, + const int J_in_time, const int rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 9c3c17b0b..712e287f0 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -548,7 +548,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri int ncell, int delta, amrex::IntVect ref_ratio, Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain, - const bool do_multi_J, + const int J_in_time, const int rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, @@ -738,7 +738,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { #ifndef WARPX_USE_PSATD - amrex::ignore_unused(lev, dt, do_multi_J); + amrex::ignore_unused(lev, dt, J_in_time, rho_in_time); # if(AMREX_SPACEDIM!=3) amrex::ignore_unused(noy_fft); # endif @@ -759,7 +759,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri spectral_solver_fp = std::make_unique<SpectralSolver>(lev, realspace_ba, dm, nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, do_multi_J, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); #endif } @@ -878,7 +878,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri spectral_solver_cp = std::make_unique<SpectralSolver>(lev, realspace_cba, cdm, nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, do_multi_J, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); #endif } } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 1a2136021..202e91643 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -564,16 +564,19 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // 4) Deposit J at relative time -dt with time step dt // (dt[0] denotes the time step on mesh refinement level 0) - auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp; - mypc->DepositCurrent(current, dt[0], -dt[0]); - // Synchronize J: filter, exchange boundary, and interpolate across levels. - // With current centering, the nodal current is deposited in 'current', - // namely 'current_fp_nodal': SyncCurrent stores the result of its centering - // into 'current_fp' and then performs both filtering, if used, and exchange - // of guard cells. - SyncCurrent(current_fp, current_cp); - // Forward FFT of J - PSATDForwardTransformJ(current_fp, current_cp); + if (J_in_time == JInTime::Linear) + { + auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp; + mypc->DepositCurrent(current, dt[0], -dt[0]); + // Synchronize J: filter, exchange boundary, and interpolate across levels. + // With current centering, the nodal current is deposited in 'current', + // namely 'current_fp_nodal': SyncCurrent stores the result of its centering + // into 'current_fp' and then performs both filtering, if used, and exchange + // of guard cells. + SyncCurrent(current_fp, current_cp); + // Forward FFT of J + PSATDForwardTransformJ(current_fp, current_cp); + } // Number of depositions for multi-J scheme const int n_depose = WarpX::do_multi_J_n_depositions; @@ -587,13 +590,21 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) for (int i_depose = 0; i_depose < n_loop; i_depose++) { // Move J deposited previously, from new to old - PSATDMoveJNewToJOld(); + if (J_in_time == JInTime::Linear) + { + PSATDMoveJNewToJOld(); + } + + const amrex::Real t_depose_current = (J_in_time == JInTime::Linear) ? + (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5_rt)*sub_dt; - const amrex::Real t_depose = (i_depose-n_depose+1)*sub_dt; + // TODO Update this when rho quadratic in time is implemented + const amrex::Real t_depose_charge = (i_depose-n_depose+1)*sub_dt; - // Deposit new J at relative time t_depose with time step dt + // Deposit new J at relative time t_depose_current with time step dt // (dt[0] denotes the time step on mesh refinement level 0) - mypc->DepositCurrent(current, dt[0], t_depose); + auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp; + mypc->DepositCurrent(current, dt[0], t_depose_current); // Synchronize J: filter, exchange boundary, and interpolate across levels. // With current centering, the nodal current is deposited in 'current', // namely 'current_fp_nodal': SyncCurrent stores the result of its centering @@ -609,8 +620,8 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // Move rho deposited previously, from new to old PSATDMoveRhoNewToRhoOld(); - // Deposit rho at relative time t_depose - mypc->DepositCharge(rho_fp, t_depose); + // Deposit rho at relative time t_depose_charge + mypc->DepositCharge(rho_fp, t_depose_charge); // Filter, exchange boundary, and interpolate across levels SyncRho(); // Forward FFT of rho_new diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index eb88b1f4e..a370d4b2d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -1,6 +1,6 @@ target_sources(WarpX PRIVATE - PsatdAlgorithm.cpp + PsatdAlgorithmJConstantInTime.cpp PsatdAlgorithmJLinearInTime.cpp PsatdAlgorithmPml.cpp SpectralBaseAlgorithm.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index 7c2bf4281..c798ffb01 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,5 +1,5 @@ CEXE_sources += SpectralBaseAlgorithm.cpp -CEXE_sources += PsatdAlgorithm.cpp +CEXE_sources += PsatdAlgorithmJConstantInTime.cpp CEXE_sources += PsatdAlgorithmJLinearInTime.cpp CEXE_sources += PsatdAlgorithmPml.cpp CEXE_sources += PsatdAlgorithmComoving.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H index dd9c6a7fd..0d2d67434 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H @@ -4,8 +4,8 @@ * * License: BSD-3-Clause-LBNL */ -#ifndef WARPX_PSATD_ALGORITHM_H_ -#define WARPX_PSATD_ALGORITHM_H_ +#ifndef WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_ +#define WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_ #include "FieldSolver/SpectralSolver/SpectralFieldData.H" #include "FieldSolver/SpectralSolver/SpectralKSpace.H" @@ -24,12 +24,12 @@ /* \brief Class that updates the field in spectral space * and stores the coefficients of the corresponding update equation. */ -class PsatdAlgorithm : public SpectralBaseAlgorithm +class PsatdAlgorithmJConstantInTime : public SpectralBaseAlgorithm { public: /** - * \brief Constructor of the class PsatdAlgorithm + * \brief Constructor of the class PsatdAlgorithmJConstantInTime * * \param[in] spectral_kspace spectral space * \param[in] dm distribution mapping @@ -45,7 +45,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm * \param[in] dive_cleaning Update F as part of the field update, so that errors in divE=rho propagate away at the speed of light * \param[in] divb_cleaning Update G as part of the field update, so that errors in divB=0 propagate away at the speed of light */ - PsatdAlgorithm ( + PsatdAlgorithmJConstantInTime ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const SpectralFieldIndex& spectral_index, @@ -142,4 +142,4 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm bool m_is_galilean; }; #endif // WARPX_USE_PSATD -#endif // WARPX_PSATD_ALGORITHM_H_ +#endif // WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp index 1cbc27f0b..8971061f6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp @@ -4,7 +4,7 @@ * * License: BSD-3-Clause-LBNL */ -#include "PsatdAlgorithm.H" +#include "PsatdAlgorithmJConstantInTime.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" @@ -27,7 +27,7 @@ using namespace amrex; -PsatdAlgorithm::PsatdAlgorithm( +PsatdAlgorithmJConstantInTime::PsatdAlgorithmJConstantInTime( const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const SpectralFieldIndex& spectral_index, @@ -110,7 +110,7 @@ PsatdAlgorithm::PsatdAlgorithm( } void -PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const +PsatdAlgorithmJConstantInTime::pushSpectralFields (SpectralFieldData& f) const { const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; @@ -340,7 +340,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const } } -void PsatdAlgorithm::InitializeSpectralCoefficients ( +void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficients ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const amrex::Real dt) @@ -542,7 +542,7 @@ void PsatdAlgorithm::InitializeSpectralCoefficients ( } } -void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging ( +void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficientsAveraging ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const amrex::Real dt) @@ -733,10 +733,10 @@ void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging ( } } -void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data) +void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_data) { // Profiling - BL_PROFILE("PsatdAlgorithm::CurrentCorrection"); + BL_PROFILE("PsatdAlgorithmJConstantInTime::CurrentCorrection"); const SpectralFieldIndex& Idx = m_spectral_index; @@ -833,10 +833,10 @@ void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data) } void -PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data) +PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data) { // Profiling - BL_PROFILE("PsatdAlgorithm::VayDeposition()"); + BL_PROFILE("PsatdAlgorithmJConstantInTime::VayDeposition()"); const SpectralFieldIndex& Idx = m_spectral_index; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index 608da5fd5..097a1a9d6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -23,7 +23,8 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ bool const nodal, amrex::Real const dt_step, bool const update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning); // Redefine functions from base class @@ -62,7 +63,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; - bool m_do_multi_J; + int m_J_in_time; bool m_dive_cleaning; bool m_divb_cleaning; SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index fec0ebc8f..55b58821c 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -6,6 +6,7 @@ */ #include "PsatdAlgorithmRZ.H" #include "Utils/TextMsg.H" +#include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" @@ -23,7 +24,8 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, bool const nodal, amrex::Real const dt, bool const update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning) // Initialize members of base class @@ -32,10 +34,11 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, m_dt(dt), m_update_with_rho(update_with_rho), m_time_averaging(time_averaging), - m_do_multi_J(do_multi_J), + m_J_in_time(J_in_time), m_dive_cleaning(dive_cleaning), m_divb_cleaning(divb_cleaning) { + amrex::ignore_unused(rho_in_time); // Allocate the arrays of coefficients amrex::BoxArray const & ba = spectral_kspace.spectralspace_ba; @@ -47,28 +50,28 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, coefficients_initialized = false; - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time == JInTime::Linear) { X5_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); X6_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); } - if (time_averaging && !do_multi_J) + if (time_averaging && J_in_time != JInTime::Linear) { amrex::Abort(Utils::TextMsg::Err( - "RZ PSATD: psatd.do_time_averaging = 1 implemented only with warpx.do_multi_J = 1")); + "RZ PSATD: psatd.do_time_averaging=1 implemented only with psatd.J_in_time=linear")); } - if (dive_cleaning && !do_multi_J) + if (dive_cleaning && J_in_time != JInTime::Linear) { amrex::Abort(Utils::TextMsg::Err( - "RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with warpx.do_multi_J = 1")); + "RZ PSATD: warpx.do_dive_cleaning=1 implemented only with psatd.J_in_time=linear")); } - if (divb_cleaning && !do_multi_J) + if (divb_cleaning && J_in_time != JInTime::Linear) { amrex::Abort(Utils::TextMsg::Err( - "RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with warpx.do_multi_J = 1")); + "RZ PSATD: warpx.do_divb_cleaning=1 implemented only with psatd.J_in_time=linear")); } } @@ -80,7 +83,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; - const bool do_multi_J = m_do_multi_J; + const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false; const bool dive_cleaning = m_dive_cleaning; const bool divb_cleaning = m_divb_cleaning; @@ -109,7 +112,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) amrex::Array4<const amrex::Real> X5_arr; amrex::Array4<const amrex::Real> X6_arr; - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time_linear) { X5_arr = X5_coef[mfi].array(); X6_arr = X6_coef[mfi].array(); @@ -235,7 +238,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) G_old = fields(i,j,k,G_m); } - if (do_multi_J) + if (J_in_time_linear) { const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode; const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode; @@ -332,7 +335,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f) { const bool time_averaging = m_time_averaging; - const bool do_multi_J = m_do_multi_J; + const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false; // Fill them with the right values: // Loop over boxes and allocate the corresponding coefficients @@ -353,7 +356,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const amrex::Array4<amrex::Real> X5; amrex::Array4<amrex::Real> X6; - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time_linear) { X5 = X5_coef[mfi].array(); X6 = X6_coef[mfi].array(); @@ -392,7 +395,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0); } - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time_linear) { constexpr amrex::Real c2 = PhysConst::c; const amrex::Real dt3 = dt * dt * dt; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 811fdd4d7..4ab88f1a3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -55,7 +55,8 @@ class SpectralFieldIndex */ SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning, const bool pml, @@ -89,8 +90,13 @@ class SpectralFieldIndex int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1; int Bx_avg = -1, By_avg = -1, Bz_avg = -1; - // Multi-J, div(E) and div(B) cleaning + // J linear in time int Jx_new = -1, Jy_new = -1, Jz_new = -1; + + // rho quadratic in time + int rho_mid = -1; + + // div(E) and div(B) cleaning int F = -1, G = -1; // PML diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 00346a2c7..56189e0ff 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -35,7 +35,8 @@ using namespace amrex; SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning, const bool pml, @@ -68,13 +69,18 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, if (divb_cleaning) G = c++; - if (do_multi_J) + if (J_in_time == JInTime::Linear) { Jx_new = c++; Jy_new = c++; Jz_new = c++; } + if (rho_in_time == RhoInTime::Quadratic) + { + rho_mid = c++; + } + if (pml_rz) { Er_pml = c++; Et_pml = c++; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 684cf9586..38b242010 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -58,9 +58,10 @@ class SpectralSolver * (no domain decomposition) * \param[in] update_with_rho whether rho is used in the field update equations * \param[in] fft_do_time_averaging whether the time averaging algorithm is used - * \param[in] do_multi_J whether the multi-J algorithm is used (hence two currents - * computed at the beginning and the end of the time interval - * instead of one current computed at half time) + * \param[in] J_in_time integer that corresponds to the time dependency of J + * (constant, linear) for the PSATD algorithm + * \param[in] rho_in_time integer that corresponds to the time dependency of rho + * (linear, quadratic) for the PSATD algorithm * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in * Gauss law (new field F in the update equations) * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in @@ -79,7 +80,8 @@ class SpectralSolver const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 75c82319c..c0d7b8941 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -8,7 +8,7 @@ #include "FieldSolver/SpectralSolver/SpectralFieldData.H" #include "SpectralAlgorithms/PsatdAlgorithmComoving.H" #include "SpectralAlgorithms/PsatdAlgorithmPml.H" -#include "SpectralAlgorithms/PsatdAlgorithm.H" +#include "SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H" #include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H" #include "SpectralKSpace.H" #include "SpectralSolver.H" @@ -30,7 +30,8 @@ SpectralSolver::SpectralSolver( const bool pml, const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning) { @@ -41,8 +42,9 @@ SpectralSolver::SpectralSolver( // as well as the value of the corresponding k coordinates) const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); - m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - do_multi_J, dive_cleaning, divb_cleaning, pml); + m_spectral_index = SpectralFieldIndex( + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, + dive_cleaning, divb_cleaning, pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space @@ -64,18 +66,18 @@ SpectralSolver::SpectralSolver( } else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J { - if (do_multi_J) + if (J_in_time == JInTime::Constant) { - algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>( + algorithm = std::make_unique<PsatdAlgorithmJConstantInTime>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, - dt, fft_do_time_averaging, dive_cleaning, divb_cleaning); + v_galilean, dt, update_with_rho, fft_do_time_averaging, + dive_cleaning, divb_cleaning); } - else // standard, Galilean, averaged Galilean + else // J linear in time { - algorithm = std::make_unique<PsatdAlgorithm>( + algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, - v_galilean, dt, update_with_rho, fft_do_time_averaging, - dive_cleaning, divb_cleaning); + dt, fft_do_time_averaging, dive_cleaning, divb_cleaning); } } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 30c922512..45b55d9d4 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -37,7 +37,8 @@ class SpectralSolverRZ bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index 3f4bddd8b..23d93fe04 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -35,7 +35,8 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning) : k_space(realspace_ba, dm, dx) @@ -47,8 +48,9 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // as well as the value of the corresponding k coordinates. const bool is_pml = false; - m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - do_multi_J, dive_cleaning, divb_cleaning, is_pml, with_pml); + m_spectral_index = SpectralFieldIndex( + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, + dive_cleaning, divb_cleaning, is_pml, with_pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space @@ -60,7 +62,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // v_galilean is 0: use standard PSATD algorithm algorithm = std::make_unique<PsatdAlgorithmRZ>( k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt, - update_with_rho, fft_do_time_averaging, do_multi_J, dive_cleaning, divb_cleaning); + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, dive_cleaning, divb_cleaning); } else { // Otherwise: use the Galilean algorithm algorithm = std::make_unique<PsatdAlgorithmGalileanRZ>( diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 7499ee140..e8f20661a 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -280,9 +280,9 @@ void WarpX::PSATDForwardTransformJ ( { Idx = spectral_solver_fp[lev]->m_spectral_index; - idx_jx = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); - idx_jy = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); - idx_jz = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); + idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); ForwardTransformVect(lev, *spectral_solver_fp[lev], J_fp[lev], idx_jx, idx_jy, idx_jz); @@ -290,9 +290,9 @@ void WarpX::PSATDForwardTransformJ ( { Idx = spectral_solver_cp[lev]->m_spectral_index; - idx_jx = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); - idx_jy = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); - idx_jz = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); + idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); ForwardTransformVect(lev, *spectral_solver_cp[lev], J_cp[lev], idx_jx, idx_jy, idx_jz); } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 6f9178af7..b28c9bd59 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -497,7 +497,7 @@ WarpX::InitPML () pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(), dt[0], nox_fft, noy_fft, noz_fft, do_nodal, do_moving_window, pml_has_particles, do_pml_in_domain, - do_multi_J, + J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), @@ -534,7 +534,7 @@ WarpX::InitPML () pml_ncell, pml_delta, refRatio(lev-1), dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, do_moving_window, pml_has_particles, do_pml_in_domain, - do_multi_J, do_pml_dive_cleaning, do_pml_divb_cleaning, + J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), v_particle_pml, diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index ba752e012..97dd3b89f 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -83,6 +83,20 @@ struct GatheringAlgo { }; }; +struct JInTime { + enum { + Constant = 0, + Linear = 1 + }; +}; + +struct RhoInTime { + enum { + Linear = 1, + Quadratic = 2 + }; +}; + /** Strategy to compute weights for use in load balance. */ struct LoadBalanceCostsUpdateAlgo { diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index d0db13e12..088d5322f 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -62,6 +62,18 @@ const std::map<std::string, int> gathering_algo_to_int = { {"default", GatheringAlgo::EnergyConserving } }; +const std::map<std::string, int> J_in_time_to_int = { + {"constant", JInTime::Constant}, + {"linear", JInTime::Linear}, + {"default", JInTime::Constant} +}; + +const std::map<std::string, int> rho_in_time_to_int = { + {"linear", RhoInTime::Linear}, + {"quadratic", RhoInTime::Quadratic}, + {"default", RhoInTime::Linear} +}; + const std::map<std::string, int> load_balance_costs_update_algo_to_int = { {"timers", LoadBalanceCostsUpdateAlgo::Timers }, {"gpuclock", LoadBalanceCostsUpdateAlgo::GpuClock }, @@ -131,6 +143,10 @@ GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ algo_to_int = charge_deposition_algo_to_int; } else if (0 == std::strcmp(pp_search_key, "field_gathering")) { algo_to_int = gathering_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "J_in_time")) { + algo_to_int = J_in_time_to_int; + } else if (0 == std::strcmp(pp_search_key, "rho_in_time")) { + algo_to_int = rho_in_time_to_int; } else if (0 == std::strcmp(pp_search_key, "load_balance_costs_update")) { algo_to_int = load_balance_costs_update_algo_to_int; } else if (0 == std::strcmp(pp_search_key, "em_solver_medium")) { diff --git a/Source/WarpX.H b/Source/WarpX.H index a0f0f22f9..5c143a161 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -186,6 +186,10 @@ public: */ static amrex::Vector<ParticleBoundaryType> particle_boundary_hi; + //! Integers that correspond to the time dependency of J (constant, linear) + //! and rho (linear, quadratic) for the PSATD algorithm + static short J_in_time; + static short rho_in_time; //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid //! using finite centering of order given by #current_centering_nox, #current_centering_noy, diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1cdd2e3e0..aa6301a6e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -121,6 +121,8 @@ short WarpX::charge_deposition_algo; short WarpX::field_gathering_algo; short WarpX::particle_pusher_algo; short WarpX::maxwell_solver_id; +short WarpX::J_in_time; +short WarpX::rho_in_time; short WarpX::load_balance_costs_update_algo; bool WarpX::do_dive_cleaning = false; bool WarpX::do_divb_cleaning = false; @@ -1154,6 +1156,11 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE(noz_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used"); } + // Integers that correspond to the time dependency of J (constant, linear) + // and rho (linear, quadratic) for the PSATD algorithm + J_in_time = GetAlgorithmInteger(pp_psatd, "J_in_time"); + rho_in_time = GetAlgorithmInteger(pp_psatd, "rho_in_time"); + // Current correction activated by default, unless a charge-conserving // current deposition (Esirkepov, Vay) or the div(E) cleaning scheme // are used @@ -1308,10 +1315,28 @@ WarpX::ReadParameters () v_galilean_is_zero, "Multi-J algorithm not implemented with Galilean PSATD" ); + } - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(update_with_rho, - "psatd.update_with_rho must be set to 1 when warpx.do_multi_J = 1" - ); + if (J_in_time == JInTime::Constant) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + rho_in_time == RhoInTime::Linear, + "psatd.J_in_time=constant supports only psatd.rho_in_time=linear"); + } + + if (J_in_time == JInTime::Linear) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + update_with_rho, + "psatd.update_with_rho must be set to 1 when psatd.J_in_time=linear"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + v_galilean_is_zero, + "psatd.J_in_time=linear not implemented with Galilean PSATD"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + v_comoving_is_zero, + "psatd.J_in_time=linear not implemented with comoving PSATD"); } for (int dir = 0; dir < AMREX_SPACEDIM; dir++) @@ -2216,7 +2241,8 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector<std::unique_ptr<SpectralSo isAnyBoundaryPML(), update_with_rho, fft_do_time_averaging, - do_multi_J, + J_in_time, + rho_in_time, do_dive_cleaning, do_divb_cleaning); spectral_solver[lev] = std::move(pss); @@ -2271,7 +2297,8 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector<std::unique_ptr<SpectralSolv fft_periodic_single_box, update_with_rho, fft_do_time_averaging, - do_multi_J, + J_in_time, + rho_in_time, do_dive_cleaning, do_divb_cleaning); spectral_solver[lev] = std::move(pss); |