ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_MakeMomSources.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_NumericalDiffusion.H"
#include "ERF_PlaneAverage.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_SrcHeaders.H"
#include "ERF_Utils.H"
Include dependency graph for ERF_MakeMomSources.cpp:

Functions

void make_mom_sources (Real time, const Vector< MultiFab > &S_data, MultiFab &z_phys_nd, MultiFab &z_phys_cc, Vector< Real > &stretched_dz_h, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &wvel, MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, const MultiFab &base_state, MultiFab *forest_drag, MultiFab *terrain_blank, MultiFab *cosPhi_mf, MultiFab *sinPhi_mf, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &, const Real *dptr_u_geos, const Real *dptr_v_geos, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const Vector< Real * > d_sponge_ptrs_at_lev, InputSoundingData &input_sounding_data, bool is_slow_step)
 

Function Documentation

◆ make_mom_sources()

void make_mom_sources ( Real  time,
const Vector< MultiFab > &  S_data,
MultiFab &  z_phys_nd,
MultiFab &  z_phys_cc,
Vector< Real > &  stretched_dz_h,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  wvel,
MultiFab &  xmom_src,
MultiFab &  ymom_src,
MultiFab &  zmom_src,
const MultiFab &  base_state,
MultiFab *  forest_drag,
MultiFab *  terrain_blank,
MultiFab *  cosPhi_mf,
MultiFab *  sinPhi_mf,
const Geometry  geom,
const SolverChoice &  solverChoice,
Vector< std::unique_ptr< MultiFab >> &  ,
const Real *  dptr_u_geos,
const Real *  dptr_v_geos,
const Real *  dptr_wbar_sub,
const Vector< Real * >  d_rayleigh_ptrs_at_lev,
const Vector< Real * >  d_sponge_ptrs_at_lev,
InputSoundingData &  input_sounding_data,
bool  is_slow_step 
)

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

