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