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 (Real time, Real, const Vector< MultiFab > &S_data, const MultiFab *z_phys_nd, const MultiFab *z_phys_cc, Vector< Real > &stretched_dz_h, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &wvel, MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, const MultiFab &base_state, MultiFab *forest_drag, MultiFab *terrain_blank, MultiFab *cosPhi_mf, MultiFab *sinPhi_mf, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &, const Real *dptr_u_geos, const Real *dptr_v_geos, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const amrex::Real *d_sinesq_at_lev, const amrex::Real *d_sinesq_stag_at_lev, const Vector< Real * > d_sponge_ptrs_at_lev, const Vector< MultiFab > *forecast_state_at_lev, InputSoundingData &input_sounding_data, bool is_slow_step)
 

Function Documentation

◆ make_mom_sources()

void make_mom_sources ( Real  time,
Real  ,
const Vector< MultiFab > &  S_data,
const MultiFab *  z_phys_nd,
const MultiFab *  z_phys_cc,
Vector< Real > &  stretched_dz_h,
const MultiFab &  xvel,
const MultiFab &  yvel,
const MultiFab &  wvel,
MultiFab &  xmom_src,
MultiFab &  ymom_src,
MultiFab &  zmom_src,
const MultiFab &  base_state,
MultiFab *  forest_drag,
MultiFab *  terrain_blank,
MultiFab *  cosPhi_mf,
MultiFab *  sinPhi_mf,
const Geometry  geom,
const SolverChoice &  solverChoice,
Vector< std::unique_ptr< MultiFab >> &  ,
const Real *  dptr_u_geos,
const Real *  dptr_v_geos,
const Real *  dptr_wbar_sub,
const Vector< Real * >  d_rayleigh_ptrs_at_lev,
const amrex::Real *  d_sinesq_at_lev,
const amrex::Real *  d_sinesq_stag_at_lev,
const Vector< Real * >  d_sponge_ptrs_at_lev,
const Vector< MultiFab > *  forecast_state_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]timecurrent time
[in]dtcurrent slow or fast timestep size
[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
[in]d_sinesq_at_levsin( (pi/2) (z-z_t)/(damping depth)) at cell centers
[in]d_sinesq_stag_at_levsin( (pi/2) (z-z_t)/(damping depth)) at z-faces
67 {
68  BL_PROFILE_REGION("erf_make_mom_sources()");
69 
70  Box domain(geom.Domain());
71  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
72 
73  // Initialize sources to zero each time we may use them
74  xmom_src.setVal(0.0);
75  ymom_src.setVal(0.0);
76  zmom_src.setVal(0.0);
77 
78  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
79 
80  // flags to apply certain source terms in substep call only
81  bool use_Rayleigh_fast_uv = ( (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit) ||
82  (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastImplicit) );
83  bool use_Rayleigh_fast_w = (solverChoice.dampingChoice.rayleigh_damping_type == RayleighDampingType::FastExplicit);
84  bool use_canopy_fast = solverChoice.forest_substep;
85  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
86 
87  // *****************************************************************************
88  // Define source term for all three components of momenta from
89  // 1. Coriolis forcing for (xmom,ymom,zmom)
90  // 2. Rayleigh damping for (xmom,ymom,zmom)
91  // 3. Constant / height-dependent geostrophic forcing
92  // 4. Subsidence
93  // 5. Nudging towards input sounding data
94  // 6. Numerical diffusion for (xmom,ymom,zmom)
95  // 7. Sponge
96  // 8. Forest canopy
97  // 9a. Immersed forcing for terrain
98  // 9b. Immersed forcing for buildings
99  // 10. Constant mass flux
100  // *****************************************************************************
101  // NOTE: buoyancy is now computed in a separate routine - it should not appear here
102  // *****************************************************************************
103  //const bool l_use_ndiff = solverChoice.use_num_diff;
104 
105  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
106  if (solverChoice.do_forest_drag) {
107  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
108  }
109  }
110 
111 
112  // *****************************************************************************
113  // Data for Coriolis forcing
114  // *****************************************************************************
115  auto use_coriolis = solverChoice.use_coriolis;
116  auto coriolis_factor = solverChoice.coriolis_factor;
117  auto cosphi = solverChoice.cosphi;
118  auto sinphi = solverChoice.sinphi;
119  auto var_coriolis = solverChoice.variable_coriolis;
120  auto has_lat_lon = solverChoice.has_lat_lon;
121 
122  // *****************************************************************************
123  // Flag for Geostrophic forcing
124  // *****************************************************************************
125  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
126  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
127 
128  // *****************************************************************************
129  // Data for Rayleigh damping
130  // *****************************************************************************
131  auto rayleigh_damp_U = solverChoice.dampingChoice.rayleigh_damp_U;
132  auto rayleigh_damp_V = solverChoice.dampingChoice.rayleigh_damp_V;
133  auto rayleigh_damp_W = solverChoice.dampingChoice.rayleigh_damp_W;
134 
135  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
136  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
137  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
138 
139  // *****************************************************************************
140  // Data for constant mass flux
141  // *****************************************************************************
142  bool enforce_massflux_x = (solverChoice.const_massflux_u != 0);
143  bool enforce_massflux_y = (solverChoice.const_massflux_v != 0);
144  Real U_target = solverChoice.const_massflux_u;
145  Real V_target = solverChoice.const_massflux_v;
146  int massflux_klo = solverChoice.massflux_klo;
147  int massflux_khi = solverChoice.massflux_khi;
148 
149  // These will be updated by integrating through the planar average profiles
150  Real rhoUA_target{0};
151  Real rhoVA_target{0};
152  Real rhoUA{0};
153  Real rhoVA{0};
154 
155  // *****************************************************************************
156  // Planar averages for subsidence, nudging, or constant mass flux
157  // *****************************************************************************
158  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
159  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
160 
161  if (is_slow_step && (dptr_wbar_sub || solverChoice.nudging_from_input_sounding ||
162  enforce_massflux_x || enforce_massflux_y))
163  {
164  const int offset = 1;
165  const int u_offset = 1;
166  const int v_offset = 1;
167 
168  //
169  // We use the alias here to control ncomp inside the PlaneAverage
170  //
171  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
172 
173  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = offset;
174  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
175  r_ave.compute_averages(ZDir(), r_ave.field());
176 
177  int ncell = r_ave.ncell_line();
178  Gpu::HostVector< Real> r_plane_h(ncell);
179  Gpu::DeviceVector< Real> r_plane_d(ncell);
180 
181  r_ave.line_average(Rho_comp, r_plane_h);
182 
183  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
184 
185  Real* dptr_r = r_plane_d.data();
186 
187  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
188  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
189 
190  dptr_r_plane = r_plane_tab.table();
191  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
192  {
193  dptr_r_plane(k-offset) = dptr_r[k];
194  });
195 
196  // U and V momentum
197  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
198  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
199 
200  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
201  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
202 
203  u_ave.compute_averages(ZDir(), u_ave.field());
204  v_ave.compute_averages(ZDir(), v_ave.field());
205 
206  int u_ncell = u_ave.ncell_line();
207  int v_ncell = v_ave.ncell_line();
208  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
209  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
210 
211  u_ave.line_average(0, u_plane_h);
212  v_ave.line_average(0, v_plane_h);
213 
214  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
215  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
216 
217  Real* dptr_u = u_plane_d.data();
218  Real* dptr_v = v_plane_d.data();
219 
220  Box udomain = domain; udomain.grow(2,ng_u[2]);
221  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
222  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
223  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
224 
225  dptr_u_plane = u_plane_tab.table();
226  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
227  {
228  dptr_u_plane(k-u_offset) = dptr_u[k];
229  });
230 
231  dptr_v_plane = v_plane_tab.table();
232  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
233  {
234  dptr_v_plane(k-v_offset) = dptr_v[k];
235  });
236 
237  // sum in z for massflux adjustment
238  if (enforce_massflux_x || enforce_massflux_y) {
239  Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
240  Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
241 
242  if (solverChoice.mesh_type == MeshType::ConstantDz) {
243  // note: massflux_khi corresponds to unstaggered indices in this case
244  rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
245  u_plane_h.begin() + u_offset + massflux_khi+1, 0.0);
246  rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
247  v_plane_h.begin() + v_offset + massflux_khi+1, 0.0);
248  rhoUA_target = std::accumulate(r_plane_h.begin() + offset + massflux_klo,
249  r_plane_h.begin() + offset + massflux_khi+1, 0.0);
250  rhoVA_target = rhoUA_target;
251 
252  rhoUA *= geom.CellSize(2) * Ly;
253  rhoVA *= geom.CellSize(2) * Lx;
254  rhoUA_target *= geom.CellSize(2) * Ly;
255  rhoVA_target *= geom.CellSize(2) * Lx;
256 
257  } else if (solverChoice.mesh_type == MeshType::StretchedDz) {
258  // note: massflux_khi corresponds to staggered indices in this case
259  for (int k=massflux_klo; k < massflux_khi; ++k) {
260  rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
261  rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
262  rhoUA_target += r_plane_h[k + offset] * stretched_dz_h[k];
263  }
264  rhoVA_target = rhoUA_target;
265 
266  rhoUA *= Ly;
267  rhoVA *= Lx;
268  rhoUA_target *= Ly;
269  rhoVA_target *= Lx;
270  }
271 
272  // at this point, this is integrated rho*dA
273  rhoUA_target *= U_target;
274  rhoVA_target *= V_target;
275 
276  Print() << "Integrated mass flux : " << rhoUA << " " << rhoVA
277  << " (target: " << rhoUA_target << " " << rhoVA_target << ")"
278  << std::endl;
279  }
280  }
281 
282  // *****************************************************************************
283  // Add all the other forcings
284  // *****************************************************************************
285  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
286  {
287  Box tbx = mfi.nodaltilebox(0);
288  Box tby = mfi.nodaltilebox(1);
289  Box tbz = mfi.nodaltilebox(2);
290  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
291 
292  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
293  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
294  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
295  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
296 
297  const Array4<const Real>& u = xvel.array(mfi);
298  const Array4<const Real>& v = yvel.array(mfi);
299  const Array4<const Real>& w = wvel.array(mfi);
300 
301  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
302  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
303  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
304 
305  const Array4<const Real>& r0 = r_hse.const_array(mfi);
306 
307  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
308  Array4<const Real>{};
309  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
310  Array4<const Real>{};
311 
312  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
313  Array4<const Real>{};
314  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
315  Array4<const Real>{};
316 
317  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
318  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
319 
320 
321  // *****************************************************************************
322  // 1. Add CORIOLIS forcing (this assumes east is +x, north is +y)
323  // *****************************************************************************
324  if (use_coriolis && is_slow_step) {
325  if(solverChoice.init_type == InitType::HindCast) {
326  const Array4<const Real>& latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
327  ParallelFor(tbx, tby, tbz,
328  [=] AMREX_GPU_DEVICE (int i, int j, int k)
329  {
330  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));
331  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));
332  Real latitude = latlon_arr(i,j,k,0);
333  Real sphi_loc = std::sin(latitude*PI/180.0);
334  Real cphi_loc = std::cos(latitude*PI/180.0);
335  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
336  },
337  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
338  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));
339  Real latitude = latlon_arr(i,j,k,0);
340  Real sphi_loc = std::sin(latitude*PI/180.0);
341  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
342  },
343  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
344  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));
345  Real latitude = latlon_arr(i,j,k,0);
346  Real cphi_loc = std::cos(latitude*PI/180.0);
347  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_loc;
348  });
349  }
350  else if (var_coriolis && has_lat_lon) {
351  ParallelFor(tbx, tby, tbz,
352  [=] AMREX_GPU_DEVICE (int i, int j, int k)
353  {
354  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));
355  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));
356  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
357  Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
358  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
359  },
360  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
361  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));
362  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
363  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
364  },
365  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
366  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));
367  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
368  });
369  } else {
370  ParallelFor(tbx, tby, tbz,
371  [=] AMREX_GPU_DEVICE (int i, int j, int k)
372  {
373  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));
374  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));
375  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
376  },
377  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
378  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));
379  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
380  },
381  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
382  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));
383  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
384  });
385  } // var_coriolis
386  } // use_coriolis
387 
388  // *****************************************************************************
389  // 2. Add RAYLEIGH damping
390  // *****************************************************************************
391  Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef;
392 
393  if ( (is_slow_step && !use_Rayleigh_fast_uv) || (!is_slow_step && use_Rayleigh_fast_uv)) {
394  if (rayleigh_damp_U) {
395  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
396  {
397  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
398  Real uu = rho_u(i,j,k) / rho_on_u_face;
399  Real sinesq = d_sinesq_at_lev[k];
400  xmom_src_arr(i, j, k) -= dampcoef*sinesq * (uu - ubar[k]) * rho_on_u_face;
401  });
402  }
403 
404  if (rayleigh_damp_V) {
405  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
406  {
407  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
408  Real vv = rho_v(i,j,k) / rho_on_v_face;
409  Real sinesq = d_sinesq_at_lev[k];
410  ymom_src_arr(i, j, k) -= dampcoef*sinesq * (vv - vbar[k]) * rho_on_v_face;
411  });
412  }
413  } // fast or slow step
414 
415  if ( (is_slow_step && !use_Rayleigh_fast_w) || (!is_slow_step && use_Rayleigh_fast_w)) {
416  if (rayleigh_damp_W) {
417  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
418  {
419  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
420  Real ww = rho_w(i,j,k) / rho_on_w_face;
421  Real sinesq = d_sinesq_stag_at_lev[k];
422  zmom_src_arr(i, j, k) -= dampcoef*sinesq * (ww - wbar[k]) * rho_on_w_face;
423  });
424  }
425  } // fast or slow step
426 
427  // *****************************************************************************
428  // 3a. Add constant GEOSTROPHIC forcing
429  // *****************************************************************************
430  if (is_slow_step) {
431  ParallelFor(tbx, tby, tbz,
432  [=] AMREX_GPU_DEVICE (int i, int j, int k)
433  {
434  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
435  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
436  },
437  [=] AMREX_GPU_DEVICE (int i, int j, int k)
438  {
439  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
440  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
441  },
442  [=] AMREX_GPU_DEVICE (int i, int j, int k)
443  {
444  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
445  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
446  });
447  }
448 
449  // *****************************************************************************
450  // 3b. Add height-dependent GEOSTROPHIC forcing
451  // *****************************************************************************
452  if (geo_wind_profile && is_slow_step) {
453  ParallelFor(tbx, tby,
454  [=] AMREX_GPU_DEVICE (int i, int j, int k)
455  {
456  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
457  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
458  },
459  [=] AMREX_GPU_DEVICE (int i, int j, int k)
460  {
461  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
462  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
463  });
464  } // geo_wind_profile
465 
466  // *****************************************************************************
467  // 4. Add custom SUBSIDENCE terms
468  // *****************************************************************************
469  if (solverChoice.custom_w_subsidence && is_slow_step) {
470  if (solverChoice.custom_forcing_prim_vars) {
471  const int nr = Rho_comp;
472  ParallelFor(tbx, tby,
473  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
474  {
475  Real dzInv = 0.5*dxInv[2];
476  if (z_nd_arr) {
477  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
478  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
479  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
480  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
481  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
482  }
483  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
484  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
485  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
486  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
487  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
488  },
489  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
490  {
491  Real dzInv = 0.5*dxInv[2];
492  if (z_nd_arr) {
493  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
494  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
495  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
496  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
497  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
498  }
499  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
500  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
501  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
502  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
503  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
504  });
505  } else {
506  ParallelFor(tbx, tby,
507  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
508  {
509  Real dzInv = 0.5*dxInv[2];
510  if (z_nd_arr) {
511  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
512  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
513  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
514  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
515  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
516  }
517  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
518  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
519  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
520  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
521  },
522  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
523  {
524  Real dzInv = 0.5*dxInv[2];
525  if (z_nd_arr) {
526  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
527  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
528  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
529  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
530  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
531  }
532  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
533  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
534  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
535  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
536  });
537  }
538  }
539 
540  // *************************************************************************************
541  // 5. Add nudging towards value specified in input sounding
542  // *************************************************************************************
543  if (solverChoice.nudging_from_input_sounding && is_slow_step)
544  {
545  int itime_n = 0;
546  int itime_np1 = 0;
547  Real coeff_n = Real(1.0);
548  Real coeff_np1 = Real(0.0);
549 
550  Real tau_inv = Real(1.0) / input_sounding_data.tau_nudging;
551 
552  int n_sounding_times = input_sounding_data.input_sounding_time.size();
553 
554  for (int nt = 1; nt < n_sounding_times; nt++) {
555  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
556  }
557  if (itime_n == n_sounding_times-1) {
558  itime_np1 = itime_n;
559  } else {
560  itime_np1 = itime_n+1;
561  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
562  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
563  coeff_n = Real(1.0) - coeff_np1;
564  }
565 
566  int nr = Rho_comp;
567 
568  const Real* u_inp_sound_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr();
569  const Real* u_inp_sound_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr();
570  const Real* v_inp_sound_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr();
571  const Real* v_inp_sound_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr();
572  ParallelFor(tbx, tby,
573  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
574  {
575  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));
576  nudge_u *= tau_inv;
577  xmom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_u;
578  },
579  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
580  {
581  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));
582  nudge_v *= tau_inv;
583  ymom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_v;
584  });
585  }
586 
587  // *****************************************************************************
588  // 6. Add NUMERICAL DIFFUSION terms
589  // *****************************************************************************
590 #if 0
591  if (l_use_ndiff) {
592  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
593  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
594  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
595  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
596  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
597  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
598  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
599  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
600  }
601 #endif
602 
603  // *****************************************************************************
604  // 7. Add SPONGING
605  // *****************************************************************************
606  if (is_slow_step) {
607  if (solverChoice.spongeChoice.sponge_type == "input_sponge")
608  {
609  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
610  z_cc_arr, xmom_src_arr, ymom_src_arr,
611  rho_u, rho_v, d_sponge_ptrs_at_lev);
612  }
613  else
614  {
615  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
616  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
617  r0, z_nd_arr, z_cc_arr);
618  }
619 
620  if(solverChoice.init_type == InitType::HindCast and solverChoice.hindcast_lateral_forcing){
621 
622  const Array4<const Real>& rho_u_forecast_state = (*forecast_state_at_lev)[IntVars::xmom].array(mfi);
623  const Array4<const Real>& rho_v_forecast_state = (*forecast_state_at_lev)[IntVars::ymom].array(mfi);
624  const Array4<const Real>& rho_w_forecast_state = (*forecast_state_at_lev)[IntVars::zmom].array(mfi);
625  const Array4<const Real>& cons_forecast_state = (*forecast_state_at_lev)[IntVars::cons].array(mfi);
626  ApplyBndryForcing_Forecast(solverChoice, geom, tbx, tby, tbz, z_nd_arr,
627  xmom_src_arr, ymom_src_arr, zmom_src_arr,
628  rho_u, rho_v, rho_w,
629  rho_u_forecast_state, rho_v_forecast_state, rho_w_forecast_state,
630  cons_forecast_state);
631  }
632  }
633 
634  // *****************************************************************************
635  // 8. Add CANOPY source terms
636  // *****************************************************************************
637  if (solverChoice.do_forest_drag &&
638  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
639  ParallelFor(tbx, tby, tbz,
640  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
641  {
642  const Real ux = u(i, j, k);
643  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
644  + v(i, j+1, k ) + v(i-1, j+1, k ) );
645  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
646  + w(i, j , k+1) + w(i-1, j , k+1) );
647  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
648  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
649  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
650  },
651  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
652  {
653  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
654  + u(i+1, j , k ) + u(i+1, j-1, k ) );
655  const Real uy = v(i, j, k);
656  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
657  + w(i , j , k+1) + w(i , j-1, k+1) );
658  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
659  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
660  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
661  },
662  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
663  {
664  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
665  + u(i , j , k-1) + u(i+1, j , k-1) );
666  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
667  + v(i , j , k-1) + v(i , j+1, k-1) );
668  const amrex::Real uz = w(i, j, k);
669  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
670  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
671  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
672  });
673  }
674  // *****************************************************************************
675  // 9a. Add immersed source terms for terrain
676  // *****************************************************************************
677  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
678  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
679  // geometric properties
680  const Real* dx_arr = geom.CellSize();
681  const Real dx_x = dx_arr[0];
682  const Real dx_y = dx_arr[1];
683  const Real dx_z = dx_arr[2];
684 
685  const Real alpha_m = solverChoice.if_Cd_momentum;
686  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
687  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
688  const Real U_s = 1.0; // unit velocity scale
689 
690  // MOST parameters
691  similarity_funs sfuns;
692  const Real ggg = CONST_GRAV;
693  const Real kappa = KAPPA;
694  const Real z0 = solverChoice.if_z0;
695  const Real tflux_in = solverChoice.if_surf_temp_flux;
696  const Real Olen_in = solverChoice.if_Olen_in;
697  const bool l_use_most = solverChoice.if_use_most;
698 
699  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
700  {
701  const Real ux = u(i, j, k);
702  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
703  + v(i, j+1, k ) + v(i-1, j+1, k ) );
704  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
705  + w(i, j , k+1) + w(i-1, j , k+1) );
706  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
707  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
708  const Real t_blank_above = 0.5 * (t_blank_arr(i, j, k+1) + t_blank_arr(i-1, j, k+1));
709  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
710  const Real rho_xface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
711 
712  if ((t_blank > 0 && (t_blank_above == 0.0)) && l_use_most) { // force to MOST value
713  // calculate tangential velocity one cell above
714  const Real ux2r = u(i, j, k+1) ;
715  const Real uy2r = 0.25 * ( v(i, j , k+1) + v(i-1, j , k+1)
716  + v(i, j+1, k+1) + v(i-1, j+1, k+1) ) ;
717  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
718 
719  // MOST
720  const Real theta_xface = (0.5 * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i-1,j,k, RhoTheta_comp))) / rho_xface;
721  const Real rho_xface_below = 0.5 * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp) );
722  const Real theta_xface_below = (0.5 * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i-1,j,k-1, RhoTheta_comp))) / rho_xface_below;
723  const Real theta_surf = theta_xface_below;
724 
725  Real psi_m = 0.0;
726  Real psi_h = 0.0;
727  Real ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
728  Real tflux = (tflux_in != 1e-8) ? tflux_in : -(theta_xface - theta_surf) * ustar * kappa / (std::log(1.5 * dx_z / z0) - psi_h);
729  Real Olen = (Olen_in != 1e-8) ? Olen_in : -ustar * ustar * ustar * theta_xface / (kappa * ggg * tflux + tiny);
730  Real zeta = 1.5 * dx_z / Olen;
731 
732  // similarity functions
733  psi_m = sfuns.calc_psi_m(zeta);
734  psi_h = sfuns.calc_psi_h(zeta);
735  ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m);
736 
737  // prevent some unphysical math
738  if (!(ustar > 0.0 && !std::isnan(ustar))) { ustar = 0.0; }
739  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
740  if (psi_m > std::log(0.5 * dx_z / z0)) { psi_m = std::log(0.5 * dx_z / z0); }
741 
742  // determine target velocity
743  const Real uTarget = ustar / kappa * (std::log(0.5 * dx_z / z0) - psi_m);
744  Real uxTarget = uTarget * ux2r / (tiny + h_windspeed2r);
745  const Real bc_forcing_x = -(uxTarget - ux); // BC forcing pushes nonrelative velocity toward target velocity
746  xmom_src_arr(i, j, k) -= (1-t_blank) * rho_xface * CdM * U_s * bc_forcing_x; // if Vf low, force more strongly to MOST. If high, less forcing.
747  } else {
748  xmom_src_arr(i, j, k) -= t_blank * rho_xface * CdM * ux * windspeed;
749  }
750  });
751  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
752  {
753  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
754  + u(i+1, j , k ) + u(i+1, j-1, k ) );
755  const Real uy = v(i, j, k);
756  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
757  + w(i , j , k+1) + w(i , j-1, k+1) );
758  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
759  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
760  const Real t_blank_above = 0.5 * (t_blank_arr(i, j, k+1) + t_blank_arr(i, j-1, k+1));
761  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
762  const Real rho_yface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
763 
764  if ((t_blank > 0 && (t_blank_above == 0.0)) && l_use_most) { // force to MOST value
765  // calculate tangential velocity one cell above
766  const Real ux2r = 0.25 * ( u(i , j , k+1) + u(i , j-1, k+1)
767  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
768  const Real uy2r = v(i, j, k+1) ;
769  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
770 
771  // MOST
772  const Real theta_yface = (0.5 * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k, RhoTheta_comp))) / rho_yface;
773  const Real rho_yface_below = 0.5 * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
774  const Real theta_yface_below = (0.5 * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i,j-1,k-1, RhoTheta_comp))) / rho_yface_below;
775  const Real theta_surf = theta_yface_below;
776 
777  Real psi_m = 0.0;
778  Real psi_h = 0.0;
779  Real ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
780  Real tflux = (tflux_in != 1e-8) ? tflux_in : -(theta_yface - theta_surf) * ustar * kappa / (std::log(1.5 * dx_z / z0) - psi_h);
781  Real Olen = (Olen_in != 1e-8) ? Olen_in : -ustar * ustar * ustar * theta_yface / (kappa * ggg * tflux + tiny);
782  Real zeta = 1.5 * dx_z / Olen;
783 
784  // similarity functions
785  psi_m = sfuns.calc_psi_m(zeta);
786  psi_h = sfuns.calc_psi_h(zeta);
787  ustar = h_windspeed2r * kappa / (std::log(1.5 * dx_z / z0) - psi_m);
788 
789  // prevent some unphysical math
790  if (!(ustar > 0.0 && !std::isnan(ustar))) { ustar = 0.0; }
791  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
792  if (psi_m > std::log(0.5 * dx_z / z0)) { psi_m = std::log(0.5 * dx_z / z0); }
793 
794  // determine target velocity
795  const Real uTarget = ustar / kappa * (std::log(0.5 * dx_z / z0) - psi_m);
796  Real uyTarget = uTarget * uy2r / (tiny + h_windspeed2r);
797  const Real bc_forcing_y = -(uyTarget - uy); // BC forcing pushes nonrelative velocity toward target velocity
798  ymom_src_arr(i, j, k) -= (1 - t_blank) * rho_yface * CdM * U_s * bc_forcing_y; // if Vf low, force more strongly to MOST. If high, less forcing.
799  } else {
800  ymom_src_arr(i, j, k) -= t_blank * rho_yface * CdM * uy * windspeed;
801  }
802  });
803  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
804  {
805  const Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
806  + u(i , j , k-1) + u(i+1, j , k-1) );
807  const Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
808  + v(i , j , k-1) + v(i , j+1, k-1) );
809  const Real uz = w(i, j, k);
810  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
811  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
812  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
813  const Real rho_zface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
814  zmom_src_arr(i, j, k) -= t_blank * rho_zface * CdM * uz * windspeed;
815  });
816  }
817 
818  // *****************************************************************************
819  // 9b. Add immersed source terms for buildings
820  // *****************************************************************************
821  if ((solverChoice.buildings_type == BuildingsType::ImmersedForcing ) &&
822  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast)))
823  {
824  // geometric properties
825  const Real* dx_arr = geom.CellSize();
826  const Real dx_x = dx_arr[0];
827  const Real dx_y = dx_arr[1];
828 
829  const Real alpha_m = solverChoice.if_Cd_momentum;
830  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
831  const Real min_t_blank = 0.005; // threshold for where immersed forcing acts
832 
833  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
834  {
835  const Real ux = u(i, j, k);
836  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
837  + v(i, j+1, k ) + v(i-1, j+1, k ) );
838  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
839  + w(i, j , k+1) + w(i-1, j , k+1) );
840  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
841 
842  Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
843  if (t_blank < min_t_blank) { t_blank = 0.0; }
844  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
845  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
846  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
847  const Real rho_xface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
848  xmom_src_arr(i, j, k) -= t_blank * rho_xface * CdM * ux * windspeed;
849  });
850  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
851  {
852  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
853  + u(i+1, j , k ) + u(i+1, j-1, k ) );
854  const Real uy = v(i, j, k);
855  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
856  + w(i , j , k+1) + w(i , j-1, k+1) );
857  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
858 
859  Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
860  if (t_blank < min_t_blank) { t_blank = 0.0; }
861  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
862  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
863  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
864  const Real rho_yface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
865  ymom_src_arr(i, j, k) -= t_blank * rho_yface * CdM * uy * windspeed;
866  });
867  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
868  {
869  const Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
870  + u(i , j , k-1) + u(i+1, j , k-1) );
871  const Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
872  + v(i , j , k-1) + v(i , j+1, k-1) );
873  const Real uz = w(i, j, k);
874  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
875 
876  Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
877  if (t_blank < min_t_blank) { t_blank = 0.0; }
878  const Real dx_z = (z_nd_arr) ? (z_nd_arr(i,j,k) - z_nd_arr(i,j,k-1)) : dx_arr[2]; // ASW double check
879  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, 1./3.);
880  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
881  const Real rho_zface = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
882  zmom_src_arr(i, j, k) -= t_blank * rho_zface * CdM * uz * windspeed;
883  });
884  }
885 
886  // *****************************************************************************
887  // 10. Enforce constant mass flux
888  // *****************************************************************************
889  if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
890  Real tau_inv = Real(1.0) / solverChoice.const_massflux_tau;
891 
892  ParallelFor(tbx, tby,
893  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
894  xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
895  },
896  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
897  ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
898  });
899  }
900 
901  } // mfi
902 }
ApplyBndryForcing_Forecast
void ApplyBndryForcing_Forecast(const SolverChoice &solverChoice, const Geometry geom, const Box &tbx, const Box &tby, const Box &tbz, const Array4< const Real > &z_phys_nd, 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 > &rho_u_initial_state, const Array4< const Real > &rho_v_initial_state, const Array4< const Real > &rho_w_initial_state, const Array4< const Real > &cons_initial_state)
Definition: ERF_ApplyBndryForcing_Forecast.cpp:8
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
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
RayleighDampingType::FastExplicit
@ FastExplicit
RayleighDampingType::FastImplicit
@ FastImplicit
ubar
@ ubar
Definition: ERF_DataStruct.H:96
wbar
@ wbar
Definition: ERF_DataStruct.H:96
vbar
@ vbar
Definition: ERF_DataStruct.H:96
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
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
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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
DampingChoice::rayleigh_damp_V
bool rayleigh_damp_V
Definition: ERF_DampingStruct.H:85
DampingChoice::rayleigh_dampcoef
amrex::Real rayleigh_dampcoef
Definition: ERF_DampingStruct.H:88
DampingChoice::rayleigh_damp_W
bool rayleigh_damp_W
Definition: ERF_DampingStruct.H:86
DampingChoice::rayleigh_damping_type
RayleighDampingType rayleigh_damping_type
Definition: ERF_DampingStruct.H:101
DampingChoice::rayleigh_damp_U
bool rayleigh_damp_U
Definition: ERF_DampingStruct.H:84
InputSoundingData::V_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Definition: ERF_InputSoundingData.H:407
InputSoundingData::input_sounding_time
amrex::Vector< amrex::Real > input_sounding_time
Definition: ERF_InputSoundingData.H:394
InputSoundingData::tau_nudging
amrex::Real tau_nudging
Definition: ERF_InputSoundingData.H:391
InputSoundingData::U_inp_sound_d
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Definition: ERF_InputSoundingData.H:407
SolverChoice::coriolis_factor
amrex::Real coriolis_factor
Definition: ERF_DataStruct.H:1066
SolverChoice::mesh_type
static MeshType mesh_type
Definition: ERF_DataStruct.H:983
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Definition: ERF_DataStruct.H:1046
SolverChoice::if_use_most
bool if_use_most
Definition: ERF_DataStruct.H:1050
SolverChoice::dampingChoice
DampingChoice dampingChoice
Definition: ERF_DataStruct.H:993
SolverChoice::const_massflux_v
amrex::Real const_massflux_v
Definition: ERF_DataStruct.H:1136
SolverChoice::if_z0
amrex::Real if_z0
Definition: ERF_DataStruct.H:1045
SolverChoice::cosphi
amrex::Real cosphi
Definition: ERF_DataStruct.H:1067
SolverChoice::abl_geo_forcing
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:1109
SolverChoice::hindcast_lateral_forcing
bool hindcast_lateral_forcing
Definition: ERF_DataStruct.H:1145
SolverChoice::massflux_klo
int massflux_klo
Definition: ERF_DataStruct.H:1140
SolverChoice::custom_w_subsidence
bool custom_w_subsidence
Definition: ERF_DataStruct.H:1073
SolverChoice::nudging_from_input_sounding
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:1079
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:1038
SolverChoice::sinphi
amrex::Real sinphi
Definition: ERF_DataStruct.H:1068
SolverChoice::have_geo_wind_profile
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:1111
SolverChoice::const_massflux_u
amrex::Real const_massflux_u
Definition: ERF_DataStruct.H:1135
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Definition: ERF_DataStruct.H:1049
SolverChoice::use_coriolis
bool use_coriolis
Definition: ERF_DataStruct.H:1032
SolverChoice::num_diff_coeff
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:1095
SolverChoice::variable_coriolis
bool variable_coriolis
Definition: ERF_DataStruct.H:1114
SolverChoice::if_Cd_momentum
amrex::Real if_Cd_momentum
Definition: ERF_DataStruct.H:1043
SolverChoice::custom_forcing_prim_vars
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:1075
SolverChoice::buildings_type
static BuildingsType buildings_type
Definition: ERF_DataStruct.H:974
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:971
SolverChoice::spongeChoice
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:994
SolverChoice::init_type
static InitType init_type
Definition: ERF_DataStruct.H:965
SolverChoice::has_lat_lon
bool has_lat_lon
Definition: ERF_DataStruct.H:1113
SolverChoice::do_forest_drag
bool do_forest_drag
Definition: ERF_DataStruct.H:1132
SolverChoice::const_massflux_tau
amrex::Real const_massflux_tau
Definition: ERF_DataStruct.H:1137
SolverChoice::massflux_khi
int massflux_khi
Definition: ERF_DataStruct.H:1141
SolverChoice::forest_substep
bool forest_substep
Definition: ERF_DataStruct.H:1039
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:1116
SpongeChoice::sponge_type
std::string sponge_type
Definition: ERF_SpongeStruct.H:58
similarity_funs
Definition: ERF_MOSTStress.H:41
similarity_funs::calc_psi_m
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:91
similarity_funs::calc_psi_h
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:105
Here is the call graph for this function: