[WIP] PSATD-JRhom with arbitrary time dependencies of J and rho

Docs/source/usage/parameters.rst

@@ -2389,7 +2389,7 @@ Maxwell solver: PSATD method
 
     where :math:`\theta=\exp(i\,\boldsymbol{k}\cdot\boldsymbol{v}_G\,\Delta{t}/2)`.
 
-    This option is currently implemented only for the standard PSATD, Galilean PSATD, and averaged Galilean PSATD schemes, while it is not yet available for the multi-J algorithm.
+    This option is currently implemented only for the standard PSATD, Galilean PSATD, and averaged Galilean PSATD schemes, while it is not yet available for the PSATD-JRhom algorithm.
 
 * ``psatd.update_with_rho`` (`0` or `1`)
     If true, the update equation for the electric field is expressed in terms of both the current density and the charge density, namely :math:`\widehat{\boldsymbol{J}}^{\,n+1/2}`, :math:`\widehat\rho^{n}`, and :math:`\widehat\rho^{n+1}`.
@@ -2458,7 +2458,7 @@ 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 PSATD-JRhom (``warpx.do_psatd_JRhom=1``).
 
     Note that the update with and without rho is also supported in RZ geometry.
 
@@ -2485,12 +2485,12 @@ 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_psatd_JRhom`` (`0` or `1`; default: `0`)
+    Whether to use the PSATD-JRhom 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_psatd_JRhom_n_depositions``. Unlike sub-cycling, field gathering is performed only once per time step, as in regular PIC cycles. When ``warpx.do_psatd_JRhom = 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 PSATD-JRhom 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.
+* ``warpx.do_psatd_JRhom_n_depositions`` (integer)
+    Number of sub-steps to use with the PSATD-JRhom algorithm, when ``warpx.do_psatd_JRhom = 1``.
+    Note that this input parameter is not optional and must always be set in all input files where ``warpx.do_psatd_JRhom = 1``. No default value is provided automatically.
 
 Maxwell solver: macroscopic media
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

@@ -209,28 +209,28 @@ def get_laser(antenna_z, profile_t_peak, fill_in=True):
 
 # Electromagnetic solver
 
-psatd_algo = "multij"
+psatd_algo = "psatd_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
+    do_psatd_JRhom = None
+    do_psatd_JRhom_n_depositions = None
     J_in_time = None
     current_correction = True
     divE_cleaning = False
-elif psatd_algo == "multij":
+elif psatd_algo == "psatd_JRhom":
     n_pass_z = 4
     galilean_velocity = None
-    do_multiJ = True
-    do_multi_J_n_depositions = 2
+    do_psatd_JRhom = True
+    do_psatd_JRhom_n_depositions = 2
     J_in_time = "linear"
     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 'psatd_JRhom' or 'galilean'."
     )
 if dim == "rz":
     stencil_order = [8, 16]
@@ -329,8 +329,8 @@ 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_do_psatd_JRhom=do_psatd_JRhom,
+    warpx_do_psatd_JRhom_n_depositions=do_psatd_JRhom_n_depositions,
     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_LL2  # name
         2  # dims
         2  # nprocs
-        inputs_test_2d_langmuir_multi_psatd_multiJ  # inputs
+        inputs_test_2d_langmuir_multi_psatd_LL2  # 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_LL2_nodal  # name
         2  # dims
         2  # nprocs
-        inputs_test_2d_langmuir_multi_psatd_multiJ_nodal  # inputs
+        inputs_test_2d_langmuir_multi_psatd_LL2_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_LL2  # name
         3  # dims
         2  # nprocs
-        inputs_test_3d_langmuir_multi_psatd_multiJ  # inputs
+        inputs_test_3d_langmuir_multi_psatd_LL2  # 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_LL2_nodal  # name
         3  # dims
         2  # nprocs
-        inputs_test_3d_langmuir_multi_psatd_multiJ_nodal  # inputs
+        inputs_test_3d_langmuir_multi_psatd_LL2_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_LL4  # name
         RZ  # dims
         2  # nprocs
-        inputs_test_rz_langmuir_multi_psatd_multiJ  # inputs
+        inputs_test_rz_langmuir_multi_psatd_LL4  # 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_LL2

@@ -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_LL2_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_LL2

@@ -4,10 +4,7 @@ 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.solution_type = first-order
+psatd.JRhom = "LL2"
 psatd.update_with_rho = 1
+warpx.cfl = 0.5773502691896258
 warpx.abort_on_warning_threshold = medium

Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_multiJ_nodal → Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_LL2_nodal

@@ -4,11 +4,8 @@ FILE = inputs_base_3d
 # test input parameters
 algo.current_deposition = direct
 algo.maxwell_solver = psatd
