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

Function Documentation

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

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]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]mapfac_mmap factor at cell centers
[in]mapfac_umap factor at x-faces
[in]mapfac_vmap factor at y-faces
[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
[in]n_qstatenumber of moisture components
71 {
72  BL_PROFILE_REGION("erf_make_mom_sources()");
73 
74  Box domain(geom.Domain());
75  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
76 
77  // Initialize sources to zero since we re-compute them ever RK stage
78  xmom_src.setVal(0.0);
79  ymom_src.setVal(0.0);
80  zmom_src.setVal(0.0);
81 
82  // *****************************************************************************
83  // Define source term for all three components of momenta from
84  // 1. buoyancy for (zmom)
85  // 2. Coriolis forcing for (xmom,ymom,zmom)
86  // 3. Rayleigh damping for (xmom,ymom,zmom)
87  // 4. Constant / height-dependent geostrophic forcing
88  // 5. subsidence
89  // 6. nudging towards input sounding data
90  // 7. numerical diffusion for (xmom,ymom,zmom)
91  // 8. sponge
92  // 9. Forest canopy
93  // 10. Immersed Forcing
94  // *****************************************************************************
95  //const bool l_use_ndiff = solverChoice.use_num_diff;
96 
97  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
98  if (solverChoice.do_forest_drag) {
99  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
100  }
101  }
102 
103 
104  // *****************************************************************************
105  // Data for Coriolis forcing
106  // *****************************************************************************
107  auto use_coriolis = solverChoice.use_coriolis;
108  auto coriolis_factor = solverChoice.coriolis_factor;
109  auto cosphi = solverChoice.cosphi;
110  auto sinphi = solverChoice.sinphi;
111  auto var_coriolis = solverChoice.variable_coriolis;
112 
113  // *****************************************************************************
114  // Flag for Geostrophic forcing
115  // *****************************************************************************
116  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
117  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
118 
119  // *****************************************************************************
120  // Data for Rayleigh damping
121  // *****************************************************************************
122  auto rayleigh_damp_U = solverChoice.rayleigh_damp_U;
123  auto rayleigh_damp_V = solverChoice.rayleigh_damp_V;
124  auto rayleigh_damp_W = solverChoice.rayleigh_damp_W;
125 
126  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
127  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
128  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
129 
130  // *****************************************************************************
131  // Planar averages for subsidence terms
132  // *****************************************************************************
133  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
134  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
135 
136  if (dptr_wbar_sub || solverChoice.nudging_from_input_sounding)
137  {
138  // Rho
139  PlaneAverage r_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true);
140  r_ave.compute_averages(ZDir(), r_ave.field());
141 
142  int ncell = r_ave.ncell_line();
143  Gpu::HostVector< Real> r_plane_h(ncell);
144  Gpu::DeviceVector< Real> r_plane_d(ncell);
145 
146  r_ave.line_average(Rho_comp, r_plane_h);
147 
148  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
149 
150  Real* dptr_r = r_plane_d.data();
151 
152  IntVect ng_c = S_data[IntVars::cons].nGrowVect();
153  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
154  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
155 
156  int offset = ng_c[2];
157  dptr_r_plane = r_plane_tab.table();
158  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
159  {
160  dptr_r_plane(k-offset) = dptr_r[k];
161  });
162 
163  // U and V momentum
164  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, true);
165  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, true);
166 
167  u_ave.compute_averages(ZDir(), u_ave.field());
168  v_ave.compute_averages(ZDir(), v_ave.field());
169 
170  int u_ncell = u_ave.ncell_line();
171  int v_ncell = v_ave.ncell_line();
172  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
173  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
174 
175  u_ave.line_average(0, u_plane_h);
176  v_ave.line_average(0, v_plane_h);
177 
178  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
179  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
180 
181  Real* dptr_u = u_plane_d.data();
182  Real* dptr_v = v_plane_d.data();
183 
184  IntVect ng_u = S_data[IntVars::xmom].nGrowVect();
185  IntVect ng_v = S_data[IntVars::ymom].nGrowVect();
186  Box udomain = domain; udomain.grow(2,ng_u[2]);
187  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
188  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
189  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
190 
191  int u_offset = ng_u[2];
192  dptr_u_plane = u_plane_tab.table();
193  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
194  {
195  dptr_u_plane(k-u_offset) = dptr_u[k];
196  });
197 
198  int v_offset = ng_v[2];
199  dptr_v_plane = v_plane_tab.table();
200  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
201  {
202  dptr_v_plane(k-v_offset) = dptr_v[k];
203  });
204  }
205 
206  // *****************************************************************************
207  // 1. Create the BUOYANCY forcing term in the z-direction
208  // *****************************************************************************
209  make_buoyancy(S_data, S_prim, zmom_src, geom, solverChoice, base_state,
210  n_qstate, solverChoice.anelastic[level]);
211 
212  // *****************************************************************************
213  // Add all the other forcings
214  // *****************************************************************************
215  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
216  {
217  Box tbx = mfi.nodaltilebox(0);
218  Box tby = mfi.nodaltilebox(1);
219  Box tbz = mfi.nodaltilebox(2);
220  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
221 
222  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
223  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
224  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
225  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
226 
227  const Array4<const Real>& u = xvel.array(mfi);
228  const Array4<const Real>& v = yvel.array(mfi);
229  const Array4<const Real>& w = wvel.array(mfi);
230 
231  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
232  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
233  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
234 
235  // Map factors
236  //const Array4<const Real>& mf_m = mapfac_m->const_array(mfi);
237  //const Array4<const Real>& mf_u = mapfac_u->const_array(mfi);
238  //const Array4<const Real>& mf_v = mapfac_v->const_array(mfi);
239 
240  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
241  Array4<const Real>{};
242  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
243  Array4<const Real>{};
244 
245  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
246  Array4<const Real>{};
247  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
248  Array4<const Real>{};
249 
250  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
251  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
252 
253  // *****************************************************************************
254  // 2. Add CORIOLIS forcing (this assumes east is +x, north is +y)
255  // *****************************************************************************
256  if (use_coriolis) {
257  if (var_coriolis) {
258  ParallelFor(tbx, tby, tbz,
259  [=] AMREX_GPU_DEVICE (int i, int j, int k)
260  {
261  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));
262  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));
263  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
264  Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
265  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
266  },
267  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
268  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));
269  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
270  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
271  },
272  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
273  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));
274  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
275  });
276  } else {
277  ParallelFor(tbx, tby, tbz,
278  [=] AMREX_GPU_DEVICE (int i, int j, int k)
279  {
280  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));
281  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));
282  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
283  },
284  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
285  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));
286  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
287  },
288  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
289  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));
290  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
291  });
292  } // var_coriolis
293  } // use_coriolis
294 
295  // *****************************************************************************
296  // 3. Add RAYLEIGH damping
297  // *****************************************************************************
298  Real zlo = geom.ProbLo(2);
299  Real dz = geom.CellSize(2);
300  Real ztop = solverChoice.rayleigh_ztop;
301  Real zdamp = solverChoice.rayleigh_zdamp;
302  Real dampcoef = solverChoice.rayleigh_dampcoef;
303 
304  if (rayleigh_damp_U) {
305  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
306  {
307  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
308  Real zfrac = 1 - (ztop - zcc) / zdamp;
309  if (zfrac > 0) {
310  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
311  Real uu = rho_u(i,j,k) / rho_on_u_face;
312  Real sinefac = std::sin(PIoTwo*zfrac);
313  xmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (uu - ubar[k]) * rho_on_u_face;
314  }
315  });
316  }
317 
318  if (rayleigh_damp_V) {
319  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
320  {
321  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
322  Real zfrac = 1 - (ztop - zcc) / zdamp;
323  if (zfrac > 0) {
324  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
325  Real vv = rho_v(i,j,k) / rho_on_v_face;
326  Real sinefac = std::sin(PIoTwo*zfrac);
327  ymom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (vv - vbar[k]) * rho_on_v_face;
328  }
329  });
330  }
331 
332  if (rayleigh_damp_W) {
333  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
334  {
335  Real zstag = (z_nd_arr) ? z_nd_arr(i,j,k) : zlo + k*dz;
336  Real zfrac = 1 - (ztop - zstag) / zdamp;
337  if (zfrac > 0) {
338  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
339  Real ww = rho_w(i,j,k) / rho_on_w_face;
340  Real sinefac = std::sin(PIoTwo*zfrac);
341  zmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (ww - wbar[k]) * rho_on_w_face;
342  }
343  });
344  }
345 
346  // *****************************************************************************
347  // 4. Add constant GEOSTROPHIC forcing
348  // *****************************************************************************
349  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
350  {
351  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
352  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
353  });
354  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
355  {
356  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
357  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
358  });
359  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
360  {
361  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
362  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
363  });
364 
365  // *****************************************************************************
366  // 4. Add height-dependent GEOSTROPHIC forcing
367  // *****************************************************************************
368  if (geo_wind_profile) {
369  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
370  {
371  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
372  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
373  });
374  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
375  {
376  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
377  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
378  });
379  } // geo_wind_profile
380 
381  // *****************************************************************************
382  // 5. Add custom SUBSIDENCE terms
383  // *****************************************************************************
384  if (solverChoice.custom_w_subsidence) {
385  if (solverChoice.custom_forcing_prim_vars) {
386  const int nr = Rho_comp;
387  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
388  {
389  Real dzInv = 0.5*dxInv[2];
390  if (z_nd_arr) {
391  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
392  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
393  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
394  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
395  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
396  }
397  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
398  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
399  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
400  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
401  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
402  });
403  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
404  {
405  Real dzInv = 0.5*dxInv[2];
406  if (z_nd_arr) {
407  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
408  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
409  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
410  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
411  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
412  }
413  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
414  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
415  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
416  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
417  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
418  });
419  } else {
420  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
421  {
422  Real dzInv = 0.5*dxInv[2];
423  if (z_nd_arr) {
424  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
425  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
426  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
427  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
428  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
429  }
430  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
431  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
432  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
433  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
434  });
435  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
436  {
437  Real dzInv = 0.5*dxInv[2];
438  if (z_nd_arr) {
439  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
440  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
441  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
442  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
443  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
444  }
445  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
446  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
447  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
448  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
449  });
450  }
451  }
452 
453  // *************************************************************************************
454  // 6. Add nudging towards value specified in input sounding
455  // *************************************************************************************
456  if (solverChoice.nudging_from_input_sounding)
457  {
458  int itime_n = 0;
459  int itime_np1 = 0;
460  Real coeff_n = Real(1.0);
461  Real coeff_np1 = Real(0.0);
462 
463  Real tau_inv = Real(1.0) / input_sounding_data.tau_nudging;
464 
465  int n_sounding_times = input_sounding_data.input_sounding_time.size();
466 
467  for (int nt = 1; nt < n_sounding_times; nt++) {
468  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
469  }
470  if (itime_n == n_sounding_times-1) {
471  itime_np1 = itime_n;
472  } else {
473  itime_np1 = itime_n+1;
474  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
475  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
476  coeff_n = Real(1.0) - coeff_np1;
477  }
478 
479  int nr = Rho_comp;
480 
481  const Real* u_inp_sound_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr();
482  const Real* u_inp_sound_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr();
483  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
484  {
485  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));
486  nudge_u *= tau_inv;
487  xmom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_u;
488  });
489 
490  const Real* v_inp_sound_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr();
491  const Real* v_inp_sound_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr();
492  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
493  {
494  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));
495  nudge_v *= tau_inv;
496  ymom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_v;
497  });
498  }
499 
500  // *****************************************************************************
501  // 7. Add NUMERICAL DIFFUSION terms
502  // *****************************************************************************
503 #if 0
504  if (l_use_ndiff) {
505  const Array4<const Real>& mf_u = mapfac_u->const_array(mfi);
506  const Array4<const Real>& mf_v = mapfac_v->const_array(mfi);
507  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
508  u, cell_data, xmom_src_arr, mf_u);
509  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
510  v, cell_data, ymom_src_arr, mf_v);
511  }
512 #endif
513 
514  // *****************************************************************************
515  // 8. Add SPONGING
516  // *****************************************************************************
517  if(solverChoice.spongeChoice.sponge_type == "input_sponge")
518  {
519  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
520  z_cc_arr, xmom_src_arr, ymom_src_arr,
521  rho_u, rho_v, d_sponge_ptrs_at_lev);
522  }
523  else
524  {
525  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
526  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
527  z_nd_arr, z_cc_arr);
528  }
529 
530  // *****************************************************************************
531  // 9. Add CANOPY source terms
532  // *****************************************************************************
533  if (solverChoice.do_forest_drag) {
534  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
535  {
536  const Real ux = u(i, j, k);
537  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
538  + v(i, j+1, k ) + v(i-1, j+1, k ) );
539  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
540  + w(i, j , k+1) + w(i-1, j , k+1) );
541  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
542  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
543  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
544  });
545  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
546  {
547  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
548  + u(i+1, j , k ) + u(i+1, j-1, k ) );
549  const Real uy = v(i, j, k);
550  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
551  + w(i , j , k+1) + w(i , j-1, k+1) );
552  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
553  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
554  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
555  });
556  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
557  {
558  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
559  + u(i , j , k-1) + u(i+1, j , k-1) );
560  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
561  + v(i , j , k-1) + v(i , j+1, k-1) );
562  const amrex::Real uz = w(i, j, k);
563  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
564  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
565  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
566  });
567  }
568  // *****************************************************************************
569  // 10. Add Immersed source terms
570  // *****************************************************************************
571  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
572  const Real drag_coefficient=10.0/dz;
573  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
574  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
575  {
576  const Real ux = u(i, j, k);
577  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
578  + v(i, j+1, k ) + v(i-1, j+1, k ) );
579  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
580  + w(i, j , k+1) + w(i-1, j , k+1) );
581  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
582  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
583  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
584  xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed;
585  });
586  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
587  {
588  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
589  + u(i+1, j , k ) + u(i+1, j-1, k ) );
590  const Real uy = v(i, j, k);
591  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
592  + w(i , j , k+1) + w(i , j-1, k+1) );
593  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
594  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
595  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
596  ymom_src_arr(i, j, k) -= t_blank * CdM * uy * windspeed;
597  });
598  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
599  {
600  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
601  + u(i , j , k-1) + u(i+1, j , k-1) );
602  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
603  + v(i , j , k-1) + v(i , j+1, k-1) );
604  const amrex::Real uz = w(i, j, k);
605  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
606  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
607  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
608  zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed;
609  });
610  }
611  } // mfi
612 }
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 > &z_phys_nd, const Array4< const Real > &z_phys_cc)
Definition: ERF_ApplySpongeZoneBCs.cpp:115
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:74
wbar
@ wbar
Definition: ERF_DataStruct.H:74
vbar
@ vbar
Definition: ERF_DataStruct.H:74
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
make_buoyancy
void make_buoyancy(Vector< MultiFab > &S_data, const MultiFab &S_prim, MultiFab &buoyancy, const amrex::Geometry geom, const SolverChoice &solverChoice, const MultiFab &base_state, const int n_qstate, const int anelastic)
Definition: ERF_MakeBuoyancy.cpp:32
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 > &mf_arr)
Definition: ERF_NumericalDiffusion.cpp:84
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 > &mf_arr)
Definition: ERF_NumericalDiffusion.cpp:149
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
PlaneAverage
Definition: ERF_PlaneAverage.H:14
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:152
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:150
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:153
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:151
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142
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:721
SolverChoice::rayleigh_damp_V
bool rayleigh_damp_V
Definition: ERF_DataStruct.H:698
SolverChoice::rayleigh_zdamp
amrex::Real rayleigh_zdamp
Definition: ERF_DataStruct.H:702
SolverChoice::cosphi
amrex::Real cosphi
Definition: ERF_DataStruct.H:722
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:766
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Definition: ERF_DataStruct.H:728
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:734
SolverChoice::rayleigh_damp_U
bool rayleigh_damp_U
Definition: ERF_DataStruct.H:697
SolverChoice::rayleigh_ztop
amrex::Real rayleigh_ztop
Definition: ERF_DataStruct.H:703
SolverChoice::sinphi
amrex::Real sinphi
Definition: ERF_DataStruct.H:723
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:768
SolverChoice::use_coriolis
bool use_coriolis
Definition: ERF_DataStruct.H:694
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:750
SolverChoice::variable_coriolis
bool variable_coriolis
Definition: ERF_DataStruct.H:771
SolverChoice::anelastic
amrex::Vector< int > anelastic
Definition: ERF_DataStruct.H:681
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:730
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:659
SolverChoice::rayleigh_damp_W
bool rayleigh_damp_W
Definition: ERF_DataStruct.H:699
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:675
SolverChoice::do_forest_drag
bool do_forest_drag
Definition: ERF_DataStruct.H:793
SolverChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Definition: ERF_DataStruct.H:701
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:773
SpongeChoice::sponge_type
std::string sponge_type
Definition: ERF_SpongeStruct.H:58
Here is the call graph for this function: