[WIP] Refactor user interface of JRhom algorithm (previously multi-J)

Docs/source/usage/parameters.rst

@@ -2449,16 +2449,10 @@ Maxwell solver: PSATD method
 
     The default value for ``psatd.update_with_rho`` is ``1`` if ``psatd.v_galilean`` is non-zero and ``0`` otherwise.
     The option ``psatd.update_with_rho=0`` is not implemented with the following algorithms:
-    comoving PSATD (``psatd.v_comoving``), time averaging (``psatd.do_time_averaging=1``), div(E) cleaning (``warpx.do_dive_cleaning=1``), and multi-J (``warpx.do_multi_J=1``).
+    comoving PSATD (``psatd.v_comoving``), time averaging (``psatd.do_time_averaging=1``), div(E) cleaning (``warpx.do_dive_cleaning=1``), and JRhom (``psatd.JRhom``).
 
     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).
@@ -2476,12 +2470,8 @@ Maxwell solver: PSATD method
 * ``psatd.do_time_averaging`` (`0` or `1`; default: 0)
     Whether to use an averaged Galilean PSATD algorithm or standard Galilean PSATD.
 
-* ``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. 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``.
-    Note that this input parameter is not optional and must always be set in all input files where ``warpx.do_multi_J = 1``. No default value is provided automatically.
+* ``psatd.JRhom`` (str)
+    Whether to use the JRhom algorithm, where current deposition and fields update are performed multiple times within each time step. This is a string of the form ``CC1``, ``CL2``, etc. The first character of the string denotes the time dependency of :math:`J` (e.g., ``C`` stands for constant in time, ``L`` stands for linear in time). The second character of the string denotes the time dependency of :math:`\rho` (e.g., ``C`` stands for constant in time, ``L`` stands for linear in time). The integer after the first two characters denotes the number of subintervals for current deposition and fields update. Unlike subcycling, 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 JRhom algorithm in combination with ``psatd.do_time_averaging = 1``.
 
 Maxwell solver: macroscopic media
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Docs/source/usage/workflows/ml_materials/run_warpx_training.py

@@ -209,28 +209,24 @@ def get_laser(antenna_z, profile_t_peak, fill_in=True):
 
 # Electromagnetic solver
 
-psatd_algo = "multij"
+psatd_algo = "JRhom"
 if psatd_algo == "galilean":
     galilean_velocity = [0.0, 0.0] if dim == "3" else [0.0]
     galilean_velocity += [-c * beta_boost]
     n_pass_z = 1
-    do_multiJ = None
-    do_multi_J_n_depositions = None
-    J_in_time = None
+    JRhom = None
     current_correction = True
     divE_cleaning = False
-elif psatd_algo == "multij":
+elif psatd_algo == "JRhom":
     n_pass_z = 4
     galilean_velocity = None
-    do_multiJ = True
-    do_multi_J_n_depositions = 2
-    J_in_time = "linear"
+    JRhom = "LL2"
     current_correction = False
     divE_cleaning = True
 else:
     raise Exception(
         f"PSATD algorithm '{psatd_algo}' is not recognized!\n"
-        "Valid options are 'multiJ' or 'galilean'."
+        "Valid options are 'JRhom' or 'galilean'."
     )
 if dim == "rz":
     stencil_order = [8, 16]
@@ -252,7 +248,6 @@ def get_laser(antenna_z, profile_t_peak, fill_in=True):
     warpx_psatd_update_with_rho=True,
     warpx_current_correction=current_correction,
     divE_cleaning=divE_cleaning,
-    warpx_psatd_J_in_time=J_in_time,
 )
 
 # Diagnostics
@@ -329,8 +324,7 @@ def get_laser(antenna_z, profile_t_peak, fill_in=True):
     warpx_particle_pusher_algo="vay",
     warpx_amrex_the_arena_is_managed=False,
     warpx_amrex_use_gpu_aware_mpi=True,
