CFL-limited adaptive timestepping for electrostatic solver

Docs/source/usage/parameters.rst

@@ -2129,14 +2129,24 @@ Time step
     The ratio between the actual timestep that is used in the simulation
     and the Courant-Friedrichs-Lewy (CFL) limit. (e.g. for `warpx.cfl=1`,
     the timestep will be exactly equal to the CFL limit.)
-    This parameter will only be used with the electromagnetic solver.
+    For some speed v and grid spacing dx, this limits the timestep to `warpx.cfl * dx / v`.
+    When used with the electromagnetic solver, `v` is the speed of light.
+    For the electrostatic solver, `v` is the maximum speed among all particles in the domain.
 
 * ``warpx.const_dt`` (`float`)
     Allows direct specification of the time step size, in units of seconds.
-    When the electrostatic solver is being used, this must be supplied.
+    When the electrostatic solver is being used, this must be supplied if not using adaptive timestepping.
     This can be used with the electromagnetic solver, overriding ``warpx.cfl``, but
     it is up to the user to ensure that the CFL condition is met.
 
+* ``warpx.dt_update_interval`` (`string`) optional (default `-1`)
+    How many iterations pass between timestep adaptations when using the electrostatic solver.
+    Must be greater than `0` to use adaptive timestepping, or else ``warpx.const_dt`` must be specified.
+
+* ``warpx.max_dt`` (`float`) optional
+    The maximum timestep permitted for the electrostatic solver, when using adaptive timestepping.
+    If supplied, also sets the initial timestep for these simulations, before the first timestep update.
+
 Filtering
 ^^^^^^^^^
 
@@ -3448,39 +3458,42 @@ Reduced Diagnostics
         For 1D-Z, :math:`x`-related and :math:`y`-related quantities are not outputted.
         RZ geometry is not supported yet.
 
-* ``DifferentialLuminosity``
-    This type computes the differential luminosity between two species, defined as:
+    * ``DifferentialLuminosity``
+        This type computes the differential luminosity between two species, defined as:
 
-    .. math::
+        .. math::
+
+            \frac{d\mathcal{L}}{d\mathcal{E}^*}(\mathcal{E}^*, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\;
+             \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(\mathcal{E}^* - \mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2))
 
-        \frac{d\mathcal{L}}{d\mathcal{E}^*}(\mathcal{E}^*, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\;
-         \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(\mathcal{E}^* - \mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2))
+        where :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2(m_1 m_2 c^4
+        \gamma_1 \gamma_2 - \boldsymbol{p}_1\cdot\boldsymbol{p}_2 c^2)}` is the energy in the center-of-mass frame,
+        and :math:`f_i` is the distribution function of species :math:`i`. Note that, if :math:`\sigma^*(\mathcal{E}^*)`
+        is the center-of-mass cross-section of a given collision process, then
+        :math:`\int d\mathcal{E}^* \frac{d\mathcal{L}}{d\mathcal{E}^*} (\mathcal{E}^*, t)\sigma^*(\mathcal{E}^*)`
+        gives the total number of collisions of that process (from the beginning of the simulation up until time :math:`t`).
 
-    where :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2(m_1 m_2 c^4
-    \gamma_1 \gamma_2 - \boldsymbol{p}_1\cdot\boldsymbol{p}_2 c^2)}` is the energy in the center-of-mass frame,
-    and :math:`f_i` is the distribution function of species :math:`i`. Note that, if :math:`\sigma^*(\mathcal{E}^*)`
-    is the center-of-mass cross-section of a given collision process, then
-    :math:`\int d\mathcal{E}^* \frac{d\mathcal{L}}{d\mathcal{E}^*} (\mathcal{E}^*, t)\sigma^*(\mathcal{E}^*)`
-    gives the total number of collisions of that process (from the beginning of the simulation up until time :math:`t`).
+         The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-1}`. For collider-relevant WarpX simulations
+        involving two crossing, high-energy beams of particles, the differential luminosity in :math:`\text{s}^{-1}.\text{m}^{-2}.\text{eV}^{-1}`
+        can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams.
 
-    The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-1}`. For collider-relevant WarpX simulations
-    involving two crossing, high-energy beams of particles, the differential luminosity in :math:`\text{s}^{-1}.\text{m}^{-2}.\text{eV}^{-1}`
-    can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams.
+        In practice, the above expression of the differential luminosity is evaluated over discrete bins in energy :math:`\mathcal{E}^*`,
+        and by summing over macroparticles.
 