Parameters
[in]S_datacurrent solution
[in]xvelx-component of velocity
[in]yvely-component of velocity
[in]xmom_srcsource terms for x-momentum
[in]ymom_srcsource terms for y-momentum
[in]zmom_srcsource terms for z-momentum
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]mapfacmap factors
[in]dptr_u_geoscustom geostrophic wind profile
[in]dptr_v_geoscustom geostrophic wind profile
[in]dptr_wbar_subsubsidence source term
[in]d_rayleigh_ptrs_at_levVector of {strength of Rayleigh damping, reference value for xvel/yvel/zvel/theta} used to define Rayleigh damping
59 {
60  BL_PROFILE_REGION("erf_make_mom_sources()");
61 
62  Box domain(geom.Domain());
63  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
64 
65  // Initialize sources to zero each time we may use them
66  xmom_src.setVal(0.0);
67  ymom_src.setVal(0.0);
68  zmom_src.setVal(0.0);
69 
70  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
71 
72  // flags to apply certain source terms in substep call only
73  bool use_Rayleigh_fast = solverChoice.rayleigh_damp_substep;
74  bool use_canopy_fast = solverChoice.forest_substep;
75  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
76 
77  // *****************************************************************************
78  // Define source term for all three components of momenta from
79  // 1. Coriolis forcing for (xmom,ymom,zmom)
80  // 2. Rayleigh damping for (xmom,ymom,zmom)
81  // 3. Constant / height-dependent geostrophic forcing
82  // 4. Subsidence
83  // 5. Nudging towards input sounding data
84  // 6. Numerical diffusion for (xmom,ymom,zmom)
85  // 7. Sponge
86  // 8. Forest canopy
87  // 9. Immersed forcing
88  // 10. Constant mass flux
89  // *****************************************************************************
90  // NOTE: buoyancy is now computed in a separate routine - it should not appear here
91  // *****************************************************************************
92  //const bool l_use_ndiff = solverChoice.use_num_diff;
93 
94  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
95  if (solverChoice.do_forest_drag) {
96  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
97  }
98  }
99 
100 
101  // *****************************************************************************
102  // Data for Coriolis forcing
103  // *****************************************************************************
104  auto use_coriolis = solverChoice.use_coriolis;
105  auto coriolis_factor = solverChoice.coriolis_factor;
106  auto cosphi = solverChoice.cosphi;
107  auto sinphi = solverChoice.sinphi;
108  auto var_coriolis = solverChoice.variable_coriolis;
109  auto has_lat_lon = solverChoice.has_lat_lon;
110 
111  // *****************************************************************************
112  // Flag for Geostrophic forcing
113  // *****************************************************************************
114  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
115  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
116 
117  // *****************************************************************************
118  // Data for Rayleigh damping
119  // *****************************************************************************
120  auto rayleigh_damp_U = solverChoice.rayleigh_damp_U;
121  auto rayleigh_damp_V = solverChoice.rayleigh_damp_V;
122  auto rayleigh_damp_W = solverChoice.rayleigh_damp_W;
123 
124  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
125  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
126  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
127 
128  // *****************************************************************************
129  // Data for constant mass flux
130  // *****************************************************************************
131  bool enforce_massflux_x = (solverChoice.const_massflux_u != 0);
132  bool enforce_massflux_y = (solverChoice.const_massflux_v != 0);
133  Real U_target = solverChoice.const_massflux_u;
134  Real V_target = solverChoice.const_massflux_v;
135  int massflux_klo = solverChoice.massflux_klo;
136  int massflux_khi = solverChoice.massflux_khi;
137 
138  // These will be updated by integrating through the planar average profiles
139  Real rhoUA_target{0};
140  Real rhoVA_target{0};
141  Real rhoUA{0};
142  Real rhoVA{0};
143 
144  // *****************************************************************************
145  // Planar averages for subsidence, nudging, or constant mass flux
146  // *****************************************************************************
147  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
148  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
149 
150  if (is_slow_step && (dptr_wbar_sub || solverChoice.nudging_from_input_sounding ||
151  enforce_massflux_x || enforce_massflux_y))
152  {
153  const int offset = 1;
154  const int u_offset = 1;
155  const int v_offset = 1;
156 
157  //
158  // We use the alias here to control ncomp inside the PlaneAverage
159  //
160  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
161 
162  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = offset;
163  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
164  r_ave.compute_averages(ZDir(), r_ave.field());
165 
166  int ncell = r_ave.ncell_line();
167  Gpu::HostVector< Real> r_plane_h(ncell);
168  Gpu::DeviceVector< Real> r_plane_d(ncell);
169 
170  r_ave.line_average(Rho_comp, r_plane_h);
171 
172  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
173 
174  Real* dptr_r = r_plane_d.data();
175 
176  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
177  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
178 
179  dptr_r_plane = r_plane_tab.table();
180  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
181  {
182  dptr_r_plane(k-offset) = dptr_r[k];
183  });
184 
185  // U and V momentum
186  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
187  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
188 
189  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
190  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
191 
192  u_ave.compute_averages(ZDir(), u_ave.field());
193  v_ave.compute_averages(ZDir(), v_ave.field());
194 
195  int u_ncell = u_ave.ncell_line();
196  int v_ncell = v_ave.ncell_line();
197  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
198  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
199 
200  u_ave.line_average(0, u_plane_h);
201  v_ave.line_average(0, v_plane_h);
202 
203  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
204  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
205 
206  Real* dptr_u = u_plane_d.data();
207  Real* dptr_v = v_plane_d.data();
208 
209  Box udomain = domain; udomain.grow(2,ng_u[2]);
210  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
211  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
212  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
213 
214  dptr_u_plane = u_plane_tab.table();
215  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
216  {
217  dptr_u_plane(k-u_offset) = dptr_u[k];
218  });
219 
220  dptr_v_plane = v_plane_tab.table();
221  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
222  {
223  dptr_v_plane(k-v_offset) = dptr_v[k];
224  });
225 
226  // sum in z for massflux adjustment
227  if (enforce_massflux_x || enforce_massflux_y) {
228  Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
229  Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
230 
231  if (solverChoice.mesh_type == MeshType::ConstantDz) {
232  // note: massflux_khi corresponds to unstaggered indices in this case
233  rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
234  u_plane_h.begin() + u_offset + massflux_khi+1, 0.0);
235  rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
236  v_plane_h.begin() + v_offset + massflux_khi+1, 0.0);
237  rhoUA_target = std::accumulate(r_plane_h.begin() + offset + massflux_klo,
238  r_plane_h.begin() + offset + massflux_khi+1, 0.0);
239  rhoVA_target = rhoUA_target;
240 
241  rhoUA *= geom.CellSize(2) * Ly;
242  rhoVA *= geom.CellSize(2) * Lx;
243  rhoUA_target *= geom.CellSize(2) * Ly;
244  rhoVA_target *= geom.CellSize(2) * Lx;
245 
246  } else if (solverChoice.mesh_type == MeshType::StretchedDz) {
247  // note: massflux_khi corresponds to staggered indices in this case
248  for (int k=massflux_klo; k < massflux_khi; ++k) {
249  rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
250  rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
251  rhoUA_target += r_plane_h[k + offset] * stretched_dz_h[k];
252  }
253  rhoVA_target = rhoUA_target;
254 
255  rhoUA *= Ly;
256  rhoVA *= Lx;
257  rhoUA_target *= Ly;
258  rhoVA_target *= Lx;
259  }
260 
261  // at this point, this is integrated rho*dA
262  rhoUA_target *= U_target;
263  rhoVA_target *= V_target;
264 
265  Print() << "Integrated mass flux : " << rhoUA << " " << rhoVA
266  << " (target: " << rhoUA_target << " " << rhoVA_target << ")"
267  << std::endl;
268  }
269  }
270 
271  // *****************************************************************************
272  // Add all the other forcings
273  // *****************************************************************************
274  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
275  {
276  Box tbx = mfi.nodaltilebox(0);
277  Box tby = mfi.nodaltilebox(1);
278  Box tbz = mfi.nodaltilebox(2);
279  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
280 
281  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
282  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
283  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
284  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
285 
286  const Array4<const Real>& u = xvel.array(mfi);
287  const Array4<const Real>& v = yvel.array(mfi);
288  const Array4<const Real>& w = wvel.array(mfi);
289 
290  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
291  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
292  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
293 
294  const Array4<const Real>& r0 = r_hse.const_array(mfi);
295 
296  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
297  Array4<const Real>{};
298  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
299  Array4<const Real>{};
300 
301  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
302  Array4<const Real>{};
303  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
304  Array4<const Real>{};
305 
306  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
307  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
308 
309  // *****************************************************************************
310  // 1. Add CORIOLIS forcing (this assumes east is +x, north is +y)
311  // *****************************************************************************
312  if (use_coriolis && is_slow_step) {
313  if (var_coriolis && has_lat_lon) {
314  ParallelFor(tbx, tby, tbz,
315  [=] AMREX_GPU_DEVICE (int i, int j, int k)
316  {
317  Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k));
318  Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k));
319  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
320  Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
321  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
322  },
323  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
324  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k));
325  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
326  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
327  },
328  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
329  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1));
330  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
331  });
332  } else {
333  ParallelFor(tbx, tby, tbz,
334  [=] AMREX_GPU_DEVICE (int i, int j, int k)
335  {
336  Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k));
337  Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k));
338  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
339  },
340  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
341  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k));
342  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
343  },
344  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
345  Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j,k-1) + rho_u(i,j,k-1));
346  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
347  });
348  } // var_coriolis
349  } // use_coriolis
350 
351  // *****************************************************************************
352  // 2. Add RAYLEIGH damping
353  // *****************************************************************************
354  Real zlo = geom.ProbLo(2);
355  Real dz = geom.CellSize(2);
356  Real ztop = solverChoice.rayleigh_ztop;
357  Real zdamp = solverChoice.rayleigh_zdamp;
358  Real dampcoef = solverChoice.rayleigh_dampcoef;
359 
360  if ((is_slow_step && !use_Rayleigh_fast) || (!is_slow_step && use_Rayleigh_fast)) {
361  if (rayleigh_damp_U) {
362  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
363  {
364  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
365  Real zfrac = 1 - (ztop - zcc) / zdamp;
366  if (zfrac > 0) {
367  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
368  Real uu = rho_u(i,j,k) / rho_on_u_face;
369  Real sinefac = std::sin(PIoTwo*zfrac);
370  xmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (uu - ubar[k]) * rho_on_u_face;
371  }
372  });
373  }
374 
375  if (rayleigh_damp_V) {
376  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
377  {
378  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
379  Real zfrac = 1 - (ztop - zcc) / zdamp;
380  if (zfrac > 0) {
381  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
382  Real vv = rho_v(i,j,k) / rho_on_v_face;
383  Real sinefac = std::sin(PIoTwo*zfrac);
384  ymom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (vv - vbar[k]) * rho_on_v_face;
385  }
386  });
387  }
388 
389  if (rayleigh_damp_W) {
390  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
391  {
392  Real zstag = (z_nd_arr) ? z_nd_arr(i,j,k) : zlo + k*dz;
393  Real zfrac = 1 - (ztop - zstag) / zdamp;
394  if (zfrac > 0) {
395  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
396  Real ww = rho_w(i,j,k) / rho_on_w_face;
397  Real sinefac = std::sin(PIoTwo*zfrac);
398  zmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (ww - wbar[k]) * rho_on_w_face;
399  }
400  });
401  }
402  } // fast or slow step
403 
404  // *****************************************************************************
405  // 3a. Add constant GEOSTROPHIC forcing
406  // *****************************************************************************
407  if (is_slow_step) {
408  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
409  {
410  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
411  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
412  });
413  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
414  {
415  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
416  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
417  });
418  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
419  {
420  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
421  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
422  });
423  }
424 
425  // *****************************************************************************
426  // 3b. Add height-dependent GEOSTROPHIC forcing
427  // *****************************************************************************
428  if (geo_wind_profile && is_slow_step) {
429  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
430  {
431  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
432  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
433  });
434  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
435  {
436  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
437  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
438  });
439  } // geo_wind_profile
440 
441  // *****************************************************************************
442  // 4. Add custom SUBSIDENCE terms
443  // *****************************************************************************
444  if (solverChoice.custom_w_subsidence && is_slow_step) {
445  if (solverChoice.custom_forcing_prim_vars) {
446  const int nr = Rho_comp;
447  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
448  {
449  Real dzInv = 0.5*dxInv[2];
450  if (z_nd_arr) {
451  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
452  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
453  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
454  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
455  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
456  }
457  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
458  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
459  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
460  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
461  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
462  });
463  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
464  {
465  Real dzInv = 0.5*dxInv[2];
466  if (z_nd_arr) {
467  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
468  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
469  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
470  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
471  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
472  }
473  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
474  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
475  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
476  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
477  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
478  });
479  } else {
480  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
481  {
482  Real dzInv = 0.5*dxInv[2];
483  if (z_nd_arr) {
484  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
485  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
486  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
487  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
488  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
489  }
490  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
491  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
492  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
493  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
494  });
495  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
496  {
497  Real dzInv = 0.5*dxInv[2];
498  if (z_nd_arr) {
499  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
500  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
501  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
502  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
503  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
504  }
505  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
506  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
507  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
508  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
509  });
510  }
511  }
512 
513  // *************************************************************************************
514  // 5. Add nudging towards value specified in input sounding
515  // *************************************************************************************
516  if (solverChoice.nudging_from_input_sounding && is_slow_step)
517  {
518  int itime_n = 0;
519  int itime_np1 = 0;
520  Real coeff_n = Real(1.0);
521  Real coeff_np1 = Real(0.0);
522 
523  Real tau_inv = Real(1.0) / input_sounding_data.tau_nudging;
524 
525  int n_sounding_times = input_sounding_data.input_sounding_time.size();
526 
527  for (int nt = 1; nt < n_sounding_times; nt++) {
528  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
529  }
530  if (itime_n == n_sounding_times-1) {
531  itime_np1 = itime_n;
532  } else {
533  itime_np1 = itime_n+1;
534  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
535  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
536  coeff_n = Real(1.0) - coeff_np1;
537  }
538 
539  int nr = Rho_comp;
540 
541  const Real* u_inp_sound_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr();
542  const Real* u_inp_sound_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr();
543  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
544  {
545  Real nudge_u = (coeff_n*u_inp_sound_n[k] + coeff_np1*u_inp_sound_np1[k]) - (dptr_u_plane(k)/dptr_r_plane(k));
546  nudge_u *= tau_inv;
547  xmom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_u;
548  });
549 
550  const Real* v_inp_sound_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr();
551  const Real* v_inp_sound_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr();
552  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
553  {
554  Real nudge_v = (coeff_n*v_inp_sound_n[k] + coeff_np1*v_inp_sound_np1[k]) - (dptr_v_plane(k)/dptr_r_plane(k));
555  nudge_v *= tau_inv;
556  ymom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_v;
557  });
558  }
559 
560  // *****************************************************************************
561  // 6. Add NUMERICAL DIFFUSION terms
562  // *****************************************************************************
563 #if 0
564  if (l_use_ndiff) {
565  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
566  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
567  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
568  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
569  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
570  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
571  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
572  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
573  }
574 #endif
575 
576  // *****************************************************************************
577  // 7. Add SPONGING
578  // *****************************************************************************
579  if (is_slow_step) {
580  if (solverChoice.spongeChoice.sponge_type == "input_sponge")
581  {
582  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
583  z_cc_arr, xmom_src_arr, ymom_src_arr,
584  rho_u, rho_v, d_sponge_ptrs_at_lev);
585  }
586  else
587  {
588  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
589  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
590  r0, z_nd_arr, z_cc_arr);
591  }
592  }
593 
594  // *****************************************************************************
595  // 8. Add CANOPY source terms
596  // *****************************************************************************
597  if (solverChoice.do_forest_drag &&
598  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
599  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
600  {
601  const Real ux = u(i, j, k);
602  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
603  + v(i, j+1, k ) + v(i-1, j+1, k ) );
604  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
605  + w(i, j , k+1) + w(i-1, j , k+1) );
606  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
607  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
608  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
609  });
610  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
611  {
612  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
613  + u(i+1, j , k ) + u(i+1, j-1, k ) );
614  const Real uy = v(i, j, k);
615  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
616  + w(i , j , k+1) + w(i , j-1, k+1) );
617  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
618  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
619  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
620  });
621  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
622  {
623  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
624  + u(i , j , k-1) + u(i+1, j , k-1) );
625  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
626  + v(i , j , k-1) + v(i , j+1, k-1) );
627  const amrex::Real uz = w(i, j, k);
628  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
629  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
630  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
631  });
632  }
633  // *****************************************************************************
634  // 9. Add Immersed source terms
635  // *****************************************************************************
636  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
637  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
638  const Real drag_coefficient=10.0/dz;
639  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
640  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
641  {
642  const Real ux = u(i, j, k);
643  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
644  + v(i, j+1, k ) + v(i-1, j+1, k ) );
645  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
646  + w(i, j , k+1) + w(i-1, j , k+1) );
647  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
648  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
649  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
650  xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed;
651  });
652  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
653  {
654  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
655  + u(i+1, j , k ) + u(i+1, j-1, k ) );
656  const Real uy = v(i, j, k);
657  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
658  + w(i , j , k+1) + w(i , j-1, k+1) );
659  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
660  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
661  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
662  ymom_src_arr(i, j, k) -= t_blank * CdM * uy * windspeed;
663  });
664  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
665  {
666  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
667  + u(i , j , k-1) + u(i+1, j , k-1) );
668  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
669  + v(i , j , k-1) + v(i , j+1, k-1) );
670  const amrex::Real uz = w(i, j, k);
671  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
672  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
673  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
674  zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed;
675  });
676  }
677 
678  // *****************************************************************************
679  // 10. Enforce constant mass flux
680  // *****************************************************************************
681  if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
682  Real tau_inv = Real(1.0) / solverChoice.const_massflux_tau;
683 
684  ParallelFor(tbx, tby,
685  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
686  xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
687  },
688  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
689  ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
690  });
691  }
692  } // mfi
693 }
ApplySpongeZoneBCsForMom
void ApplySpongeZoneBCsForMom(const SpongeChoice &spongeChoice, const Geometry geom, const Box &tbx, const Box &tby, const Box &tbz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &rho_w, const Array4< const Real > &r0, const Array4< const Real > &z_phys_nd, const Array4< const Real > &z_phys_cc)
Definition: ERF_ApplySpongeZoneBCs.cpp:118
ApplySpongeZoneBCsForMom_ReadFromFile
void ApplySpongeZoneBCsForMom_ReadFromFile(const SpongeChoice &spongeChoice, const Geometry geom, const Box &tbx, const Box &tby, const Array4< const Real > &cell_data, const Array4< const Real > &z_phys_cc, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Vector< Real * > d_sponge_ptrs_at_lev)
Definition: ERF_ApplySpongeZoneBCs_ReadFromFile.cpp:8
PIoTwo
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:7
ubar
@ ubar
Definition: ERF_DataStruct.H:87
wbar
@ wbar
Definition: ERF_DataStruct.H:87
vbar
@ vbar
Definition: ERF_DataStruct.H:87
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
NumericalDiffusion_Ymom
void NumericalDiffusion_Ymom(const Box &bx, const Real 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:151
NumericalDiffusion_Xmom
void NumericalDiffusion_Xmom(const Box &bx, const Real 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:85
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
PlaneAverage
Definition: ERF_PlaneAverage.H:14
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:63
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:160
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:158
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:161
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:159
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:45
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
InputSoundingData::V_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Definition: ERF_InputSoundingData.H:326
InputSoundingData::input_sounding_time
amrex::Vector< amrex::Real > input_sounding_time
Definition: ERF_InputSoundingData.H:315
InputSoundingData::tau_nudging
amrex::Real tau_nudging
Definition: ERF_InputSoundingData.H:312
InputSoundingData::U_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Definition: ERF_InputSoundingData.H:326
SolverChoice::coriolis_factor
amrex::Real coriolis_factor
Definition: ERF_DataStruct.H:769
SolverChoice::mesh_type
static MeshType mesh_type
Definition: ERF_DataStruct.H:708
SolverChoice::rayleigh_damp_V
bool rayleigh_damp_V
Definition: ERF_DataStruct.H:741
SolverChoice::rayleigh_damp_substep
bool rayleigh_damp_substep
Definition: ERF_DataStruct.H:747
SolverChoice::rayleigh_zdamp
amrex::Real rayleigh_zdamp
Definition: ERF_DataStruct.H:745
SolverChoice::const_massflux_v
amrex::Real const_massflux_v
Definition: ERF_DataStruct.H:844
SolverChoice::cosphi
amrex::Real cosphi
Definition: ERF_DataStruct.H:770
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:812
SolverChoice::massflux_klo
int massflux_klo
Definition: ERF_DataStruct.H:848
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Definition: ERF_DataStruct.H:776
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:782
SolverChoice::rayleigh_damp_U
bool rayleigh_damp_U
Definition: ERF_DataStruct.H:740
SolverChoice::rayleigh_ztop
amrex::Real rayleigh_ztop
Definition: ERF_DataStruct.H:746
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:750
SolverChoice::sinphi
amrex::Real sinphi
Definition: ERF_DataStruct.H:771
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:814
SolverChoice::const_massflux_u
amrex::Real const_massflux_u
Definition: ERF_DataStruct.H:843
SolverChoice::use_coriolis
bool use_coriolis
Definition: ERF_DataStruct.H:737
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:795
SolverChoice::variable_coriolis
bool variable_coriolis
Definition: ERF_DataStruct.H:817
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:778
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:702
SolverChoice::rayleigh_damp_W
bool rayleigh_damp_W
Definition: ERF_DataStruct.H:742
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:718
SolverChoice::has_lat_lon
bool has_lat_lon
Definition: ERF_DataStruct.H:816
SolverChoice::do_forest_drag
bool do_forest_drag
Definition: ERF_DataStruct.H:840
SolverChoice::const_massflux_tau
amrex::Real const_massflux_tau
Definition: ERF_DataStruct.H:845
SolverChoice::massflux_khi
int massflux_khi
Definition: ERF_DataStruct.H:849
SolverChoice::forest_substep
bool forest_substep
Definition: ERF_DataStruct.H:751
SolverChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Definition: ERF_DataStruct.H:744
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:819
SpongeChoice::sponge_type
std::string sponge_type
Definition: ERF_SpongeStruct.H:58
Here is the call graph for this function: