diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms')
4 files changed, 82 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H index 9c52327e3..bb480364a 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H @@ -8,6 +8,8 @@ #ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_ #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_ +#include "Utils/WarpXConst.H" + #include <AMReX_REAL.H> #include <AMReX_Array4.H> #include <AMReX_Gpu.H> @@ -15,6 +17,7 @@ #include <algorithm> #include <array> + /** * This struct contains only static functions to initialize the stencil coefficients * and to compute finite-difference derivatives for the Cartesian CKC algorithm. @@ -91,6 +94,20 @@ struct CartesianCKCAlgorithm { } /** + * Compute the maximum timestep, for which the scheme remains stable + * (Courant-Friedrichs-Levy limit) */ + static amrex::Real ComputeMaxDt ( amrex::Real const * const dx ) { +#if (defined WARPX_DIM_XZ) + // - In Cartesian 2D geometry: determined by the minimum cell size in all direction + amrex::Real const delta_t = std::min( dx[0], dx[1] )/PhysConst::c; +#else + // - In Cartesian 3D geometry: determined by the minimum cell size in all direction + amrex::Real const delta_t = std::min( dx[0], std::min( dx[1], dx[2] ) ) / PhysConst::c; +#endif + return delta_t; + } + + /** * Perform derivative along x on a cell-centered grid, from a nodal field `F` */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real UpwardDx ( diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H index e5f89ff98..99510cc52 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H @@ -8,11 +8,15 @@ #ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_NODAL_H_ #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_NODAL_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 @@ -38,6 +42,19 @@ struct CartesianNodalAlgorithm { } /** + * Compute the maximum timestep, for which the scheme remains stable + * (Courant-Friedrichs-Levy limit) */ + static amrex::Real ComputeMaxDt ( amrex::Real const * const dx ) { + using namespace amrex::literals; + amrex::Real const delta_t = 1._rt / ( std::sqrt( AMREX_D_TERM( + 1._rt/(dx[0]*dx[0]), + + 1._rt/(dx[1]*dx[1]), + + 1._rt/(dx[2]*dx[2]) + ) ) * PhysConst::c ); + return delta_t; + } + + /** * Perform derivative along x * (For a solver on a staggered grid, `UpwardDx` and `DownwardDx` take into * account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H index 3a98bea0e..1c8609dfa 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -8,11 +8,15 @@ #ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_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 @@ -37,6 +41,19 @@ struct CartesianYeeAlgorithm { } /** + * Compute the maximum timestep, for which the scheme remains stable + * (Courant-Friedrichs-Levy limit) */ + static amrex::Real ComputeMaxDt ( amrex::Real const * const dx ) { + using namespace amrex::literals; + amrex::Real const delta_t = 1._rt / ( std::sqrt( AMREX_D_TERM( + 1._rt / (dx[0]*dx[0]), + + 1._rt / (dx[1]*dx[1]), + + 1._rt / (dx[2]*dx[2]) + ) ) * PhysConst::c ); + return delta_t; + } + + /** * Perform derivative along x on a cell-centered grid, from a nodal field `F`*/ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real UpwardDx ( 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. |