ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MakeSources.cpp File Reference
#include <AMReX_MultiFab.H>
#include <AMReX_ArrayLim.H>
#include <AMReX_BCRec.H>
#include <AMReX_TableData.H>
#include <AMReX_GpuContainers.H>
#include <ERF_ImmersedForcing.H>
#include <ERF_NumericalDiffusion.H>
#include <ERF_PlaneAverage.H>
#include <ERF_SrcHeaders.H>
#include <ERF_TI_slow_headers.H>
#include <ERF_MOSTStress.H>
Include dependency graph for ERF_MakeSources.cpp:

Functions

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)
 

Function Documentation

◆ 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 
)

Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum.

Parameters
[in]levellevel of resolution
[in]nrkwhich RK stage
[in]dtslow time step
[in]S_datacurrent solution
[in]S_primprimitive variables (i.e. conserved variables divided by density)
[in]sourcesource terms for conserved variables
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]mapfacmap factors
[in]dptr_rhotheta_srccustom temperature source term
[in]dptr_rhoqt_srccustom moisture source term
[in]dptr_wbar_subsubsidence source term
[in]d_rayleigh_ptrs_at_levVector of {strength of Rayleigh damping, reference value of theta} used to define Rayleigh damping
[in]d_sinesq_at_levsin( (pi/2) (z-z_t)/(damping depth)) at cell centers
61 {
62  BL_PROFILE_REGION("erf_make_sources()");
63 
64  Real time = static_cast<Real>(time_d);
65 
66  // *****************************************************************************
67  // Initialize source to zero since we re-compute it every RK stage
68  // *****************************************************************************
69  if (is_slow_step) {
70  source.setVal(0.);
71  } else {
72  source.setVal(0.0,Rho_comp,2);
73  }
74 
75  const bool l_use_ndiff = solverChoice.use_num_diff;
76 
77  TurbChoice tc = solverChoice.turbChoice[level];
78  const bool l_use_KE = tc.use_tke;
79  const bool l_diff_KE = tc.diffuse_tke_3D;
80 
81  const Box& domain = geom.Domain();
82 
83  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
84  const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
85 
86  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
87  MultiFab th_hse (base_state, make_alias, BaseState::th0_comp , 1);
88  MultiFab qv_hse (base_state, make_alias, BaseState::qv0_comp , 1);
89 
90  Real* thetabar = d_rayleigh_ptrs_at_lev[Rayleigh::thetabar];
91 
92  // flags to apply certain source terms in substep call only
93  bool use_Rayleigh_fast = ( (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit) ||
94  (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastImplicit) );
95  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
96 
97  // flag for a moisture model
98  bool has_moisture = (solverChoice.moisture_type != MoistureType::None);
99 
100  // *****************************************************************************
101  // Planar averages for subsidence terms and immersed forcing
102  // *****************************************************************************
103  Table1D<Real> dptr_r_plane, dptr_t_plane, dptr_qv_plane, dptr_qc_plane;
104  Table1D<Real> dptr_r_plane_if, dptr_t_plane_if; // Separate tables for immersed forcing
105  TableData<Real, 1> r_plane_tab, t_plane_tab, qv_plane_tab, qc_plane_tab;
106 
107  bool use_immersed_forcing = (solverChoice.terrain_type == TerrainType::ImmersedForcing ||
108  solverChoice.buildings_type == BuildingsType::ImmersedForcing);
109 
110  // Use passed-in planar averages for immersed forcing (computed once in ERF_AdvanceDycore.cpp)
111  AMREX_ALWAYS_ASSERT(!use_immersed_forcing || (r_plane_avg && t_plane_avg));
112  if (use_immersed_forcing && r_plane_avg && t_plane_avg) {
113  // Use pre-computed values from persistent storage (already in Table1D format)
114  dptr_r_plane_if = r_plane_avg;
115  dptr_t_plane_if = t_plane_avg;
116  }
117 
118  // Compute planar averages if needed for subsidence
119  bool compute_averages = (is_slow_step && dptr_wbar_sub);
120 
121  if (compute_averages)
122  {
123  // The plane averaging operates at fixed k-index (not fixed height), which is
124  // acceptable for immersed forcing but not correct for subsidence with variable dz
125  AMREX_ALWAYS_ASSERT(solverChoice.mesh_type != MeshType::VariableDz);
126 
127  //
128  // The call to "compute_averages" currently does all the components in one call
129  // We can then extract each component separately with the "line_average" call
130  //
131  // We need just one ghost cell in the vertical
132  //
133  IntVect ng_c(S_data[IntVars::cons].nGrowVect()); ng_c[2] = 1;
134  //
135  // With no moisture we only (rho) and (rho theta); with moisture we also do qv and qc
136  // We use the alias here to control ncomp inside the PlaneAverage
137  //
138  int ncomp = (!has_moisture) ? 2 : RhoQ2_comp+1;
139  MultiFab cons(S_data[IntVars::cons], make_alias, 0, ncomp);
140 
141  PlaneAverage cons_ave(&cons, geom, solverChoice.ave_plane, ng_c);
142  cons_ave.compute_averages(ZDir(), cons_ave.field());
143 
144  int ncell = cons_ave.ncell_line();
145 
146  Gpu::HostVector< Real> r_plane_h(ncell);
147  Gpu::DeviceVector< Real> r_plane_d(ncell);
148 
149  Gpu::HostVector< Real> t_plane_h(ncell);
150  Gpu::DeviceVector< Real> t_plane_d(ncell);
151 
152  cons_ave.line_average(Rho_comp , r_plane_h);
153  cons_ave.line_average(RhoTheta_comp, t_plane_h);
154 
155  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
156  Gpu::copyAsync(Gpu::hostToDevice, t_plane_h.begin(), t_plane_h.end(), t_plane_d.begin());
157 
158  Real* dptr_r = r_plane_d.data();
159  Real* dptr_t = t_plane_d.data();
160 
161  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
162  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
163  t_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
164 
165  int offset = ng_c[2];
166 
167  dptr_r_plane = r_plane_tab.table();
168  dptr_t_plane = t_plane_tab.table();
169  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
170  {
171  dptr_r_plane(k-offset) = dptr_r[k];
172  dptr_t_plane(k-offset) = dptr_t[k];
173  });
174 
175  if (has_moisture)
176  {
177  Gpu::HostVector< Real> qv_plane_h(ncell), qc_plane_h(ncell);
178  Gpu::DeviceVector<Real> qv_plane_d(ncell), qc_plane_d(ncell);
179 
180  // Water vapor
181  cons_ave.line_average(RhoQ1_comp, qv_plane_h);
182  Gpu::copyAsync(Gpu::hostToDevice, qv_plane_h.begin(), qv_plane_h.end(), qv_plane_d.begin());
183 
184  // Cloud water
185  cons_ave.line_average(RhoQ2_comp, qc_plane_h);
186  Gpu::copyAsync(Gpu::hostToDevice, qc_plane_h.begin(), qc_plane_h.end(), qc_plane_d.begin());
187 
188  Real* dptr_qv = qv_plane_d.data();
189  Real* dptr_qc = qc_plane_d.data();
190 
191  qv_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
192  qc_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
193 
194  dptr_qv_plane = qv_plane_tab.table();
195  dptr_qc_plane = qc_plane_tab.table();
196  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
197  {
198  dptr_qv_plane(k-offset) = dptr_qv[k];
199  dptr_qc_plane(k-offset) = dptr_qc[k];
200  });
201  }
202  }
203 
204  // *****************************************************************************
205  // Vertical extent used by the four stream radiation approximation
206  // *****************************************************************************
207  // NOTE: The fluxes live on w-faces
208  //
209  // NOTE: There is deliberately no scratch storage for the Q integral or the
210  // radiative flux here. Those used to be two nk-length device vectors
211  // shared by every (i,j) thread of the 2-D ParallelFor below, so all the
212  // columns raced on the same buffers (issue 3714). The kernel now keeps
213  // the running integral and the lower-face flux in thread-private scalars.
214  int klo = domain.smallEnd(2);
215  int khi = domain.bigEnd(2);
216 
217  // *****************************************************************************
218  // Define source term for cell-centered conserved variables, from
219  // 1. user-defined source terms for (rho theta) and (rho q_t)
220  // 2. radiation for (rho theta)
221  // 3. Rayleigh damping for (rho theta)
222  // 4. custom forcing for (rho theta) and (rho Q1)
223  // 5. custom subsidence for (rho theta) and (rho Q1)
224  // 6. numerical diffusion for (rho theta)
225  // 7. sponging
226  // 8. turbulent perturbation
227  // 9. nudging towards input sounding values (only for theta)
228  // 10a. Immersed forcing for terrain
229  // 10b. Immersed forcing for buildings
230  // Real(11.) Four stream radiation source for (rho theta)
231  // *****************************************************************************
232 
233  // ***********************************************************************************************
234  // Add remaining source terms
235  // ***********************************************************************************************
236 #ifdef _OPENMP
237 #pragma omp parallel if (Gpu::notInLaunchRegion())
238 #endif
239  {
240  for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
241  {
242  Box bx = mfi.tilebox();
243 
244  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
245  const Array4<const Real>& cell_prim = S_prim.array(mfi);
246  const Array4<Real> & cell_src = source.array(mfi);
247 
248  const Array4<const Real>& r0 = r_hse.const_array(mfi);
249  const Array4<const Real>& th0 = th_hse.const_array(mfi);
250  const Array4<const Real>& qv0 = qv_hse.const_array(mfi);
251 
252  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
253 
254  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
255  Array4<const Real>{};
256 
257 
258  // *************************************************************************************
259  // 2. Add radiation source terms to (rho theta)
260  // *************************************************************************************
261  if (solverChoice.rad_type != RadiationType::None && is_slow_step) {
262  auto const& qheating_arr = qheating_rates->const_array(mfi);
263  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
264  {
265  // Short-wavelength and long-wavelength radiation source terms
266  cell_src(i,j,k,RhoTheta_comp) += cell_data(i,j,k,Rho_comp) * ( qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1) );
267  });
268  }
269 
270 
271  // *************************************************************************************
272  // 3. Add Rayleigh damping for (rho theta)
273  // *************************************************************************************
274  Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef;
275 
276  if (solverChoice.dampingChoice.rayleigh_damp_T) {
277  if ((is_slow_step && !use_Rayleigh_fast) || (!is_slow_step && use_Rayleigh_fast)) {
278  int n = RhoTheta_comp;
279  int nr = Rho_comp;
280  int np = PrimTheta_comp;
281 
282  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
283  {
284  Real theta = cell_prim(i,j,k,np);
285  Real sinesq = d_sinesq_at_lev[k];
286  cell_src(i, j, k, n) -= dampcoef*sinesq * (theta - thetabar[k]) * cell_data(i,j,k,nr);
287  });
288  }
289  }
290 
291  // *************************************************************************************
292  // 4. Add custom forcing for (rho theta)
293  // *************************************************************************************
294  if (solverChoice.custom_rhotheta_forcing && is_slow_step) {
295  const int n = RhoTheta_comp;
296  auto const& rhotheta_src_arr = rhotheta_src->const_array(mfi);
297  if (solverChoice.spatial_rhotheta_forcing)
298  {
299  if (solverChoice.custom_forcing_prim_vars) {
300  const int nr = Rho_comp;
301  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
302  {
303  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhotheta_src_arr(i, j, k);
304  });
305  } else {
306  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
307  {
308  cell_src(i, j, k, n) += rhotheta_src_arr(i, j, k);
309  });
310  }
311  } else {
312  if (solverChoice.custom_forcing_prim_vars) {
313  const int nr = Rho_comp;
314  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
315  {
316  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhotheta_src_arr(0, 0, k);
317  });
318  } else {
319  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
320  {
321  cell_src(i, j, k, n) += rhotheta_src_arr(0, 0, k);
322  });
323  }
324  }
325  }
326 
327  // *************************************************************************************
328  // 4. Add custom forcing for RhoQ1
329  // *************************************************************************************
330  if (solverChoice.custom_moisture_forcing && is_slow_step) {
331  const int n = RhoQ1_comp;
332  auto const& rhoqt_src_arr = rhoqt_src->const_array(mfi);
333  if (solverChoice.spatial_moisture_forcing)
334  {
335  if (solverChoice.custom_forcing_prim_vars) {
336  const int nr = Rho_comp;
337  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
338  {
339  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhoqt_src_arr(i, j, k);
340  });
341  } else {
342  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
343  {
344  cell_src(i, j, k, n) += rhoqt_src_arr(i, j, k);
345  });
346  }
347  } else {
348  if (solverChoice.custom_forcing_prim_vars) {
349  const int nr = Rho_comp;
350  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
351  {
352  cell_src(i, j, k, n) += cell_data(i,j,k,nr) * rhoqt_src_arr(0, 0, k);
353  });
354  } else {
355  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
356  {
357  cell_src(i, j, k, n) += rhoqt_src_arr(0, 0, k);
358  });
359  }
360  }
361  }
362 
363  // *************************************************************************************
364  // 5. Add custom subsidence for (rho theta)
365  // *************************************************************************************
366  if (solverChoice.custom_w_subsidence && is_slow_step && solverChoice.do_theta_advection) {
367  const int n = RhoTheta_comp;
368  if (solverChoice.custom_forcing_prim_vars) {
369  const int nr = Rho_comp;
370  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
371  {
372  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
373  Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
374  Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
375  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
376  cell_src(i, j, k, n) -= cell_data(i,j,k,nr) * wbar_cc * (T_hi - T_lo) * dzInv;
377  });
378  } else {
379  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
380  {
381  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
382  Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
383  Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
384  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
385  cell_src(i, j, k, n) -= wbar_cc * (T_hi - T_lo) * dzInv;
386  });
387  }
388  }
389 
390  // *************************************************************************************
391  // 5. Add custom subsidence for RhoQ1 and RhoQ2
392  // *************************************************************************************
393  if (solverChoice.custom_w_subsidence && (solverChoice.moisture_type != MoistureType::None) && is_slow_step) {
394  const int nv = RhoQ1_comp;
395  if (solverChoice.custom_forcing_prim_vars) {
396  const int nr = Rho_comp;
397  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
398  {
399  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
400  Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
401  Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
402  Real Qc_hi = dptr_qc_plane(k+1) / dptr_r_plane(k+1);
403  Real Qc_lo = dptr_qc_plane(k-1) / dptr_r_plane(k-1);
404  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
405  cell_src(i, j, k, nv ) -= cell_data(i,j,k,nr) * wbar_cc * (Qv_hi - Qv_lo) * dzInv;
406  cell_src(i, j, k, nv+1) -= cell_data(i,j,k,nr) * wbar_cc * (Qc_hi - Qc_lo) * dzInv;
407  });
408  } else {
409  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
410  {
411  Real dzInv = (z_cc_arr) ? one/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : myhalf*dxInv[2];
412  Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
413  Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
414  Real Qc_hi = dptr_qc_plane(k+1) / dptr_r_plane(k+1);
415  Real Qc_lo = dptr_qc_plane(k-1) / dptr_r_plane(k-1);
416  Real wbar_cc = myhalf * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
417  cell_src(i, j, k, nv ) -= wbar_cc * (Qv_hi - Qv_lo) * dzInv;
418  cell_src(i, j, k, nv+1) -= wbar_cc * (Qc_hi - Qc_lo) * dzInv;
419  });
420  }
421  }
422 
423  // *************************************************************************************
424  // 6. Add numerical diffusion for rho and (rho theta)
425  // *************************************************************************************
426  if (l_use_ndiff && is_slow_step)
427  {
428  const Array4<const Real>& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi);
429  const Array4<const Real>& mf_my = mapfac[MapFacType::m_y]->const_array(mfi);
430 
431  // Rho is a special case
432  NumericalDiffusion_Scal(bx, 0, 1, dt, solverChoice.num_diff_coeff,
433  cell_data, cell_data, cell_src, mf_mx, mf_my);
434 
435  // Other scalars proceed as normal
436  NumericalDiffusion_Scal(bx, 1, 1, dt, solverChoice.num_diff_coeff,
437  cell_prim, cell_data, cell_src, mf_mx, mf_my);
438 
439 
440  if (l_use_KE && l_diff_KE) {
441  NumericalDiffusion_Scal(bx, RhoKE_comp, 1, dt, solverChoice.num_diff_coeff,
442  cell_prim, cell_data, cell_src, mf_mx, mf_my);
443  }
444 
445  NumericalDiffusion_Scal(bx, RhoScalar_comp, NSCALARS, dt, solverChoice.num_diff_coeff,
446  cell_prim, cell_data, cell_src, mf_mx, mf_my);
447  }
448 
449  // *************************************************************************************
450  // 7. Add sponging
451  // *************************************************************************************
452  if ( (solverChoice.spongeChoice.sponge_type != SpongeType::None) && is_slow_step) {
453  const int n_qstate = S_data[IntVars::cons].nComp() - (NDRY + NSCALARS);
454  ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data, r0, th0, qv0, z_cc_arr, n_qstate);
455  }
456 
457  // *************************************************************************************
458  // 8. Add perturbation
459  // *************************************************************************************
460  if (solverChoice.use_source_perturbation(level) && is_slow_step) {
461  auto m_ixtype = S_data[IntVars::cons].boxArray().ixType(); // Conserved term
462  const amrex::Array4<const amrex::Real>& pert_cell = turbPert.pb_cell[level].const_array(mfi);
463  turbPert.apply_tpi(level, bx, RhoTheta_comp, m_ixtype, cell_src, pert_cell); // Applied as source term
464  }
465 
466  // *************************************************************************************
467  // 10a. Add immersed source terms for terrain
468  // *************************************************************************************
469  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
470  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
471 
472  const Array4<const Real>& u = xvel.array(mfi);
473  const Array4<const Real>& v = yvel.array(mfi);
474 
475  AMREX_ALWAYS_ASSERT(dptr_r_plane_if && dptr_t_plane_if); // Tables must be filled before immersed forcing
476  ImmersedForcingTerrain_Scalar(bx, u, v, cell_data, t_blank_arr, z_cc_arr,
477  cell_src, geom, solverChoice, dptr_r_plane_if, dptr_t_plane_if, time);
478  }
479 
480  // *************************************************************************************
481  // 10b. Add immersed source terms for buildings
482  // *************************************************************************************
483  // geometric properties
484  const Real* dx_arr = geom.CellSize();
485  const Real delta_xy = std::sqrt(dx_arr[0] * dx_arr[1]);
486 
487  if ((solverChoice.buildings_type == BuildingsType::ImmersedForcing) &&
488  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)) &&
489  (delta_xy <= 50.0)) { // only apply immersed forcing when grid spacing is less than 50m
490 
491  const Array4<const Real>& u = xvel.array(mfi);
492  const Array4<const Real>& v = yvel.array(mfi);
493  const Array4<const Real>& w = zvel.array(mfi);
494 
495  AMREX_ALWAYS_ASSERT(dptr_r_plane_if && dptr_t_plane_if); // Tables must be filled before immersed forcing
496  ImmersedForcingBuildings_Scalar(bx, u, v, w, cell_data, t_blank_arr, z_cc_arr,
497  cell_src, geom, solverChoice, dptr_r_plane_if, dptr_t_plane_if, time);
498  }
499 
500  // *************************************************************************************
501  // Real(11.) Add 4 stream radiation src to RhoTheta
502  // *************************************************************************************
503  if (solverChoice.four_stream_radiation && has_moisture && is_slow_step)
504  {
505  AMREX_ALWAYS_ASSERT((bx.smallEnd(2) == klo) && (bx.bigEnd(2) == khi));
506  Real D = Real(3.75e-6); // [s^-1]
507  Real F0 = 70; // [W/m^2]
508  Real F1 = 22; // [W/m^2]
509  Real krad = 85; // [m^2 kg^-1]
510  Real qt_i = Real(0.008);
511 
512  Box xybx = makeSlab(bx,2,klo);
513  ParallelFor(xybx, [=]
514  AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
515  {
516  // Pass 1: total Q integral through the column (also find "i" values).
517  // NOTE: the running integral is kept in a thread-private scalar rather
518  // than a shared nk-length buffer, and is simply re-formed in the
519  // second pass below. The accumulation order is unchanged, so the
520  // values are identical to the old per-level array.
521  Real q_int_inf = zero;
522  Real zi = myhalf * (z_cc_arr(i,j,khi) + z_cc_arr(i,j,khi-1));
523  Real rhoi = myhalf * (cell_data(i,j,khi,Rho_comp) + cell_data(i,j,khi-1,Rho_comp));
524  for (int k(klo+1); k<=khi+1; ++k) {
525  // Average to w-faces when looping w-faces
526  Real dz = (z_cc_arr) ? myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-2)) : dx[2];
527  q_int_inf += krad * cell_data(i,j,k-1,Rho_comp) * cell_data(i,j,k-1,RhoQ2_comp) * dz;
528  Real qt_hi = cell_data(i,j,k ,RhoQ1_comp) + cell_data(i,j,k ,RhoQ2_comp);
529  Real qt_lo = cell_data(i,j,k-1,RhoQ1_comp) + cell_data(i,j,k-1,RhoQ2_comp);
530  if ( (qt_lo > qt_i) && (qt_hi < qt_i) ) {
531  zi = myhalf * (z_cc_arr(i,j,k) + z_cc_arr(i,j,k-1));
532  rhoi = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp));
533  }
534  }
535 
536  // Pass 2: re-form the integral level by level, decompose it into the
537  // w-face flux, and difference the flux to get the heating. The flux on
538  // the lower face is carried across iterations so each face is evaluated
539  // once, exactly as when both faces were read from the shared array.
540  Real q_int = zero; // face klo
541  Real z_lo = myhalf * (z_cc_arr(i,j,klo) + z_cc_arr(i,j,klo-1));
542  Real flux_lo = F1*std::exp(-q_int) + F0*std::exp(-(q_int_inf - q_int));
543  if (z_lo > zi) {
544  flux_lo += rhoi * Cp_d * D * ( std::pow(z_lo-zi,Real(4.)/three)/Real(4.) + zi*std::pow(z_lo-zi,one/three) ) ;
545  }
546 
547  for (int k(klo); k<=khi; ++k) {
548  // Advance the integral from face k to face k+1
549  Real dz = (z_cc_arr) ? myhalf * (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : dx[2];
550  q_int += krad * cell_data(i,j,k,Rho_comp) * cell_data(i,j,k,RhoQ2_comp) * dz;
551 
552  Real z_hi = myhalf * (z_cc_arr(i,j,k+1) + z_cc_arr(i,j,k));
553  Real flux_hi = F1*std::exp(-q_int) + F0*std::exp(-(q_int_inf - q_int));
554  if (z_hi > zi) {
555  flux_hi += rhoi * Cp_d * D * ( std::pow(z_hi-zi,Real(4.)/three)/Real(4.) + zi*std::pow(z_hi-zi,one/three) ) ;
556  }
557 
558  // Average to w-faces when looping CC
559  Real dzInv = (z_cc_arr) ? one/ (myhalf * (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1))) : dxInv[2];
560  // NOTE: Fnet = Up - Dn (all fluxes are up here)
561  // dT/dt = dF/dz * (1/(-rho*Cp))
562  Real dTdt = (flux_hi - flux_lo) * dzInv / (-cell_data(i,j,k,Rho_comp)*Cp_d);
563  Real qv = cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp);
564  Real iexner = one/getExnergivenRTh(cell_data(i,j,k,RhoTheta_comp), RdoCp, qv);
565  // Convert dT/dt to dTheta/dt and multiply rho
566  cell_src(i,j,k,RhoTheta_comp) += cell_data(i,j,k,Rho_comp) * dTdt * iexner;
567 
568  flux_lo = flux_hi;
569  }
570  });
571  }
572  } // mfi
573  } // OMP
574 }
ApplySpongeZoneBCsForCC
void ApplySpongeZoneBCsForCC(const SpongeChoice &spongeChoice, const Geometry geom, const Box &bx, const Array4< Real > &cell_rhs, const Array4< const Real > &cell_data, const Array4< const Real > &r0, const Array4< const Real > &th0, const Array4< const Real > &qv0, const Array4< const Real > &z_phys_cc, int n_qstate)
Apply sponge zone damping to cell-centered state variables.
Definition: ERF_ApplySpongeZoneBCs.cpp:21
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
Cp_d
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:49
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
RdoCp
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:54
RayleighDampingType::FastExplicit
@ FastExplicit
RayleighDampingType::FastImplicit
@ FastImplicit
thetabar
@ thetabar
Definition: ERF_DataStruct.H:153
m_y
@ m_y
Definition: ERF_DataStruct.H:29
m_x
@ m_x
Definition: ERF_DataStruct.H:28
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:55
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
ImmersedForcingTerrain_Scalar
void ImmersedForcingTerrain_Scalar(const Box &bx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &cell_src, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real time)
Definition: ERF_ImmersedForcing.cpp:840
ImmersedForcingBuildings_Scalar
void ImmersedForcingBuildings_Scalar(const Box &bx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &cell_src, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real time)
Definition: ERF_ImmersedForcing.cpp:972
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:43
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
NDRY
#define NDRY
Definition: ERF_IndexDefines.H:13
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:46
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:16
PrimTheta_comp
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:58
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:41
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
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);})
NumericalDiffusion_Scal
void NumericalDiffusion_Scal(const Box &bx, const int start_comp, const int num_comp, const double dt, const Real num_diff_coeff, const Array4< const Real > &prim_data, const Array4< const Real > &cell_data, const Array4< Real > &rhs, const Array4< const Real > &mfx_arr, const Array4< const Real > &mfy_arr)
Definition: ERF_NumericalDiffusion.cpp:18
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:31
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
PlaneAverage
Definition: ERF_PlaneAverage.H:14
BaseState::qv0_comp
@ qv0_comp
Definition: ERF_IndexDefines.H:80
BaseState::th0_comp
@ th0_comp
Definition: ERF_IndexDefines.H:79
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:76
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:232
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:20
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:30
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:46
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:217
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:270
WDM6SedNodeScratch::zi
@ zi
Definition: ERF_AdvanceWDM6.cpp:274
module_mp_morr_two_moment::rhoi
real(c_double), private rhoi
Definition: ERF_module_mp_morr_two_moment.F90:188
mp_wdm6::r0
real(kind=kind_phys), parameter, private r0
Definition: ERF_module_mp_wdm6.F90:75
DampingChoice::rayleigh_damping_type
RayleighDampingType rayleigh_damping_type
Selected Rayleigh damping time-integration treatment.
Definition: ERF_DampingStruct.H:111
DampingChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Rayleigh damping inverse time scale [1/s].
Definition: ERF_DampingStruct.H:98
DampingChoice::rayleigh_damp_T
bool rayleigh_damp_T
Whether Rayleigh damping is applied to potential temperature.
Definition: ERF_DampingStruct.H:97
SolverChoice::turbChoice
amrex::Vector< TurbChoice > turbChoice
Turbulence options for each AMR level.
Definition: ERF_DataStruct.H:1864
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:2124
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2141
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Numerical diffusion coefficient after input scaling.
Definition: ERF_DataStruct.H:2121
SolverChoice::do_theta_advection
bool do_theta_advection
Whether custom vertical subsidence is applied to rho-theta.
Definition: ERF_DataStruct.H:1964
SolverChoice::spatial_moisture_forcing
bool spatial_moisture_forcing
Whether spatially varying moisture forcing is enabled.
Definition: ERF_DataStruct.H:1969
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:1848
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Sponge-layer options.
Definition: ERF_DataStruct.H:1863
SolverChoice::four_stream_radiation
bool four_stream_radiation
Whether the four-stream radiation approximation is enabled.
Definition: ERF_DataStruct.H:1916
SolverChoice::dampingChoice
DampingChoice dampingChoice
Damping-related options.
Definition: ERF_DataStruct.H:1862
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:1839
SolverChoice::buildings_type
static BuildingsType buildings_type
Building representation.
Definition: ERF_DataStruct.H:1842
SolverChoice::custom_rhotheta_forcing
bool custom_rhotheta_forcing
Whether custom rho-theta forcing is enabled.
Definition: ERF_DataStruct.H:1961
SolverChoice::spatial_rhotheta_forcing
bool spatial_rhotheta_forcing
Whether spatially varying rho-theta forcing is enabled.
Definition: ERF_DataStruct.H:1968
SolverChoice::use_source_perturbation
bool use_source_perturbation(int lev) const
Query whether source-term turbulent perturbations are enabled on a level.
Definition: ERF_DataStruct.H:2044
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Whether custom vertical subsidence is enabled.
Definition: ERF_DataStruct.H:1963
SolverChoice::custom_moisture_forcing
bool custom_moisture_forcing
Whether custom moisture forcing is enabled.
Definition: ERF_DataStruct.H:1962
SolverChoice::rad_type
RadiationType rad_type
Radiation model.
Definition: ERF_DataStruct.H:2128
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:1919
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Whether custom forcing operates on primitive variables.
Definition: ERF_DataStruct.H:1967
SolverChoice::use_num_diff
bool use_num_diff
Whether sixth-order numerical diffusion is enabled.
Definition: ERF_DataStruct.H:2120
SpongeChoice::sponge_type
static SpongeType sponge_type
Selected sponge damping model.
Definition: ERF_SpongeStruct.H:100
TurbChoice
Definition: ERF_TurbStruct.H:114
TurbChoice::diffuse_tke_3D
bool diffuse_tke_3D
Whether three-dimensional numerical diffusion is applied to TKE/QKE.
Definition: ERF_TurbStruct.H:736
TurbChoice::use_tke
bool use_tke
Whether any TKE or QKE closure is active.
Definition: ERF_TurbStruct.H:671
TurbulentPerturbation::pb_cell
amrex::Vector< amrex::MultiFab > pb_cell
Per-cell perturbation amplitude storage.
Definition: ERF_TurbPertStruct.H:762
TurbulentPerturbation::apply_tpi
void apply_tpi(const int &lev, const amrex::Box &vbx, const int &comp, const amrex::IndexType &m_ixtype, const amrex::Array4< amrex::Real > &src_arr, const amrex::Array4< amrex::Real const > &pert_cell)
Apply stored turbulent perturbations to a source or state array.
Definition: ERF_TurbPertStruct.H:408
Here is the call graph for this function: