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