ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_TI_substep_fun.H
Go to the documentation of this file.
1 #include "ERF_SrcHeaders.H"
2 
3 /**
4  * Wrapper for calling the routine that creates the fast RHS
5  */
6 auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
7  Vector<MultiFab>& S_slow_rhs,
8  const Vector<MultiFab>& S_old,
9  Vector<MultiFab>& S_stage,
10  Vector<MultiFab>& S_data,
11  Vector<MultiFab>& S_scratch,
12  const Real dtau,
13  const Real inv_fac,
14  const Real old_substep_time,
15  const Real new_substep_time)
16  {
17  BL_PROFILE("fast_rhs_fun");
18  if (verbose) amrex::Print() << "Fast time integration at level " << level
19  << std::setprecision(timeprecision)
20  << " from " << old_substep_time << " to " << new_substep_time
21  << " with dt = " << dtau << std::endl;
22 
23  // Define beta_s here so that it is consistent between where we make the fast coefficients
24  // and where we use them
25  // Per p2902 of Klemp-Skamarock-Dudhia-2007
26  // beta_s = -1.0 : fully explicit
27  // beta_s = 1.0 : fully implicit
28  Real beta_s;
29  if (solverChoice.substepping_type[level] == SubsteppingType::Implicit) {
30  beta_s = 0.1;
31  } else { // Fully explicit
32  beta_s = -1.0;
33  }
34 
35  // Canopy data for mom sources
36  MultiFab* forest_drag = (solverChoice.do_forest_drag) ?
37  m_forest_drag[level]->get_drag_field() : nullptr;
38 
39  // Immersed Forcing
40  MultiFab* terrain_blank = (solverChoice.terrain_type == TerrainType::ImmersedForcing) ?
41  terrain_blanking[level].get() : nullptr;
42 
43  Vector<MultiFab> Svec_to_use;
44  if (fast_step == 0) {
45  // If this is the first substep we pass in S_old as the previous step's solution
46  Svec_to_use.push_back(MultiFab(S_old[IntVars::cons],make_alias,0,S_old[IntVars::cons].nComp()));
47  Svec_to_use.push_back(MultiFab(S_old[IntVars::xmom],make_alias,0,1));
48  Svec_to_use.push_back(MultiFab(S_old[IntVars::ymom],make_alias,0,1));
49  Svec_to_use.push_back(MultiFab(S_old[IntVars::zmom],make_alias,0,1));
50  } else {
51  // If this is not the first substep we pass in S_data as the previous step's solution
52  Svec_to_use.push_back(MultiFab(S_data[IntVars::cons],make_alias,0,S_data[IntVars::cons].nComp()));
53  Svec_to_use.push_back(MultiFab(S_data[IntVars::xmom],make_alias,0,1));
54  Svec_to_use.push_back(MultiFab(S_data[IntVars::ymom],make_alias,0,1));
55  Svec_to_use.push_back(MultiFab(S_data[IntVars::zmom],make_alias,0,1));
56  }
57 
58  make_sources(level, nrk, dtau, old_substep_time,
59  Svec_to_use, S_prim, cc_src, base_state[level], z_phys_cc[level].get(),
60  qheating_rates[level].get(),
61  terrain_blank, fine_geom, solverChoice,
62  mapfac[level],
63  dptr_rhotheta_src, dptr_rhoqt_src,
64  dptr_wbar_sub, d_rayleigh_ptrs_at_lev,
65  input_sounding_data, turbPert, false);
66 
67  // *************************************************************************
68  // Set up flux registers if using two_way coupling
69  // *************************************************************************
70  const bool l_reflux = ( (solverChoice.coupling_type == CouplingType::TwoWay) && (nrk == 2) && (finest_level > 0) );
71 
72  YAFluxRegister* fr_as_crse = nullptr;
73  YAFluxRegister* fr_as_fine = nullptr;
74  if (l_reflux) {
75  if (level < finest_level) {
76  fr_as_crse = getAdvFluxReg(level+1);
77  }
78  if (level > 0) {
79  fr_as_fine = getAdvFluxReg(level);
80  }
81  }
82 
83  MultiFab base_to_use;
84  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh ) {
85  base_to_use = MultiFab(base_state_new[level], make_alias, 0, BaseState::num_comps);
86  } else {
87  base_to_use = MultiFab(base_state[level], make_alias, 0, BaseState::num_comps);
88  }
89 
90  make_mom_sources(old_substep_time, Svec_to_use,
91  *z_phys_nd[level].get(), *z_phys_cc[level].get(), stretched_dz_h[level],
92  xvel_new, yvel_new, zvel_new,
93  xmom_src, ymom_src, zmom_src,
94  base_to_use, forest_drag, terrain_blank,
95  cosPhi_m[level].get(), sinPhi_m[level].get(), fine_geom, solverChoice,
96  mapfac[level],
97  (solverChoice.have_geo_wind_profile) ? d_u_geos[level].data() : nullptr,
98  (solverChoice.have_geo_wind_profile) ? d_v_geos[level].data() : nullptr,
99  dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev,
100  input_sounding_data, false);
101 
102  // Moving terrain
103  std::unique_ptr<MultiFab> z_t_pert;
104  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
105  {
106  z_t_pert = std::make_unique<MultiFab>(S_data[IntVars::zmom].boxArray(), S_data[IntVars::zmom].DistributionMap(), 1, 1);
107  update_terrain_substep(level, old_substep_time, new_substep_time, dtau, S_data, z_t_pert);
108  }
109 
110  if ( (fast_step == 0) || (solverChoice.terrain_type == TerrainType::MovingFittedMesh) )
111  {
112  make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom,
113  l_use_moisture, SolverChoice::mesh_type, solverChoice.gravity, solverChoice.c_p,
114  detJ_cc[level], r0, pi0, dtau, beta_s, phys_bc_type);
115  }
116 
117  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
118  {
119  erf_fast_rhs_MT(fast_step, nrk, level, finest_level,
120  S_slow_rhs, Svec_to_use, S_stage, S_prim, qt, pi_stage, fast_coeffs,
121  S_data, S_scratch,
122  cc_src, xmom_src, ymom_src, zmom_src,
123  fine_geom,
124  solverChoice.gravity, solverChoice.use_lagged_delta_rt,
125  z_t_rk[level], z_t_pert.get(),
126  z_phys_nd[level], z_phys_nd_new[level], z_phys_nd_src[level],
127  detJ_cc[level], detJ_cc_new[level], detJ_cc_src[level],
128  dtau, beta_s, inv_fac, mapfac[level],
129  fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
130 
131  } else if (solverChoice.mesh_type != MeshType::ConstantDz) {
132 
133  erf_fast_rhs_T(fast_step, nrk, level, finest_level,
134  S_slow_rhs, Svec_to_use, S_stage, S_prim, qt, pi_stage, fast_coeffs,
135  S_data, S_scratch,
136  cc_src, xmom_src, ymom_src, zmom_src,
137  fine_geom, solverChoice.gravity,
138  z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac,
139  mapfac[level],
140  fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
141 
142  } else {
143 
144  erf_fast_rhs_N(fast_step, nrk, level, finest_level,
145  S_slow_rhs, Svec_to_use, S_stage, S_prim, qt, pi_stage, fast_coeffs,
146  S_data, S_scratch,
147  cc_src, xmom_src, ymom_src, zmom_src,
148  fine_geom, solverChoice.gravity,
149  dtau, beta_s, inv_fac, mapfac[level],
150  fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
151  }
152 
153  // Even if we update all the conserved variables we don't need
154  // to fillpatch the slow ones every acoustic substep
155 
156  // NOTE: Numerical diffusion is tested on in FillPatchIntermediate and dictates the size of the
157  // box over which VelocityToMomentum is computed. V2M requires one more ghost cell be
158  // filled for rho than velocity. This logical condition ensures we fill enough ghost cells
159  // when use_num_diff is true.
160  int ng_cons = S_data[IntVars::cons].nGrowVect().max() - 1;
161  int ng_vel = S_data[IntVars::xmom].nGrowVect().max();
162  if (!solverChoice.use_num_diff) {
163  ng_cons = 1;
164  ng_vel = 1;
165  }
166  apply_bcs(S_data, new_substep_time, ng_cons, ng_vel, fast_only=true, vel_and_mom_synced=false);
167  };
void erf_fast_rhs_MT(int step, int, int level, int finest_level, Vector< MultiFab > &S_slow_rhs, const Vector< MultiFab > &S_prev, Vector< MultiFab > &S_stg_data, const MultiFab &S_stg_prim, const MultiFab &qt, const MultiFab &pi_stage, const MultiFab &fast_coeffs, Vector< MultiFab > &S_data, Vector< MultiFab > &S_scratch, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const Geometry geom, const Real gravity, const bool use_lagged_delta_rt, std::unique_ptr< MultiFab > &z_t_rk, const MultiFab *z_t_pert, std::unique_ptr< MultiFab > &z_phys_nd_old, std::unique_ptr< MultiFab > &z_phys_nd_new, std::unique_ptr< MultiFab > &z_phys_nd_stg, std::unique_ptr< MultiFab > &detJ_cc_old, std::unique_ptr< MultiFab > &detJ_cc_new, std::unique_ptr< MultiFab > &detJ_cc_stg, const Real dtau, const Real beta_s, const Real facinv, Vector< std::unique_ptr< MultiFab >> &mapfac, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine, bool l_use_moisture, bool l_reflux, bool)
Definition: ERF_FastRhs_MT.cpp:47
void erf_fast_rhs_N(int step, int nrk, int level, int finest_level, Vector< MultiFab > &S_slow_rhs, const Vector< MultiFab > &S_prev, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &qt, const MultiFab &pi_stage, const MultiFab &fast_coeffs, Vector< MultiFab > &S_data, Vector< MultiFab > &S_scratch, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const Geometry geom, const Real gravity, const Real dtau, const Real beta_s, const Real facinv, Vector< std::unique_ptr< MultiFab >> &mapfac, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine, bool l_use_moisture, bool l_reflux, bool l_implicit_substepping)
Definition: ERF_FastRhs_N.cpp:38
void erf_fast_rhs_T(int step, int, int level, int finest_level, Vector< MultiFab > &S_slow_rhs, const Vector< MultiFab > &S_prev, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &qt, const MultiFab &pi_stage, const MultiFab &fast_coeffs, Vector< MultiFab > &S_data, Vector< MultiFab > &S_scratch, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const Geometry geom, const Real gravity, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &detJ_cc, const Real dtau, const Real beta_s, const Real facinv, Vector< std::unique_ptr< MultiFab >> &mapfac, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine, bool l_use_moisture, bool l_reflux, bool)
Definition: ERF_FastRhs_T.cpp:40
void make_fast_coeffs(int, MultiFab &fast_coeffs, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &pi_stage, const amrex::Geometry geom, bool l_use_moisture, MeshType mesh_type, Real gravity, Real c_p, std::unique_ptr< MultiFab > &detJ_cc, const MultiFab *r0, const MultiFab *pi0, Real dtau, Real beta_s, amrex::GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type)
Definition: ERF_MakeFastCoeffs.cpp:26
void make_mom_sources(Real time, const Vector< MultiFab > &S_data, MultiFab &z_phys_nd, MultiFab &z_phys_cc, Vector< Real > &stretched_dz_h, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &wvel, MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, const MultiFab &base_state, MultiFab *forest_drag, MultiFab *terrain_blank, MultiFab *cosPhi_mf, MultiFab *sinPhi_mf, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &, const Real *dptr_u_geos, const Real *dptr_v_geos, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const Vector< Real * > d_sponge_ptrs_at_lev, InputSoundingData &input_sounding_data, bool is_slow_step)
Definition: ERF_MakeMomSources.cpp:33
void make_sources(int level, int, Real dt, Real time, const Vector< MultiFab > &S_data, const MultiFab &S_prim, MultiFab &source, const MultiFab &base_state, const MultiFab *z_phys_cc, const MultiFab *qheating_rates, MultiFab *terrain_blank, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &mapfac, const Real *dptr_rhotheta_src, const Real *dptr_rhoqt_src, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, InputSoundingData &input_sounding_data, TurbulentPerturbation &turbPert, bool is_slow_step)
Definition: ERF_MakeSources.cpp:31
auto fast_rhs_fun
Definition: ERF_TI_substep_fun.H:6
auto apply_bcs
Definition: ERF_TI_utils.H:71
auto update_terrain_substep
Definition: ERF_TI_utils.H:236
@ num_comps
Definition: ERF_IndexDefines.H:68
@ ymom
Definition: ERF_IndexDefines.H:160
@ cons
Definition: ERF_IndexDefines.H:158
@ zmom
Definition: ERF_IndexDefines.H:161
@ xmom
Definition: ERF_IndexDefines.H:159
@ qt
Definition: ERF_Kessler.H:27
static MeshType mesh_type
Definition: ERF_DataStruct.H:755