ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_TI_slow_rhs_pre.H
Go to the documentation of this file.
1 #include "ERF_SrcHeaders.H"
2 
3 /**
4  * Wrapper for calling the routine that creates the slow RHS
5  */
6  auto slow_rhs_fun_pre = [&,zero_d=zero,one_d=one](Vector<MultiFab>& S_rhs,
7  Vector<MultiFab>& S_old,
8  Vector<MultiFab>& S_data,
9  const double old_step_time,
10  const double old_stage_time,
11  const double new_stage_time,
12  const int nrk)
13  {
14  BL_PROFILE("slow_rhs_fun_pre");
15  //
16  // Define primitive variables for all later RK stages
17  // (We have already done this for the first RK step)
18  // Note that it is essential this happen before the call to make_mom_sources
19  // because some of the buoyancy routines use the primitive variables
20  //
21  if (nrk > 0) {
22  int ng_cons = S_data[IntVars::cons].nGrow();
23  cons_to_prim(S_data[IntVars::cons], S_prim, ng_cons);
24  make_pi_stage(S_data[IntVars::cons]);
25  }
26 
27  if (verbose) Print() << std::setprecision(timeprecision)
28  << "Making slow rhs at time " << old_stage_time
29  << " for fast variables advancing from " << old_step_time
30  << " to " << new_stage_time << std::endl;
31 
32  double slow_dt = new_stage_time - old_step_time;
33 
34  const GpuArray<Real, AMREX_SPACEDIM> dxInv = fine_geom.InvCellSizeArray();
35 
36  // *************************************************************************
37  // Set up flux registers if using two_way coupling
38  // *************************************************************************
39  YAFluxRegister* fr_as_crse = nullptr;
40  YAFluxRegister* fr_as_fine = nullptr;
41  if (solverChoice.coupling_type == CouplingType::TwoWay && finest_level > 0) {
42  if (level < finest_level) {
43  fr_as_crse = getAdvFluxReg(level+1);
44  fr_as_crse->reset();
45  }
46  if (level > 0) {
47  fr_as_fine = getAdvFluxReg(level);
48  }
49  }
50 
51  // *************************************************************************
52  // Get multifab pointers
53  // *************************************************************************
54 
55  // Canopy data for mom sources
56  MultiFab* forest_drag = (solverChoice.do_forest_drag) ?
57  m_forest_drag[level]->get_drag_field() : nullptr;
58  MultiFab* frontal_area = (solverChoice.do_forest_drag) ?
59  m_forest_drag[level]->get_frontal_area() : nullptr;
60 
61  // Immersed Forcing
62  MultiFab* terrain_blank = (solverChoice.terrain_type == TerrainType::ImmersedForcing ||
63  solverChoice.buildings_type == BuildingsType::ImmersedForcing) ?
64  terrain_blanking[level].get() : nullptr;
65  MultiFab* terrain_blank_xface = (solverChoice.terrain_type == TerrainType::ImmersedForcing ||
66  solverChoice.buildings_type == BuildingsType::ImmersedForcing) ?
67  terrain_blanking_xface[level].get() : nullptr;
68  MultiFab* terrain_blank_yface = (solverChoice.terrain_type == TerrainType::ImmersedForcing ||
69  solverChoice.buildings_type == BuildingsType::ImmersedForcing) ?
70  terrain_blanking_yface[level].get() : nullptr;
71  MultiFab* terrain_blank_zface = (solverChoice.terrain_type == TerrainType::ImmersedForcing ||
72  solverChoice.buildings_type == BuildingsType::ImmersedForcing) ?
73  terrain_blanking_zface[level].get() : nullptr;
74 
75  // Update the total moisture variable *before* computing sources since this is used in
76  // the buoyancy calculation
77  if (solverChoice.moisture_type != MoistureType::None) {
78  int n_qstate_into_total = micro->Get_Qstate_Moist_Size() - micro->Get_Qstate_Moist_NumConc_Size();
79  make_qt(S_data[IntVars::cons], qt, n_qstate_into_total);
80  }
81 
82  MultiFab p0_to_use, base_to_use;
83  MultiFab *zpn_to_use, *zpc_to_use, *ax_to_use, *ay_to_use, *az_to_use, *dJ_to_use;
84 
85  // Moving terrain
86  std::unique_ptr<MultiFab> z_t_pert;
87  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
88  {
89  z_t_pert = std::make_unique<MultiFab>(S_data[IntVars::zmom].boxArray(), S_data[IntVars::zmom].DistributionMap(), 1, 1);
90  update_terrain_stage(level, old_step_time, old_stage_time, new_stage_time, slow_dt);
91 
92  p0_to_use = MultiFab(base_state_new[level], make_alias, BaseState::p0_comp, 1);
93  base_to_use = MultiFab(base_state_new[level], make_alias, 0, BaseState::num_comps);
94  zpn_to_use = z_phys_nd_src[level].get();
95  zpc_to_use = z_phys_cc_src[level].get();
96  ax_to_use = ax_src[level].get();
97  ay_to_use = ay_src[level].get();
98  az_to_use = az_src[level].get();
99  dJ_to_use = detJ_cc_src[level].get();
100 
101  } else {
102  p0_to_use = MultiFab(base_state[level], make_alias, BaseState::p0_comp, 1);
103  base_to_use = MultiFab(base_state[level], make_alias, 0, BaseState::num_comps);
104  zpn_to_use = z_phys_nd[level].get();
105  zpc_to_use = z_phys_cc[level].get();
106  ax_to_use = ax[level].get();
107  ay_to_use = ay[level].get();
108  az_to_use = az[level].get();
109  dJ_to_use = detJ_cc[level].get();
110  }
111 
112  // Get planar averages from persistent storage for immersed forcing
113  bool l_use_IF = (solverChoice.terrain_type == TerrainType::ImmersedForcing ||
114  solverChoice.buildings_type == BuildingsType::ImmersedForcing);
115  Table1D<Real> r_avg_to_pass = l_use_IF ? r_plane_avg[level].table() : Table1D<Real>();
116  Table1D<Real> t_avg_to_pass = l_use_IF ? t_plane_avg[level].table() : Table1D<Real>();
117 
118  // *****************************************************************************
119  // Construct the source terms for the cell-centered (conserved) variables
120  // *****************************************************************************
121  make_sources(level, nrk, slow_dt, old_stage_time,
122  S_data, S_prim, cc_src, base_state[level], zpc_to_use,
123  xvel_new, yvel_new, zvel_new,
124  qheating_rates[level].get(),
125  terrain_blank, fine_geom, solverChoice,
126  mapfac[level],
127  rhotheta_src[level].get(), rhoqt_src[level].get(),
128  dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sinesq_at_lev,
129  turbPert,
130  r_avg_to_pass, t_avg_to_pass,
131  true);
132 
133  // *****************************************************************************
134  // Add nudging terms on moist variables if we have a non-zero relaxation region
135  // *****************************************************************************
136 #if defined(ERF_USE_NETCDF)
137  if ( solverChoice.use_real_bcs && (level==0) &&
138  (solverChoice.moisture_type != MoistureType::None) )
139  {
140  Real bdy_factor = solverChoice.bdy_nudge_factor;
141  Real l_rdOcp = solverChoice.rdOcp;
142  Real l_c_p = solverChoice.c_p;
143  int moist_nudge_type = solverChoice.bdy_moist_nudge_type;
144  int num_q = micro->Get_Qstate_Moist_Size() - micro->Get_Qstate_Moist_NumConc_Size();
145  AMREX_ALWAYS_ASSERT( (solverChoice.bdy_moist_nudge_type != 2) ||
146  (solverChoice.moisture_type != MoistureType::Morrison_NoIce &&
147  solverChoice.moisture_type != MoistureType::SAM_NoIce) );
148  add_moist_nudging_terms(S_data[IntVars::cons], cc_src, num_q, slow_dt,
149  start_time+old_stage_time,
150  start_bdy_time, final_bdy_time, bdy_time_interval,
151  bdy_factor, real_width, geom[level],
152  bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi,
153  m_r2d, l_c_p, l_rdOcp,
154  solverChoice.moisture_indices,
155  solverChoice.use_wrf_bdy_density, moist_nudge_type);
156  }
157 #endif
158 
159  // Accumulate fixed-leaf-temperature canopy exchange in the staged
160  // source storage. The existing slow-pre and slow-post paths then
161  // apply rho-theta and rho-qv with their normal RK/anelastic weighting.
162  if (solverChoice.do_forest_drag && solverChoice.forest_biophysics &&
163  solverChoice.forest_biophysics_heat && frontal_area) {
164  AddCanopyBiophysicsHeatSources(cc_src,
165  S_data[IntVars::cons],
166  xvel_new, yvel_new, zvel_new,
167  frontal_area, base_to_use,
168  solverChoice);
169  }
170 
171  // *****************************************************************************
172  // Define the pressure gradient
173  // *****************************************************************************
174  make_gradp_pert(level, solverChoice, fine_geom, S_data,
175  p0_to_use, *zpn_to_use, *zpc_to_use,
176  mapfac[level],
177  get_eb(level), gradp[level]);
178 
179  // *****************************************************************************
180  // Define the buoyancy forcing term in the z-direction
181  // *****************************************************************************
182  int num_q = micro->Get_Qstate_Moist_Size() - micro->Get_Qstate_Moist_NumConc_Size();
183  make_buoyancy(level, S_data, S_prim, qt, buoyancy, fine_geom, solverChoice, base_to_use,
184  num_q, get_eb(level), solverChoice.anelastic[level]);
185 
186  // *****************************************************************************
187  // Make remaining (not gradp or buoyancy) momentum sources
188  // *****************************************************************************
189  make_mom_sources(old_stage_time, slow_dt,
190  S_data, zpn_to_use, zpc_to_use, stretched_dz_h[level],
191  xvel_new, yvel_new, zvel_new,
192  xmom_src, ymom_src, zmom_src,
193  base_to_use, forest_drag, terrain_blank,
194  terrain_blank_xface, terrain_blank_yface, terrain_blank_zface,
195  cosPhi_m[level].get(), sinPhi_m[level].get(), fine_geom, solverChoice,
196  mapfac[level],
197  (solverChoice.have_geo_wind_profile) ? d_u_geos[level].data(): nullptr,
198  (solverChoice.have_geo_wind_profile) ? d_v_geos[level].data(): nullptr,
199  dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sinesq_at_lev, d_sinesq_stag_at_lev,
200  d_sponge_ptrs_at_lev,
201  (solverChoice.hindcast_lateral_forcing? &forecast_state_interp[level] : nullptr),
202  input_sounding_data, lsf, lsf_data[level],
203  get_eb(level), true);
204 
205  // *****************************************************************************
206  // Add body sources if doing flow around a body
207  // *****************************************************************************
208  add_thin_body_sources(xmom_src, ymom_src, zmom_src,
209  xflux_imask[level], yflux_imask[level], zflux_imask[level],
210  thin_xforce[level], thin_yforce[level], thin_zforce[level]);
211 
212  // *****************************************************************************
213  // Define RHS for rho, rho_theta and momenta
214  // *****************************************************************************
215  erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data,
216  S_prim, qt, avg_xmom[level], avg_ymom[level], avg_zmom[level],
217  xvel_new, yvel_new, zvel_new,
218  z_t_rk[level], cc_src, xmom_src, ymom_src, zmom_src, buoyancy,
219  (level > 0) ? &zmom_crse_rhs[level] : nullptr,
220  Tau[level], Tau_corr[level], Tau_EB[level],
221  SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss, Hfx3_EB,
222  fine_geom, solverChoice, m_SurfaceLayer, domain_bcs_type_d, domain_bcs_type,
223  *zpn_to_use, *zpc_to_use, *ax_to_use, *ay_to_use, *az_to_use, *dJ_to_use,
224  stretched_dz_d[level], gradp[level],
225  mapfac[level], get_eb(level),
226 #ifdef ERF_USE_EAMXX_SHOC
227  eamxx_shoc_interface[level].get(),
228 #endif
229  native_shoc_driver[level].get(),
230  fr_as_crse, fr_as_fine, &base_to_use,
231  cloud_chamber_config.active ? &cloud_chamber_config : nullptr,
232  cloud_chamber_budget.get());
233 
234  if ((solverChoice.vert_implicit_fac[level][nrk] > zero_d) && solverChoice.implicit_before_substep) {
235  const Real stage_dt = static_cast<Real>(slow_dt);
236  /**
237  * @brief Temporary MultiFab for conserved variable updates.
238  */
239  MultiFab scratch(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(), 2,
240  S_data[IntVars::cons].nGrowVect());
241  MultiFab::Copy(scratch, S_old[IntVars::cons], 0, 0, 2, S_data[IntVars::cons].nGrowVect()); // scratch := S_old (for rho, rhotheta)
242  MultiFab::Saxpy(scratch, stage_dt, S_rhs[IntVars::cons], 0, 0, 2, 0); // scratch := S_old + stage_dt*Src (for rho, rhotheta)
243  scratch.FillBoundary(geom[level].periodicity());
244 
245  /**
246  * @brief Temporary MultiFab for x-momentum updates.
247  */
248  MultiFab scratch_xmom(S_data[IntVars::xmom].boxArray(),
249  S_data[IntVars::xmom].DistributionMap(), 1,
250  S_data[IntVars::xmom].nGrowVect());
251  /**
252  * @brief Temporary MultiFab for y-momentum updates.
253  */
254  MultiFab scratch_ymom(S_data[IntVars::ymom].boxArray(),
255  S_data[IntVars::ymom].DistributionMap(), 1,
256  S_data[IntVars::ymom].nGrowVect());
257 #ifdef ERF_IMPLICIT_W
258  /**
259  * @brief Temporary MultiFab for z-momentum updates.
260  */
261  MultiFab scratch_zmom(S_data[IntVars::zmom].boxArray(),
262  S_data[IntVars::zmom].DistributionMap(), 1,
263  S_data[IntVars::zmom].nGrowVect());
264 #endif
265  if (solverChoice.implicit_momentum_diffusion) {
266  MultiFab::Copy(scratch_xmom, S_old[IntVars::xmom], 0, 0, 1, S_data[IntVars::xmom].nGrowVect()); // scratch := S_old
267  MultiFab::Saxpy(scratch_xmom, stage_dt, S_rhs[IntVars::xmom], 0, 0, 1, 0); // scratch := S_old + stage_dt*Src
268  scratch_xmom.FillBoundary(geom[level].periodicity());
269 
270  MultiFab::Copy(scratch_ymom, S_old[IntVars::ymom], 0, 0, 1, S_data[IntVars::ymom].nGrowVect()); // scratch := S_old
271  MultiFab::Saxpy(scratch_ymom, stage_dt, S_rhs[IntVars::ymom], 0, 0, 1, 0); // scratch := S_old + stage_dt*Src
272  scratch_ymom.FillBoundary(geom[level].periodicity());
273 #ifdef ERF_IMPLICIT_W
274  MultiFab::Copy(scratch_zmom, S_old[IntVars::zmom], 0, 0, 1, S_data[IntVars::zmom].nGrowVect()); // scratch := S_old
275  MultiFab::Saxpy(scratch_zmom, stage_dt, S_rhs[IntVars::zmom], 0, 0, 1, 0); // scratch := S_old + stage_dt*Src
276  scratch_zmom.FillBoundary(geom[level].periodicity());
277 #endif
278  }
279 
280 #include "ERF_ImplicitPre.H"
281 
282  MultiFab::Saxpy(scratch, -one_d, S_old[IntVars::cons], 1, 1, 1, 0); // scratch := (S_new - S_old) (for rhotheta only)
283  scratch.mult(one_d / stage_dt); // scratch := (S_new - S_old) / stage_dt
284  MultiFab::Copy(S_rhs[IntVars::cons], scratch, 1, 1, 1, 0); // slow_rhs := (S_new - S_old) / stage_dt (for rhotheta only)
285 
286  if (solverChoice.implicit_momentum_diffusion) {
287  MultiFab::Saxpy(scratch_xmom, -one_d, S_old[IntVars::xmom], 0, 0, 1, 0); // scratch := (S_new - S_old)
288  scratch_xmom.mult(one_d / stage_dt); // scratch := (S_new - S_old) / stage_dt
289  MultiFab::Copy(S_rhs[IntVars::xmom], scratch_xmom, 0, 0, 1, 0); // slow_rhs := (S_new - S_old) / stage_dt
290 
291  MultiFab::Saxpy(scratch_ymom, -one_d, S_old[IntVars::ymom], 0, 0, 1, 0); // scratch := (S_new - S_old)
292  scratch_ymom.mult(one_d / stage_dt); // scratch := (S_new - S_old) / stage_dt
293  MultiFab::Copy(S_rhs[IntVars::ymom], scratch_ymom, 0, 0, 1, 0); // slow_rhs := (S_new - S_old) / stage_dt
294 #ifdef ERF_IMPLICIT_W
295  MultiFab::Saxpy(scratch_zmom, -one_d, S_old[IntVars::zmom], 0, 0, 1, 0); // scratch := (S_new - S_old)
296  scratch_zmom.mult(one_d / stage_dt); // scratch := (S_new - S_old) / stage_dt
297  MultiFab::Copy(S_rhs[IntVars::zmom], scratch_zmom, 0, 0, 1, 0); // slow_rhs := (S_new - S_old) / stage_dt
298 #endif
299  }
300  }
301 
302 #ifdef ERF_USE_EAMXX_SHOC
303  if (solverChoice.turbChoice[level].uses_eamxx_shoc() && eamxx_shoc_interface[level]) {
304  eamxx_shoc_interface[level]->add_fast_tend(S_rhs);
305  }
306 #endif
307  if (solverChoice.turbChoice[level].uses_native_shoc() && native_shoc_driver[level]) {
308  // Native SHOC applies its increment directly to state; ERF no
309  // longer adds a separate fast RHS contribution.
310  }
311 
312  // *****************************************************************************
313  // Update for moving terrain
314  // *****************************************************************************
315  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
316  {
317  MultiFab r_hse_new (base_state_new[level], make_alias, BaseState::r0_comp, 1);
318  MultiFab p_hse_new (base_state_new[level], make_alias, BaseState::p0_comp, 1);
319  MultiFab pi_hse_new (base_state_new[level], make_alias, BaseState::pi0_comp, 1);
320  MultiFab th_hse_new (base_state_new[level], make_alias, BaseState::th0_comp, 1);
321 
322  MultiFab* r0_new = &r_hse_new;
323  MultiFab* p0_new = &p_hse_new;
324  MultiFab* pi0_new = &pi_hse_new;
325  MultiFab* th0_new = &th_hse_new;
326 
327  // We define and evolve (rho theta)_0 in order to re-create p_0 in a way that is consistent
328  // with our update of (rho theta) but does NOT maintain dp_0 / dz = -rho_0 g. This is why
329  // we no longer discretize the vertical pressure gradient in perturbational form.
330  MultiFab rt0(p0->boxArray(),p0->DistributionMap(),1,1);
331  MultiFab rt0_new(p0->boxArray(),p0->DistributionMap(),1,1);
332  MultiFab r0_temp(p0->boxArray(),p0->DistributionMap(),1,1);
333 
334  // Remember this does NOT maintain dp_0 / dz = -rho_0 g, so we can no longer
335  // discretize the vertical pressure gradient in perturbational form.
336  AMREX_ALWAYS_ASSERT(solverChoice.advChoice.dycore_horiz_adv_type == AdvType::Centered_2nd);
337  AMREX_ALWAYS_ASSERT(solverChoice.advChoice.dycore_vert_adv_type == AdvType::Centered_2nd);
338 
339  double dt_base = new_stage_time - old_step_time;
340 
341  const Real l_rdOcp = solverChoice.rdOcp;
342 
343 #ifdef _OPENMP
344 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
345 #endif
346  for ( MFIter mfi(*p0,TilingIfNotGPU()); mfi.isValid(); ++mfi)
347  {
348  const Array4<Real > rt0_arr = rt0.array(mfi);
349  const Array4<Real > rt0_tmp_arr = rt0_new.array(mfi);
350 
351  const Array4<Real const> r0_arr = r0->const_array(mfi);
352  const Array4<Real > r0_new_arr = r0_new->array(mfi);
353  const Array4<Real > r0_tmp_arr = r0_temp.array(mfi);
354 
355  const Array4<Real const> p0_arr = p0->const_array(mfi);
356  const Array4<Real > p0_new_arr = p0_new->array(mfi);
357  const Array4<Real > pi0_new_arr = pi0_new->array(mfi);
358  const Array4<Real > th0_new_arr = th0_new->array(mfi);
359 
360  const Array4<Real >& z_t_arr = z_t_rk[level]->array(mfi);
361 
362  const Array4<Real const>& dJ_old_arr = detJ_cc[level]->const_array(mfi);
363  const Array4<Real const>& dJ_new_arr = detJ_cc_new[level]->const_array(mfi);
364  const Array4<Real const>& dJ_src_arr = detJ_cc_src[level]->const_array(mfi);
365 
366  Box gbx = mfi.growntilebox({1,1,1});
367  amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
368  {
369  rt0_arr(i,j,k) = getRhoThetagivenP(p0_arr(i,j,k));
370  rt0_tmp_arr(i,j,k) = getRhoThetagivenP(p0_new_arr(i,j,k));
371  r0_tmp_arr(i,j,k) = r0_new_arr(i,j,k);
372  });
373 
374  Box gbx2 = mfi.growntilebox({1,1,0});
375  amrex::ParallelFor(gbx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
376  {
377  Real zflux_r_lo = -z_t_arr(i,j,k ) * myhalf * (r0_tmp_arr(i,j,k) + r0_tmp_arr(i,j,k-1));
378  Real zflux_r_hi = -z_t_arr(i,j,k+1) * myhalf * (r0_tmp_arr(i,j,k) + r0_tmp_arr(i,j,k+1));
379 
380  Real zflux_rt_lo = zflux_r_lo * myhalf * (rt0_tmp_arr(i,j,k)/r0_tmp_arr(i,j,k) + rt0_tmp_arr(i,j,k-1)/r0_tmp_arr(i,j,k-1));
381  Real zflux_rt_hi = zflux_r_hi * myhalf * (rt0_tmp_arr(i,j,k)/r0_tmp_arr(i,j,k) + rt0_tmp_arr(i,j,k+1)/r0_tmp_arr(i,j,k+1));
382 
383  Real invdetJ = one / dJ_src_arr(i,j,k);
384 
385  Real src_r = - invdetJ * ( zflux_r_hi - zflux_r_lo ) * dxInv[2];
386  Real src_rt = - invdetJ * ( zflux_rt_hi - zflux_rt_lo ) * dxInv[2];
387 
388  Real dt_b = static_cast<Real>(dt_base);
389  Real rho0_new = dJ_old_arr(i,j,k) * r0_arr(i,j,k) + dt_b * dJ_src_arr(i,j,k) * src_r;
390  Real rt0_tmp_new = dJ_old_arr(i,j,k) * rt0_arr(i,j,k) + dt_b * dJ_src_arr(i,j,k) * src_rt;
391 
392  r0_new_arr(i,j,k) = rho0_new / dJ_new_arr(i,j,k);
393  rt0_tmp_new /= dJ_new_arr(i,j,k);
394 
395  p0_new_arr(i,j,k) = getPgivenRTh(rt0_tmp_new);
396  pi0_new_arr(i,j,k) = getExnergivenRTh(rt0_tmp_new, l_rdOcp);
397  th0_new_arr(i,j,k) = rt0_tmp_new / r0_new_arr(i,j,k);
398  });
399  } // MFIter
400  r0_new->FillBoundary(fine_geom.periodicity());
401  p0_new->FillBoundary(fine_geom.periodicity());
402  th0_new->FillBoundary(fine_geom.periodicity());
403  }
404 
405 #ifdef ERF_USE_NETCDF
406  // Populate RHS for relaxation zones if using real bcs
407  if (solverChoice.use_real_bcs && (level == 0)) {
408  const Real bdy_factor = solverChoice.bdy_nudge_factor;
409  if (real_width>0) {
410  //
411  // Note that old_stage_time is elapsed time, but (start_time+old_stage_time) is total time
412  // start_bdy_time and final_bdy_time are total time
413  //
414  double total_time = start_time + old_stage_time;
415  Real l_rdOcp = solverChoice.rdOcp;
416  Real l_c_p = solverChoice.c_p;
417  realbdy_compute_interior_ghost_rhs(total_time, slow_dt,
418  start_bdy_time, final_bdy_time, bdy_time_interval,
419  bdy_factor, real_width, fine_geom,
420  S_rhs, S_data,
421  bdy_data_xlo, bdy_data_xhi,
422  bdy_data_ylo, bdy_data_yhi,
423  m_r2d, l_c_p, l_rdOcp,
424  solverChoice.use_wrf_bdy_density,
425  solverChoice.bdy_rho_nudge_factor);
426  }
427  }
428 #endif
429  }; // end slow_rhs_fun_pre
add_thin_body_sources
void add_thin_body_sources(MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, std::unique_ptr< iMultiFab > &xflux_imask_lev, std::unique_ptr< iMultiFab > &yflux_imask_lev, std::unique_ptr< iMultiFab > &zflux_imask_lev, std::unique_ptr< MultiFab > &thin_xforce_lev, std::unique_ptr< MultiFab > &thin_yforce_lev, std::unique_ptr< MultiFab > &thin_zforce_lev)
Definition: ERF_AddThinBodySources.cpp:27
AddCanopyBiophysicsHeatSources
AMREX_INLINE void AddCanopyBiophysicsHeatSources(amrex::MultiFab &cell_source, const amrex::MultiFab &S_data, const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &zvel, const amrex::MultiFab *frontal_area, const amrex::MultiFab &base_state, const SolverChoice &solver_choice)
Definition: ERF_CanopyBiophysics.H:28
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
getRhoThetagivenP
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:172
getExnergivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:156
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
ERF_ImplicitPre.H
AdvType::Centered_2nd
@ Centered_2nd
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
get
pp get("wavelength", wavelength)
realbdy_compute_interior_ghost_rhs
void realbdy_compute_interior_ghost_rhs(const double &time, const double &delta_t_d, const double &start_bdy_time, const double &final_bdy_time, const double &bdy_time_interval, const Real &nudge_factor, int width, const Geometry &geom, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_cur_data, Vector< Vector< FArrayBox >> &bdy_data_xlo, Vector< Vector< FArrayBox >> &bdy_data_xhi, Vector< Vector< FArrayBox >> &bdy_data_ylo, Vector< Vector< FArrayBox >> &bdy_data_yhi, std::unique_ptr< ReadBndryPlanes > &m_r2d, const Real &c_p, const Real &rdOcp, const bool use_wrf_bdy_density, const Real &bdy_rho_nudge_factor)
Definition: ERF_InteriorGhostCells.cpp:161
make_buoyancy
void make_buoyancy(int lev, const Vector< MultiFab > &S_data, const MultiFab &S_prim, const MultiFab &qt, MultiFab &buoyancy, const Geometry geom, const SolverChoice &solverChoice, const MultiFab &base_state, const int n_qstate, const eb_ &ebfact, const int anelastic)
Definition: ERF_MakeBuoyancy.cpp:32
make_gradp_pert
void make_gradp_pert(int level, const SolverChoice &solverChoice, const Geometry &geom, Vector< MultiFab > &S_data, const MultiFab &p0, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, Vector< MultiFab > &gradp)
Definition: ERF_MakeGradP.cpp:28
make_mom_sources
void make_mom_sources(double time_d, double dt, const Vector< MultiFab > &S_data, const MultiFab *z_phys_nd, const 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 *terrain_blank_xface, MultiFab *terrain_blank_yface, MultiFab *terrain_blank_zface, 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 amrex::Real *d_sinesq_at_lev, const amrex::Real *d_sinesq_stag_at_lev, const Vector< Real * > d_sponge_ptrs_at_lev, const Vector< MultiFab > *forecast_state_at_lev, InputSoundingData &input_sounding_data, LargeScaleForcingData &lsf_data, std::unique_ptr< amrex::MultiFab > &lsf_tendencies, const eb_ &ebfact, bool is_slow_step)
Definition: ERF_MakeMomSources.cpp:38
make_sources
void make_sources(int level, int, double dt, double time_d, const Vector< MultiFab > &S_data, const MultiFab &S_prim, MultiFab &source, const MultiFab &base_state, const MultiFab *z_phys_cc, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, const MultiFab *qheating_rates, MultiFab *terrain_blank, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &mapfac, const MultiFab *rhotheta_src, const MultiFab *rhoqt_src, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const Real *d_sinesq_at_lev, TurbulentPerturbation &turbPert, const Table1D< Real > r_plane_avg, const Table1D< Real > t_plane_avg, bool is_slow_step)
Definition: ERF_MakeSources.cpp:35
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
erf_slow_rhs_pre
void erf_slow_rhs_pre(int level, int finest_level, int nrk, double dt, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_data, const MultiFab &S_prim, const MultiFab &qt, MultiFab &avg_xmom, MultiFab &avg_ymom, MultiFab &avg_zmom, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, std::unique_ptr< MultiFab > &z_t_mf, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const MultiFab &buoyancy, const MultiFab *zmom_crse_rhs, Vector< std::unique_ptr< MultiFab >> &Tau_lev, Vector< std::unique_ptr< MultiFab >> &Tau_corr_lev, Vector< Vector< std::unique_ptr< MultiFab >>> &Tau_EB, MultiFab *SmnSmn, MultiFab *eddyDiffs, MultiFab *Hfx1, MultiFab *Hfx2, MultiFab *Hfx3, MultiFab *Q1fx1, MultiFab *Q1fx2, MultiFab *Q1fx3, MultiFab *Q2fx3, MultiFab *Diss, MultiFab *Hfx3_EB, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, const Gpu::DeviceVector< BCRec > &domain_bcs_type_d, const Vector< BCRec > &domain_bcs_type_h, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, const MultiFab &ax, const MultiFab &ay, const MultiFab &az, const MultiFab &detJ, Gpu::DeviceVector< Real > &stretched_dz_d, Vector< MultiFab > &gradp, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, ShocDriver *native_shoc_lev, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine, const MultiFab *cloud_chamber_base_state, const erf_cloud_chamber::Config *cloud_chamber_config, CloudChamberBudget *cloud_budget)
Definition: ERF_SlowRhsPre.cpp:68
ERF_SrcHeaders.H
slow_rhs_fun_pre
auto slow_rhs_fun_pre
Definition: ERF_TI_slow_rhs_pre.H:6
make_pi_stage
auto make_pi_stage
Definition: ERF_TI_utils.H:4
update_terrain_stage
auto update_terrain_stage
Definition: ERF_TI_utils.H:125
cons_to_prim
void cons_to_prim(const MultiFab &cons_state, MultiFab &S_prim, int ng)
Definition: ERF_Utils.cpp:13
make_qt
void make_qt(const MultiFab &cons_state, MultiFab &qt, int n_qstate_into_total)
Definition: ERF_Utils.cpp:108
BaseState::num_comps
@ num_comps
Definition: ERF_IndexDefines.H:81
BaseState::pi0_comp
@ pi0_comp
Definition: ERF_IndexDefines.H:78
BaseState::p0_comp
@ p0_comp
Definition: ERF_IndexDefines.H:77
BaseState::th0_comp
@ th0_comp
Definition: ERF_IndexDefines.H:79
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:76
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:234
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:232
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:235
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:233
MicVar_Kess::qt
@ qt
Definition: ERF_Kessler.H:29
module_model_constants::p0
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40
mp_wdm6::r0
real(kind=kind_phys), parameter, private r0
Definition: ERF_module_mp_wdm6.F90:75