ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
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 dt, const Vector< MultiFab > &S_data, MultiFab &z_phys_nd, 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 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  dt,
const Vector< MultiFab > &  S_data,
MultiFab &  z_phys_nd,
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 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
63 {
64  BL_PROFILE_REGION("erf_make_mom_sources()");
65 
66  Box domain(geom.Domain());
67  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
68 
69  // Initialize sources to zero each time we may use them
70  xmom_src.setVal(0.0);
71  ymom_src.setVal(0.0);
72  zmom_src.setVal(0.0);
73 
74  MultiFab r_hse (base_state, make_alias, BaseState::r0_comp , 1);
75 
76  // flags to apply certain source terms in substep call only
77  bool use_Rayleigh_fast = solverChoice.rayleigh_damp_substep;
78  bool use_canopy_fast = solverChoice.forest_substep;
79  bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
80 
81  // *****************************************************************************
82  // Define source term for all three components of momenta from
83  // 1. Coriolis forcing for (xmom,ymom,zmom)
84  // 2. Rayleigh damping for (xmom,ymom,zmom)
85  // 3. Constant / height-dependent geostrophic forcing
86  // 4. Subsidence
87  // 5. Nudging towards input sounding data
88  // 6. Numerical diffusion for (xmom,ymom,zmom)
89  // 7. Sponge
90  // 8. Forest canopy
91  // 9. Immersed forcing
92  // 10. Constant mass flux
93  // 11. Vertical-velocity damping
94  // *****************************************************************************
95  // NOTE: buoyancy is now computed in a separate routine - it should not appear here
96  // *****************************************************************************
97  //const bool l_use_ndiff = solverChoice.use_num_diff;
98 
99  if (solverChoice.terrain_type == TerrainType::ImmersedForcing) {
100  if (solverChoice.do_forest_drag) {
101  amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
102  }
103  }
104 
105 
106  // *****************************************************************************
107  // Data for Coriolis forcing
108  // *****************************************************************************
109  auto use_coriolis = solverChoice.use_coriolis;
110  auto coriolis_factor = solverChoice.coriolis_factor;
111  auto cosphi = solverChoice.cosphi;
112  auto sinphi = solverChoice.sinphi;
113  auto var_coriolis = solverChoice.variable_coriolis;
114  auto has_lat_lon = solverChoice.has_lat_lon;
115 
116  // *****************************************************************************
117  // Flag for Geostrophic forcing
118  // *****************************************************************************
119  auto abl_geo_forcing = solverChoice.abl_geo_forcing;
120  auto geo_wind_profile = solverChoice.have_geo_wind_profile;
121 
122  // *****************************************************************************
123  // Data for Rayleigh damping
124  // *****************************************************************************
125  auto rayleigh_damp_U = solverChoice.rayleigh_damp_U;
126  auto rayleigh_damp_V = solverChoice.rayleigh_damp_V;
127  auto rayleigh_damp_W = solverChoice.rayleigh_damp_W;
128 
129  Real* ubar = d_rayleigh_ptrs_at_lev[Rayleigh::ubar];
130  Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar];
131  Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar];
132 
133  // *****************************************************************************
134  // Data for constant mass flux
135  // *****************************************************************************
136  bool enforce_massflux_x = (solverChoice.const_massflux_u != 0);
137  bool enforce_massflux_y = (solverChoice.const_massflux_v != 0);
138  Real U_target = solverChoice.const_massflux_u;
139  Real V_target = solverChoice.const_massflux_v;
140  int massflux_klo = solverChoice.massflux_klo;
141  int massflux_khi = solverChoice.massflux_khi;
142 
143  // These will be updated by integrating through the planar average profiles
144  Real rhoUA_target{0};
145  Real rhoVA_target{0};
146  Real rhoUA{0};
147  Real rhoVA{0};
148 
149  // *****************************************************************************
150  // Planar averages for subsidence, nudging, or constant mass flux
151  // *****************************************************************************
152  Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
153  TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
154 
155  if (is_slow_step && (dptr_wbar_sub || solverChoice.nudging_from_input_sounding ||
156  enforce_massflux_x || enforce_massflux_y))
157  {
158  const int offset = 1;
159  const int u_offset = 1;
160  const int v_offset = 1;
161 
162  //
163  // We use the alias here to control ncomp inside the PlaneAverage
164  //
165  MultiFab cons(S_data[IntVars::cons], make_alias, 0, 1);
166 
167  IntVect ng_c = S_data[IntVars::cons].nGrowVect(); ng_c[2] = offset;
168  PlaneAverage r_ave(&cons, geom, solverChoice.ave_plane, ng_c);
169  r_ave.compute_averages(ZDir(), r_ave.field());
170 
171  int ncell = r_ave.ncell_line();
172  Gpu::HostVector< Real> r_plane_h(ncell);
173  Gpu::DeviceVector< Real> r_plane_d(ncell);
174 
175  r_ave.line_average(Rho_comp, r_plane_h);
176 
177  Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
178 
179  Real* dptr_r = r_plane_d.data();
180 
181  Box tdomain = domain; tdomain.grow(2,ng_c[2]);
182  r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
183 
184  dptr_r_plane = r_plane_tab.table();
185  ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
186  {
187  dptr_r_plane(k-offset) = dptr_r[k];
188  });
189 
190  // U and V momentum
191  IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
192  PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, ng_u);
193 
194  IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
195  PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, ng_v);
196 
197  u_ave.compute_averages(ZDir(), u_ave.field());
198  v_ave.compute_averages(ZDir(), v_ave.field());
199 
200  int u_ncell = u_ave.ncell_line();
201  int v_ncell = v_ave.ncell_line();
202  Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
203  Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
204 
205  u_ave.line_average(0, u_plane_h);
206  v_ave.line_average(0, v_plane_h);
207 
208  Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
209  Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
210 
211  Real* dptr_u = u_plane_d.data();
212  Real* dptr_v = v_plane_d.data();
213 
214  Box udomain = domain; udomain.grow(2,ng_u[2]);
215  Box vdomain = domain; vdomain.grow(2,ng_v[2]);
216  u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
217  v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
218 
219  dptr_u_plane = u_plane_tab.table();
220  ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
221  {
222  dptr_u_plane(k-u_offset) = dptr_u[k];
223  });
224 
225  dptr_v_plane = v_plane_tab.table();
226  ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
227  {
228  dptr_v_plane(k-v_offset) = dptr_v[k];
229  });
230 
231  // sum in z for massflux adjustment
232  if (enforce_massflux_x || enforce_massflux_y) {
233  Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
234  Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
235 
236  if (solverChoice.mesh_type == MeshType::ConstantDz) {
237  // note: massflux_khi corresponds to unstaggered indices in this case
238  rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
239  u_plane_h.begin() + u_offset + massflux_khi+1, 0.0);
240  rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
241  v_plane_h.begin() + v_offset + massflux_khi+1, 0.0);
242  rhoUA_target = std::accumulate(r_plane_h.begin() + offset + massflux_klo,
243  r_plane_h.begin() + offset + massflux_khi+1, 0.0);
244  rhoVA_target = rhoUA_target;
245 
246  rhoUA *= geom.CellSize(2) * Ly;
247  rhoVA *= geom.CellSize(2) * Lx;
248  rhoUA_target *= geom.CellSize(2) * Ly;
249  rhoVA_target *= geom.CellSize(2) * Lx;
250 
251  } else if (solverChoice.mesh_type == MeshType::StretchedDz) {
252  // note: massflux_khi corresponds to staggered indices in this case
253  for (int k=massflux_klo; k < massflux_khi; ++k) {
254  rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
255  rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
256  rhoUA_target += r_plane_h[k + offset] * stretched_dz_h[k];
257  }
258  rhoVA_target = rhoUA_target;
259 
260  rhoUA *= Ly;
261  rhoVA *= Lx;
262  rhoUA_target *= Ly;
263  rhoVA_target *= Lx;
264  }
265 
266  // at this point, this is integrated rho*dA
267  rhoUA_target *= U_target;
268  rhoVA_target *= V_target;
269 
270  Print() << "Integrated mass flux : " << rhoUA << " " << rhoVA
271  << " (target: " << rhoUA_target << " " << rhoVA_target << ")"
272  << std::endl;
273  }
274  }
275 
276  // *****************************************************************************
277  // Add all the other forcings
278  // *****************************************************************************
279  for ( MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi)
280  {
281  Box tbx = mfi.nodaltilebox(0);
282  Box tby = mfi.nodaltilebox(1);
283  Box tbz = mfi.nodaltilebox(2);
284  if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
285 
286  const Array4<const Real>& cell_data = S_data[IntVars::cons].array(mfi);
287  const Array4<const Real>& rho_u = S_data[IntVars::xmom].array(mfi);
288  const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
289  const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);
290 
291  const Array4<const Real>& u = xvel.array(mfi);
292  const Array4<const Real>& v = yvel.array(mfi);
293  const Array4<const Real>& w = wvel.array(mfi);
294 
295  const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
296  const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
297  const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
298 
299  const Array4<const Real>& r0 = r_hse.const_array(mfi);
300 
301  const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
302  Array4<const Real>{};
303  const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
304  Array4<const Real>{};
305 
306  const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
307  Array4<const Real>{};
308  const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
309  Array4<const Real>{};
310 
311  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
312  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
313 
314 
315  // *****************************************************************************
316  // 1. Add CORIOLIS forcing (this assumes east is +x, north is +y)
317  // *****************************************************************************
318  if (use_coriolis && is_slow_step) {
319  if(solverChoice.init_type == InitType::HindCast) {
320  const Array4<const Real>& latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
321  ParallelFor(tbx, tby, tbz,
322  [=] AMREX_GPU_DEVICE (int i, int j, int k)
323  {
324  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));
325  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));
326  Real latitude = latlon_arr(i,j,k,0);
327  Real sphi_loc = std::sin(latitude*PI/180.0);
328  Real cphi_loc = std::cos(latitude*PI/180.0);
329  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
330  },
331  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
332  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));
333  Real latitude = latlon_arr(i,j,k,0);
334  Real sphi_loc = std::sin(latitude*PI/180.0);
335  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_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,k-1) + rho_u(i,j,k-1));
339  Real latitude = latlon_arr(i,j,k,0);
340  Real cphi_loc = std::cos(latitude*PI/180.0);
341  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_loc;
342  });
343  }
344  else if (var_coriolis && has_lat_lon) {
345  ParallelFor(tbx, tby, tbz,
346  [=] AMREX_GPU_DEVICE (int i, int j, int k)
347  {
348  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));
349  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));
350  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
351  Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
352  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
353  },
354  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
355  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));
356  Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
357  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
358  },
359  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
360  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));
361  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
362  });
363  } else {
364  ParallelFor(tbx, tby, tbz,
365  [=] AMREX_GPU_DEVICE (int i, int j, int k)
366  {
367  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));
368  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));
369  xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
370  },
371  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
372  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));
373  ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
374  },
375  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
376  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));
377  zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
378  });
379  } // var_coriolis
380  } // use_coriolis
381 
382  // *****************************************************************************
383  // 2. Add RAYLEIGH damping
384  // *****************************************************************************
385  Real zlo = geom.ProbLo(2);
386  Real dz = geom.CellSize(2);
387  Real ztop = solverChoice.rayleigh_ztop;
388  Real zdamp = solverChoice.rayleigh_zdamp;
389  Real dampcoef = solverChoice.rayleigh_dampcoef;
390 
391  if ((is_slow_step && !use_Rayleigh_fast) || (!is_slow_step && use_Rayleigh_fast)) {
392  if (rayleigh_damp_U) {
393  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
394  {
395  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
396  Real zfrac = 1 - (ztop - zcc) / zdamp;
397  if (zfrac > 0) {
398  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
399  Real uu = rho_u(i,j,k) / rho_on_u_face;
400  Real sinefac = std::sin(PIoTwo*zfrac);
401  xmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (uu - ubar[k]) * rho_on_u_face;
402  }
403  });
404  }
405 
406  if (rayleigh_damp_V) {
407  ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
408  {
409  Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
410  Real zfrac = 1 - (ztop - zcc) / zdamp;
411  if (zfrac > 0) {
412  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
413  Real vv = rho_v(i,j,k) / rho_on_v_face;
414  Real sinefac = std::sin(PIoTwo*zfrac);
415  ymom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (vv - vbar[k]) * rho_on_v_face;
416  }
417  });
418  }
419 
420  if (rayleigh_damp_W) {
421  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
422  {
423  Real zstag = (z_nd_arr) ? z_nd_arr(i,j,k) : zlo + k*dz;
424  Real zfrac = 1 - (ztop - zstag) / zdamp;
425  if (zfrac > 0) {
426  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
427  Real ww = rho_w(i,j,k) / rho_on_w_face;
428  Real sinefac = std::sin(PIoTwo*zfrac);
429  zmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (ww - wbar[k]) * rho_on_w_face;
430  }
431  });
432  }
433  } // fast or slow step
434 
435  // *****************************************************************************
436  // 3a. Add constant GEOSTROPHIC forcing
437  // *****************************************************************************
438  if (is_slow_step) {
439  ParallelFor(tbx, tby, tbz,
440  [=] AMREX_GPU_DEVICE (int i, int j, int k)
441  {
442  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
443  xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
444  },
445  [=] AMREX_GPU_DEVICE (int i, int j, int k)
446  {
447  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
448  ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
449  },
450  [=] AMREX_GPU_DEVICE (int i, int j, int k)
451  {
452  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
453  zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
454  });
455  }
456 
457  // *****************************************************************************
458  // 3b. Add height-dependent GEOSTROPHIC forcing
459  // *****************************************************************************
460  if (geo_wind_profile && is_slow_step) {
461  ParallelFor(tbx, tby,
462  [=] AMREX_GPU_DEVICE (int i, int j, int k)
463  {
464  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
465  xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
466  },
467  [=] AMREX_GPU_DEVICE (int i, int j, int k)
468  {
469  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
470  ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
471  });
472  } // geo_wind_profile
473 
474  // *****************************************************************************
475  // 4. Add custom SUBSIDENCE terms
476  // *****************************************************************************
477  if (solverChoice.custom_w_subsidence && is_slow_step) {
478  if (solverChoice.custom_forcing_prim_vars) {
479  const int nr = Rho_comp;
480  ParallelFor(tbx, tby,
481  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
482  {
483  Real dzInv = 0.5*dxInv[2];
484  if (z_nd_arr) {
485  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
486  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
487  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
488  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
489  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
490  }
491  Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
492  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
493  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
494  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
495  xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
496  },
497  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
498  {
499  Real dzInv = 0.5*dxInv[2];
500  if (z_nd_arr) {
501  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
502  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
503  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
504  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
505  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
506  }
507  Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
508  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
509  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
510  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
511  ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
512  });
513  } else {
514  ParallelFor(tbx, tby,
515  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
516  {
517  Real dzInv = 0.5*dxInv[2];
518  if (z_nd_arr) {
519  Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
520  + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
521  Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
522  + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
523  dzInv = 1.0 / (z_xf_hi - z_xf_lo);
524  }
525  Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
526  Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
527  Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
528  xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
529  },
530  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
531  {
532  Real dzInv = 0.5*dxInv[2];
533  if (z_nd_arr) {
534  Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
535  + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
536  Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
537  + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
538  dzInv = 1.0 / (z_yf_hi - z_yf_lo);
539  }
540  Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
541  Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
542  Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
543  ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
544  });
545  }
546  }
547 
548  // *************************************************************************************
549  // 5. Add nudging towards value specified in input sounding
550  // *************************************************************************************
551  if (solverChoice.nudging_from_input_sounding && is_slow_step)
552  {
553  int itime_n = 0;
554  int itime_np1 = 0;
555  Real coeff_n = Real(1.0);
556  Real coeff_np1 = Real(0.0);
557 
558  Real tau_inv = Real(1.0) / input_sounding_data.tau_nudging;
559 
560  int n_sounding_times = input_sounding_data.input_sounding_time.size();
561 
562  for (int nt = 1; nt < n_sounding_times; nt++) {
563  if (time > input_sounding_data.input_sounding_time[nt]) itime_n = nt;
564  }
565  if (itime_n == n_sounding_times-1) {
566  itime_np1 = itime_n;
567  } else {
568  itime_np1 = itime_n+1;
569  coeff_np1 = (time - input_sounding_data.input_sounding_time[itime_n]) /
570  (input_sounding_data.input_sounding_time[itime_np1] - input_sounding_data.input_sounding_time[itime_n]);
571  coeff_n = Real(1.0) - coeff_np1;
572  }
573 
574  int nr = Rho_comp;
575 
576  const Real* u_inp_sound_n = input_sounding_data.U_inp_sound_d[itime_n].dataPtr();
577  const Real* u_inp_sound_np1 = input_sounding_data.U_inp_sound_d[itime_np1].dataPtr();
578  const Real* v_inp_sound_n = input_sounding_data.V_inp_sound_d[itime_n].dataPtr();
579  const Real* v_inp_sound_np1 = input_sounding_data.V_inp_sound_d[itime_np1].dataPtr();
580  ParallelFor(tbx, tby,
581  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
582  {
583  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));
584  nudge_u *= tau_inv;
585  xmom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_u;
586  },
587  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
588  {
589  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));
590  nudge_v *= tau_inv;
591  ymom_src_arr(i, j, k) += cell_data(i, j, k, nr) * nudge_v;
592  });
593  }
594 
595  // *****************************************************************************
596  // 6. Add NUMERICAL DIFFUSION terms
597  // *****************************************************************************
598 #if 0
599  if (l_use_ndiff) {
600  const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
601  const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
602  const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
603  const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
604  NumericalDiffusion_Xmom(tbx, dt, solverChoice.num_diff_coeff,
605  u, cell_data, xmom_src_arr, mf_ux, mf_uy);
606  NumericalDiffusion_Ymom(tby, dt, solverChoice.num_diff_coeff,
607  v, cell_data, ymom_src_arr, mf_vx, mf_vy);
608  }
609 #endif
610 
611  // *****************************************************************************
612  // 7. Add SPONGING
613  // *****************************************************************************
614  if (is_slow_step) {
615  if (solverChoice.spongeChoice.sponge_type == "input_sponge")
616  {
617  ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
618  z_cc_arr, xmom_src_arr, ymom_src_arr,
619  rho_u, rho_v, d_sponge_ptrs_at_lev);
620  }
621  else
622  {
623  ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
624  xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
625  r0, z_nd_arr, z_cc_arr);
626  }
627 
628  if(solverChoice.init_type == InitType::HindCast and solverChoice.hindcast_lateral_forcing){
629 
630  const Array4<const Real>& rho_u_forecast_state = (*forecast_state_at_lev)[IntVars::xmom].array(mfi);
631  const Array4<const Real>& rho_v_forecast_state = (*forecast_state_at_lev)[IntVars::ymom].array(mfi);
632  const Array4<const Real>& rho_w_forecast_state = (*forecast_state_at_lev)[IntVars::zmom].array(mfi);
633  const Array4<const Real>& cons_forecast_state = (*forecast_state_at_lev)[IntVars::cons].array(mfi);
634  ApplyBndryForcing_Forecast(solverChoice, geom, tbx, tby, tbz, z_nd_arr,
635  xmom_src_arr, ymom_src_arr, zmom_src_arr,
636  rho_u, rho_v, rho_w,
637  rho_u_forecast_state, rho_v_forecast_state, rho_w_forecast_state,
638  cons_forecast_state);
639  }
640  }
641 
642  // *****************************************************************************
643  // 8. Add CANOPY source terms
644  // *****************************************************************************
645  if (solverChoice.do_forest_drag &&
646  ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
647  ParallelFor(tbx, tby, tbz,
648  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
649  {
650  const Real ux = u(i, j, k);
651  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
652  + v(i, j+1, k ) + v(i-1, j+1, k ) );
653  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
654  + w(i, j , k+1) + w(i-1, j , k+1) );
655  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
656  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
657  xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
658  },
659  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
660  {
661  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
662  + u(i+1, j , k ) + u(i+1, j-1, k ) );
663  const Real uy = v(i, j, k);
664  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
665  + w(i , j , k+1) + w(i , j-1, k+1) );
666  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
667  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
668  ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
669  },
670  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
671  {
672  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
673  + u(i , j , k-1) + u(i+1, j , k-1) );
674  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
675  + v(i , j , k-1) + v(i , j+1, k-1) );
676  const amrex::Real uz = w(i, j, k);
677  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
678  const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
679  zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
680  });
681  }
682  // *****************************************************************************
683  // 9. Add Immersed source terms
684  // *****************************************************************************
685  if (solverChoice.terrain_type == TerrainType::ImmersedForcing &&
686  ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
687  const Real drag_coefficient=10.0/dz;
689  ParallelFor(tbx, tby, tbz,
690  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
691  {
692  const Real ux = u(i, j, k);
693  const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
694  + v(i, j+1, k ) + v(i-1, j+1, k ) );
695  const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
696  + w(i, j , k+1) + w(i-1, j , k+1) );
697  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
698  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
699  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
700  xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed;
701  },
702  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
703  {
704  const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
705  + u(i+1, j , k ) + u(i+1, j-1, k ) );
706  const Real uy = v(i, j, k);
707  const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
708  + w(i , j , k+1) + w(i , j-1, k+1) );
709  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
710  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
711  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
712  ymom_src_arr(i, j, k) -= t_blank * CdM * uy * windspeed;
713  },
714  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
715  {
716  const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
717  + u(i , j , k-1) + u(i+1, j , k-1) );
718  const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
719  + v(i , j , k-1) + v(i , j+1, k-1) );
720  const amrex::Real uz = w(i, j, k);
721  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
722  const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
723  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
724  zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed;
725  });
726  }
727 
728  // *****************************************************************************
729  // 10. Enforce constant mass flux
730  // *****************************************************************************
731  if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
732  Real tau_inv = Real(1.0) / solverChoice.const_massflux_tau;
733 
734  ParallelFor(tbx, tby,
735  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
736  xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
737  },
738  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
739  ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
740  });
741  }
742 
743  // *****************************************************************************
744  // 11. Add w-damping
745  // *****************************************************************************
746  bool w_damping = solverChoice.w_damping;
747  amrex::Real w_damping_coeff = solverChoice.w_damping_coeff;
748  amrex::Real cflw_lim = solverChoice.w_damping_cfl;
749 
750  if (w_damping && is_slow_step) {
751  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
752  {
753  Real dzInv;
754  if (z_nd_arr) {
755  // average z_nd to face then diff
756  if (k == domain.smallEnd(2)) {
757  dzInv = 4.0 / (
758  z_nd_arr(i ,j ,k+1) - z_nd_arr(i ,j ,k)
759  + z_nd_arr(i+1,j ,k+1) - z_nd_arr(i+1,j ,k)
760  + z_nd_arr(i ,j+1,k+1) - z_nd_arr(i ,j+1,k)
761  + z_nd_arr(i+1,j+1,k+1) - z_nd_arr(i+1,j+1,k) );
762  } else if (k == domain.bigEnd(2)+1) {
763  dzInv = 4.0 / (
764  z_nd_arr(i ,j ,k) - z_nd_arr(i ,j ,k-1)
765  + z_nd_arr(i+1,j ,k) - z_nd_arr(i+1,j ,k-1)
766  + z_nd_arr(i ,j+1,k) - z_nd_arr(i ,j+1,k-1)
767  + z_nd_arr(i+1,j+1,k) - z_nd_arr(i+1,j+1,k-1) );
768  } else {
769  // dz = 0.5 * (dz[i,j,k] + dz[i,j,k+1])
770  // = 0.5 * (z_nd_face[i,j,k+1] - z_nd_face[i,j,k-1])
771  dzInv = 8.0 / (
772  z_nd_arr(i ,j ,k+1) - z_nd_arr(i ,j ,k-1)
773  + z_nd_arr(i+1,j ,k+1) - z_nd_arr(i+1,j ,k-1)
774  + z_nd_arr(i ,j+1,k+1) - z_nd_arr(i ,j+1,k-1)
775  + z_nd_arr(i+1,j+1,k+1) - z_nd_arr(i+1,j+1,k-1) );
776  }
777  } else {
778  dzInv = dxInv[2];
779  }
780 
781  Real rho_on_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
782  Real wmag = std::abs(rho_w(i,j,k) / rho_on_w_face);
783 
784  Real cflw = wmag * dt * dzInv;
785  if (cflw > cflw_lim) {
786 #ifdef AMREX_USE_GPU
787  AMREX_DEVICE_PRINTF("w-damping applied at (%d,%d,%d) for w-CFL = %f > %f\n",
788  i,j,k, cflw, cflw_lim);
789 #else
790  printf("w-damping applied at (%d,%d,%d) for w-CFL = %f > %f\n",
791  i,j,k, cflw, cflw_lim);
792 #endif
793  Real sgn_w = (rho_w(i,j,k) > 0) ? 1.0 : -1.0;
794  zmom_src_arr(i, j, k) -= rho_on_w_face * sgn_w * w_damping_coeff * (cflw - cflw_lim);
795  }
796  });
797  }
798  } // mfi
799 }
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
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
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
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:7
@ ubar
Definition: ERF_DataStruct.H:91
@ wbar
Definition: ERF_DataStruct.H:91
@ vbar
Definition: ERF_DataStruct.H:91
if(l_use_mynn &&start_comp<=RhoKE_comp &&end_comp >=RhoKE_comp)
Definition: ERF_DiffQKEAdjustment.H:2
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
#define Rho_comp
Definition: ERF_IndexDefines.H:36
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
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
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Definition: ERF_PlaneAverage.H:14
@ r0_comp
Definition: ERF_IndexDefines.H:63
@ ymom
Definition: ERF_IndexDefines.H:160
@ cons
Definition: ERF_IndexDefines.H:158
@ zmom
Definition: ERF_IndexDefines.H:161
@ xmom
Definition: ERF_IndexDefines.H:159
@ nr
Definition: ERF_Morrison.H:45
@ xvel
Definition: ERF_IndexDefines.H:141
@ cons
Definition: ERF_IndexDefines.H:140
@ yvel
Definition: ERF_IndexDefines.H:142
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Definition: ERF_InputSoundingData.H:410
amrex::Vector< amrex::Real > input_sounding_time
Definition: ERF_InputSoundingData.H:397
amrex::Real tau_nudging
Definition: ERF_InputSoundingData.H:394
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Definition: ERF_InputSoundingData.H:410
amrex::Real w_damping_cfl
Definition: ERF_DataStruct.H:900
amrex::Real coriolis_factor
Definition: ERF_DataStruct.H:923
static MeshType mesh_type
Definition: ERF_DataStruct.H:852
bool rayleigh_damp_V
Definition: ERF_DataStruct.H:891
bool rayleigh_damp_substep
Definition: ERF_DataStruct.H:897
amrex::Real rayleigh_zdamp
Definition: ERF_DataStruct.H:895
amrex::Real const_massflux_v
Definition: ERF_DataStruct.H:993
amrex::Real cosphi
Definition: ERF_DataStruct.H:924
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:966
bool hindcast_lateral_forcing
Definition: ERF_DataStruct.H:1001
int massflux_klo
Definition: ERF_DataStruct.H:997
bool custom_w_subsidence
Definition: ERF_DataStruct.H:930
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:936
bool rayleigh_damp_U
Definition: ERF_DataStruct.H:890
bool w_damping
Definition: ERF_DataStruct.H:899
amrex::Real rayleigh_ztop
Definition: ERF_DataStruct.H:896
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:904
amrex::Real sinphi
Definition: ERF_DataStruct.H:925
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:968
amrex::Real const_massflux_u
Definition: ERF_DataStruct.H:992
amrex::Real w_damping_coeff
Definition: ERF_DataStruct.H:901
bool use_coriolis
Definition: ERF_DataStruct.H:887
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:952
bool variable_coriolis
Definition: ERF_DataStruct.H:971
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:932
static TerrainType terrain_type
Definition: ERF_DataStruct.H:843
bool rayleigh_damp_W
Definition: ERF_DataStruct.H:892
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:862
static InitType init_type
Definition: ERF_DataStruct.H:837
bool has_lat_lon
Definition: ERF_DataStruct.H:970
bool do_forest_drag
Definition: ERF_DataStruct.H:989
amrex::Real const_massflux_tau
Definition: ERF_DataStruct.H:994
int massflux_khi
Definition: ERF_DataStruct.H:998
bool forest_substep
Definition: ERF_DataStruct.H:905
amrex::Real rayleigh_dampcoef
Definition: ERF_DataStruct.H:894
int ave_plane
Definition: ERF_DataStruct.H:973
std::string sponge_type
Definition: ERF_SpongeStruct.H:58
Here is the call graph for this function: