ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_TI_slow_rhs_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 slow RHS
5  */
6  auto slow_rhs_fun_pre = [&](Vector<MultiFab>& S_rhs,
7  Vector<MultiFab>& S_old,
8  Vector<MultiFab>& S_data,
9  Vector<MultiFab>& S_scratch,
10  const Real old_step_time,
11  const Real old_stage_time,
12  const Real new_stage_time,
13  const int nrk)
14  {
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], ng_cons);
24  }
25 
26  BL_PROFILE("slow_rhs_fun_pre");
27  if (verbose) Print() << "Making slow rhs at time " << old_stage_time << " for fast variables advancing from " <<
28  old_step_time << " to " << new_stage_time << std::endl;
29 
30  Real slow_dt = new_stage_time - old_step_time;
31 
32  int n_qstate = micro->Get_Qstate_Size();
33 
34  // *************************************************************************
35  // Set up flux registers if using two_way coupling
36  // *************************************************************************
37  YAFluxRegister* fr_as_crse = nullptr;
38  YAFluxRegister* fr_as_fine = nullptr;
39  if (solverChoice.coupling_type == CouplingType::TwoWay) {
40  if (level < finest_level) {
41  fr_as_crse = getAdvFluxReg(level+1);
42  fr_as_crse->reset();
43  }
44  if (level > 0) {
45  fr_as_fine = getAdvFluxReg(level);
46  }
47  }
48 
49  Real* dptr_u_geos = solverChoice.have_geo_wind_profile ? d_u_geos[level].data(): nullptr;
50  Real* dptr_v_geos = solverChoice.have_geo_wind_profile ? d_v_geos[level].data(): nullptr;
51 
52  // Construct the source terms for the cell-centered (conserved) variables
53  make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level],
54 #if defined(ERF_USE_RRTMGP)
55  qheating_rates[level].get(),
56 #endif
57  fine_geom, solverChoice,
58  mapfac_u[level], mapfac_v[level], mapfac_m[level],
59  dptr_rhotheta_src, dptr_rhoqt_src,
60  dptr_wbar_sub, d_rayleigh_ptrs_at_lev,
61  input_sounding_data, turbPert);
62 
63  // Canopy data for mom sources
64  MultiFab* forest_drag = nullptr;
65  if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); }
66  // Immersed Forcing
67  MultiFab* terrain_blank = nullptr;
68  if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); }
69 
70  // Moving terrain
71  if ( solverChoice.terrain_type == TerrainType::Moving )
72  {
73  // Note that the "old" and "new" metric terms correspond to
74  // t^n and the RK stage (either t^*, t^** or t^{n+1} that this source
75  // will be used to advance to
76 
77  // The "src" metric terms correspond to the time at which we are evaluating the source here,
78  // aka old_stage_time
79 
80  if (verbose) Print() << "Re-making old geometry at old time : " << old_step_time << std::endl;
81  prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_step_time);
82  init_terrain_grid (level,fine_geom,*z_phys_nd[level], zlevels_stag[level], phys_bc_type);
83  make_J (fine_geom,*z_phys_nd[level], *detJ_cc[level]);
84  make_areas (fine_geom,*z_phys_nd[level], *ax[level], *ay[level], *az[level]);
85 
86  if (verbose) Print() << "Making src geometry at old_stage_time: " << old_stage_time << std::endl;
87  prob->init_custom_terrain(fine_geom,*z_phys_nd_src[level],old_stage_time);
88  init_terrain_grid (level,fine_geom,*z_phys_nd_src[level], zlevels_stag[level], phys_bc_type);
89  make_J (fine_geom,*z_phys_nd_src[level], *detJ_cc_src[level]);
90  make_areas (fine_geom,*z_phys_nd_src[level], *ax_src[level], *ay_src[level], *az_src[level]);
91 
92  if (verbose) Print() << "Making new geometry at new_stage_time: " << new_stage_time << std::endl;
93  prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_stage_time);
94  init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag[level], phys_bc_type);
95  make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]);
96  make_areas (fine_geom,*z_phys_nd_new[level], *ax_new[level], *ay_new[level], *az_new[level]);
97 
98  Real inv_dt = 1./slow_dt;
99 
100 #ifdef _OPENMP
101 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
102 #endif
103  for (MFIter mfi(*z_t_rk[level],TilingIfNotGPU()); mfi.isValid(); ++mfi)
104  {
105  Box gbx = mfi.growntilebox(IntVect(1,1,0));
106 
107  const Array4<Real >& z_t_arr = z_t_rk[level]->array(mfi);
108  const Array4<Real const>& z_nd_new_arr = z_phys_nd_new[level]->const_array(mfi);
109  const Array4<Real const>& z_nd_old_arr = z_phys_nd[level]->const_array(mfi);
110 
111  // Loop over horizontal plane
112  amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
113  {
114  // Evaluate between RK stages assuming the geometry is linear between old and new time
115  z_t_arr(i,j,k) = 0.25 * inv_dt * (z_nd_new_arr(i+1,j+1,k) - z_nd_old_arr(i+1,j+1,k)
116  +z_nd_new_arr(i ,j+1,k) - z_nd_old_arr( i,j+1,k)
117  +z_nd_new_arr(i+1,j ,k) - z_nd_old_arr(i+1,j ,k)
118  +z_nd_new_arr(i ,j ,k) - z_nd_old_arr(i ,j ,k));
119  });
120 
121  } // mfi
122 
123  MultiFab r_hse_new (base_state_new[level], make_alias, BaseState::r0_comp, 1);
124  MultiFab p_hse_new (base_state_new[level], make_alias, BaseState::p0_comp, 1);
125  MultiFab pi_hse_new (base_state_new[level], make_alias, BaseState::pi0_comp, 1);
126  MultiFab th_hse_new (base_state_new[level], make_alias, BaseState::th0_comp, 1);
127 
128  MultiFab* r0_new = &r_hse_new;
129  MultiFab* p0_new = &p_hse_new;
130  MultiFab* pi0_new = &pi_hse_new;
131  MultiFab* th0_new = &th_hse_new;
132 
133  make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim,
134  z_phys_nd[level], z_phys_cc[level],
135  xvel_new, yvel_new, zvel_new,
136  xmom_src, ymom_src, zmom_src,
137  base_state_new[level], forest_drag, terrain_blank, fine_geom, solverChoice,
138  mapfac_m[level], mapfac_u[level], mapfac_v[level],
139  dptr_u_geos, dptr_v_geos, dptr_wbar_sub,
140  d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev,
141  input_sounding_data, n_qstate);
142 
143  erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch,
144  xvel_new, yvel_new, zvel_new,
145  z_t_rk[level], cc_src, xmom_src, ymom_src, zmom_src,
146  (level > 0) ? &zmom_crse_rhs[level] : nullptr,
147  Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(),
148  Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(),
149  Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
150  fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
151  z_phys_nd_src[level], ax_src[level], ay_src[level], az_src[level], detJ_cc_src[level], p0_new,
152  pp_inc[level],
153  mapfac_m[level], mapfac_u[level], mapfac_v[level],
154 #ifdef ERF_USE_EB
155  EBFactory(level),
156 #endif
157  fr_as_crse, fr_as_fine);
158 
159  add_thin_body_sources(xmom_src, ymom_src, zmom_src,
160  xflux_imask[level], yflux_imask[level], zflux_imask[level],
161  thin_xforce[level], thin_yforce[level], thin_zforce[level]);
162 
163  // We define and evolve (rho theta)_0 in order to re-create p_0 in a way that is consistent
164  // with our update of (rho theta) but does NOT maintain dp_0 / dz = -rho_0 g. This is why
165  // we no longer discretize the vertical pressure gradient in perturbational form.
166  MultiFab rt0(p0->boxArray(),p0->DistributionMap(),1,1);
167  MultiFab rt0_new(p0->boxArray(),p0->DistributionMap(),1,1);
168  MultiFab r0_temp(p0->boxArray(),p0->DistributionMap(),1,1);
169 
170  // Remember this does NOT maintain dp_0 / dz = -rho_0 g, so we can no longer
171  // discretize the vertical pressure gradient in perturbational form.
172  AMREX_ALWAYS_ASSERT(solverChoice.advChoice.dycore_horiz_adv_type == AdvType::Centered_2nd);
173  AMREX_ALWAYS_ASSERT(solverChoice.advChoice.dycore_vert_adv_type == AdvType::Centered_2nd);
174 
175  Real dt_base = (new_stage_time - old_step_time);
176 
177  const GpuArray<Real, AMREX_SPACEDIM> dxInv = fine_geom.InvCellSizeArray();
178 
179  const Real l_rdOcp = solverChoice.rdOcp;
180 
181 #ifdef _OPENMP
182 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
183 #endif
184  for ( MFIter mfi(*p0,TilingIfNotGPU()); mfi.isValid(); ++mfi)
185  {
186  const Array4<Real > rt0_arr = rt0.array(mfi);
187  const Array4<Real > rt0_tmp_arr = rt0_new.array(mfi);
188 
189  const Array4<Real const> r0_arr = r0->const_array(mfi);
190  const Array4<Real > r0_new_arr = r0_new->array(mfi);
191  const Array4<Real > r0_tmp_arr = r0_temp.array(mfi);
192 
193  const Array4<Real const> p0_arr = p0->const_array(mfi);
194  const Array4<Real > p0_new_arr = p0_new->array(mfi);
195  const Array4<Real > pi0_new_arr = pi0_new->array(mfi);
196  const Array4<Real > th0_new_arr = th0_new->array(mfi);
197 
198  const Array4<Real >& z_t_arr = z_t_rk[level]->array(mfi);
199 
200  const Array4<Real const>& dJ_old_arr = detJ_cc[level]->const_array(mfi);
201  const Array4<Real const>& dJ_new_arr = detJ_cc_new[level]->const_array(mfi);
202  const Array4<Real const>& dJ_src_arr = detJ_cc_src[level]->const_array(mfi);
203 
204  Box gbx = mfi.growntilebox({1,1,1});
205  amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
206  {
207  rt0_arr(i,j,k) = getRhoThetagivenP(p0_arr(i,j,k));
208  rt0_tmp_arr(i,j,k) = getRhoThetagivenP(p0_new_arr(i,j,k));
209  r0_tmp_arr(i,j,k) = r0_new_arr(i,j,k);
210  });
211 
212  Box gbx2 = mfi.growntilebox({1,1,0});
213  amrex::ParallelFor(gbx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
214  {
215  Real zflux_r_lo = -z_t_arr(i,j,k ) * 0.5 * (r0_tmp_arr(i,j,k) + r0_tmp_arr(i,j,k-1));
216  Real zflux_r_hi = -z_t_arr(i,j,k+1) * 0.5 * (r0_tmp_arr(i,j,k) + r0_tmp_arr(i,j,k+1));
217 
218  Real zflux_rt_lo = zflux_r_lo * 0.5 * (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));
219  Real zflux_rt_hi = zflux_r_hi * 0.5 * (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));
220 
221  Real invdetJ = 1.0 / dJ_src_arr(i,j,k);
222 
223  Real src_r = - invdetJ * ( zflux_r_hi - zflux_r_lo ) * dxInv[2];
224  Real src_rt = - invdetJ * ( zflux_rt_hi - zflux_rt_lo ) * dxInv[2];
225 
226  Real rho0_new = dJ_old_arr(i,j,k) * r0_arr(i,j,k) + dt_base * dJ_src_arr(i,j,k) * src_r;
227  Real rt0_tmp_new = dJ_old_arr(i,j,k) * rt0_arr(i,j,k) + dt_base * dJ_src_arr(i,j,k) * src_rt;
228 
229  r0_new_arr(i,j,k) = rho0_new / dJ_new_arr(i,j,k);
230  rt0_tmp_new /= dJ_new_arr(i,j,k);
231 
232  p0_new_arr(i,j,k) = getPgivenRTh(rt0_tmp_new);
233  pi0_new_arr(i,j,k) = getExnergivenRTh(rt0_tmp_new, l_rdOcp);
234  th0_new_arr(i,j,k) = rt0_tmp_new / r0_new_arr(i,j,k);
235  });
236  } // MFIter
237  r0_new->FillBoundary(fine_geom.periodicity());
238  p0_new->FillBoundary(fine_geom.periodicity());
239 
240  } else { // If not moving_terrain
241 
242  make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim,
243  z_phys_nd[level], z_phys_cc[level],
244  xvel_new, yvel_new, zvel_new,
245  xmom_src, ymom_src, zmom_src,
246  base_state[level], forest_drag, terrain_blank, fine_geom, solverChoice,
247  mapfac_m[level], mapfac_u[level], mapfac_v[level],
248  dptr_u_geos, dptr_v_geos, dptr_wbar_sub,
249  d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev,
250  input_sounding_data, n_qstate);
251 
252  erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch,
253  xvel_new, yvel_new, zvel_new,
254  z_t_rk[level], cc_src, xmom_src, ymom_src, zmom_src,
255  (level > 0) ? &zmom_crse_rhs[level] : nullptr,
256  Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(),
257  Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(),
258  Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3,Q2fx3, Diss,
259  fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
260  z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0,
261  pp_inc[level],
262  mapfac_m[level], mapfac_u[level], mapfac_v[level],
263 #ifdef ERF_USE_EB
264  EBFactory(level),
265 #endif
266  fr_as_crse, fr_as_fine);
267 
268  add_thin_body_sources(xmom_src, ymom_src, zmom_src,
269  xflux_imask[level], yflux_imask[level], zflux_imask[level],
270  thin_xforce[level], thin_yforce[level], thin_zforce[level]);
271  }
272 
273 #ifdef ERF_USE_NETCDF
274  // Populate RHS for relaxation zones if using real bcs
275  if (use_real_bcs && (level == 0)) {
276  if (real_width>0) {
277  realbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt,
278  real_width, real_set_width, fine_geom,
279  S_rhs, S_old, S_data,
280  bdy_data_xlo, bdy_data_xhi,
281  bdy_data_ylo, bdy_data_yhi);
282  }
283  }
284 #endif
285 
286 #if 0
287  // HACK -- NO RELAXATION INSIDE FINE GRIDS
288  // Compute RHS for fine interior ghost
289  if (level > 0 && cf_width > 0) {
290  fine_compute_interior_ghost_rhs(new_stage_time, slow_dt,
291  cf_width, cf_set_width, fine_geom,
292  &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1],
293  domain_bcs_type, S_rhs, S_data);
294  }
295 #endif
296  }; // end slow_rhs_fun_pre
297 
298  // *************************************************************
299  // The "slow" integrator for MRI and the only integrator for SRI
300  // *************************************************************
301  auto slow_rhs_fun_post = [&](Vector<MultiFab>& S_rhs,
302  Vector<MultiFab>& S_old,
303  Vector<MultiFab>& S_new,
304  Vector<MultiFab>& S_data,
305  Vector<MultiFab>& S_scratch,
306  const Real old_step_time,
307  const Real old_stage_time,
308  const Real new_stage_time,
309  const int nrk)
310  {
311  amrex::ignore_unused(nrk);
312 
313  // Note that the "old" and "new" metric terms correspond to
314  // t^n and the RK stage (either t^*, t^** or t^{n+1} that this source
315  // will be used to advance to
316  Real slow_dt = new_stage_time - old_step_time;
317 
318  if (verbose) amrex::Print() << "Time integration of scalars at level " << level
319  << " from " << old_step_time << " to " << new_stage_time
320  << " with dt = " << slow_dt
321  << " using RHS created at " << old_stage_time << std::endl;
322 
323  int n_qstate = micro->Get_Qstate_Size();
324 
325 #if defined(ERF_USE_NETCDF)
326  bool moist_set_rhs = false;
327  if ( use_real_bcs &&
328  (level==0) &&
329  (real_set_width > 0) &&
330  (solverChoice.moisture_type != MoistureType::None) )
331  {
332  moist_set_rhs = true;
333  }
334 #endif
335 
336  // *************************************************************************
337  // Set up flux registers if using two_way coupling
338  // *************************************************************************
339  YAFluxRegister* fr_as_crse = nullptr;
340  YAFluxRegister* fr_as_fine = nullptr;
341  if (solverChoice.coupling_type == CouplingType::TwoWay)
342  {
343  if (level < finest_level) {
344  fr_as_crse = getAdvFluxReg(level+1);
345  }
346  if (level > 0) {
347  fr_as_fine = getAdvFluxReg(level);
348  }
349  }
350 
351  // Moving terrain
352  if ( solverChoice.terrain_type == TerrainType::Moving ) {
353  erf_slow_rhs_post(level, finest_level, nrk, slow_dt, n_qstate,
354  S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
355  xvel_new, yvel_new, zvel_new, cc_src, SmnSmn, eddyDiffs,
356  Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
357  fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
358  z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc_new[level],
359  mapfac_m[level], mapfac_u[level], mapfac_v[level],
360 #ifdef ERF_USE_EB
361  EBFactory(level),
362 #endif
363 #if defined(ERF_USE_NETCDF)
364  moist_set_rhs, bdy_time_interval, start_bdy_time, new_stage_time,
365  real_width, real_set_width,
366  bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi,
367 #endif
368  fr_as_crse, fr_as_fine);
369  } else {
370  erf_slow_rhs_post(level, finest_level, nrk, slow_dt, n_qstate,
371  S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
372  xvel_new, yvel_new, zvel_new, cc_src, SmnSmn, eddyDiffs,
373  Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
374  fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
375  z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], detJ_cc[level],
376  mapfac_m[level], mapfac_u[level], mapfac_v[level],
377 #ifdef ERF_USE_EB
378  EBFactory(level),
379 #endif
380 #if defined(ERF_USE_NETCDF)
381  moist_set_rhs, bdy_time_interval, start_bdy_time, new_stage_time,
382  real_width, real_set_width,
383  bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi,
384 #endif
385  fr_as_crse, fr_as_fine);
386  }
387 
388  // Apply boundary conditions on all the state variables that have been updated
389  // in both the fast and slow integrators
390  apply_bcs(S_new, new_stage_time, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow(),
391  fast_only=false, vel_and_mom_synced=false);
392 
393  }; // end slow_rhs_fun_post
394 
395  auto slow_rhs_fun_inc = [&](Vector<MultiFab>& S_rhs,
396  Vector<MultiFab>& S_old,
397  Vector<MultiFab>& S_data,
398  Vector<MultiFab>& S_scratch,
399  const Real old_step_time,
400  const Real old_stage_time,
401  const Real new_stage_time,
402  const int nrk)
403  {
404  BL_PROFILE("slow_rhs_fun_inc");
405  if (verbose) Print() << "Making slow rhs at time " << old_stage_time << " for fast variables advancing from " <<
406  old_step_time << " to " << new_stage_time << std::endl;
407  //
408  // Define primitive variables for all later RK stages
409  // (We have already done this for the first RK step)
410  //
411  if (nrk > 0) {
412  int ng_cons = S_data[IntVars::cons].nGrow();
413  cons_to_prim(S_data[IntVars::cons], ng_cons);
414  }
415 
416  Real slow_dt = new_stage_time - old_step_time;
417 
418  // *************************************************************************
419  // Set up flux registers if using two_way coupling
420  // *************************************************************************
421  YAFluxRegister* fr_as_crse = nullptr;
422  YAFluxRegister* fr_as_fine = nullptr;
423  if (solverChoice.coupling_type == CouplingType::TwoWay) {
424  if (level < finest_level) {
425  fr_as_crse = getAdvFluxReg(level+1);
426  fr_as_crse->reset();
427  }
428  if (level > 0) {
429  fr_as_fine = getAdvFluxReg(level);
430  }
431  }
432 
433  Real* dptr_u_geos = solverChoice.have_geo_wind_profile ? d_u_geos[level].data(): nullptr;
434  Real* dptr_v_geos = solverChoice.have_geo_wind_profile ? d_v_geos[level].data(): nullptr;
435 
436  // Canopy data for mom sources
437  MultiFab* forest_drag = nullptr;
438  if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); }
439  // Immersed Forcing
440  MultiFab* terrain_blank = nullptr;
441  if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); }
442 
443  make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level],
444 #if defined(ERF_USE_RRTMGP)
445  qheating_rates[level],
446 #endif
447  fine_geom, solverChoice,
448  mapfac_u[level], mapfac_v[level], mapfac_m[level],
449  dptr_rhotheta_src, dptr_rhoqt_src,
450  dptr_wbar_sub, d_rayleigh_ptrs_at_lev,
451  input_sounding_data, turbPert);
452 
453  int n_qstate = micro->Get_Qstate_Size();
454  make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim,
455  z_phys_nd[level], z_phys_cc[level],
456  xvel_new, yvel_new, zvel_new,
457  xmom_src, ymom_src, zmom_src,
458  base_state[level], forest_drag, terrain_blank, fine_geom, solverChoice,
459  mapfac_m[level], mapfac_u[level], mapfac_v[level],
460  dptr_u_geos, dptr_v_geos, dptr_wbar_sub,
461  d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev,
462  input_sounding_data, n_qstate);
463 
464  erf_slow_rhs_pre(level, finest_level, nrk, slow_dt,
465  S_rhs, S_old, S_data, S_prim, S_scratch,
466  xvel_new, yvel_new, zvel_new,
467  z_t_rk[level], cc_src, xmom_src, ymom_src, zmom_src,
468  (level > 0) ? &zmom_crse_rhs[level] : nullptr,
469  Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(),
470  Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(),
471  Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
472  fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
473  z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0,
474  pp_inc[level],
475  mapfac_m[level], mapfac_u[level], mapfac_v[level],
476 #ifdef ERF_USE_EB
477  EBFactory(level),
478 #endif
479  fr_as_crse, fr_as_fine);
480 
481  add_thin_body_sources(xmom_src, ymom_src, zmom_src,
482  xflux_imask[level], yflux_imask[level], zflux_imask[level],
483  thin_xforce[level], thin_yforce[level], thin_zforce[level]);
484 
485 #ifdef ERF_USE_NETCDF
486  // Populate RHS for relaxation zones if using real bcs
487  if (use_real_bcs && (level == 0)) {
488  if (real_width>0) {
489  realbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt,
490  real_width, real_set_width, fine_geom,
491  S_rhs, S_old, S_data,
492  bdy_data_xlo, bdy_data_xhi,
493  bdy_data_ylo, bdy_data_yhi);
494  }
495  }
496 #endif
497  }; // end slow_rhs_fun_inc
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:29
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:84
getExnergivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:159
getRhoThetagivenP
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:175
AdvType::Centered_2nd
@ Centered_2nd
realbdy_compute_interior_ghost_rhs
void realbdy_compute_interior_ghost_rhs(const Real &bdy_time_interval, const Real &start_bdy_time, const Real &time, const Real &delta_t, int width, int set_width, const Geometry &geom, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old_data, 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)
Definition: ERF_InteriorGhostCells.cpp:107
fine_compute_interior_ghost_rhs
void fine_compute_interior_ghost_rhs(const Real &time, const Real &delta_t, const int &width, const int &set_width, const Geometry &geom, ERFFillPatcher *FPr_c, ERFFillPatcher *FPr_u, ERFFillPatcher *FPr_v, ERFFillPatcher *FPr_w, Vector< BCRec > &domain_bcs_type, Vector< MultiFab > &S_rhs_f, Vector< MultiFab > &S_data_f)
Definition: ERF_InteriorGhostCells.cpp:520
make_mom_sources
void make_mom_sources(int level, int, Real, Real time, Vector< MultiFab > &S_data, const MultiFab &S_prim, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &z_phys_cc, 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, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< MultiFab > &, std::unique_ptr< MultiFab > &, 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, int n_qstate)
Definition: ERF_MakeMomSources.cpp:40
make_sources
void make_sources(int level, int, Real dt, Real time, Vector< MultiFab > &S_data, const MultiFab &S_prim, MultiFab &source, std::unique_ptr< MultiFab > &z_phys_cc, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< MultiFab > &, std::unique_ptr< MultiFab > &, std::unique_ptr< MultiFab > &mapfac_m, 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)
Definition: ERF_MakeSources.cpp:32
erf_slow_rhs_post
void erf_slow_rhs_post(int level, int finest_level, int nrk, Real dt, int n_qstate, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_new, Vector< MultiFab > &S_data, const MultiFab &S_prim, Vector< MultiFab > &S_scratch, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &, const MultiFab &source, const MultiFab *SmnSmn, const MultiFab *eddyDiffs, MultiFab *Hfx1, MultiFab *Hfx2, MultiFab *Hfx3, MultiFab *Q1fx1, MultiFab *Q1fx2, MultiFab *Q1fx3, MultiFab *Q2fx3, MultiFab *Diss, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< ABLMost > &most, const Gpu::DeviceVector< BCRec > &domain_bcs_type_d, const Vector< BCRec > &domain_bcs_type_h, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &ax, std::unique_ptr< MultiFab > &ay, std::unique_ptr< MultiFab > &az, std::unique_ptr< MultiFab > &detJ, std::unique_ptr< MultiFab > &detJ_new, std::unique_ptr< MultiFab > &mapfac_m, std::unique_ptr< MultiFab > &mapfac_u, std::unique_ptr< MultiFab > &mapfac_v, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine)
Definition: ERF_SlowRhsPost.cpp:45
erf_slow_rhs_pre
void erf_slow_rhs_pre(int level, int finest_level, int nrk, Real dt, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_data, const MultiFab &S_prim, Vector< MultiFab > &S_scratch, 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 *zmom_crse_rhs, MultiFab *Tau11, MultiFab *Tau22, MultiFab *Tau33, MultiFab *Tau12, MultiFab *Tau13, MultiFab *Tau21, MultiFab *Tau23, MultiFab *Tau31, MultiFab *Tau32, MultiFab *SmnSmn, MultiFab *eddyDiffs, MultiFab *Hfx1, MultiFab *Hfx2, MultiFab *Hfx3, MultiFab *Q1fx1, MultiFab *Q1fx2, MultiFab *Q1fx3, MultiFab *Q2fx3, MultiFab *Diss, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< ABLMost > &most, const Gpu::DeviceVector< BCRec > &domain_bcs_type_d, const Vector< BCRec > &domain_bcs_type_h, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &ax, std::unique_ptr< MultiFab > &ay, std::unique_ptr< MultiFab > &az, std::unique_ptr< MultiFab > &detJ, const MultiFab *p0, const MultiFab &pp_inc, std::unique_ptr< MultiFab > &mapfac_m, std::unique_ptr< MultiFab > &mapfac_u, std::unique_ptr< MultiFab > &mapfac_v, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine)
Definition: ERF_SlowRhsPre.cpp:67
ERF_SrcHeaders.H
slow_rhs_fun_inc
auto slow_rhs_fun_inc
Definition: ERF_TI_slow_rhs_fun.H:395
slow_rhs_fun_pre
auto slow_rhs_fun_pre
Definition: ERF_TI_slow_rhs_fun.H:6
slow_rhs_fun_post
auto slow_rhs_fun_post
Definition: ERF_TI_slow_rhs_fun.H:301
apply_bcs
auto apply_bcs
Definition: ERF_TI_utils.H:50
cons_to_prim
auto cons_to_prim
Definition: ERF_TI_utils.H:4
init_terrain_grid
void init_terrain_grid(int lev, const Geometry &geom, MultiFab &z_phys_nd, Vector< Real > const &z_levels_h, GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type)
Definition: ERF_TerrainMetrics.cpp:118
make_areas
void make_areas(const Geometry &geom, MultiFab &z_phys_nd, MultiFab &ax, MultiFab &ay, MultiFab &az)
Definition: ERF_TerrainMetrics.cpp:651
make_J
void make_J(const Geometry &geom, MultiFab &z_phys_nd, MultiFab &detJ_cc)
Definition: ERF_TerrainMetrics.cpp:613
BaseState::pi0_comp
@ pi0_comp
Definition: ERF_IndexDefines.H:65
BaseState::p0_comp
@ p0_comp
Definition: ERF_IndexDefines.H:64
BaseState::th0_comp
@ th0_comp
Definition: ERF_IndexDefines.H:66
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:63
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:139
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:140