port some routines to C++

Source/driver/Castro.H

@@ -574,6 +574,17 @@ public:
 ///
     amrex::Real estTimeStep (amrex::Real dt_old);
 
+
+///
+/// Compute the CFL timestep
+///
+    amrex::Real estdt_cfl(const amrex::Real time);
+
+///
+/// Diffusion-limited timestep
+///
+    amrex::Real estdt_temp_diffusion();
+
 ///
 /// Compute initial time step.
 ///

Source/driver/Castro.cpp

@@ -1287,7 +1287,7 @@ Castro::estTimeStep (Real dt_old)
     Real estdt = max_dt;
 
     const MultiFab& stateMF = get_new_data(State_Type);
-
+    Real time = state[State_Type].curTime();
     const Real* dx = geom.CellSize();
 
     std::string limiter = "castro.max_dt";
@@ -1337,24 +1337,7 @@ Castro::estTimeStep (Real dt_old)
         {
 #endif
 
-#ifdef _OPENMP
-#pragma omp parallel reduction(min:estdt_hydro)
-#endif
-            {
-                Real dt = max_dt / cfl;
-
-                for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-                {
-                    const Box& box = mfi.tilebox();
-
-#pragma gpu box(box)
-                    ca_estdt(AMREX_INT_ANYD(box.loVect()), AMREX_INT_ANYD(box.hiVect()),
-                             BL_TO_FORTRAN_ANYD(stateMF[mfi]),
-                             AMREX_REAL_ANYD(dx),
-                             AMREX_MFITER_REDUCE_MIN(&dt));
-                }
-                estdt_hydro = std::min(estdt_hydro, dt);
-            }
+          estdt_hydro = estdt_cfl(time);
 
 #ifdef RADIATION
         }
@@ -1384,23 +1367,7 @@ Castro::estTimeStep (Real dt_old)
 
     if (diffuse_temp)
     {
-#ifdef _OPENMP
-#pragma omp parallel reduction(min:estdt_diffusion)
-#endif
-        {
-            Real dt = max_dt / cfl;
-
-            for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-            {
-                const Box& box = mfi.tilebox();
-
-#pragma gpu box(box)
-                ca_estdt_temp_diffusion(AMREX_INT_ANYD(box.loVect()), AMREX_INT_ANYD(box.hiVect()),
-                                        BL_TO_FORTRAN_ANYD(stateMF[mfi]),
-                                        AMREX_REAL_ANYD(dx), AMREX_MFITER_REDUCE_MIN(&dt));
-            }
-            estdt_diffusion = std::min(estdt_diffusion, dt);
-        }
+      estdt_diffusion = estdt_temp_diffusion();
     }
 
     ParallelDescriptor::ReduceRealMin(estdt_diffusion);

Source/driver/Castro_F.H

@@ -160,17 +160,6 @@ extern "C"
     (const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(S_new),
      const int verbose);
 
-  void ca_estdt
-    (const int* lo, const int* hi,
-     const BL_FORT_FAB_ARG_3D(state),
-     const amrex::Real* dx, amrex::Real* dt);
-
-#ifdef DIFFUSION
-  void ca_estdt_temp_diffusion
-    (const int* lo, const int* hi,
-     const BL_FORT_FAB_ARG_3D(state),
-     const amrex::Real* dx, amrex::Real* dt);
-#endif
 
 #ifdef RADIATION
   void ca_estdt_rad

Source/driver/Make.package

@@ -62,3 +62,4 @@ endif
 ca_F90EXE_sources += prob_params_nd.F90
 ca_F90EXE_sources += Tagging_nd.F90
 ca_F90EXE_sources += timestep_nd.F90
+CEXE_sources += timestep.cpp

Source/driver/timestep.cpp

