add option for explicit substepping

Docs/sphinx_doc/TimeAdvance.rst

@@ -135,8 +135,14 @@ Then the acoustic substepping evolves the equations in the form
           - \frac{\partial (\beta_1 W^{\prime \prime, \tau} + \beta_2 W^{\prime \prime, \tau + \delta \tau})}{\partial z} +  R^t_{\rho}
             \right)
 
-where :math:`\beta_1 = 0.5 (1 - \beta_s)` and :math:`\beta_2 = 0.5 (1 + \beta_s)` with :math:`\beta_s = 0.1`.
-:math:`\beta_s` is the acoustic step off-centering coefficient and 0.1 is the typical WRF value. This off-centering is intended to provide damping of both horizontally and vertically propagating sound waves by biasing the time average toward the future time step.
+where :math:`\beta_1 = 0.5 (1 - \beta_s)` and :math:`\beta_2 = 0.5 (1 + \beta_s)`.
+
+:math:`\beta_s` is the acoustic step off-centering coefficient.  When we do implicit substepping, we use
+the typical WRF value of 0.1. This off-centering is intended to provide damping of both horizontally
+and vertically propagating sound waves by biasing the time average toward the future time step.
+
+When we do fully explicit substepping, we set :math:`\beta_s = -1.0`, which sets
+:math:`\beta_1 = 1` and :math:`\beta_2 = 0`.
 
 To solve the coupled system, we first evolve the equations for :math:`U^{\prime \prime, \tau + \delta \tau}`  and
 :math:`V^{\prime \prime, \tau + \delta \tau}` explicitly using :math:`\Theta^{\prime \prime, \tau}` which is already known.
@@ -149,10 +155,10 @@ to control horizontally propagating sound waves.
 .. math::
 
    p^{\prime\prime,\tau*} = p^{\prime\prime,\tau}
-     + \beta_d \left( p^{\prime\prime,\tau} + p^{\prime\prime,\tau-\delta\tau} \right)
+     + \beta_d \left( p^{\prime\prime,\tau} - p^{\prime\prime,\tau-\delta\tau} \right)
 
 where :math:`\tau*` is the forward projected value used in RHS of the acoustic
 substepping equations for horizontal momentum. According to Skamarock et al,
-This is equivalent to including a horizontal diffusion term in the continuity
+this is equivalent to including a horizontal diffusion term in the continuity
 equation. A typical damping coefficient of :math:`\beta_d = 0.1` is used, as in
 WRF.

Exec/CMakeLists.txt

@@ -20,7 +20,6 @@ elseif (ERF_ENABLE_REGRESSION_TESTS_ONLY)
 else ()
   add_subdirectory(ABL)
   add_subdirectory(SuperCell)
-  add_subdirectory(Radiation)
   add_subdirectory(SquallLine_2D)
   add_subdirectory(RegTests/Bubble)
   add_subdirectory(RegTests/Couette_Poiseuille)
@@ -35,9 +34,10 @@ else ()
   add_subdirectory(RegTests/WPS_Test)
   add_subdirectory(RegTests/Bomex)
   add_subdirectory(RegTests/TurbulentInflow)
-  add_subdirectory(DevTests/MovingTerrain)
-  add_subdirectory(DevTests/MetGrid)
   add_subdirectory(DevTests/LandSurfaceModel)
+  add_subdirectory(DevTests/MetGrid)
+  add_subdirectory(DevTests/MovingTerrain)
+  add_subdirectory(DevTests/Radiation)
   add_subdirectory(DevTests/TemperatureSource)
   add_subdirectory(DevTests/TropicalCyclone)
 endif()

Source/TimeIntegration/ERF_MRI.H

@@ -273,6 +273,10 @@ public:
                     nsubsteps =             1;   dtau =     timestep;
                 } else {
                     nsubsteps = substep_ratio;   dtau = sub_timestep;
+
+                    // STRT HACK -- this hack can be used to approximate the no-substepping algorithm
+                    // nsubsteps = 1;   dtau = timestep;
+                    // END HACK
                 }
                 time_stage = time + timestep;
             }

Source/TimeIntegration/ERF_TI_fast_headers.H

