diff options
21 files changed, 90 insertions, 109 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index c967064ba..176a2b6ed 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1563,7 +1563,7 @@ Numerics and algorithms https://ieeexplore.ieee.org/document/8659392. * ``warpx.do_multi_J`` (`0` or `1`; default: `0`) - Whether to use the multi-J algorithm, where current deposition and field update are performed multiple times within each time step. The number of sub-steps is determined by the input parameter ``warpx.do_multi_J_n_depositions``. Unlike sub-cycling, field gathering is performed only once per time step, as in regular PIC cycles. For simulations with strong numerical Cherenkov instability (NCI), it is recommended to use the multi-J algorithm in combination with ``psatd.do_time_averaging = 1``. + Whether to use the multi-J algorithm, where current deposition and field update are performed multiple times within each time step. The number of sub-steps is determined by the input parameter ``warpx.do_multi_J_n_depositions``. Unlike sub-cycling, field gathering is performed only once per time step, as in regular PIC cycles. When ``warpx.do_multi_J = 1``, we perform linear interpolation of two distinct currents deposited at the beginning and the end of the time step, instead of using one single current deposited at half time. For simulations with strong numerical Cherenkov instability (NCI), it is recommended to use the multi-J algorithm in combination with ``psatd.do_time_averaging = 1``. * ``warpx.do_multi_J_n_depositions`` (integer) Number of sub-steps to use with the multi-J algorithm, when ``warpx.do_multi_J = 1``. @@ -1697,9 +1697,6 @@ Numerics and algorithms * ``psatd.do_time_averaging`` (`0` or `1`; default: 0) Whether to use an averaged Galilean PSATD algorithm or standard Galilean PSATD. -* ``psatd.J_linear_in_time`` (`0` or `1`; default: `0`) - Whether to perform linear interpolation of two distinct currents deposited at the beginning and the end of the time step (``psatd.J_linear_in_time = 1``), instead of using one single current deposited at half time (``psatd.J_linear_in_time = 0``), for the field update in Fourier space. Currently requires ``psatd.update_with_rho = 1``, ``warpx.do_dive_cleaning = 1``, and ``warpx.do_divb_cleaning = 1``. - * ``warpx.override_sync_intervals`` (`string`) optional (default `1`) Using the `Intervals parser`_ syntax, this string defines the timesteps at which synchronization of sources (`rho` and `J`) and fields (`E` and `B`) on grid nodes at box diff --git a/Examples/Tests/multi_J/inputs_2d b/Examples/Tests/multi_J/inputs_2d index d85ff29f8..e09615503 100644 --- a/Examples/Tests/multi_J/inputs_2d +++ b/Examples/Tests/multi_J/inputs_2d @@ -39,7 +39,6 @@ warpx.moving_window_v = 1. # Spectral solver psatd.do_time_averaging = 1 -psatd.J_linear_in_time = 1 psatd.update_with_rho = 1 # Multi-J scheme diff --git a/Examples/Tests/multi_J/inputs_2d_pml b/Examples/Tests/multi_J/inputs_2d_pml index 9098d55cd..6cde7b9d1 100644 --- a/Examples/Tests/multi_J/inputs_2d_pml +++ b/Examples/Tests/multi_J/inputs_2d_pml @@ -39,7 +39,6 @@ warpx.moving_window_v = 1. # Spectral solver psatd.do_time_averaging = 0 -psatd.J_linear_in_time = 1 psatd.update_with_rho = 1 # Multi-J scheme diff --git a/Examples/Tests/multi_J/inputs_rz b/Examples/Tests/multi_J/inputs_rz index 704ed7fa7..061b1dcea 100644 --- a/Examples/Tests/multi_J/inputs_rz +++ b/Examples/Tests/multi_J/inputs_rz @@ -37,7 +37,6 @@ warpx.do_dive_cleaning = 1 warpx.do_divb_cleaning = 1 warpx.do_multi_J = 1 warpx.do_multi_J_n_depositions = 2 -psatd.J_linear_in_time = 1 psatd.do_time_averaging = 1 # PSATD diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 2783a6bf6..b39f29578 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -115,7 +115,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 J_linear_in_time, + const bool do_multi_J, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, int max_guard_EB, const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 837f66c2e..e6af648cc 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -477,7 +477,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 J_linear_in_time, + const bool do_multi_J, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, int max_guard_EB, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) @@ -669,7 +669,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, J_linear_in_time); + amrex::ignore_unused(lev, dt, do_multi_J); # if(AMREX_SPACEDIM!=3) amrex::ignore_unused(noy_fft); # endif @@ -690,7 +690,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, WarpX::fill_guards, v_galilean_zero, v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, J_linear_in_time, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, do_multi_J, m_dive_cleaning, m_divb_cleaning); #endif } @@ -808,7 +808,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, WarpX::fill_guards, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, J_linear_in_time, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, do_multi_J, m_dive_cleaning, m_divb_cleaning); #endif } } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index ad3ec8483..547f632e9 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -545,18 +545,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) PSATDForwardTransformRho(0, 1); } - // 4) Deposit J if needed - if (WarpX::J_linear_in_time) - { - // 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]); - // Filter, exchange boundary, and interpolate across levels - SyncCurrent(); - // Forward FFT of J - PSATDForwardTransformJ(); - } + // 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]); + // Filter, exchange boundary, and interpolate across levels + SyncCurrent(); + // Forward FFT of J + PSATDForwardTransformJ(); // Number of depositions for multi-J scheme const int n_depose = WarpX::do_multi_J_n_depositions; @@ -569,19 +565,13 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // Loop over multi-J depositions for (int i_depose = 0; i_depose < n_loop; i_depose++) { - // Move rho deposited previously, from new to old - PSATDMoveRhoNewToRhoOld(); - // Move J deposited previously, from new to old - // (when using assumption of J linear in time) - if (WarpX::J_linear_in_time) PSATDMoveJNewToJOld(); + PSATDMoveJNewToJOld(); - const amrex::Real t_depose = (WarpX::J_linear_in_time) ? - (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5)*sub_dt; + const amrex::Real t_depose = (i_depose-n_depose+1)*sub_dt; // Deposit new J at relative time t_depose 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], t_depose); // Filter, exchange boundary, and interpolate across levels SyncCurrent(); @@ -591,8 +581,11 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // Deposit new rho if (WarpX::update_with_rho) { - // Deposit rho at relative time (i_depose-n_depose+1)*sub_dt - mypc->DepositCharge(rho_fp, (i_depose-n_depose+1)*sub_dt); + // Move rho deposited previously, from new to old + PSATDMoveRhoNewToRhoOld(); + + // Deposit rho at relative time t_depose + mypc->DepositCharge(rho_fp, t_depose); // Filter, exchange boundary, and interpolate across levels SyncRho(); // Forward FFT of rho_new diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 41053fc9e..5549eda07 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -43,8 +43,9 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm * \param[in] dt time step of the simulation * \param[in] update_with_rho whether the update equation for E uses rho or not * \param[in] time_averaging whether to use time averaging for large time steps - * \param[in] J_linear_in_time whether to use two currents computed at the beginning and the end - * of the time interval (instead of using one current computed at half time) + * \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] 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 */ @@ -61,7 +62,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm const amrex::Real dt, const bool update_with_rho, const bool time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning); @@ -171,7 +172,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; - bool m_J_linear_in_time; + bool m_do_multi_J; bool m_dive_cleaning; bool m_divb_cleaning; bool m_is_galilean; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index fe9562dba..2b1e09c2d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -39,7 +39,7 @@ PsatdAlgorithm::PsatdAlgorithm( const amrex::Real dt, const bool update_with_rho, const bool time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning) // Initializer list @@ -59,7 +59,7 @@ PsatdAlgorithm::PsatdAlgorithm( m_dt(dt), m_update_with_rho(update_with_rho), m_time_averaging(time_averaging), - m_J_linear_in_time(J_linear_in_time), + m_do_multi_J(do_multi_J), m_dive_cleaning(dive_cleaning), m_divb_cleaning(divb_cleaning) { @@ -84,7 +84,7 @@ PsatdAlgorithm::PsatdAlgorithm( InitializeSpectralCoefficients(spectral_kspace, dm, dt); // Allocate these coefficients only with time averaging - if (time_averaging && !J_linear_in_time) + if (time_averaging && !do_multi_J) { Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); @@ -95,8 +95,8 @@ PsatdAlgorithm::PsatdAlgorithm( InitializeSpectralCoefficientsAveraging(spectral_kspace, dm, dt); } // Allocate these coefficients only with time averaging - // and with the assumption that J is linear in time - else if (time_averaging && J_linear_in_time) + // and with the assumption that J is linear in time (always with multi-J algorithm) + else if (time_averaging && do_multi_J) { X5_coef = SpectralComplexCoefficients(ba, dm, 1, 0); X6_coef = SpectralComplexCoefficients(ba, dm, 1, 0); @@ -123,7 +123,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const { const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; - const bool J_linear_in_time = m_J_linear_in_time; + const bool do_multi_J = m_do_multi_J; const bool dive_cleaning = m_dive_cleaning; const bool divb_cleaning = m_divb_cleaning; const bool is_galilean = m_is_galilean; @@ -163,7 +163,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const amrex::Array4<const Complex> Y3_arr; amrex::Array4<const Complex> Y4_arr; - if (time_averaging && !J_linear_in_time) + if (time_averaging && !do_multi_J) { Psi1_arr = Psi1_coef[mfi].array(); Psi2_arr = Psi2_coef[mfi].array(); @@ -175,7 +175,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const Array4<const Complex> X5_arr; Array4<const Complex> X6_arr; - if (time_averaging && J_linear_in_time) + if (time_averaging && do_multi_J) { X5_arr = X5_coef[mfi].array(); X6_arr = X6_coef[mfi].array(); @@ -320,7 +320,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B; } - if (J_linear_in_time) + if (do_multi_J) { const Complex Jx_new = fields(i,j,k,Idx.Jx_new); const Complex Jy_new = fields(i,j,k,Idx.Jy_new); @@ -391,7 +391,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const } // Additional update equations for averaged Galilean algorithm - if (time_averaging && !J_linear_in_time) + if (time_averaging && !do_multi_J) { // These coefficients are initialized in the function InitializeSpectralCoefficients below const Complex Psi1 = Psi1_arr(i,j,k); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index a23afbf4c..fba8f0881 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -23,7 +23,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ bool const nodal, amrex::Real const dt_step, bool const update_with_rho, const bool time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning); // Redefine functions from base class @@ -74,7 +74,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; - bool m_J_linear_in_time; + bool m_do_multi_J; 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 60a30e6f6..efc9337af 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -22,7 +22,7 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, bool const nodal, amrex::Real const dt, bool const update_with_rho, const bool time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning) // Initialize members of base class @@ -31,7 +31,7 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, m_dt(dt), m_update_with_rho(update_with_rho), m_time_averaging(time_averaging), - m_J_linear_in_time(J_linear_in_time), + m_do_multi_J(do_multi_J), m_dive_cleaning(dive_cleaning), m_divb_cleaning(divb_cleaning) { @@ -46,25 +46,25 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, coefficients_initialized = false; - if (time_averaging && J_linear_in_time) + if (time_averaging && do_multi_J) { X5_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); X6_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); } - if (time_averaging && !J_linear_in_time) + if (time_averaging && !do_multi_J) { - amrex::Abort("RZ PSATD: psatd.do_time_averaging = 1 implemented only with psatd.J_linear_in_time = 1"); + amrex::Abort("RZ PSATD: psatd.do_time_averaging = 1 implemented only with warpx.do_multi_J = 1"); } - if (dive_cleaning && !J_linear_in_time) + if (dive_cleaning && !do_multi_J) { - amrex::Abort("RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with psatd.J_linear_in_time = 1"); + amrex::Abort("RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with warpx.do_multi_J = 1"); } - if (divb_cleaning && !J_linear_in_time) + if (divb_cleaning && !do_multi_J) { - amrex::Abort("RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with psatd.J_linear_in_time = 1"); + amrex::Abort("RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with warpx.do_multi_J = 1"); } } @@ -76,7 +76,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; - const bool J_linear_in_time = m_J_linear_in_time; + const bool do_multi_J = m_do_multi_J; const bool dive_cleaning = m_dive_cleaning; const bool divb_cleaning = m_divb_cleaning; @@ -105,7 +105,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) amrex::Array4<const amrex::Real> X5_arr; amrex::Array4<const amrex::Real> X6_arr; - if (time_averaging && J_linear_in_time) + if (time_averaging && do_multi_J) { X5_arr = X5_coef[mfi].array(); X6_arr = X6_coef[mfi].array(); @@ -231,7 +231,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) G_old = fields(i,j,k,G_m); } - if (J_linear_in_time) + if (do_multi_J) { const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode; const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode; @@ -328,7 +328,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f) { const bool time_averaging = m_time_averaging; - const bool J_linear_in_time = m_J_linear_in_time; + const bool do_multi_J = m_do_multi_J; // Fill them with the right values: // Loop over boxes and allocate the corresponding coefficients @@ -349,7 +349,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const amrex::Array4<amrex::Real> X5; amrex::Array4<amrex::Real> X6; - if (time_averaging && J_linear_in_time) + if (time_averaging && do_multi_J) { X5 = X5_coef[mfi].array(); X6 = X6_coef[mfi].array(); @@ -388,7 +388,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0); } - if (time_averaging && J_linear_in_time) + if (time_averaging && do_multi_J) { 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 16f8e179c..9b748a048 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -39,21 +39,21 @@ class SpectralFieldIndex * Set integer indices to access data in spectral space * and total number of fields to be stored. * - * \param[in] update_with_rho whether rho is used in the field update equations - * \param[in] time_averaging whether the time averaging algorithm is used - * \param[in] J_linear_in_time whether to use two currents computed at the beginning and - * the end of the time interval (instead of using one current - * computed at half time) - * \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 - * div(B) = 0 law (new field G in the update equations) - * \param[in] pml whether the indices are used to access spectral data - * for the PML spectral solver + * \param[in] update_with_rho whether rho is used in the field update equations + * \param[in] 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] 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 + * div(B) = 0 law (new field G in the update equations) + * \param[in] pml whether the indices are used to access spectral data + * for the PML spectral solver */ SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning, const bool pml); @@ -86,7 +86,7 @@ class SpectralFieldIndex int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1; int Bx_avg = -1, By_avg = -1, Bz_avg = -1; - // J linear in time + // Multi-J, div(E) and div(B) cleaning int Jx_new = -1, Jy_new = -1, Jz_new = -1; int F = -1, G = -1; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 9d8563b6b..86fe8f7bd 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -34,7 +34,7 @@ using namespace amrex; SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning, const bool pml) @@ -66,9 +66,11 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, if (divb_cleaning) G = c++; - if (J_linear_in_time) + if (do_multi_J) { - Jx_new = c++; Jy_new = c++; Jz_new = c++; + Jx_new = c++; + Jy_new = c++; + Jz_new = c++; } } else // PML diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 5c3036838..f5cfb814b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -59,9 +59,9 @@ 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] J_linear_in_time whether to use two currents computed at the beginning and - * the end of the time interval (instead of using one current - * computed at half time) + * \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] 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 @@ -81,7 +81,7 @@ class SpectralSolver const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index a6b2f0aa4..70800f732 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -30,7 +30,7 @@ SpectralSolver::SpectralSolver( const bool pml, const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning) { @@ -42,7 +42,7 @@ SpectralSolver::SpectralSolver( const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - J_linear_in_time, dive_cleaning, divb_cleaning, pml); + do_multi_J, dive_cleaning, divb_cleaning, pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space @@ -63,7 +63,7 @@ SpectralSolver::SpectralSolver( else { algorithm = std::make_unique<PsatdAlgorithm>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, - v_galilean, dt, update_with_rho, fft_do_time_averaging, J_linear_in_time, + v_galilean, dt, update_with_rho, fft_do_time_averaging, do_multi_J, dive_cleaning, divb_cleaning); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index b97487401..827d226c2 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -36,7 +36,7 @@ class SpectralSolverRZ amrex::RealVect const dx, amrex::Real const dt, bool const update_with_rho, const bool fft_do_time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index b47dfa4ad..e187b9d27 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -34,7 +34,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, amrex::RealVect const dx, amrex::Real const dt, bool const update_with_rho, const bool fft_do_time_averaging, - const bool J_linear_in_time, + const bool do_multi_J, const bool dive_cleaning, const bool divb_cleaning) : k_space(realspace_ba, dm, dx) @@ -47,7 +47,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, const bool pml = false; m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - J_linear_in_time, dive_cleaning, divb_cleaning, pml); + do_multi_J, dive_cleaning, divb_cleaning, pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space @@ -56,7 +56,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, J_linear_in_time, dive_cleaning, divb_cleaning); + update_with_rho, fft_do_time_averaging, do_multi_J, dive_cleaning, divb_cleaning); } else { // Otherwise: use the Galilean algorithm algorithm = std::make_unique<GalileanPsatdAlgorithmRZ>( diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index ea112bbff..a096401bc 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -236,12 +236,9 @@ WarpX::PSATDForwardTransformJ () { const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; - const int idx_jx = (WarpX::J_linear_in_time) ? static_cast<int>(Idx.Jx_new) - : static_cast<int>(Idx.Jx); - const int idx_jy = (WarpX::J_linear_in_time) ? static_cast<int>(Idx.Jy_new) - : static_cast<int>(Idx.Jy); - const int idx_jz = (WarpX::J_linear_in_time) ? static_cast<int>(Idx.Jz_new) - : static_cast<int>(Idx.Jz); + const int idx_jx = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); + const int idx_jy = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); + const int idx_jz = (WarpX::do_multi_J) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); for (int lev = 0; lev <= finest_level; ++lev) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index d92345888..be810d21e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -243,7 +243,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, - J_linear_in_time, + do_multi_J, do_pml_dive_cleaning, do_pml_divb_cleaning, guard_cells.ng_FieldSolver.max(), do_pml_Lo_corrected, do_pml_Hi); @@ -274,7 +274,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, - J_linear_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, + do_multi_J, do_pml_dive_cleaning, do_pml_divb_cleaning, guard_cells.ng_FieldSolver.max(), do_pml_Lo_MR, do_pml_Hi_MR); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 59fe14eba..8160a428b 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -270,7 +270,6 @@ public: static int do_subcycling; static int do_multi_J; static int do_multi_J_n_depositions; - static int J_linear_in_time; static bool do_device_synchronize; static bool safe_guard_cells; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c30eb79d1..b2f6480a7 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -183,7 +183,6 @@ int WarpX::self_fields_verbosity = 2; int WarpX::do_subcycling = 0; int WarpX::do_multi_J = 0; int WarpX::do_multi_J_n_depositions; -int WarpX::J_linear_in_time = 0; bool WarpX::safe_guard_cells = 0; IntVect WarpX::filter_npass_each_dir(1); @@ -1088,7 +1087,6 @@ WarpX::ReadParameters () pp_psatd.query("current_correction", current_correction); pp_psatd.query("do_time_averaging", fft_do_time_averaging); - pp_psatd.query("J_linear_in_time", J_linear_in_time); if (!fft_periodic_single_box && current_correction) amrex::Abort( @@ -1200,12 +1198,9 @@ WarpX::ReadParameters () { amrex::Abort("Multi-J algorithm not implemented with Galilean PSATD"); } - } - if (J_linear_in_time) - { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(update_with_rho, - "psatd.update_with_rho must be set to 1 when psatd.J_linear_in_time = 1"); + "psatd.update_with_rho must be set to 1 when warpx.do_multi_J = 1"); } for (int dir = 0; dir < AMREX_SPACEDIM; dir++) @@ -2016,7 +2011,7 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector<std::unique_ptr<SpectralSo solver_dt, update_with_rho, fft_do_time_averaging, - J_linear_in_time, + do_multi_J, do_dive_cleaning, do_divb_cleaning); spectral_solver[lev] = std::move(pss); @@ -2072,7 +2067,7 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector<std::unique_ptr<SpectralSolv fft_periodic_single_box, update_with_rho, fft_do_time_averaging, - J_linear_in_time, + do_multi_J, do_dive_cleaning, do_divb_cleaning); spectral_solver[lev] = std::move(pss); |