CFL constraint shouldn't see z-direction when using substepping (#1169)

* CFL constraint shouldn't see z-direction when using substepping

* clean up how some parameters are stored

* fix trailing white space

Source/DataStruct.H

@@ -189,6 +189,18 @@ struct SolverChoice {
  rhoAlpha_T = rho0_trans * alpha_T;
  rhoAlpha_C = rho0_trans * alpha_C;
 
+ // Turn off acoustic substepping?
+ pp.query("no_substepping", no_substepping);
+
+ pp.query("force_stage1_single_substep", force_stage1_single_substep);
+ pp.query("incompressible", incompressible);
+
+ // If this is set, it must be even
+ if (incompressible != 0 && no_substepping == 0)
+ {
+ amrex::Abort("If you specify incompressible, you must specific no_substepping");
+ }
+
  // Order of spatial discretization
  pp.query("horiz_spatial_order", horiz_spatial_order);
  pp.query("vert_spatial_order", vert_spatial_order);
@@ -275,24 +287,27 @@ struct SolverChoice {
  void display()
  {
  amrex::Print() << "SOLVER CHOICE: " << std::endl;
- amrex::Print() << "use_coriolis : " << use_coriolis << std::endl;
- amrex::Print() << "use_rayleigh_damping : " << use_rayleigh_damping << std::endl;
- amrex::Print() << "use_gravity : " << use_gravity << std::endl;
- amrex::Print() << "rho0_trans : " << rho0_trans << std::endl;
- amrex::Print() << "alpha_T : " << alpha_T << std::endl;
- amrex::Print() << "alpha_C : " << alpha_C << std::endl;
- amrex::Print() << "dynamicViscosity : " << dynamicViscosity << std::endl;
- amrex::Print() << "Cs : " << Cs << std::endl;
- amrex::Print() << "CI : " << CI << std::endl;
- amrex::Print() << "Ce : " << Ce << std::endl;
- amrex::Print() << "Ce at wall : " << Ce_wall << std::endl;
- amrex::Print() << "Ck : " << Ck << std::endl;
- amrex::Print() << "reference theta : " << theta_ref << std::endl;
- amrex::Print() << "sigma_k : " << sigma_k << std::endl;
- amrex::Print() << "Pr_t : " << Pr_t << std::endl;
- amrex::Print() << "Sc_t : " << Sc_t << std::endl;
- amrex::Print() << "horiz spatial_order : " << horiz_spatial_order << std::endl;
- amrex::Print() << "vert spatial_order : " << vert_spatial_order << std::endl;
+ amrex::Print() << "no_substepping : " << no_substepping << std::endl;
+ amrex::Print() << "force_stage1_single_substep : " << force_stage1_single_substep << std::endl;
+ amrex::Print() << "incompressible : " << incompressible << std::endl;
+ amrex::Print() << "use_coriolis : " << use_coriolis << std::endl;
+ amrex::Print() << "use_rayleigh_damping : " << use_rayleigh_damping << std::endl;
+ amrex::Print() << "use_gravity : " << use_gravity << std::endl;
+ amrex::Print() << "rho0_trans : " << rho0_trans << std::endl;
+ amrex::Print() << "alpha_T : " << alpha_T << std::endl;
+ amrex::Print() << "alpha_C : " << alpha_C << std::endl;
+ amrex::Print() << "dynamicViscosity : " << dynamicViscosity << std::endl;
+ amrex::Print() << "Cs : " << Cs << std::endl;
+ amrex::Print() << "CI : " << CI << std::endl;
+ amrex::Print() << "Ce : " << Ce << std::endl;
+ amrex::Print() << "Ce at wall : " << Ce_wall << std::endl;
+ amrex::Print() << "Ck : " << Ck << std::endl;
+ amrex::Print() << "reference theta : " << theta_ref << std::endl;
+ amrex::Print() << "sigma_k : " << sigma_k << std::endl;
+ amrex::Print() << "Pr_t : " << Pr_t << std::endl;
+ amrex::Print() << "Sc_t : " << Sc_t << std::endl;
+ amrex::Print() << "horiz spatial_order : " << horiz_spatial_order << std::endl;
+ amrex::Print() << "vert spatial_order : " << vert_spatial_order << std::endl;
 
  if (abl_driver_type == ABLDriverType::None) {
  amrex::Print() << "ABL Driver Type: " << "None" << std::endl;
@@ -364,6 +379,10 @@ struct SolverChoice {
  // Default prefix
  std::string pp_prefix {"erf"};
 
+ int no_substepping = 0;
+ int force_stage1_single_substep = 1;
+ int incompressible = 0;
+
  bool use_terrain = false;
  bool test_mapfactor = false;
  int terrain_type = 0;

Source/ERF.H

@@ -641,11 +641,6 @@ private:
  static SolverChoice solverChoice;
 
  static int verbose;
- static int use_native_mri;
- static int no_substepping;
- static int force_stage1_single_substep;
-
- static int incompressible;
 
  // mesh refinement
  static std::string coupling_type;

Source/ERF.cpp

@@ -43,12 +43,6 @@ std::string ERF::coupling_type = "OneWay";
 // Dictate verbosity in screen output
 int ERF::verbose = 0;
 
-// Use the native ERF MRI integrator
-int ERF::use_native_mri = 1;
-int ERF::no_substepping = 0;
-int ERF::force_stage1_single_substep = 1;
-int ERF::incompressible = 0;
-
 // Frequency of diagnostic output
 int ERF::sum_interval = -1;
 amrex::Real ERF::sum_per = -1.0;
@@ -788,17 +782,6 @@ ERF::ReadParameters ()
  // Verbosity
  pp.query("v", verbose);
 
- // Use the native ERF MRI integrator
- pp.query("use_native_mri", use_native_mri);
- pp.query("no_substepping", no_substepping);
- pp.query("force_stage1_single_substep", force_stage1_single_substep);
- pp.query("incompressible", incompressible);
-
- // If this is set, it must be even
- if (incompressible != 0 && no_substepping == 0)
- {
- amrex::Abort("If you specify incompressible, you must specific no_substepping");
- }
 
  // Frequency of diagnostic output
  pp.query("sum_interval", sum_interval);

Source/ERF_make_new_level.cpp

@@ -486,9 +486,9 @@ ERF::initialize_integrator(int lev, MultiFab& cons_mf, MultiFab& vel_mf)
  }
 
  mri_integrator_mem[lev] = std::make_unique<MRISplitIntegrator<amrex::Vector<amrex::MultiFab> > >(int_state);
- mri_integrator_mem[lev]->setNoSubstepping(no_substepping);
- mri_integrator_mem[lev]->setIncompressible(incompressible);
- mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(force_stage1_single_substep);
+ mri_integrator_mem[lev]->setNoSubstepping(solverChoice.no_substepping);
+ mri_integrator_mem[lev]->setIncompressible(solverChoice.incompressible);
+ mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep);
 
  physbcs[lev] = std::make_unique<ERFPhysBCFunct> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
  solverChoice.terrain_type, m_bc_extdir_vals, m_bc_neumann_vals,

Source/TimeIntegration/ERF_ComputeTimestep.cpp

@@ -64,6 +64,8 @@ ERF::estTimeStep(int level, long& dt_fast_ratio) const
  &vars_new[level][Vars::yvel],
  &vars_new[level][Vars::zvel]});
 
+ int l_no_substepping = solverChoice.no_substepping;
+
  Real estdt_comp_inv = amrex::ReduceMax(S_new, ccvel, 0,
  [=] AMREX_GPU_HOST_DEVICE (Box const& b,
  Array4<Real const> const& s,
@@ -81,9 +83,19 @@ ERF::estTimeStep(int level, long& dt_fast_ratio) const
  amrex::Real pressure = getPgivenRTh(rhotheta);
  amrex::Real c = std::sqrt(Gamma * pressure / rho);
 
- new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
- ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]),
- ((amrex::Math::abs(u(i,j,k,2))+c)*dxinv[2]), new_comp_dt);
+ // If we are not doing the acoustic substepping, then the z-direction contributes
+ // to the computation of the time step
+ if (l_no_substepping) {
+ new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
+ ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]),
+ ((amrex::Math::abs(u(i,j,k,2))+c)*dxinv[2]), new_comp_dt);
+
+ // If we are doing the acoustic substepping, then the z-direction does not contribute
+ // to the computation of the time step
+ } else {
+ new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]),
+ ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt);
+ }
  });
  return new_comp_dt;
  });
