From fb34ee127216bb5b44fa51d415150482d1e9454f Mon Sep 17 00:00:00 2001
From: Zhi Chen <62574124+zhichen3@users.noreply.github.com>
Date: Thu, 19 Sep 2024 11:49:50 -0400
Subject: [PATCH] address timestep and cfl violation for Spherical2D (#2962)

Change timestep and cfl_violation so we use rdtheta, rdx[1], which is the physical cell length instead of just dx[1] for cell length in theta direction for spherical 2d.
---
 Source/driver/timestep.cpp    | 38 ++++++++++++++++++++++++++++++++++-
 Source/hydro/Castro_hydro.cpp |  8 ++++++++
 2 files changed, 45 insertions(+), 1 deletion(-)

diff --git a/Source/driver/timestep.cpp b/Source/driver/timestep.cpp
index a7b5936e4e..fd7e029a4b 100644
--- a/Source/driver/timestep.cpp
+++ b/Source/driver/timestep.cpp
@@ -35,6 +35,9 @@ Castro::estdt_cfl (int is_new)
   // Courant-condition limited timestep
 
   const auto dx = geom.CellSizeArray();
+  const auto problo = geom.ProbLoArray();
+  const auto coord = geom.Coord();
+  amrex::ignore_unused(problo, coord);
 
   const MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type);
 
@@ -82,6 +85,11 @@ Castro::estdt_cfl (int is_new)
       Real dt2;
 #if AMREX_SPACEDIM >= 2
       dt2 = dx[1]/(c + std::abs(uy));
+      if (coord == 2) {
+          // dx[1] in Spherical2D is just dtheta, need rdtheta for physical length
+          // so just multiply by the smallest r
+          dt2 *= problo[0] + 0.5_rt * dx[0];
+      }
 #else
       dt2 = dt1;
 #endif
@@ -127,6 +135,9 @@ Castro::estdt_mhd (int is_new)
 
   // MHD timestep limiter
   const auto dx = geom.CellSizeArray();
+  const auto problo = geom.ProbLoArray();
+  const auto coord = geom.Coord();
+  amrex::ignore_unused(problo, coord);
 
   const MultiFab& U_state = is_new ? get_new_data(State_Type) : get_old_data(State_Type);
 
@@ -206,6 +217,9 @@ Castro::estdt_mhd (int is_new)
       Real dt2;
 #if AMREX_SPACEDIM >= 2
       dt2 = dx[1]/(cy + std::abs(uy));
+      if (coord == 2) {
+          dt2 *= problo[0] + 0.5_rt * dx[0];
+      }
 #else
       dt2 = dt1;
 #endif
@@ -238,6 +252,9 @@ Castro::estdt_temp_diffusion (int is_new)
   // where D = k/(rho c_v), and k is the conductivity
 
   const auto dx = geom.CellSizeArray();
+  const auto problo = geom.ProbLoArray();
+  const auto coord = geom.Coord();
+  amrex::ignore_unused(problo, coord);
 
   const MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type);
 
@@ -287,6 +304,10 @@ Castro::estdt_temp_diffusion (int is_new)
           Real dt2;
 #if AMREX_SPACEDIM >= 2
           dt2 = 0.5_rt * dx[1]*dx[1] / D;
+          if (coord == 2) {
+              Real r = problo[0] + 0.5_rt * dx[0];
+              dt2 *= r * r;
+          }
 #else
           dt2 = dt1;
 #endif
@@ -320,6 +341,9 @@ Castro::estdt_burning (int is_new)
     }
 
     const auto dx = geom.CellSizeArray();
+    const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    amrex::ignore_unused(problo, coord);
 
     MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type);
 
@@ -368,7 +392,13 @@ Castro::estdt_burning (int is_new)
 #if AMREX_SPACEDIM == 1
         burn_state.dx = dx[0];
 #else
-        burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx[1], dx[2]));
+        Real dx1 = dx[1];
+#if AMREX_SPACEDIM >= 2
+        if (coord == 2) {
+            dx1 *= problo[0] + 0.5_rt * dx[0];
+        }
+#endif
+        burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx1, dx[2]));
 #endif
 
         burn_state.rho = S(i,j,k,URHO);
@@ -464,6 +494,9 @@ Real
 Castro::estdt_rad (int is_new)
 {
     auto dx = geom.CellSizeArray();
+    const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    amrex::ignore_unused(problo, coord);
 
     const MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type);
     const MultiFab& radMF = is_new ? get_new_data(Rad_Type) : get_old_data(Rad_Type);
@@ -523,6 +556,9 @@ Castro::estdt_rad (int is_new)
             Real dt1 = dx[0] / (c + std::abs(ux));
 #if AMREX_SPACEDIM >= 2
             Real dt2 = dx[1] / (c + std::abs(uy));
+            if (coord == 2) {
+                dt2 *= problo[0] + 0.5_rt * dx[0];
+            }
 #else
             Real dt2 = std::numeric_limits<Real>::max();
 #endif
diff --git a/Source/hydro/Castro_hydro.cpp b/Source/hydro/Castro_hydro.cpp
index 0dedbc17aa..0e2250ed20 100644
--- a/Source/hydro/Castro_hydro.cpp
+++ b/Source/hydro/Castro_hydro.cpp
@@ -237,12 +237,20 @@ Castro::check_for_cfl_violation(const MultiFab& State, const Real dt)
     int cfl_violation = 0;
 
     auto dx = geom.CellSizeArray();
+    const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    amrex::ignore_unused(problo, coord);
 
     Real dtdx = dt / dx[0];
 
     Real dtdy = 0.0_rt;
     if (AMREX_SPACEDIM >= 2) {
       dtdy = dt / dx[1];
+      if (coord == 2) {
+          // dx[1] in Spherical2D is just rdtheta, need rdtheta for physical length
+          // Just choose to divide by the smallest r
+          dtdy /= problo[0] + 0.5_rt * dx[0];
+      }
     }
 
     Real dtdz = 0.0_rt;