Use RK4 to integrate the B-field in time in the hybrid-PIC algorithm (BLAST-WarpX#4461)

* add external current support to the hybrid-PIC solver

* add RZ support for `FiniteDifferenceSolver::CalculateCurrentAmpere`

* allow an initial Bz field to be set in RZ

* code cleanup and addition of CI test

* revert unwanted changes

* WIP transition to use RK4 in B-field time integration in hybrid-PIC algorithm

Co-authored-by: Avigdor Veksler <[email protected]>

* fix some clang-tidy issues

* Complete BfieldEvolveRK (#2)

* continued transition to RK4 in B-field time integration in hybrid-PIC algorithm

* RK-integration over all levels

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* update some of the CI tests

* more CI test updates

* Various code cleanups

* more updates to CI tests and checksum benchmarks

* more updates to CI tests and checksum benchmarks; remove commented code

* update documentation

* fix 1 for RZ CI test after merging of BLAST-WarpX#4464

* Avoid using `const` with `Real` passed by value

Co-authored-by: Axel Huebl <[email protected]>

* reduce default number of substeps to 10

* formatting fix

---------

Co-authored-by: Avigdor Veksler <[email protected]>
Co-authored-by: Avigdor Veksler <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Axel Huebl <[email protected]>

Author: 4 people

Docs/source/theory/kinetic_fluid_hybrid_model.rst

@@ -116,8 +116,7 @@ be interpolated to the correct time, using :math:`\vec{J}_i^n = 1/2(\vec{J}_i^{n
 The electron pressure is simply calculated using :math:`\rho^n` and the B-field is also already
 known at the correct time since it was calculated for :math:`t=t_n` at the end of the last step.
 Once :math:`\vec{E}^n` is calculated, it is used to push :math:`\vec{B}^n` forward in time
-(using the Maxwell-Faraday equation, i.e. the same as in the regular PIC routine with ``WarpX::EvolveB()``)
-to :math:`\vec{B}^{n+1/2}`.
+(using the Maxwell-Faraday equation) to :math:`\vec{B}^{n+1/2}`.
 
 Second half step
 """"""""""""""""
@@ -147,10 +146,11 @@ Sub-stepping
 ^^^^^^^^^^^^
 
 It is also well known that hybrid PIC routines require the B-field to be
-updated with a smaller timestep than needed for the particles. The update steps
-as outlined above are therefore wrapped in loops that enable the B-field to be
-sub-stepped. The exact number of sub-steps used can be specified by the user
-through a runtime simulation parameter (see :ref:`input parameters section <running-cpp-parameters-hybrid-model>`).
+updated with a smaller timestep than needed for the particles. A 4th order
+Runge-Kutta scheme is used to update the B-field. The RK scheme is repeated a
+number of times during each half-step outlined above. The number of sub-steps
+used can be specified by the user through a runtime simulation parameter
+(see :ref:`input parameters section <running-cpp-parameters-hybrid-model>`).
 
 .. _theory-hybrid-model-elec-temp:
 

Docs/source/usage/parameters.rst

@@ -2256,7 +2256,7 @@ Maxwell solver: kinetic-fluid hybrid
 * ``hybrid_pic_model.n_floor`` (`float`) optional (default ``1``)
     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the plasma density floor, in :math:`m^{-3}`, which is useful since the generalized Ohm's law used to calculate the E-field includes a :math:`1/n` term.
 
-* ``hybrid_pic_model.substeps`` (`int`) optional (default ``100``)
+* ``hybrid_pic_model.substeps`` (`int`) optional (default ``10``)
     If ``algo.maxwell_solver`` is set to ``hybrid``, this sets the number of sub-steps to take during the B-field update.
 
 .. note::

Examples/Tests/ohm_solver_EM_modes/PICMI_inputs.py

@@ -58,7 +58,7 @@ class EMModes(object):
     # Plasma resistivity - used to dampen the mode excitation
     eta = [[1e-7, 1e-7], [1e-7, 1e-5], [1e-7, 1e-4]]
     # Number of substeps used to update B
-    substeps = [[250, 500], [250, 750], [250, 1000]]
+    substeps = 20
 
     def __init__(self, test, dim, B_dir, verbose):
         """Get input parameters for the specific case desired."""
@@ -149,7 +149,6 @@ def get_simulation_parameters(self):
 
         self.NPPC = self.NPPC[self.dim-1]
         self.eta = self.eta[self.dim-1][idx]
-        self.substeps = self.substeps[self.dim-1][idx]
 
     def get_plasma_quantities(self):
         """Calculate various plasma parameters based on the simulation input."""

Examples/Tests/ohm_solver_EM_modes/PICMI_inputs_rz.py

@@ -53,7 +53,7 @@ class CylindricalNormalModes(object):
     # Plasma resistivity - used to dampen the mode excitation
     eta = 5e-4
     # Number of substeps used to update B
-    substeps = 250
+    substeps = 20
 
     def __init__(self, test, verbose):
         """Get input parameters for the specific case desired."""

Examples/Tests/ohm_solver_EM_modes/analysis_rz.py

@@ -154,7 +154,7 @@ def process(it):
     amps = np.abs(F_kw[2, 1, len(kz)//2-2:len(kz)//2+2])
     print("Amplitude sample: ", amps)
     assert np.allclose(
-        amps, np.array([61.4170941, 19.39380715, 101.08640009, 11.09261815])
+        amps, np.array([61.50945919, 19.74831134, 101.01820349, 10.8974811])
     )
 
 if sim.test:

Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py

@@ -56,7 +56,7 @@ class IonLandauDamping(object):
     # Plasma resistivity - used to dampen the mode excitation
     eta = 1e-7
     # Number of substeps used to update B
-    substeps = 100
+    substeps = 10
 
 
     def __init__(self, test, dim, m, T_ratio, verbose):

Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py

@@ -54,7 +54,7 @@ class HybridPICBeamInstability(object):
     # Plasma resistivity - used to dampen the mode excitation
     eta = 1e-7
     # Number of substeps used to update B
-    substeps = 400
+    substeps = 10
 
     # Beam parameters
     n_beam = [0.02, 0.1]
@@ -95,11 +95,9 @@ def __init__(self, test, dim, resonant, verbose):
         diag_period = 1 / 4.0 # Output interval (ion cyclotron periods)
         self.diag_steps = int(diag_period / self.DT)
 
-        # if this is a test case run for only 25 cyclotron periods and use
-        # fewer substeps to speed up the simulation
+        # if this is a test case run for only 25 cyclotron periods
         if self.test:
             self.LT = 25.0
-            self.substeps = 100
 
         self.total_steps = int(np.ceil(self.LT / self.DT))
 

Examples/Tests/ohm_solver_ion_beam_instability/analysis.py

@@ -190,11 +190,11 @@
     # assert's fail, the full benchmark should be rerun (same as the test but
     # without the `--test` argument) and the growth rates (up to saturation)
     # compared to the theoretical ones to determine if the physics test passes.
-    # At creation, the full test had the following errors when ran on 8 procs:
-    # m4_rms_error = 4.476; m5_rms_error = 9.211; m6_rms_error = 3.252
-    assert m4_rms_error < 1.55
-    assert m5_rms_error < 0.75
-    assert m6_rms_error < 0.40
+    # At creation, the full test (3d) had the following errors (ran on 1 V100):
+    # m4_rms_error = 3.329; m5_rms_error = 1.052; m6_rms_error = 2.583
+    assert np.isclose(m4_rms_error, 1.515, atol=0.01)
+    assert np.isclose(m5_rms_error, 0.718, atol=0.01)
+    assert np.isclose(m6_rms_error, 0.357, atol=0.01)
 
     # checksum check
     import os

Examples/Tests/ohm_solver_magnetic_reconnection/PICMI_inputs.py

@@ -59,7 +59,7 @@ class ForceFreeSheetReconnection(object):
     # Plasma resistivity - used to dampen the mode excitation
     eta = 6e-3  # normalized resistivity
     # Number of substeps used to update B
-    substeps = 750
+    substeps = 20
 
     def __init__(self, test, verbose):
 

Regression/Checksum/benchmarks_json/Python_ohms_law_solver_EM_modes_1d.json

@@ -1,17 +1,17 @@
 {
   "lev=0": {
-    "Bx": 0.0846232979356012,
-    "By": 0.0748927888009307,
+    "Bx": 0.08460920626188952,
+    "By": 0.07485331147996604,
     "Bz": 256.0,
-    "Ex": 1960.2420822749025,
-    "Ey": 2368.794327568534,
-    "Ez": 4873.496044339846
+    "Ex": 1954.3267519602405,
+    "Ey": 2363.4281756166347,
+    "Ez": 4873.508158589938
   },
   "ions": {
-    "particle_momentum_x": 1.615113301362171e-19,
-    "particle_momentum_y": 1.6152332353795267e-19,
-    "particle_momentum_z": 1.6134581585733078e-19,
-    "particle_position_x": 3678.484650863704,
+    "particle_momentum_x": 1.6151135948675135e-19,
+    "particle_momentum_y": 1.6152336151551518e-19,
+    "particle_momentum_z": 1.6134581268392543e-19,
+    "particle_position_x": 3678.484650899751,
     "particle_weight": 4.220251350277737e+21
   }
 }