1 files changed, 31 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
index cdc693a78..4e01041b2 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
@@ -8,11 +8,15 @@
 #ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_
 #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_
 
+#include "Utils/WarpXConst.H"
+
 #include <AMReX_REAL.H>
 #include <AMReX_Array4.H>
 #include <AMReX_Gpu.H>
 
 #include <array>
+#include <cmath>
+
 
 /**
  * This struct contains only static functions to initialize the stencil coefficients
@@ -33,6 +37,33 @@ struct CylindricalYeeAlgorithm {
         stencil_coefs_z[0] = 1._rt/cell_size[2];  // 1./dz
     }
 
+    /** Compute the maximum, CFL-stable timestep
+     *
+     * Compute the maximum timestep, for which the scheme remains stable
+     * under the Courant-Friedrichs-Levy limit.
+     */
+    static amrex::Real ComputeMaxDt ( amrex::Real const * const dx,
+                                      int const n_rz_azimuthal_modes ) {
+        using namespace amrex::literals;
+        // In the rz case, the Courant limit has been evaluated
+        // semi-analytically by R. Lehe, and resulted in the following
+        // coefficients.
+        std::array< amrex::Real, 6 > const multimode_coeffs = {{ 0.2105, 1.0, 3.5234, 8.5104, 15.5059, 24.5037 }};
+        amrex::Real multimode_alpha;
+        if (n_rz_azimuthal_modes < 7) {
+          // Use the table of the coefficients
+          multimode_alpha = multimode_coeffs[n_rz_azimuthal_modes-1];
+        } else {
+          // Use a realistic extrapolation
+          multimode_alpha = (n_rz_azimuthal_modes - 1._rt)*(n_rz_azimuthal_modes - 1._rt) - 0.4_rt;
+        }
+        amrex::Real delta_t = 1._rt / ( std::sqrt(
+                                     (1._rt + multimode_alpha) / (dx[0]*dx[0])
+                                    + 1._rt / (dx[1]*dx[1])
+                              ) * PhysConst::c );
+        return delta_t;
+    }
+
     /** Applies the differential operator `1/r * d(rF)/dr`,
      * where `F` is on a *nodal* grid in `r`
      * and the differential operator is evaluated on a *cell-centered* grid.