@@ -37,7 +37,8 @@ void erf_fast_rhs_N (int step, int nrk, int level, int finest_level,
                      std::unique_ptr<amrex::MultiFab>& mapfac_v,
                      amrex::YAFluxRegister* fr_as_crse,
                      amrex::YAFluxRegister* fr_as_fine,
-                     bool l_use_moisture, bool l_reflux);
+                     bool l_use_moisture, bool l_reflux,
+                     bool l_implicit_substepping);
 
 /**
  * Function for computing the fast RHS with fixed terrain
@@ -64,7 +65,8 @@ void erf_fast_rhs_T (int step, int nrk, int level, int finest_level,
                      std::unique_ptr<amrex::MultiFab>& mapfac_v,
                      amrex::YAFluxRegister* fr_as_crse,
                      amrex::YAFluxRegister* fr_as_fine,
-                     bool l_use_moisture, bool l_reflux);
+                     bool l_use_moisture, bool l_reflux,
+                     bool l_implicit_substepping);
 
 /**
  * Function for computing the fast RHS with moving terrain
@@ -98,7 +100,8 @@ void erf_fast_rhs_MT (int step, int nrk, int level, int finest_level,
                       std::unique_ptr<amrex::MultiFab>& mapfac_v,
                       amrex::YAFluxRegister* fr_as_crse,
                       amrex::YAFluxRegister* fr_as_fine,
-                      bool l_use_moisture, bool l_reflux);
+                      bool l_use_moisture, bool l_reflux,
+                      bool l_implicit_substepping);
 
 /**
  * Function for computing the coefficients for the tridiagonal solver used in the fast

Source/TimeIntegration/ERF_TI_fast_rhs_fun.H

@@ -20,7 +20,13 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
         // Per p2902 of Klemp-Skamarock-Dudhia-2007
         // beta_s = -1.0 : fully explicit
         // beta_s =  1.0 : fully implicit
-        Real beta_s = 0.1;
+        Real beta_s;
+        if (solverChoice.substepping_type[level] == SubsteppingType::Implicit)
+        {
+            beta_s = 0.1;
+        } else { // Fully explicit
+            beta_s = -1.0;
+        }
 
         // *************************************************************************
         // Set up flux registers if using two_way coupling
@@ -99,7 +105,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                                   detJ_cc[level],   detJ_cc_new[level],   detJ_cc_src[level],
                                 dtau, beta_s, inv_fac,
                                 mapfac_m[level], mapfac_u[level], mapfac_v[level],
-                                fr_as_crse, fr_as_fine, l_use_moisture, l_reflux);
+                                fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             } else {
                 // If this is not the first substep we pass in S_data as the previous step's solution
                 erf_fast_rhs_MT(fast_step, nrk, level, finest_level,
@@ -111,7 +117,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                                   detJ_cc[level],   detJ_cc_new[level],   detJ_cc_src[level],
                                 dtau, beta_s, inv_fac,
                                 mapfac_m[level], mapfac_u[level], mapfac_v[level],
-                                fr_as_crse, fr_as_fine, l_use_moisture, l_reflux);
+                                fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             }
         } else if (solverChoice.use_terrain && solverChoice.terrain_type == TerrainType::Static) {
             if (fast_step == 0) {
@@ -127,15 +133,15 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                                S_data, S_scratch, fine_geom, solverChoice.gravity, Omega,
                                z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac,
                                mapfac_m[level], mapfac_u[level], mapfac_v[level],
-                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux);
+                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             } else {
                 // If this is not the first substep we pass in S_data as the previous step's solution
                 erf_fast_rhs_T(fast_step, nrk, level, finest_level,
                                S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs,
                                S_data, S_scratch, fine_geom, solverChoice.gravity, Omega,
                                z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac,
                                mapfac_m[level], mapfac_u[level], mapfac_v[level],
-                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux);
+                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             }
         } else {
             if (fast_step == 0) {
@@ -151,15 +157,15 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
                                S_data, S_scratch, fine_geom, solverChoice.gravity,
                                dtau, beta_s, inv_fac,
                                mapfac_m[level], mapfac_u[level], mapfac_v[level],
-                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux);
+                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             } else {
                 // If this is not the first substep we pass in S_data as the previous step's solution
                 erf_fast_rhs_N(fast_step, nrk, level, finest_level,
                                S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs,
                                S_data, S_scratch, fine_geom, solverChoice.gravity,
                                dtau, beta_s, inv_fac,
                                mapfac_m[level], mapfac_u[level], mapfac_v[level],
-                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux);
+                               fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
             }
         }
 

Source/TimeIntegration/ERF_advance_dycore.cpp

@@ -90,6 +90,7 @@ void ERF::advance_dycore(int level,
     bool l_use_kturb   = ( (tc.les_type != LESType::None)   ||
                            (tc.pbl_type != PBLType::None) );
     bool l_use_moisture = ( solverChoice.moisture_type != MoistureType::None );
+    bool l_implicit_substepping = ( solverChoice.substepping_type[level] == SubsteppingType::Implicit );
 
     const bool use_most = (m_most != nullptr);
     const bool exp_most = (solverChoice.use_explicit_most);

Source/TimeIntegration/ERF_fast_rhs_MT.cpp

@@ -7,38 +7,40 @@ using namespace amrex;
  * Function for computing the fast RHS with moving terrain
  *
  * @param[in]    step  which fast time step within each Runge-Kutta step
- * @param[in]    nrk   which Runge-Kutta step
- * @param[in]    level level of resolution
- * @param[in]    finest_level finest level of resolution
- * @param[in]    S_slow_rhs slow RHS computed in erf_slow_rhs_pre
- * @param[in]    S_prev previous solution
- * @param[in]    S_stg_data solution            at previous RK stage
- * @param[in]    S_stg_prim primitive variables at previous RK stage
- * @param[in]    pi_stage   Exner function      at previous RK stage
- * @param[in]    fast_coeffs coefficients for the tridiagonal solve used in the fast integrator
- * @param[out]   S_data current solution
- * @param[in]    S_scratch scratch space
- * @param[in]    geom container for geometric information
- * @param[in]    gravity Magnitude of gravity
- * @param[in]    use_lagged_delta_rt define lagged_delta_rt for our next step
- * @param[in]    Omega component of the momentum normal to the z-coordinate surface
- * @param[in]    z_t_rk rate of change of grid height -- only relevant for moving terrain
- * @param[in]    z_t_pert rate of change of grid height -- interpolated between RK stages
- * @param[in]    z_phys_nd_old height coordinate at nodes at old time
- * @param[in]    z_phys_nd_new height coordinate at nodes at new time
- * @param[in]    z_phys_nd_stg height coordinate at nodes at previous stage
- * @param[in]    detJ_cc_old Jacobian of the metric transformation at old time
- * @param[in]    detJ_cc_new Jacobian of the metric transformation at new time
- * @param[in]    detJ_cc_stg Jacobian of the metric transformation at previous stage
- * @param[in]    dtau fast time step
- * @param[in]    beta_s  Coefficient which determines how implicit vs explicit the solve is
- * @param[in]    facinv inverse factor for time-averaging the momenta
- * @param[in]    mapfac_m map factor at cell centers
- * @param[in]    mapfac_u map factor at x-faces
- * @param[in]    mapfac_v map factor at y-faces
+ * @param[in   ] nrk   which Runge-Kutta step
+ * @param[in   ] level level of resolution
+ * @param[in   ] finest_level finest level of resolution
+ * @param[in   ] S_slow_rhs slow RHS computed in erf_slow_rhs_pre
+ * @param[in   ] S_prev previous solution
+ * @param[in   ] S_stg_data solution            at previous RK stage
+ * @param[in   ] S_stg_prim primitive variables at previous RK stage
+ * @param[in   ] pi_stage   Exner function      at previous RK stage
+ * @param[in   ] fast_coeffs coefficients for the tridiagonal solve used in the fast integrator
+ * @param[  out] S_data current solution
+ * @param[in   ] S_scratch scratch space
+ * @param[in   ] geom container for geometric information
+ * @param[in   ] gravity Magnitude of gravity
+ * @param[in   ] use_lagged_delta_rt define lagged_delta_rt for our next step
+ * @param[in   ] Omega component of the momentum normal to the z-coordinate surface
+ * @param[in   ] z_t_rk rate of change of grid height -- only relevant for moving terrain
+ * @param[in   ] z_t_pert rate of change of grid height -- interpolated between RK stages
+ * @param[in   ] z_phys_nd_old height coordinate at nodes at old time
+ * @param[in   ] z_phys_nd_new height coordinate at nodes at new time
+ * @param[in   ] z_phys_nd_stg height coordinate at nodes at previous stage
+ * @param[in   ] detJ_cc_old Jacobian of the metric transformation at old time
+ * @param[in   ] detJ_cc_new Jacobian of the metric transformation at new time
+ * @param[in   ] detJ_cc_stg Jacobian of the metric transformation at previous stage
+ * @param[in   ] dtau fast time step
+ * @param[in   ] beta_s  Coefficient which determines how implicit vs explicit the solve is
+ * @param[in   ] facinv inverse factor for time-averaging the momenta
+ * @param[in   ] mapfac_m map factor at cell centers
+ * @param[in   ] mapfac_u map factor at x-faces
+ * @param[in   ] mapfac_v map factor at y-faces
  * @param[inout] fr_as_crse YAFluxRegister at level l at level l   / l+1 interface
  * @param[inout] fr_as_fine YAFluxRegister at level l at level l-1 / l   interface
- * @param[in]    l_reflux should we add fluxes to the FluxRegisters?
+ * @param[in   ]  l_use_moisture
+ * @param[in   ]  l_reflux should we add fluxes to the FluxRegisters?
+ * @param[in   ]  l_implicit_substepping
  */
 
 void erf_fast_rhs_MT (int step, int nrk,
@@ -71,7 +73,8 @@ void erf_fast_rhs_MT (int step, int nrk,
                       YAFluxRegister* fr_as_crse,
                       YAFluxRegister* fr_as_fine,
                       bool l_use_moisture,
-                      bool l_reflux)
+                      bool l_reflux,
+                      bool /*l_implicit_substepping*/)
 {
     BL_PROFILE_REGION("erf_fast_rhs_MT()");