-    warpx_do_multi_J=do_multiJ,
-    warpx_do_multi_J_n_depositions=do_multi_J_n_depositions,
+    warpx_JRhom=JRhom,
     warpx_grid_type=grid_type,
     # default: 2 for staggered grids, 8 for hybrid grids
     warpx_field_centering_order=[16, 16, 16],

Examples/Tests/langmuir/CMakeLists.txt

@@ -133,10 +133,10 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_2d_langmuir_multi_psatd_multiJ  # name
+        test_2d_langmuir_multi_psatd_JRhom  # name
         2  # dims
         2  # nprocs
-        inputs_test_2d_langmuir_multi_psatd_multiJ  # inputs
+        inputs_test_2d_langmuir_multi_psatd_JRhom  # inputs
         "analysis_2d.py diags/diag1000080"  # analysis
         "analysis_default_regression.py --path diags/diag1000080"  # checksum
         OFF  # dependency
@@ -145,10 +145,10 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_2d_langmuir_multi_psatd_multiJ_nodal  # name
+        test_2d_langmuir_multi_psatd_JRhom_nodal  # name
         2  # dims
         2  # nprocs
-        inputs_test_2d_langmuir_multi_psatd_multiJ_nodal  # inputs
+        inputs_test_2d_langmuir_multi_psatd_JRhom_nodal  # inputs
         "analysis_2d.py diags/diag1000080"  # analysis
         "analysis_default_regression.py --path diags/diag1000080"  # checksum
         OFF  # dependency
@@ -295,10 +295,10 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_3d_langmuir_multi_psatd_multiJ  # name
+        test_3d_langmuir_multi_psatd_JRhom  # name
         3  # dims
         2  # nprocs
-        inputs_test_3d_langmuir_multi_psatd_multiJ  # inputs
+        inputs_test_3d_langmuir_multi_psatd_JRhom  # inputs
         "analysis_3d.py diags/diag1000040"  # analysis
         "analysis_default_regression.py --path diags/diag1000040"  # checksum
         OFF  # dependency
@@ -307,10 +307,10 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_3d_langmuir_multi_psatd_multiJ_nodal  # name
+        test_3d_langmuir_multi_psatd_JRhom_nodal  # name
         3  # dims
         2  # nprocs
-        inputs_test_3d_langmuir_multi_psatd_multiJ_nodal  # inputs
+        inputs_test_3d_langmuir_multi_psatd_JRhom_nodal  # inputs
         "analysis_3d.py diags/diag1000040"  # analysis
         "analysis_default_regression.py --path diags/diag1000040"  # checksum
         OFF  # dependency
@@ -400,10 +400,10 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_rz_langmuir_multi_psatd_multiJ  # name
+        test_rz_langmuir_multi_psatd_JRhom  # name
         RZ  # dims
         2  # nprocs
-        inputs_test_rz_langmuir_multi_psatd_multiJ  # inputs
+        inputs_test_rz_langmuir_multi_psatd_JRhom  # inputs
         "analysis_rz.py diags/diag1000080"  # analysis
         "analysis_default_regression.py --path diags/diag1000080"  # checksum
         OFF  # dependency

Examples/Tests/langmuir/inputs_test_2d_langmuir_multi_psatd_multiJ → Examples/Tests/langmuir/inputs_test_2d_langmuir_multi_psatd_JRhom

@@ -3,10 +3,8 @@ FILE = inputs_base_2d
 
 # test input parameters
 algo.maxwell_solver = psatd
-psatd.J_in_time = linear
+psatd.JRhom = "LL2"
 psatd.solution_type = first-order
 psatd.update_with_rho = 1
 warpx.abort_on_warning_threshold = medium
 warpx.cfl = 0.7071067811865475
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 2

Examples/Tests/langmuir/inputs_test_2d_langmuir_multi_psatd_multiJ_nodal → Examples/Tests/langmuir/inputs_test_2d_langmuir_multi_psatd_JRhom_nodal

