ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_MakeMomSources.cpp File Reference
#include <AMReX_MultiFab.H>
#include <AMReX_ArrayLim.H>
#include <AMReX_BCRec.H>
#include <AMReX_TableData.H>
#include <AMReX_GpuContainers.H>
#include "ERF_NumericalDiffusion.H"
#include "ERF_PlaneAverage.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_SrcHeaders.H"
#include "ERF_Utils.H"
Include dependency graph for ERF_MakeMomSources.cpp:

Functions

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

Function Documentation

◆ make_mom_sources()

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

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

Parameters
[in]S_datacurrent solution
[in]xvelx-component of velocity
[in]yvely-component of velocity
[in]xmom_srcsource terms for x-momentum
[in]ymom_srcsource terms for y-momentum
[in]zmom_srcsource terms for z-momentum
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]mapfacmap factors
[in]dptr_u_geoscustom geostrophic wind profile
[in]dptr_v_geoscustom geostrophic wind profile
[in]dptr_wbar_subsubsidence source term
[in]d_rayleigh_ptrs_at_levVector of {strength of Rayleigh damping, reference value for xvel/yvel/zvel/theta} used to define Rayleigh damping
58 {
59  BL_PROFILE_REGION("erf_make_mom_sources()");
60 
61  Box domain(geom.Domain());
62  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
63 
64  // Initialize sources to zero each time we may use them
65  xmom_src.setVal(0.0);
66  ymom_src.setVal(0.0);
67  zmom_src.setVal(0.0);
68 
69  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
70 
71  // flags to apply certain source terms in substep call only
72  bool use_Rayleigh_fast = solverChoice.rayleigh_damp_substep;
73  bool use_canopy_fast = solverChoice.forest_substep;
74  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
75 
76  // *****************************************************************************
77  // Define source term for all three components of momenta from
78  // 1. Coriolis forcing for (xmom,ymom,zmom)
79  // 2. Rayleigh damping for (xmom,ymom,zmom)
80  // 3. Constant / height-dependent geostrophic forcing
81  // 4. subsidence
82  // 5. nudging towards input sounding data
83  // 6. numerical diffusion for (xmom,ymom,zmom)
84  // 7. sponge
85  // 8. Forest canopy
86  // 9. Immersed Forcing
87  // *****************************************************************************
88  // NOTE: buoyancy is now computed in a separate routine - it should not appear here
89  // *****************************************************************************
90  //const bool l_use_ndiff = solverChoice.use_num_diff;
91 
92  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
93  if (solverChoice.do_forest_drag) {
94  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
95  }
96  }
97 
98 
99  // *****************************************************************************
100  // Data for Coriolis forcing
101  // *****************************************************************************
102  auto use_coriolis = solverChoice.use_coriolis;
103  auto coriolis_factor = solverChoice.coriolis_factor;
104  auto cosphi = solverChoice.cosphi;
105  auto sinphi = solverChoice.sinphi;
106  auto var_coriolis = solverChoice.variable_coriolis;
107  auto has_lat_lon = solverChoice.has_lat_lon;
108 
109  // *****************************************************************************
110  // Flag for Geostrophic forcing
111  // *****************************************************************************
112  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
113  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
114 
115  // *****************************************************************************
116  // Data for Rayleigh damping
117  // *****************************************************************************
118  auto rayleigh_damp_U = solverChoice.rayleigh_damp_U;
119  auto rayleigh_damp_V = solverChoice.rayleigh_damp_V;
120  auto rayleigh_damp_W = solverChoice.rayleigh_damp_W;
121 
122  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
123  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
124  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
125 
126  // *****************************************************************************
127  // Planar averages for subsidence terms
128  // *****************************************************************************
129  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
130  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
131 
132  if (is_slow_step && (dptr_wbar_sub || solverChoice.nudging_from_input_sounding))
133  {
134  //
135  // We use the alias here to control ncomp inside the PlaneAverage
136  //
137  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
138 
139  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = 1;
140  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
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  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  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = 1;
165  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
166 
167  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = 1;
168  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
169 
170  u_ave.compute_averages(ZDir(), u_ave.field());
171  v_ave.compute_averages(ZDir(), v_ave.field());
172 
173  int u_ncell = u_ave.ncell_line();
174  int v_ncell = v_ave.ncell_line();
175  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
176  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
177 
178  u_ave.line_average(0, u_plane_h);
179  v_ave.line_average(0, v_plane_h);
180 
181  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
182  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
183 
184  Real* dptr_u = u_plane_d.data();
185  Real* dptr_v = v_plane_d.data();
186 
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  // Add all the other forcings
209  // *****************************************************************************
210  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
211  {
212  Box tbx = mfi.nodaltilebox(0);
213  Box tby = mfi.nodaltilebox(1);
214  Box tbz = mfi.nodaltilebox(2);
215  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
216 
217  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
218  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
219  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
220  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
221 
222  const Array4<const Real>& u = xvel.array(mfi);
223  const Array4<const Real>& v = yvel.array(mfi);
224  const Array4<const Real>& w = wvel.array(mfi);
225 
226  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
227  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
228  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
229 
230  const Array4<const Real>& r0 = r_hse.const_array(mfi);
231 
232  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
233  Array4<const Real>{};
234  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
235  Array4<const Real>{};
236 
237  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
238  Array4<const Real>{};
239  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
240  Array4<const Real>{};
241 
242  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
243  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
244 
245  // *****************************************************************************
246  // 2. Add CORIOLIS forcing (this assumes east is +x, north is +y)
247  // *****************************************************************************
248  if (use_coriolis && is_slow_step) {
249  if (var_coriolis && has_lat_lon) {
250  ParallelFor(tbx, tby, tbz,
251  [=] AMREX_GPU_DEVICE (int i, int j, int k)
252  {
253  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));
254  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));
255  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
256  Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
257  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
258  },
259  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
260  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));
261  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
262  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
263  },
264  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
265  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));
266  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
267  });
268  } else {
269  ParallelFor(tbx, tby, tbz,
270  [=] AMREX_GPU_DEVICE (int i, int j, int k)
271  {
272  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));
273  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));
274  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
275  },
276  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
277  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));
278  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
279  },
280  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
281  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));
282  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
283  });
284  } // var_coriolis
285  } // use_coriolis
286 
287  // *****************************************************************************
288  // 3. Add RAYLEIGH damping
289  // *****************************************************************************
290  Real zlo = geom.ProbLo(2);
291  Real dz = geom.CellSize(2);
292  Real ztop = solverChoice.rayleigh_ztop;
293  Real zdamp = solverChoice.rayleigh_zdamp;
294  Real dampcoef = solverChoice.rayleigh_dampcoef;
295 
296  if ((is_slow_step && !use_Rayleigh_fast) || (!is_slow_step && use_Rayleigh_fast)) {
297  if (rayleigh_damp_U) {
298  ParallelFor(tbx, [=] 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_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
304  Real uu = rho_u(i,j,k) / rho_on_u_face;
305  Real sinefac = std::sin(PIoTwo*zfrac);
306  xmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (uu - ubar[k]) * rho_on_u_face;
307  }
308  });
309  }
310 
311  if (rayleigh_damp_V) {
312  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
313  {
314  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
315  Real zfrac = 1 - (ztop - zcc) / zdamp;
316  if (zfrac > 0) {
317  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
318  Real vv = rho_v(i,j,k) / rho_on_v_face;
319  Real sinefac = std::sin(PIoTwo*zfrac);
320  ymom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (vv - vbar[k]) * rho_on_v_face;
321  }
322  });
323  }
324 
325  if (rayleigh_damp_W) {
326  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
327  {
328  Real zstag = (z_nd_arr) ? z_nd_arr(i,j,k) : zlo + k*dz;
329  Real zfrac = 1 - (ztop - zstag) / zdamp;
330  if (zfrac > 0) {
331  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
332  Real ww = rho_w(i,j,k) / rho_on_w_face;
333  Real sinefac = std::sin(PIoTwo*zfrac);
334  zmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (ww - wbar[k]) * rho_on_w_face;
335  }
336  });
337  }
338  } // fast or slow step
339 
340  // *****************************************************************************
341  // 4. Add constant GEOSTROPHIC forcing
342  // *****************************************************************************
343  if (is_slow_step) {
344  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
345  {
346  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
347  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
348  });
349  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
350  {
351  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
352  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
353  });
354  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
355  {
356  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
357  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
358  });
359  }
360 
361  // *****************************************************************************
362  // 4. Add height-dependent GEOSTROPHIC forcing
363  // *****************************************************************************
364  if (geo_wind_profile && is_slow_step) {
365  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
366  {
367  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
368  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
369  });
370  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
371  {
372  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
373  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
374  });
375  } // geo_wind_profile
376 
377  // *****************************************************************************
378  // 5. Add custom SUBSIDENCE terms
379  // *****************************************************************************
380  if (solverChoice.custom_w_subsidence && is_slow_step) {
381  if (solverChoice.custom_forcing_prim_vars) {
382  const int nr = Rho_comp;
383  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
384  {
385  Real dzInv = 0.5*dxInv[2];
386  if (z_nd_arr) {
387  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
388  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
389  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
390  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
391  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
392  }
393  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
394  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
395  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
396  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
397  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
398  });
399  ParallelFor(tby, [=] 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_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
404  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
405  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
406  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
407  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
408  }
409  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
410  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
411  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
412  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
413  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
414  });
415  } else {
416  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
417  {
418  Real dzInv = 0.5*dxInv[2];
419  if (z_nd_arr) {
420  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
421  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
422  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
423  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
424  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
425  }
426  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
427  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
428  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
429  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
430  });
431  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
432  {
433  Real dzInv = 0.5*dxInv[2];
434  if (z_nd_arr) {
435  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
436  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
437  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
438  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
439  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
440  }
441  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
442  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
443  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
444  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
445  });
446  }
447  }
448 
449  // *************************************************************************************
450  // 6. Add nudging towards value specified in input sounding
451  // *************************************************************************************
452  if (solverChoice.nudging_from_input_sounding && is_slow_step)
453  {
454  int itime_n = 0;
455  int itime_np1 = 0;
456  Real coeff_n = Real(1.0);
457  Real coeff_np1 = Real(0.0);
458 
459  Real tau_inv = Real(1.0) / input_sounding_data.tau_nudging;
460 
461  int n_sounding_times = input_sounding_data.input_sounding_time.size();
462 
463  for (int nt = 1; nt < n_sounding_times; nt++) {
464  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
465  }
466  if (itime_n == n_sounding_times-1) {
467  itime_np1 = itime_n;
468  } else {
469  itime_np1 = itime_n+1;
470  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
471  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
472  coeff_n = Real(1.0) - coeff_np1;
473  }
474 
475  int nr = Rho_comp;
476 
477  const Real* u_inp_sound_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr();
478  const Real* u_inp_sound_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr();
479  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
480  {
481  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));
482  nudge_u *= tau_inv;
483  xmom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_u;
484  });
485 
486  const Real* v_inp_sound_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr();
487  const Real* v_inp_sound_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr();
488  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
489  {
490  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));
491  nudge_v *= tau_inv;
492  ymom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_v;
493  });
494  }
495 
496  // *****************************************************************************
497  // 7. Add NUMERICAL DIFFUSION terms
498  // *****************************************************************************
499 #if 0
500  if (l_use_ndiff) {
501  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
502  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
503  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
504  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
505  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
506  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
507  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
508  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
509  }
510 #endif
511 
512  // *****************************************************************************
513  // 8. Add SPONGING
514  // *****************************************************************************
515  if (is_slow_step) {
516  if (solverChoice.spongeChoice.sponge_type == "input_sponge")
517  {
518  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
519  z_cc_arr, xmom_src_arr, ymom_src_arr,
520  rho_u, rho_v, d_sponge_ptrs_at_lev);
521  }
522  else
523  {
524  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
525  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
526  r0, z_nd_arr, z_cc_arr);
527  }
528  }
529 
530  // *****************************************************************************
531  // 9. Add CANOPY source terms
532  // *****************************************************************************
533  if (solverChoice.do_forest_drag &&
534  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
535  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
536  {
537  const Real ux = u(i, j, k);
538  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
539  + v(i, j+1, k ) + v(i-1, j+1, k ) );
540  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
541  + w(i, j , k+1) + w(i-1, j , k+1) );
542  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
543  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
544  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
545  });
546  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
547  {
548  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
549  + u(i+1, j , k ) + u(i+1, j-1, k ) );
550  const Real uy = v(i, j, k);
551  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
552  + w(i , j , k+1) + w(i , j-1, k+1) );
553  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
554  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
555  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
556  });
557  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
558  {
559  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
560  + u(i , j , k-1) + u(i+1, j , k-1) );
561  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
562  + v(i , j , k-1) + v(i , j+1, k-1) );
563  const amrex::Real uz = w(i, j, k);
564  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
565  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
566  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
567  });
568  }
569  // *****************************************************************************
570  // 10. Add Immersed source terms
571  // *****************************************************************************
572  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
573  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
574  const Real drag_coefficient=10.0/dz;
575  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
576  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
577  {
578  const Real ux = u(i, j, k);
579  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
580  + v(i, j+1, k ) + v(i-1, j+1, k ) );
581  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
582  + w(i, j , k+1) + w(i-1, j , k+1) );
583  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
584  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
585  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
586  xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed;
587  });
588  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
589  {
590  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
591  + u(i+1, j , k ) + u(i+1, j-1, k ) );
592  const Real uy = v(i, j, k);
593  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
594  + w(i , j , k+1) + w(i , j-1, k+1) );
595  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
596  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
597  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
598  ymom_src_arr(i, j, k) -= t_blank * CdM * uy * windspeed;
599  });
600  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
601  {
602  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
603  + u(i , j , k-1) + u(i+1, j , k-1) );
604  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
605  + v(i , j , k-1) + v(i , j+1, k-1) );
606  const amrex::Real uz = w(i, j, k);
607  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
608  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
609  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
610  zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed;
611  });
612  }
613  } // mfi
614 }
ApplySpongeZoneBCsForMom
void ApplySpongeZoneBCsForMom(const SpongeChoice &spongeChoice, const Geometry geom, const Box &tbx, const Box &tby, const Box &tbz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Array4< const Real > &rho_w, const Array4< const Real > &r0, const Array4< const Real > &z_phys_nd, const Array4< const Real > &z_phys_cc)
Definition: ERF_ApplySpongeZoneBCs.cpp:118
ApplySpongeZoneBCsForMom_ReadFromFile
void ApplySpongeZoneBCsForMom_ReadFromFile(const SpongeChoice &spongeChoice, const Geometry geom, const Box &tbx, const Box &tby, const Array4< const Real > &cell_data, const Array4< const Real > &z_phys_cc, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< const Real > &rho_u, const Array4< const Real > &rho_v, const Vector< Real * > d_sponge_ptrs_at_lev)
Definition: ERF_ApplySpongeZoneBCs_ReadFromFile.cpp:8
PIoTwo
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:7
ubar
@ ubar
Definition: ERF_DataStruct.H:79
wbar
@ wbar
Definition: ERF_DataStruct.H:79
vbar
@ vbar
Definition: ERF_DataStruct.H:79
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
NumericalDiffusion_Ymom
void NumericalDiffusion_Ymom(const Box &bx, const Real dt, const Real num_diff_coeff, const Array4< const Real > &prim_data, const Array4< const Real > &cell_data, const Array4< Real > &rhs, const Array4< const Real > &mfx_arr, const Array4< const Real > &mfy_arr)
Definition: ERF_NumericalDiffusion.cpp:151
NumericalDiffusion_Xmom
void NumericalDiffusion_Xmom(const Box &bx, const Real dt, const Real num_diff_coeff, const Array4< const Real > &prim_data, const Array4< const Real > &cell_data, const Array4< Real > &rhs, const Array4< const Real > &mfx_arr, const Array4< const Real > &mfy_arr)
Definition: ERF_NumericalDiffusion.cpp:85
offset
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
PlaneAverage
Definition: ERF_PlaneAverage.H:14
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:63
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:160
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:158
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:161
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:159
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:45
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
InputSoundingData::V_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Definition: ERF_InputSoundingData.H:326
InputSoundingData::input_sounding_time
amrex::Vector< amrex::Real > input_sounding_time
Definition: ERF_InputSoundingData.H:315
InputSoundingData::tau_nudging
amrex::Real tau_nudging
Definition: ERF_InputSoundingData.H:312
InputSoundingData::U_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Definition: ERF_InputSoundingData.H:326
SolverChoice::coriolis_factor
amrex::Real coriolis_factor
Definition: ERF_DataStruct.H:726
SolverChoice::rayleigh_damp_V
bool rayleigh_damp_V
Definition: ERF_DataStruct.H:698
SolverChoice::rayleigh_damp_substep
bool rayleigh_damp_substep
Definition: ERF_DataStruct.H:704
SolverChoice::rayleigh_zdamp
amrex::Real rayleigh_zdamp
Definition: ERF_DataStruct.H:702
SolverChoice::cosphi
amrex::Real cosphi
Definition: ERF_DataStruct.H:727
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:768
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Definition: ERF_DataStruct.H:733
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:739
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::immersed_forcing_substep
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:707
SolverChoice::sinphi
amrex::Real sinphi
Definition: ERF_DataStruct.H:728
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:770
SolverChoice::use_coriolis
bool use_coriolis
Definition: ERF_DataStruct.H:694
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:752
SolverChoice::variable_coriolis
bool variable_coriolis
Definition: ERF_DataStruct.H:773
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:735
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::has_lat_lon
bool has_lat_lon
Definition: ERF_DataStruct.H:772
SolverChoice::do_forest_drag
bool do_forest_drag
Definition: ERF_DataStruct.H:796
SolverChoice::forest_substep
bool forest_substep
Definition: ERF_DataStruct.H:708
SolverChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Definition: ERF_DataStruct.H:701
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:775
SpongeChoice::sponge_type
std::string sponge_type
Definition: ERF_SpongeStruct.H:58
Here is the call graph for this function: