Add new time step size limiters (#6840)

This add time step limiters based on the plasma frequency and cyclotron
frequency. This is in addition to the existing limiter based on the
particle speed and controlled by the parameter
`warpx.dt_update_interval`. This adds two new input parameters
`warpx.max_omegap_dt` and `warpx.max_omegac_dt` for the plasma and
cyclotron frequency limits.

For this PR, when calculating the cyclotron frequency, only the B-field
MultiFabs and the constant external fields are used when finding the
maximum B-field. This ignores the rest of the external fields due to the
complication or cost required to include them. If required, a later PR
could include them, perhaps by fetching the B-field for the particles
(which could be expensive) or by calculating the external B-fields at
the grid locations, which would be complicated.

Also, currently, for the initial time steps the `warpx.max_dt` or
`warpx.cfl` with speed `c` is used. It would be better to apply the
limiters initially rather that requiring the user to specify the initial
step. This can be done in this PR or a separate PR.

Todo:
- [x] Add CI test
- [x] Add documentation

Author: dpgrote

Docs/source/usage/parameters.rst

@@ -3184,17 +3184,40 @@ Time step
 
 .. pp:param:: warpx.dt_update_interval
     :type: ``string``
-    :default: ``-1``
     :optional:
 
-    How many iterations pass between timestep adaptations when using the electrostatic solver.
-    Must be greater than ``0`` to use adaptive timestepping, or else :pp:param:`warpx.const_dt` must be specified.
+    This controls adaptive timestepping, where the time step size is updated based on the conditions of the simulation, and only applies when using the explicit electrostatic or theta-implicit solvers.
+    This specifies time step intervals when the time step size is updated.
+    The value must be greater than ``0``.
+    When specified, :pp:param:`warpx.const_dt` must not also be specified.
+    The time step size is updated using the limits specified by :pp:param:`warpx.cfl`, :pp:param:`warpx.max_omegap_dt`, and :pp:param:`warpx.max_omegac_dt`.
+
+.. pp:param:: warpx.dt_update_diagnostic_file
+    :type: ``string``
+    :optional:
+
+    When adaptive timestepping is activated, information about the new time step and the simulation conditions are output to the file specified by this parameter.
+
+.. pp:param:: warpx.max_omegap_dt
+    :type: ``float``
+    :optional:
+
+    With adaptive timestepping, the time step size is limited to be less than or equal to the value specified divided by the global plasma frequency.
+    The application of this limit is controlled by :pp:param:`warpx.dt_update_interval`, and is only applied when using the explicit electrostatic or theta-implicit solver..
+
+.. pp:param:: warpx.max_omegac_dt
+    :type: ``float``
+    :optional:
+
+    With adaptive timestepping, the time step size is limited to be less than or equal to the value specified divided by the maximum cyclotron frequency.
+    Note that the maximum B-field is calculated from using only the constant applied B field (as set by :pp:param:`particles.B_external_particle`) and the B-field grid data.
+    The application of this limit is controlled by :pp:param:`warpx.dt_update_interval`, and is only applied when using the explicit electrostatic or theta-implicit solver..
 
 .. pp:param:: warpx.max_dt
     :type: ``float``
     :optional:
 
-    The maximum timestep permitted for the electrostatic solver, when using adaptive timestepping.
+    The maximum timestep permitted when using adaptive timestepping.
     If supplied, also sets the initial timestep for these simulations, before the first timestep update.
 
 Filtering

Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive

@@ -1,6 +1,9 @@
 stop_time = 60e-6
 warpx.cfl = 0.2
 warpx.dt_update_interval = 10
+warpx.dt_update_diagnostic_file = diags/reducedfiles/dt_updates.txt
+warpx.max_omegap_dt = 0.2
+warpx.max_omegac_dt = 0.2
 warpx.max_dt = 1.5e-6
 amr.n_cell = 64 64 64
 amr.max_level = 0

Source/Evolve/WarpXComputeDt.cpp

@@ -15,20 +15,26 @@
 #   include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
 #   include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
 #endif
+#include "Fields.H"
 #include "Particles/MultiParticleContainer.H"
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXAlgorithmSelection.H"
 #include "Utils/WarpXConst.H"
+#include "Utils/Parser/ParserUtils.H"
+
+#include <ablastr/coarsen/sample.H>
 
 #include <AMReX.H>
 #include <AMReX_Geometry.H>
 #include <AMReX_IntVect.H>
+#include <AMReX_MFIter.H>
 #include <AMReX_Print.H>
 #include <AMReX_REAL.H>
 #include <AMReX_Vector.H>
 
 #include <algorithm>
 #include <memory>
+#include <filesystem>
 
 /**
  * Compute the minimum of array x, where x has dimension AMREX_SPACEDIM
@@ -108,35 +114,229 @@ WarpX::ComputeDt ()
 }
 
 /**
- * Determine the simulation timestep from the maximum speed of all particles
- * Sets timestep so that a particle can only cross cfl*dx cells per timestep.
+ * Used to determine the simulation timestep from the maximum speed of all particles
+ * Timestep will be set so that a particle can cross at most cfl*dx cells per timestep.
  */
-void
-WarpX::UpdateDtFromParticleSpeeds ()
+amrex::Real
+WarpX::ParticleGridSpeedMax ()
 {
     const amrex::Real* dx = geom[max_level].CellSize();
     const amrex::Real dx_min = minDim(dx);
 
     const amrex::ParticleReal max_v = mypc->maxParticleVelocity();
-    amrex::Real deltat_new = 0.;
 
-    // Protections from overly-large timesteps
-    if (max_v == 0) {
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(), "Particles at rest and no constant or maximum timestep specified. Aborting.");
-        deltat_new = m_max_dt.value();
-    } else {
-        deltat_new = cfl * dx_min / max_v;
+    return max_v/dx_min;
+}
+
+amrex::Real
+WarpX::GlobalPlasmaFrequencyMax ()
+{
+    const std::unique_ptr<amrex::MultiFab> global_plasma_frequency = mypc->GetGlobalPlasmaFrequency(0);
+    const amrex::Real global_plasma_frequency_max = global_plasma_frequency->max(0);
+    return global_plasma_frequency_max;
+}
+
+amrex::Real
+WarpX::GlobalCyclotronFrequencyMax ()
+{
+    using ablastr::fields::Direction;
+    using warpx::fields::FieldType;
+
+    const amrex::ParmParse pp_particles("particles");
+    amrex::Vector<amrex::Real> B_external_particle(3, 0.);
+    utils::parser::queryArrWithParser(pp_particles, "B_external_particle", B_external_particle);
+
+    const amrex::Real Bx_external = B_external_particle[0];
+    const amrex::Real By_external = B_external_particle[1];
+    const amrex::Real Bz_external = B_external_particle[2];
+
+    amrex::Real B_max = 0.;
+
+    // loop over refinement levels
+    for (int lev = 0; lev <= finestLevel(); ++lev)
+    {
+        // get MultiFab data at lev
+        const amrex::MultiFab & Bx = *m_fields.get(FieldType::Bfield_aux, Direction{0}, lev);
+        const amrex::MultiFab & By = *m_fields.get(FieldType::Bfield_aux, Direction{1}, lev);
+        const amrex::MultiFab & Bz = *m_fields.get(FieldType::Bfield_aux, Direction{2}, lev);
+
+        // Prepare interpolation of field components to cell center
+        // The arrays below store the index type (staggering) of each MultiFab, with the third
+        // component set to zero in the two-dimensional case.
+        auto Bxtype = amrex::GpuArray<int,3>{0, 0, 0};
+        auto Bytype = amrex::GpuArray<int,3>{0, 0, 0};
+        auto Bztype = amrex::GpuArray<int,3>{0, 0, 0};
+        for (int i = 0; i < AMREX_SPACEDIM; ++i){
+            Bxtype[i] = Bx.ixType()[i];
+            Bytype[i] = By.ixType()[i];
+            Bztype[i] = Bz.ixType()[i];
+        }
+
+        // General preparation of interpolation and reduction operations
+        const amrex::GpuArray<int,3> cellCenteredtype{0,0,0};
+        const amrex::GpuArray<int,3> reduction_coarsening_ratio{1,1,1};
+        constexpr int reduction_comp = 0;
+
+        amrex::ReduceOps<amrex::ReduceOpMax> reduce_op;
+        amrex::ReduceData<amrex::Real> reduce_data(reduce_op);
+        using ReduceTuple = typename decltype(reduce_data)::Type;
+
+        // MFIter loop to interpolate fields to cell center and get maximum value
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (amrex::MFIter mfi(Bx, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            // Make the box cell centered in preparation for the interpolation (and to avoid
+            // including ghost cells in the calculation)
+            const amrex::Box & box = enclosedCells(mfi.nodaltilebox());
+            const auto& arrBx = Bx[mfi].array();
+            const auto& arrBy = By[mfi].array();
+            const auto& arrBz = Bz[mfi].array();
+
+            reduce_op.eval(box, reduce_data,
+            [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
+            {
+                const amrex::Real Bx_interp = ablastr::coarsen::sample::Interp(arrBx, Bxtype, cellCenteredtype,
+                                                                               reduction_coarsening_ratio, i, j, k, reduction_comp);
+                const amrex::Real By_interp = ablastr::coarsen::sample::Interp(arrBy, Bytype, cellCenteredtype,
+                                                                               reduction_coarsening_ratio, i, j, k, reduction_comp);
+                const amrex::Real Bz_interp = ablastr::coarsen::sample::Interp(arrBz, Bztype, cellCenteredtype,
+                                                                               reduction_coarsening_ratio, i, j, k, reduction_comp);
+                return {amrex::Math::powi<2>(Bx_interp + Bx_external) +
+                        amrex::Math::powi<2>(By_interp + By_external) +
+                        amrex::Math::powi<2>(Bz_interp + Bz_external)};
+            });
+        }
+
+        const amrex::Real hv_Bsq = amrex::get<0>(reduce_data.value()); // highest value of |B|**2
+
+        B_max = std::max(B_max, std::sqrt(hv_Bsq));
+
+    }
+
+    // MPI reduce
+    amrex::ParallelDescriptor::ReduceRealMax({B_max});
+
+    amrex::Real omegac_max = 0.;
+
+    const int n_containers = mypc->nContainers();
+    for (int i = 0; i < n_containers; i++)
+    {
+        const WarpXParticleContainer& pc = mypc->GetParticleContainer(i);
+        if (pc.getMass() > 0.) {
+            const amrex::Real pc_omegac = pc.getCharge()*B_max/pc.getMass();
+            omegac_max = std::max(omegac_max, pc_omegac);
+        }
+    }
+
+    return omegac_max;
+}
+
+void
+WarpX::ApplyDtLimiters ()
+{
+    using namespace amrex::literals;
+
+    // Calculate limiting values from the simulation conditions
+    const amrex::Real vmax_o_dx = ParticleGridSpeedMax();
+    const amrex::Real omegap_max = m_max_omegap_dt.has_value() ? GlobalPlasmaFrequencyMax() : 0._rt;
+    const amrex::Real omegac_max = m_max_omegac_dt.has_value() ? GlobalCyclotronFrequencyMax() : 0._rt;
+
+    // Ensure that a valid time step value exists, either from the simulation conditions or from max_dt
+    if (vmax_o_dx == 0._rt &&
+        (!m_max_omegap_dt.has_value() || omegap_max == 0._rt) &&
+        (!m_max_omegac_dt.has_value() || omegac_max == 0._rt)) {
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(),
+                                         "No valid time step size limit found, warpx.max_dt must be specified");
+    }
+
+    amrex::Real dt_new = std::numeric_limits<amrex::Real>::max();
+
+    if (vmax_o_dx > 0._rt) {
+        dt_new = std::min(dt_new, cfl/vmax_o_dx);
+    }
+    if (m_max_omegap_dt.has_value() && omegap_max > 0._rt) {
+        dt_new = std::min(dt_new, m_max_omegap_dt.value()/omegap_max);
+    }
+    if (m_max_omegac_dt.has_value() && omegac_max > 0._rt) {
+        dt_new = std::min(dt_new, m_max_omegac_dt.value()/omegac_max);
     }
 
-    // Restrict to be less than user-specified maximum timestep, if present
     if (m_max_dt.has_value()) {
-        deltat_new = std::min(deltat_new, m_max_dt.value());
+        dt_new = std::min(dt_new, m_max_dt.value());
     }
 
     // Update dt
-    dt[max_level] = deltat_new;
+    dt[max_level] = dt_new;
 
     for (int lev = max_level-1; lev >= 0; --lev) {
         dt[lev] = dt[lev+1] * refRatio(lev)[0];
     }
+
+    // Write diagnostics if requested
+    if (amrex::ParallelDescriptor::IOProcessor()
+        && !m_dt_update_diagnostic_file.empty()
+        && !amrex::FileExists(m_dt_update_diagnostic_file)) {
+
+        std::filesystem::path const diagnostic_path(m_dt_update_diagnostic_file);
+        std::filesystem::path const diagnostic_dir = diagnostic_path.parent_path();
+        if (!diagnostic_dir.empty()) {
+            std::filesystem::create_directories(diagnostic_dir);
+        }
+
+        std::ofstream diagnostic_file{m_dt_update_diagnostic_file, std::ofstream::out | std::ofstream::trunc};
+        if (!diagnostic_file.is_open()) {
+            amrex::Abort("Failed to open file: " + m_dt_update_diagnostic_file);
+        }
+
+        int c = 0;
+        diagnostic_file << "#";
+        diagnostic_file << "[" << c++ << "]step()";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]time(s)";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]new_dt";
+        diagnostic_file << " ";
+        diagnostic_file << "[" << c++ << "]vmax_dt";
+        if (m_max_omegap_dt.has_value()) {
+            diagnostic_file << " ";
+            diagnostic_file << "[" << c++ << "]omegap_dt";
+        }
+        if (m_max_omegac_dt.has_value()) {
+            diagnostic_file << " ";
+            diagnostic_file << "[" << c++ << "]omegac_dt";
+        }
+        diagnostic_file << "\n";
+        diagnostic_file.close();
+    }
+
+    if (amrex::ParallelDescriptor::IOProcessor()
+        && !m_dt_update_diagnostic_file.empty()) {
+
+        std::ofstream diagnostic_file{m_dt_update_diagnostic_file, std::ofstream::out | std::ofstream::app};
+        if (!diagnostic_file.is_open()) {
+            amrex::Abort("Failed to open file: " + m_dt_update_diagnostic_file);
+        }
+
+        diagnostic_file << std::setprecision(14);
+        diagnostic_file << istep[0] + 1;
+        diagnostic_file << " ";
+        diagnostic_file << t_new[0];
+        diagnostic_file << " ";
+        diagnostic_file << dt_new;
+        diagnostic_file << " ";
+        diagnostic_file << vmax_o_dx*dt_new;
+
+        if (m_max_omegap_dt.has_value()) {
+            diagnostic_file << " ";
+            diagnostic_file << omegap_max*dt_new;
+        }
+        if (m_max_omegac_dt.has_value()) {
+            diagnostic_file << " ";
+            diagnostic_file << omegac_max*dt_new;
+        }
+        diagnostic_file << "\n";
+        diagnostic_file.close();
+    }
 }

Source/Evolve/WarpXEvolve.cpp

@@ -195,11 +195,11 @@ WarpX::Evolve (int numsteps)
         // Update the timestep for solvers that support adaptive timestepping
         // (electrostatic and theta-implicit EM), provided const_dt is not specified.
         if (m_dt_update_interval.contains(step+1)) {
+            SynchronizeVelocityWithPosition();
+            ApplyDtLimiters();
             if (verbose_step) {
-                amrex::Print() << Utils::TextMsg::Info("updating timestep");
+                amrex::Print() << Utils::TextMsg::Info("updating timestep to dt = " + std::to_string(dt[0]));
             }
-            SynchronizeVelocityWithPosition();
-            UpdateDtFromParticleSpeeds();
         }
 
         // If position and velocity are synchronized, push velocity backward one half step

Source/Particles/MultiParticleContainer.H

@@ -196,6 +196,7 @@ public:
 
     std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false);
 
+    std::unique_ptr<amrex::MultiFab> GetGlobalPlasmaFrequency (int lev);
     void GenerateGlobalDebyeLength ();
 
     void doFieldIonization (int lev,

Source/Particles/MultiParticleContainer.cpp

@@ -693,6 +693,47 @@ MultiParticleContainer::GetChargeDensity (int lev, bool local)
     return rho;
 }
 
+std::unique_ptr<amrex::MultiFab>
+MultiParticleContainer::GetGlobalPlasmaFrequency (int lev)
+{
+    const WarpX & warpx = WarpX::GetInstance();
+
+    amrex::BoxArray const & ba = warpx.boxArray(lev);
+    amrex::DistributionMapping const & dmap = warpx.DistributionMap(lev);
+    int const ncomps = 1;
+    const amrex::IntVect ng = amrex::IntVect::TheZeroVector();
+    auto global_plasma_frequency = std::make_unique<amrex::MultiFab>(ba, dmap, ncomps, ng);
+    global_plasma_frequency->setVal(amrex::Real(0.0));
+
+    for (auto& pc : allcontainers) {
+
+        if (pc->getMass() == 0. || pc->getCharge() == 0.) {
+            continue;
+        }
+
+        std::unique_ptr<amrex::MultiFab> plasma_frequency = pc->GetPlasmaFrequency(lev);
+
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (amrex::MFIter mfi(*global_plasma_frequency, TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            const amrex::Box box = mfi.tilebox();
+
+            amrex::Array4<amrex::Real> const& omegap_array = plasma_frequency->array(mfi);
+            amrex::Array4<amrex::Real> const& global_omegap_array = global_plasma_frequency->array(mfi);
+
+            amrex::ParallelFor(box,
+                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    amrex::Real const omegap = omegap_array(i,j,k);
+                    amrex::Real const global_omegap = global_omegap_array(i,j,k);
+                    global_omegap_array(i,j,k) = std::sqrt(global_omegap*global_omegap + omegap*omegap);
+                });
+        }
+    }
+    return global_plasma_frequency;
+}
+
 void
 MultiParticleContainer::GenerateGlobalDebyeLength ()
 {