@@ -3,11 +3,9 @@ FILE = inputs_base_2d
 
 # test input parameters
 algo.maxwell_solver = psatd
-psatd.J_in_time = linear
+psatd.JRhom = "LL2"
 psatd.solution_type = first-order
 psatd.update_with_rho = 1
 warpx.abort_on_warning_threshold = medium
 warpx.cfl = 0.7071067811865475
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 2
 warpx.grid_type = collocated

Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_multiJ → Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_JRhom

@@ -4,10 +4,8 @@ FILE = inputs_base_3d
 # test input parameters
 algo.current_deposition = direct
 algo.maxwell_solver = psatd
-warpx.cfl = 0.5773502691896258
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 2
-psatd.J_in_time = linear
+psatd.JRhom = "LL2"
 psatd.solution_type = first-order
 psatd.update_with_rho = 1
 warpx.abort_on_warning_threshold = medium
+warpx.cfl = 0.5773502691896258

Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_multiJ_nodal → Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_JRhom_nodal

@@ -4,11 +4,9 @@ FILE = inputs_base_3d
 # test input parameters
 algo.current_deposition = direct
 algo.maxwell_solver = psatd
-psatd.J_in_time = linear
+psatd.JRhom = "LL2"
 psatd.solution_type = first-order
 psatd.update_with_rho = 1
 warpx.abort_on_warning_threshold = medium
 warpx.cfl = 0.5773502691896258
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 2
 warpx.grid_type = collocated

Examples/Tests/langmuir/inputs_test_rz_langmuir_multi_psatd_multiJ → Examples/Tests/langmuir/inputs_test_rz_langmuir_multi_psatd_JRhom

@@ -13,10 +13,9 @@ electrons.random_theta = 0
 ions.num_particles_per_cell_each_dim = 2 4 2
 ions.random_theta = 0
 psatd.current_correction = 0
+psatd.JRhom = "CL4"
 psatd.update_with_rho = 1
 warpx.abort_on_warning_threshold = medium
 warpx.do_dive_cleaning = 0
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 4
 warpx.n_rz_azimuthal_modes = 2
 warpx.use_filter = 1

Examples/Tests/nci_psatd_stability/CMakeLists.txt

@@ -147,11 +147,11 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_3d_uniform_plasma_multiJ  # name
+        test_3d_uniform_plasma_JRhom  # name
         3  # dims
         2  # nprocs
-        inputs_test_3d_uniform_plasma_multiJ  # inputs
-        "analysis_multiJ.py diags/diag1000300"  # analysis
+        inputs_test_3d_uniform_plasma_JRhom  # inputs
+        "analysis_JRhom.py diags/diag1000300"  # analysis
         "analysis_default_regression.py --path diags/diag1000300"  # checksum
         OFF  # dependency
     )
@@ -195,13 +195,13 @@ endif()
 
 if(WarpX_FFT)
     add_warpx_test(
-        test_rz_multiJ_psatd  # name
+        test_rz_JRhom_psatd  # name
         RZ  # dims
         2  # nprocs
-        inputs_test_rz_multiJ_psatd  # inputs
+        inputs_test_rz_JRhom_psatd  # inputs
         OFF  # analysis
         "analysis_default_regression.py --path diags/diag1000025"  # checksum
         OFF  # dependency
     )
-    label_warpx_test(test_rz_multiJ_psatd slow)
+    label_warpx_test(test_rz_JRhom_psatd slow)
 endif()

Examples/Tests/nci_psatd_stability/inputs_test_3d_uniform_plasma_multiJ → Examples/Tests/nci_psatd_stability/inputs_test_3d_uniform_plasma_JRhom

@@ -3,11 +3,8 @@ FILE = inputs_base_3d
 
 # test input parameters
 diag1.fields_to_plot = Bx By Bz divE Ex Ey Ez F G jx jy jz rho
-psatd.J_in_time = constant
-psatd.rho_in_time = constant
+psatd.JRhom = "CC1"
 psatd.solution_type = first-order
 warpx.abort_on_warning_threshold = medium
 warpx.do_divb_cleaning = 1
 warpx.do_dive_cleaning = 1
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 1