@@ -0,0 +1,225 @@
+#include "Castro.H"
+#include "Castro_F.H"
+
+#ifdef DIFFUSION
+#include "conductivity.H"
+#endif
+
+using namespace amrex;
+
+Real
+Castro::estdt_cfl(const Real time)
+{
+
+  // Courant-condition limited timestep
+
+  GpuArray<Real, 3> center;
+  ca_get_center(center.begin());
+
+#ifdef ROTATION
+  GpuArray<Real, 3> omega;
+  get_omega(time, omega.begin());
+#endif
+
+  const auto dx = geom.CellSizeArray();
+
+  ReduceOps<ReduceOpMin> reduce_op;
+  ReduceData<Real> reduce_data(reduce_op);
+  using ReduceTuple = typename decltype(reduce_data)::Type;
+
+  const MultiFab& stateMF = get_new_data(State_Type);
+
+  const int ltime_integration_method = time_integration_method;
+#ifdef ROTATION
+  const int ldo_rotation = do_rotation;
+  const int lstate_in_rotating_frame = state_in_rotating_frame;
+#endif
+
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+  for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+    const Box& box = mfi.tilebox();
+
+    auto u = stateMF.array(mfi);
+
+    reduce_op.eval(box, reduce_data,
+    [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept -> ReduceTuple
+    {
+
+      Real rhoInv = 1.0_rt / u(i,j,k,URHO);
+
+      eos_t eos_state;
+      eos_state.rho = u(i,j,k,URHO);
+      eos_state.T = u(i,j,k,UTEMP);
+      eos_state.e = u(i,j,k,UEINT) * rhoInv;
+      for (int n = 0; n < NumSpec; n++) {
+        eos_state.xn[n] = u(i,j,k,UFS+n) * rhoInv;
+      }
+      for (int n = 0; n < NumAux; n++) {
+        eos_state.aux[n] = u(i,j,k,UFX+n) * rhoInv;
+      }
+
+      eos(eos_input_re, eos_state);
+
+      // Compute velocity and then calculate CFL timestep.
+
+      Real ux = u(i,j,k,UMX) * rhoInv;
+      Real uy = u(i,j,k,UMY) * rhoInv;
+      Real uz = u(i,j,k,UMZ) * rhoInv;
+
+#ifdef ROTATION
+      if (ldo_rotation == 1 && lstate_in_rotating_frame != 1) {
+        Real vel[3];
+        vel[0] = ux;
+        vel[1] = uy;
+        vel[2] = uz;
+
+        GeometryData geomdata = geom.data();
+
+        inertial_to_rotational_velocity_c(i, j, k, geomdata,
+                                          center.begin(), omega.begin(), time, vel);
+
+        ux = vel[0];
+        uy = vel[1];
+        uz = vel[2];
+      }
+#endif
+
+      Real c = eos_state.cs;
+
+      Real dt1 = dx[0]/(c + std::abs(ux));
+
+      Real dt2;
+#if AMREX_SPACEDIM >= 2
+      dt2 = dx[1]/(c + std::abs(uy));
+#else
+      dt2 = dt1;
+#endif
+
+      Real dt3;
+#if AMREX_SPACEDIM == 3
+      dt3 = dx[2]/(c + std::abs(uz));
+#else
+      dt3 = dt1;
+#endif
+
+      // The CTU method has a less restrictive timestep than MOL-based
+      // schemes (including the true SDC).  Since the simplified SDC
+      // solver is based on CTU, we can use its timestep.
+      if (ltime_integration_method == 0 || ltime_integration_method == 3) {
+        return {amrex::min(dt1, dt2, dt3)};
+
+      } else {
+        // method of lines-style constraint is tougher
+        Real dt_tmp = 1.0_rt/dt1;
+#if AMREX_SPACEIM >= 2
+        dt_tmp += 1.0_rt/dt2;
+#endif
+#if AMREX_SPACEDIM == 3
+        dt_tmp += 1.0_rt/dt3;
+#endif
+
+        return 1.0_rt/dt_tmp;
+      }
+
+    });
+
+  }
+
+  ReduceTuple hv = reduce_data.value();
+  Real estdt_hydro = amrex::get<0>(hv);
+
+  return estdt_hydro;
+
+}
+
+#ifdef DIFFUSION
+
+Real
+Castro::estdt_temp_diffusion(void)
+{
+
+  // Diffusion-limited timestep
+  //
+  // dt < 0.5 dx**2 / D
+  // where D = k/(rho c_v), and k is the conductivity
+
+  const auto dx = geom.CellSizeArray();
+
+  ReduceOps<ReduceOpMin> reduce_op;
+  ReduceData<Real> reduce_data(reduce_op);
+  using ReduceTuple = typename decltype(reduce_data)::Type;
+
+  const MultiFab& stateMF = get_new_data(State_Type);
+
+  const Real ldiffuse_cutoff_density = diffuse_cutoff_density;
+  const Real lmax_dt = max_dt;
+  const Real lcfl = cfl;
+
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+  for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+    const Box& box = mfi.tilebox();
+
+    auto ustate = stateMF.array(mfi);
+
+    reduce_op.eval(box, reduce_data,
+                   [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept -> ReduceTuple
+                   {
+
+                     if (ustate(i,j,k,URHO) > ldiffuse_cutoff_density) {
+
+                       Real rho_inv = 1.0_rt/ustate(i,j,k,URHO);
+
+                       // we need cv
+                       eos_t eos_state;
+                       eos_state.rho = ustate(i,j,k,URHO);
+                       eos_state.T = ustate(i,j,k,UTEMP);
+                       eos_state.e = ustate(i,j,k,UEINT) * rho_inv;
+                       for (int n = 0; n < NumSpec; n++) {
+                         eos_state.xn[n]  = ustate(i,j,k,UFS+n) * rho_inv;
+                       }
+                       for (int n = 0; n < NumAux; n++) {
+                         eos_state.aux[n] = ustate(i,j,k,UFX+n) * rho_inv;
+                       }
+
+                       eos(eos_input_re, eos_state);
+
+                       // we also need the conductivity
+                       conductivity(eos_state);
+
+                       // maybe we should check (and take action) on negative cv here?
+                       Real D = eos_state.conductivity * rho_inv / eos_state.cv;
+
+                       Real dt1 = 0.5_rt * dx[0]*dx[0] / D;
+
+                       Real dt2;
+#if AMREX_SPACEDIM >= 2
+                       dt2 = 0.5_rt * dx[1]*dx[1] / D;
+#else
+                       dt2 = dt1;
+#endif
+
+                       Real dt3;
+#if AMREX_SPACEDIM >= 3
+                       dt3 = 0.5_rt * dx[2]*dx[2] / D;
+#else
+                       dt3 = dt1;
+#endif
+
+                       return {amrex::min(dt1, dt2, dt3)};
+
+                     } else {
+                       return lmax_dt/lcfl;
+      }
+    });
+  }
+
+  ReduceTuple hv = reduce_data.value();
+  Real estdt_diff = amrex::get<0>(hv);
+
+  return estdt_diff;
+}
+#endif