-psatd.J_in_time = linear
-psatd.solution_type = first-order
+psatd.JRhom = "LL2"
 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_LL4

@@ -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_psatd_CC1  # name
         3  # dims
         2  # nprocs
-        inputs_test_3d_uniform_plasma_multiJ  # inputs
-        "analysis_multiJ.py diags/diag1000300"  # analysis
+        inputs_test_3d_uniform_plasma_psatd_CC1  # inputs
+        "analysis_psatd_CC1.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_psatd_LL2  # name
         RZ  # dims
         2  # nprocs
-        inputs_test_rz_multiJ_psatd  # inputs
+        inputs_test_rz_psatd_LL2  # 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_test_rz_psatd_LL2 slow)
 endif()

Examples/Tests/nci_psatd_stability/analysis_multiJ.py → Examples/Tests/nci_psatd_stability/analysis_psatd_CC1.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 """
-This script is used to test the results of the multi-J PSATD
-first-order equations, with one J deposition. It compares the
+This script is used to test the results of the PSATD-JRhom
+first-order and second-order equations, with one J deposition. It compares the
 energy of the electric field with a given reference energy.
 
 The reference energy is computed by running the same test with J constant

Examples/Tests/nci_psatd_stability/inputs_test_3d_uniform_plasma_multiJ → Examples/Tests/nci_psatd_stability/inputs_test_3d_uniform_plasma_psatd_CC1

@@ -3,11 +3,7 @@ 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.solution_type = first-order
+psatd.JRhom = "CC1"
 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_psatd_LL2

@@ -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.update_with_rho = 1
 #psatd.v_galilean = 0. 0. -0.9373391857121336
-psatd.J_in_time = linear
+psatd.JRhom = "LL2"
 
 # Particles
 

Python/pywarpx/picmi.py

@@ -2799,14 +2799,14 @@ 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_do_psatd_JRhom: bool, default=0
+        Whether to use the PSATD-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``.
+    warpx_do_psatd_JRhom_n_depositions: integer
+        Number of sub-steps to use with the PSATD-JRhom algorithm, when ``warpx_do_psatd_JRhom=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.
+        input files where ``warpx.do_psatd_JRhom=1``. No default value is provided automatically.
 
     warpx_grid_type: {'collocated', 'staggered', 'hybrid'}, default='staggered'
         Whether to use a collocated grid (all fields defined at the cell nodes),
@@ -2975,8 +2975,10 @@ 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.do_psatd_JRhom = kw.pop("warpx_do_psatd_JRhom", None)
+        self.do_psatd_JRhom_n_depositions = kw.pop(
+            "warpx_do_psatd_JRhom_n_depositions", 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)
@@ -3082,8 +3084,8 @@ 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.do_psatd_JRhom = self.do_psatd_JRhom
+        pywarpx.warpx.do_psatd_JRhom_n_depositions = self.do_psatd_JRhom_n_depositions
         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

Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_CC2_nodal.json

@@ -0,0 +1,29 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 5.634812498911814,
+    "Bz": 0.0,
+    "Ex": 3753948071065.9404,
+    "Ey": 0.0,
+    "Ez": 3753948071065.936,
+    "jx": 1.008689565719933e+16,
+    "jy": 0.0,
+    "jz": 1.0086895657199332e+16
+  },
+  "positrons": {
+    "particle_momentum_x": 5.667364968811065e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 5.667364968811068e-20,
+    "particle_position_x": 0.6553599999190434,
+    "particle_position_y": 0.6553599999190433,
+    "particle_weight": 3200000000000000.5
+  },
+  "electrons": {
+    "particle_momentum_x": 5.667364968811043e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 5.667364968811046e-20,
+    "particle_position_x": 0.6553599999190443,
+    "particle_position_y": 0.6553599999190443,
+    "particle_weight": 3200000000000000.5
+  }
+}

Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_LL1.json

@@ -0,0 +1,29 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 5.718744708950617,
+    "Bz": 0.0,
+    "Ex": 3760487299466.9355,
+    "Ey": 0.0,
+    "Ez": 3760487299466.922,
+    "jx": 1.0110520336288102e+16,
+    "jy": 0.0,
+    "jz": 1.0110520336288108e+16
+  },
+  "positrons": {
+    "particle_momentum_x": 5.673851134639854e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 5.673851134639858e-20,
+    "particle_position_x": 0.65536,
+    "particle_position_y": 0.65536,
+    "particle_weight": 3200000000000000.5
+  },
+  "electrons": {
+    "particle_momentum_x": 5.673851134639966e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 5.673851134639968e-20,
+    "particle_position_x": 0.65536,
+    "particle_position_y": 0.6553599999999999,
+    "particle_weight": 3200000000000000.5
+  }
+}