@@ -119,6 +131,7 @@ ERF::estTimeStep(int level, long& dt_fast_ratio) const
  amrex::Print() << "Slow dt at level " << level << ": undefined " << std::endl;
  }
  }
+
  if (fixed_dt > 0.0) {
  amrex::Print() << "Based on cfl of 1.0 " << std::endl;
  amrex::Print() << "Fast dt at level " << level << " would be: " << estdt_comp/cfl << std::endl;
@@ -143,7 +156,7 @@ ERF::estTimeStep(int level, long& dt_fast_ratio) const
  }
 
  // Force time step ratio to be an even value
- if (force_stage1_single_substep) {
+ if (solverChoice.force_stage1_single_substep) {
  if ( dt_fast_ratio%2 != 0) dt_fast_ratio += 1;
  } else {
  if ( dt_fast_ratio%6 != 0) {
@@ -159,7 +172,7 @@ ERF::estTimeStep(int level, long& dt_fast_ratio) const
  if (fixed_dt > 0.0) {
  return fixed_dt;
  } else {
- if (no_substepping) {
+ if (l_no_substepping) {
  return estdt_comp;
  } else {
  return estdt_lowM;

Source/TimeIntegration/ERF_advance_dycore.cpp

@@ -295,7 +295,7 @@ void ERF::advance_dycore(int level,
  mri_integrator.set_post_update(post_update_fun);
 
 #ifdef ERF_USE_POISSON_SOLVE
- if (incompressible) {
+ if (solverChoice.incompressible) {
  mri_integrator.set_slow_rhs_inc(slow_rhs_fun_inc);
  }
 #endif

Source/TimeIntegration/ERF_slow_rhs_post.cpp

@@ -42,7 +42,6 @@ using namespace amrex;
  * @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] incompressible are we running the incompressible algorithm
  */
 
 void erf_slow_rhs_post (int /*level*/,
@@ -69,9 +68,9 @@ void erf_slow_rhs_post (int /*level*/,
  std::unique_ptr<MultiFab>& detJ_new,
  std::unique_ptr<MultiFab>& mapfac_m,
  std::unique_ptr<MultiFab>& mapfac_u,
- std::unique_ptr<MultiFab>& mapfac_v,
+ std::unique_ptr<MultiFab>& mapfac_v
 #if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP))
-  const bool& moist_relax,
+ ,const bool& moist_relax,
  const bool& moist_zero,
  const Real& bdy_time_interval,
  const Real& start_bdy_time,
@@ -81,9 +80,9 @@ void erf_slow_rhs_post (int /*level*/,
  Vector<Vector<FArrayBox>>& bdy_data_xlo,
  Vector<Vector<FArrayBox>>& bdy_data_xhi,
  Vector<Vector<FArrayBox>>& bdy_data_ylo,
- Vector<Vector<FArrayBox>>& bdy_data_yhi,
+ Vector<Vector<FArrayBox>>& bdy_data_yhi
 #endif
- int incompressible)
+ )
 {
  BL_PROFILE_REGION("erf_slow_rhs_post()");
 
@@ -216,7 +215,7 @@ void erf_slow_rhs_post (int /*level*/,
  // We have projected the velocities stored in S_data but we will use
  // the velocities stored in S_scratch to update the scalars, so
  // we need to copy from S_data (projected) into S_scratch
- if (incompressible) {
+ if (solverChoice.incompressible) {
  Box tbx_inc = mfi.nodaltilebox(0);
  Box tby_inc = mfi.nodaltilebox(1);
  Box tbz_inc = mfi.nodaltilebox(2);

Source/TimeIntegration/TI_headers.H

@@ -88,9 +88,9 @@ void erf_slow_rhs_post (int level,
  std::unique_ptr<amrex::MultiFab>& dJ_new,
  std::unique_ptr<amrex::MultiFab>& mapfac_m,
  std::unique_ptr<amrex::MultiFab>& mapfac_u,
- std::unique_ptr<amrex::MultiFab>& mapfac_v,
+ std::unique_ptr<amrex::MultiFab>& mapfac_v
 #if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP))
-  const bool& moist_relax,
+ ,const bool& moist_relax,
  const bool& moist_zero,
  const amrex::Real& bdy_time_interval,
  const amrex::Real& start_bdy_time,
@@ -100,9 +100,9 @@ void erf_slow_rhs_post (int level,
  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_xlo,
  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_xhi,
  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_ylo,
- amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_yhi,
+ amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_yhi
 #endif
- const int incompressible);
+ );
 
 /**
  * Function for computing the fast RHS with no terrain

Source/TimeIntegration/TI_slow_rhs_fun.H

@@ -274,13 +274,13 @@
  Hfx3, Diss,
  fine_geom, solverChoice, m_most, domain_bcs_type_d,
  z_phys_nd_src[level], detJ_cc[level], detJ_cc_new[level],
- mapfac_m[level], mapfac_u[level], mapfac_v[level],
+ mapfac_m[level], mapfac_u[level], mapfac_v[level]
 #if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP))
-  moist_relax, moist_zero, bdy_time_interval, start_bdy_time, new_stage_time,
+ ,moist_relax, moist_zero, bdy_time_interval, start_bdy_time, new_stage_time,
  wrfbdy_width-1, wrfbdy_set_width,
- bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi,
+ bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi
 #endif
- incompressible);
+ );
  } else {
  erf_slow_rhs_post(level, slow_dt,
  S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
@@ -289,13 +289,13 @@
  Hfx3, Diss,
  fine_geom, solverChoice, m_most, domain_bcs_type_d,
  z_phys_nd[level], detJ_cc[level], detJ_cc[level],
- mapfac_m[level], mapfac_u[level], mapfac_v[level],
+ mapfac_m[level], mapfac_u[level], mapfac_v[level]
 #if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP))
-  moist_relax, moist_zero, bdy_time_interval, start_bdy_time, new_stage_time,
+ ,moist_relax, moist_zero, bdy_time_interval, start_bdy_time, new_stage_time,
  wrfbdy_width-1, wrfbdy_set_width,
- bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi,
+ bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi
 #endif
- incompressible);
+ );
  }
  }; // end slow_rhs_fun_post