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