-    In practice, the above expression of the differential luminosity is evaluated over discrete bins in energy :math:`\mathcal{E}^*`,
-    and by summing over macroparticles.
+        * ``<reduced_diags_name>.species`` (`list of two strings`)
+            The names of the two species for which the differential luminosity is computed.
 
-    * ``<reduced_diags_name>.species`` (`list of two strings`)
-        The names of the two species for which the differential luminosity is computed.
+        * ``<reduced_diags_name>.bin_number`` (`int` > 0)
+            The number of bins in energy :math:`\mathcal{E}^*`
 
-    * ``<reduced_diags_name>.bin_number`` (`int` > 0)
-        The number of bins in energy :math:`\mathcal{E}^*`
+        * ``<reduced_diags_name>.bin_max`` (`float`, in eV)
+            The minimum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
 
-    * ``<reduced_diags_name>.bin_max`` (`float`, in eV)
-        The minimum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
+        * ``<reduced_diags_name>.bin_min`` (`float`, in eV)
+            The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
 
-    * ``<reduced_diags_name>.bin_min`` (`float`, in eV)
-        The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
+    * ``Timestep``
+        This type outputs the simulation's physical timestep (in seconds) at each mesh refinement level.
 
 * ``<reduced_diags_name>.intervals`` (`string`)
     Using the `Intervals Parser`_ syntax, this string defines the timesteps at which reduced

Examples/Tests/electrostatic_sphere/CMakeLists.txt

@@ -41,6 +41,16 @@ add_warpx_test(
     OFF  # dependency
 )
 
+add_warpx_test(
+    test_3d_electrostatic_sphere_adaptive  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_electrostatic_sphere_adaptive  # inputs
+    analysis_electrostatic_sphere.py  # analysis
+    diags/diag1000054  # output
+    OFF  # dependency
+)
+
 add_warpx_test(
     test_rz_electrostatic_sphere  # name
     RZ  # dims

Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive

@@ -0,0 +1,47 @@
+stop_time = 60e-6
+warpx.cfl = 0.2
+warpx.dt_update_interval = 10
+warpx.max_dt = 1.5e-6
+amr.n_cell = 64 64 64
+amr.max_level = 0
+amr.blocking_factor = 8
+amr.max_grid_size = 128
+geometry.dims = 3
+geometry.prob_lo = -0.5 -0.5 -0.5
+geometry.prob_hi =  0.5  0.5  0.5
+boundary.field_lo = pec pec pec
+boundary.field_hi = pec pec pec
+warpx.do_electrostatic = relativistic
+
+particles.species_names = electron
+
+algo.field_gathering = momentum-conserving
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+my_constants.n0 = 1.49e6
+my_constants.R0 = 0.1
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "NUniformPerCell"
+electron.num_particles_per_cell_each_dim = 2 2 2
+electron.profile = parse_density_function
+electron.density_function(x,y,z) = "(x*x + y*y + z*z < R0*R0)*n0"
+electron.momentum_distribution_type = at_rest
+
+diagnostics.diags_names = diag1 diag2
+
+diag1.intervals = 30
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez rho
+
+diag2.intervals = 30
+diag2.diag_type = Full
+diag2.fields_to_plot = none
+diag2.format = openpmd
+
+warpx.reduced_diags_names = timestep
+timestep.intervals = 1
+timestep.type = Timestep

Python/pywarpx/picmi.py

@@ -1522,8 +1522,7 @@ def solver_initialize_inputs(self):
         # --- Same method names are used, though mapped to lower case.
         pywarpx.algo.maxwell_solver = self.method
 
-        if self.cfl is not None:
-            pywarpx.warpx.cfl = self.cfl
+        pywarpx.warpx.cfl = self.cfl
 
         if self.source_smoother is not None:
             self.source_smoother.smoother_initialize_inputs(self)
@@ -1880,20 +1879,38 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver):
 
     warpx_self_fields_verbosity: integer, default=2
         Level of verbosity for the lab frame solver
