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