Examples/Tests/nci_psatd_stability/inputs_test_rz_multiJ_psatd → Examples/Tests/nci_psatd_stability/inputs_test_rz_JRhom_psatd

@@ -36,14 +36,12 @@ warpx.cfl = 1.
 
 warpx.do_dive_cleaning = 1
 warpx.do_divb_cleaning = 1
-warpx.do_multi_J = 1
-warpx.do_multi_J_n_depositions = 2
 psatd.do_time_averaging = 1
 
 # PSATD
+psatd.JRhom = "LL2"
 psatd.update_with_rho = 1
 #psatd.v_galilean = 0. 0. -0.9373391857121336
-psatd.J_in_time = linear
 
 # Particles
 

Python/pywarpx/picmi.py

@@ -1512,15 +1512,6 @@ class ElectromagneticSolver(picmistandard.PICMI_ElectromagneticSolver):
     warpx_psatd_do_time_averaging: bool, optional
         Whether to do the time averaging for the spectral solver
 
-    warpx_psatd_J_in_time: {'constant', '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.
-
-    warpx_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.
-
     warpx_do_pml_in_domain: bool, default=False
         Whether to do the PML boundaries within the domain (versus
         in the guard cells)
@@ -1549,8 +1540,6 @@ def init(self, kw):
             self.psatd_current_correction = kw.pop("warpx_current_correction", None)
             self.psatd_update_with_rho = kw.pop("warpx_psatd_update_with_rho", None)
             self.psatd_do_time_averaging = kw.pop("warpx_psatd_do_time_averaging", None)
-            self.psatd_J_in_time = kw.pop("warpx_psatd_J_in_time", None)
-            self.psatd_rho_in_time = kw.pop("warpx_psatd_rho_in_time", None)
 
         self.do_pml_in_domain = kw.pop("warpx_do_pml_in_domain", None)
         self.pml_has_particles = kw.pop("warpx_pml_has_particles", None)
@@ -1566,8 +1555,6 @@ def solver_initialize_inputs(self):
             pywarpx.psatd.current_correction = self.psatd_current_correction
             pywarpx.psatd.update_with_rho = self.psatd_update_with_rho
             pywarpx.psatd.do_time_averaging = self.psatd_do_time_averaging
-            pywarpx.psatd.J_in_time = self.psatd_J_in_time
-            pywarpx.psatd.rho_in_time = self.psatd_rho_in_time
 
             if self.grid.guard_cells is not None:
                 pywarpx.psatd.nx_guard = self.grid.guard_cells[0]
@@ -2706,14 +2693,10 @@ class Simulation(picmistandard.PICMI_Simulation):
     warpx_use_filter: bool, optional
         Whether to use filtering. The default depends on the conditions.
 
-    warpx_do_multi_J: bool, default=0
-        Whether to use the multi-J algorithm, where current deposition and
+    warpx_JRhom: string
+        Whether to use the JRhom algorithm, where current deposition and
         field update are performed multiple times within each time step.
-
-    warpx_do_multi_J_n_depositions: integer
-        Number of sub-steps to use with the multi-J algorithm, when ``warpx_do_multi_J=1``.
-        Note that this input parameter is not optional and must always be set in all
-        input files where ``warpx.do_multi_J=1``. No default value is provided automatically.
+        This is an empty string by default.
 
     warpx_grid_type: {'collocated', 'staggered', 'hybrid'}, default='staggered'
         Whether to use a collocated grid (all fields defined at the cell nodes),
@@ -2877,8 +2860,7 @@ def init(self, kw):
         self.field_gathering_algo = kw.pop("warpx_field_gathering_algo", None)
         self.particle_pusher_algo = kw.pop("warpx_particle_pusher_algo", None)
         self.use_filter = kw.pop("warpx_use_filter", None)
-        self.do_multi_J = kw.pop("warpx_do_multi_J", None)
-        self.do_multi_J_n_depositions = kw.pop("warpx_do_multi_J_n_depositions", None)
+        self.JRhom = kw.pop("warpx_JRhom", None)
         self.grid_type = kw.pop("warpx_grid_type", None)
         self.do_current_centering = kw.pop("warpx_do_current_centering", None)
         self.field_centering_order = kw.pop("warpx_field_centering_order", None)
@@ -2982,8 +2964,7 @@ def initialize_inputs(self):
         pywarpx.warpx.grid_type = self.grid_type
         pywarpx.warpx.do_current_centering = self.do_current_centering
         pywarpx.warpx.use_filter = self.use_filter
-        pywarpx.warpx.do_multi_J = self.do_multi_J
-        pywarpx.warpx.do_multi_J_n_depositions = self.do_multi_J_n_depositions
+        pywarpx.warpx.JRhom = self.JRhom
         pywarpx.warpx.serialize_initial_conditions = self.serialize_initial_conditions
         pywarpx.warpx.random_seed = self.random_seed
         pywarpx.warpx.used_inputs_file = self.used_inputs_file

Source/BoundaryConditions/PML.H

@@ -146,7 +146,7 @@ public:
          ablastr::utils::enums::GridType grid_type,
          int do_moving_window, int pml_has_particles, int do_pml_in_domain,
          PSATDSolutionType psatd_solution_type,
-         JInTime J_in_time, RhoInTime rho_in_time,
+         TimeDependencyJ time_dependency_J, TimeDependencyRho time_dependency_Rho,
          bool do_pml_dive_cleaning, bool do_pml_divb_cleaning,
          const amrex::IntVect& fill_guards_fields,
          const amrex::IntVect& fill_guards_current,

Source/BoundaryConditions/PML.cpp

@@ -556,7 +556,7 @@ PML::PML (const int lev, const BoxArray& grid_ba,
           ablastr::utils::enums::GridType grid_type,
           int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain,
           const PSATDSolutionType psatd_solution_type,
-          const JInTime J_in_time, const RhoInTime rho_in_time,
+          const TimeDependencyJ time_dependency_J, const TimeDependencyRho time_dependency_Rho,
           const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning,
           const amrex::IntVect& fill_guards_fields,
           const amrex::IntVect& fill_guards_current,
@@ -772,7 +772,7 @@ PML::PML (const int lev, const BoxArray& grid_ba,
 
     if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
 #ifndef WARPX_USE_FFT
-        amrex::ignore_unused(lev, dt, psatd_solution_type, J_in_time, rho_in_time);
+        amrex::ignore_unused(lev, dt, psatd_solution_type, time_dependency_J, time_dependency_Rho);
 #   if(AMREX_SPACEDIM!=3)
         amrex::ignore_unused(noy_fft);
 #   endif
@@ -793,7 +793,7 @@ PML::PML (const int lev, const BoxArray& grid_ba,
         spectral_solver_fp = std::make_unique<SpectralSolver>(lev, realspace_ba, dm,
             nox_fft, noy_fft, noz_fft, grid_type, v_galilean,
             v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho,
-            fft_do_time_averaging, psatd_solution_type, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning);
+            fft_do_time_averaging, psatd_solution_type, time_dependency_J, time_dependency_Rho, m_dive_cleaning, m_divb_cleaning);
 #endif
     }
 
@@ -905,7 +905,7 @@ PML::PML (const int lev, const BoxArray& grid_ba,
             spectral_solver_cp = std::make_unique<SpectralSolver>(lev, realspace_cba, cdm,
                 nox_fft, noy_fft, noz_fft, grid_type, v_galilean,
                 v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho,
-                fft_do_time_averaging, psatd_solution_type, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning);
+                fft_do_time_averaging, psatd_solution_type, time_dependency_J, time_dependency_Rho, m_dive_cleaning, m_divb_cleaning);
 #endif
         }
     }