+
+    warpx_dt_update_interval: string, optional (default = -1)
+        How frequently the timestep is updated. Adaptive timestepping is disabled when this is <= 0.
+
+    warpx_cfl: float, optional
+        Fraction of the CFL condition for particle velocity vs grid size, used to set the timestep when `dt_update_interval > 0`.
+
+    warpx_max_dt: float, optional
+        The maximum allowable timestep when `dt_update_interval > 0`.
+
     """
 
     def init(self, kw):
         self.relativistic = kw.pop("warpx_relativistic", False)
         self.absolute_tolerance = kw.pop("warpx_absolute_tolerance", None)
         self.self_fields_verbosity = kw.pop("warpx_self_fields_verbosity", None)
         self.magnetostatic = kw.pop("warpx_magnetostatic", False)
+        self.cfl = kw.pop("warpx_cfl", None)
+        self.dt_update_interval = kw.pop("dt_update_interval", None)
+        self.max_dt = kw.pop("warpx_max_dt", None)
 
     def solver_initialize_inputs(self):
         # Open BC means FieldBoundaryType::Open for electrostatic sims, rather than perfectly-matched layer
         BC_map["open"] = "open"
 
         self.grid.grid_initialize_inputs()
 
+        # set adaptive timestepping parameters
+        pywarpx.warpx.cfl = self.cfl
+        pywarpx.warpx.dt_update_interval = self.dt_update_interval
+        pywarpx.warpx.max_dt = self.max_dt
+
         if self.relativistic:
             pywarpx.warpx.do_electrostatic = "relativistic"
         else:
@@ -3890,6 +3907,7 @@ def __init__(
             "ParticleNumber",
             "LoadBalanceCosts",
             "LoadBalanceEfficiency",
+            "Timestep",
         ]
         # The species diagnostics require a species to be provided
         self._species_reduced_diagnostics = [

Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json

@@ -0,0 +1,17 @@
+{
+  "lev=0": {
+    "Ex": 5.177444767224255,
+    "Ey": 5.177444767224254,
+    "Ez": 5.177444767224256,
+    "rho": 2.6092568008333797e-10
+  },
+  "electron": {
+    "particle_momentum_x": 1.3215019655285216e-23,
+    "particle_momentum_y": 1.3215019655285214e-23,
+    "particle_momentum_z": 1.3215019655285217e-23,
+    "particle_position_x": 912.2310003741203,
+    "particle_position_y": 912.2310003741203,
+    "particle_position_z": 912.2310003741202,
+    "particle_weight": 6212.501525878906
+  }
+}

Source/Diagnostics/ReducedDiags/CMakeLists.txt

@@ -3,26 +3,27 @@ foreach(D IN LISTS WarpX_DIMS)
     target_sources(lib_${SD}
       PRIVATE
         BeamRelevant.cpp
+        ChargeOnEB.cpp
         ColliderRelevant.cpp
         DifferentialLuminosity.cpp
         FieldEnergy.cpp
+        FieldMaximum.cpp
+        FieldMomentum.cpp
         FieldProbe.cpp
         FieldProbeParticleContainer.cpp
-        FieldMomentum.cpp
+        FieldReduction.cpp
+        FieldProbe.cpp
         LoadBalanceCosts.cpp
         LoadBalanceEfficiency.cpp
         MultiReducedDiags.cpp
         ParticleEnergy.cpp
-        ParticleMomentum.cpp
+        ParticleExtrema.cpp
         ParticleHistogram.cpp
         ParticleHistogram2D.cpp
+        ParticleMomentum.cpp
+        ParticleNumber.cpp
         ReducedDiags.cpp
-        FieldMaximum.cpp
-        ParticleExtrema.cpp
         RhoMaximum.cpp
-        ParticleNumber.cpp
-        FieldReduction.cpp
-        FieldProbe.cpp
-        ChargeOnEB.cpp
+        Timestep.cpp
     )
 endforeach()

Source/Diagnostics/ReducedDiags/Make.package

@@ -1,24 +1,24 @@
 CEXE_sources += MultiReducedDiags.cpp
 CEXE_sources += ReducedDiags.cpp
-CEXE_sources += ParticleEnergy.cpp
-CEXE_sources += ParticleMomentum.cpp
-CEXE_sources += FieldEnergy.cpp
-CEXE_sources += FieldProbe.cpp
-CEXE_sources += FieldProbeParticleContainer.cpp
-CEXE_sources += FieldMomentum.cpp
 CEXE_sources += BeamRelevant.cpp
+CEXE_sources += ChargeOnEB.cpp
 CEXE_sources += ColliderRelevant.cpp
 CEXE_sources += DifferentialLuminosity.cpp
+CEXE_sources += FieldEnergy.cpp
+CEXE_sources += FieldMaximum.cpp
+CEXE_sources += FieldMomentum.cpp
+CEXE_sources += FieldProbe.cpp
+CEXE_sources += FieldProbeParticleContainer.cpp
+CEXE_sources += FieldReduction.cpp
 CEXE_sources += LoadBalanceCosts.cpp
 CEXE_sources += LoadBalanceEfficiency.cpp
+CEXE_sources += ParticleEnergy.cpp
+CEXE_sources += ParticleExtrema.cpp
 CEXE_sources += ParticleHistogram.cpp
 CEXE_sources += ParticleHistogram2D.cpp
-CEXE_sources += FieldMaximum.cpp
-CEXE_sources += FieldProbe.cpp
-CEXE_sources += ParticleExtrema.cpp
-CEXE_sources += RhoMaximum.cpp
+CEXE_sources += ParticleMomentum.cpp
 CEXE_sources += ParticleNumber.cpp
-CEXE_sources += FieldReduction.cpp
-CEXE_sources += ChargeOnEB.cpp
+CEXE_sources += RhoMaximum.cpp
+CEXE_sources += Timestep.cpp
 
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags

Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp

@@ -12,8 +12,8 @@
 #include "DifferentialLuminosity.H"
 #include "FieldEnergy.H"
 #include "FieldMaximum.H"
-#include "FieldProbe.H"
 #include "FieldMomentum.H"
+#include "FieldProbe.H"
 #include "FieldReduction.H"
 #include "LoadBalanceCosts.H"
 #include "LoadBalanceEfficiency.H"
@@ -24,6 +24,7 @@
 #include "ParticleMomentum.H"
 #include "ParticleNumber.H"
 #include "RhoMaximum.H"
+#include "Timestep.H"
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXProfilerWrapper.H"
 
@@ -52,24 +53,25 @@ MultiReducedDiags::MultiReducedDiags ()
     using CS = const std::string& ;
     const auto reduced_diags_dictionary =
         std::map<std::string, std::function<std::unique_ptr<ReducedDiags>(CS)>>{
+            {"BeamRelevant",          [](CS s){return std::make_unique<BeamRelevant>(s);}},
+            {"ChargeOnEB",            [](CS s){return std::make_unique<ChargeOnEB>(s);}},
+            {"ColliderRelevant",      [](CS s){return std::make_unique<ColliderRelevant>(s);}},
+            {"DifferentialLuminosity",[](CS s){return std::make_unique<DifferentialLuminosity>(s);}},
             {"ParticleEnergy",        [](CS s){return std::make_unique<ParticleEnergy>(s);}},
+            {"ParticleExtrema",       [](CS s){return std::make_unique<ParticleExtrema>(s);}},
+            {"ParticleHistogram",     [](CS s){return std::make_unique<ParticleHistogram>(s);}},
+            {"ParticleHistogram2D",   [](CS s){return std::make_unique<ParticleHistogram2D>(s);}},
             {"ParticleMomentum",      [](CS s){return std::make_unique<ParticleMomentum>(s);}},
+            {"ParticleNumber",        [](CS s){return std::make_unique<ParticleNumber>(s);}},
             {"FieldEnergy",           [](CS s){return std::make_unique<FieldEnergy>(s);}},
-            {"FieldMomentum",         [](CS s){return std::make_unique<FieldMomentum>(s);}},
             {"FieldMaximum",          [](CS s){return std::make_unique<FieldMaximum>(s);}},
+            {"FieldMomentum",         [](CS s){return std::make_unique<FieldMomentum>(s);}},
             {"FieldProbe",            [](CS s){return std::make_unique<FieldProbe>(s);}},
             {"FieldReduction",        [](CS s){return std::make_unique<FieldReduction>(s);}},
-            {"RhoMaximum",            [](CS s){return std::make_unique<RhoMaximum>(s);}},
-            {"BeamRelevant",          [](CS s){return std::make_unique<BeamRelevant>(s);}},
-            {"ColliderRelevant",      [](CS s){return std::make_unique<ColliderRelevant>(s);}},
-            {"DifferentialLuminosity",[](CS s){return std::make_unique<DifferentialLuminosity>(s);}},
             {"LoadBalanceCosts",      [](CS s){return std::make_unique<LoadBalanceCosts>(s);}},
             {"LoadBalanceEfficiency", [](CS s){return std::make_unique<LoadBalanceEfficiency>(s);}},
-            {"ParticleHistogram",     [](CS s){return std::make_unique<ParticleHistogram>(s);}},
-            {"ParticleHistogram2D",   [](CS s){return std::make_unique<ParticleHistogram2D>(s);}},
-            {"ParticleNumber",        [](CS s){return std::make_unique<ParticleNumber>(s);}},
-            {"ParticleExtrema",       [](CS s){return std::make_unique<ParticleExtrema>(s);}},
-            {"ChargeOnEB",  [](CS s){return std::make_unique<ChargeOnEB>(s);}}
+            {"RhoMaximum",            [](CS s){return std::make_unique<RhoMaximum>(s);}},
+            {"Timestep",              [](CS s){return std::make_unique<Timestep>(s);}}
     };
     // loop over all reduced diags and fill m_multi_rd with requested reduced diags
     std::transform(m_rd_names.begin(), m_rd_names.end(), std::back_inserter(m_multi_rd),

Source/Diagnostics/ReducedDiags/Timestep.H

@@ -0,0 +1,35 @@
+/* Copyright 2024 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Thomas Marks
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_TIMESTEP_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_TIMESTEP_H_
+
+#include "ReducedDiags.H"
+#include <string>
+
+/**
+ * This class contains a function for retrieving the current simulation timestep as a diagnostic.
+ * Useful mainly for simulations using adaptive timestepping.
+ */
+class Timestep: public ReducedDiags {
+public:
+    /**
+     * constructor
+     * @param[in] rd_name reduced diags name
+     */
+    Timestep (const std::string& rd_name);
+
+    /**
+     * This function gets the current physical timestep of the simulation at all refinement levels.
+     * @param[in] step current time step
+     */
+    void ComputeDiags (int step) final;
+};
+
+#endif //WARPX_DIAGNOSTICS_REDUCEDDIAGS_TIMESTEP_H_