Merge branch 'development' into implicit_solvers_pc

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

Examples/Tests/open_bc_poisson_solver/CMakeLists.txt

@@ -12,3 +12,15 @@ if(WarpX_FFT)
         OFF  # dependency
     )
 endif()
+
+if(WarpX_HEFFTE)
+    add_warpx_test(
+        test_3d_open_bc_poisson_solver_heffte # name
+        3  # dims
+        2  # nprocs
+        inputs_test_3d_open_bc_poisson_solver_heffte  # inputs
+        analysis.py  # analysis
+        diags/diag1000001  # output
+        OFF  # dependency
+    )
+endif()

Examples/Tests/open_bc_poisson_solver/inputs_test_3d_open_bc_poisson_solver_heffte

@@ -0,0 +1 @@
+FILE = inputs_test_3d_open_bc_poisson_solver

GNUmakefile

@@ -38,6 +38,7 @@ USE_OPENPMD = FALSE
 WarpxBinDir = Bin
 
 USE_FFT = FALSE
+USE_HEFFTE = FALSE
 USE_RZ = FALSE
 
 USE_EB = FALSE

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
+  }
+}

Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver.json

@@ -1,19 +1,19 @@
 {
   "lev=0": {
-    "Bx": 100915975.15792876,
-    "By": 157610677.31483692,
-    "Bz": 2.404060922276648e-13,
-    "Ex": 4.725066923361703e+16,
-    "Ey": 3.0253961494347724e+16,
-    "Ez": 3276584.4383433666,
+    "Bx": 100915975.15403552,
+    "By": 157610677.3147734,
+    "Bz": 1.2276713711194638e-13,
+    "Ex": 4.725066923359797e+16,
+    "Ey": 3.025396149317578e+16,
+    "Ez": 3276584.4383433824,
     "rho": 10994013582437.197
   },
   "electron": {
-    "particle_momentum_x": 5.701279599504008e-19,
-    "particle_momentum_y": 3.650453172860547e-19,
+    "particle_momentum_x": 5.701279599509506e-19,
+    "particle_momentum_y": 3.650453172383178e-19,
     "particle_momentum_z": 1.145432768297242e-10,
     "particle_position_x": 17.31408691249785,
-    "particle_position_y": 0.2583691267187801,
+    "particle_position_y": 0.25836912671878015,
     "particle_position_z": 10066.329600000008,
     "particle_weight": 19969036501.910976
   }

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