Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum.
64 BL_PROFILE_REGION(
"erf_make_mom_sources()");
66 Box domain(geom.Domain());
67 const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
99 if (solverChoice.
terrain_type == TerrainType::ImmersedForcing) {
101 amrex::Error(
" Currently forest canopy cannot be used with immersed forcing");
111 auto cosphi = solverChoice.
cosphi;
112 auto sinphi = solverChoice.
sinphi;
144 Real rhoUA_target{0};
145 Real rhoVA_target{0};
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;
156 enforce_massflux_x || enforce_massflux_y))
159 const int u_offset = 1;
160 const int v_offset = 1;
169 r_ave.compute_averages(
ZDir(), r_ave.field());
171 int ncell = r_ave.ncell_line();
172 Gpu::HostVector< Real> r_plane_h(ncell);
173 Gpu::DeviceVector< Real> r_plane_d(ncell);
175 r_ave.line_average(
Rho_comp, r_plane_h);
177 Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
179 Real* dptr_r = r_plane_d.data();
181 Box tdomain = domain; tdomain.grow(2,ng_c[2]);
182 r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
184 dptr_r_plane = r_plane_tab.table();
185 ParallelFor(ncell, [=] AMREX_GPU_DEVICE (
int k) noexcept
187 dptr_r_plane(k-
offset) = dptr_r[k];
191 IntVect ng_u = S_data[
IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
194 IntVect ng_v = S_data[
IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
197 u_ave.compute_averages(
ZDir(), u_ave.field());
198 v_ave.compute_averages(
ZDir(), v_ave.field());
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);
205 u_ave.line_average(0, u_plane_h);
206 v_ave.line_average(0, v_plane_h);
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());
211 Real* dptr_u = u_plane_d.data();
212 Real* dptr_v = v_plane_d.data();
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)});
219 dptr_u_plane = u_plane_tab.table();
220 ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (
int k) noexcept
222 dptr_u_plane(k-u_offset) = dptr_u[k];
225 dptr_v_plane = v_plane_tab.table();
226 ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (
int k) noexcept
228 dptr_v_plane(k-v_offset) = dptr_v[k];
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);
236 if (solverChoice.
mesh_type == MeshType::ConstantDz) {
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;
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;
251 }
else if (solverChoice.
mesh_type == MeshType::StretchedDz) {
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];
258 rhoVA_target = rhoUA_target;
267 rhoUA_target *= U_target;
268 rhoVA_target *= V_target;
270 Print() <<
"Integrated mass flux : " << rhoUA <<
" " << rhoVA
271 <<
" (target: " << rhoUA_target <<
" " << rhoVA_target <<
")"
279 for ( MFIter mfi(S_data[
IntVars::cons]); mfi.isValid(); ++mfi)
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);
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);
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);
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);
299 const Array4<const Real>& r0 = r_hse.const_array(mfi);
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>{};
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>{};
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);
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)
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);
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;
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;
344 else if (var_coriolis && has_lat_lon) {
345 ParallelFor(tbx, tby, tbz,
346 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
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);
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;
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);
364 ParallelFor(tbx, tby, tbz,
365 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
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);
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;
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;
385 Real zlo = geom.ProbLo(2);
386 Real dz = geom.CellSize(2);
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)
395 Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
396 Real zfrac = 1 - (ztop - zcc) / zdamp;
399 Real uu = rho_u(i,j,k) / rho_on_u_face;
401 xmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (uu -
ubar[k]) * rho_on_u_face;
406 if (rayleigh_damp_V) {
407 ParallelFor(tby, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
409 Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
410 Real zfrac = 1 - (ztop - zcc) / zdamp;
413 Real vv = rho_v(i,j,k) / rho_on_v_face;
415 ymom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (vv -
vbar[k]) * rho_on_v_face;
420 if (rayleigh_damp_W) {
421 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
423 Real zstag = (z_nd_arr) ? z_nd_arr(i,j,k) : zlo + k*dz;
424 Real zfrac = 1 - (ztop - zstag) / zdamp;
427 Real ww = rho_w(i,j,k) / rho_on_w_face;
429 zmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (ww -
wbar[k]) * rho_on_w_face;
439 ParallelFor(tbx, tby, tbz,
440 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
443 xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
445 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
448 ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
450 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
453 zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
460 if (geo_wind_profile && is_slow_step) {
461 ParallelFor(tbx, tby,
462 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
465 xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
467 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
470 ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
480 ParallelFor(tbx, tby,
481 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
483 Real dzInv = 0.5*dxInv[2];
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);
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;
497 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
499 Real dzInv = 0.5*dxInv[2];
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);
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;
514 ParallelFor(tbx, tby,
515 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
517 Real dzInv = 0.5*dxInv[2];
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);
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;
530 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
532 Real dzInv = 0.5*dxInv[2];
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);
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;
562 for (
int nt = 1; nt < n_sounding_times; nt++) {
565 if (itime_n == n_sounding_times-1) {
568 itime_np1 = itime_n+1;
571 coeff_n =
Real(1.0) - coeff_np1;
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
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));
585 xmom_src_arr(i, j, k) += cell_data(i, j, k,
nr) * nudge_u;
587 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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));
591 ymom_src_arr(i, j, k) += cell_data(i, j, k,
nr) * nudge_v;
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);
605 u, cell_data, xmom_src_arr, mf_ux, mf_uy);
607 v, cell_data, ymom_src_arr, mf_vx, mf_vy);
618 z_cc_arr, xmom_src_arr, ymom_src_arr,
619 rho_u, rho_v, d_sponge_ptrs_at_lev);
624 xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
625 r0, z_nd_arr, z_cc_arr);
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);
635 xmom_src_arr, ymom_src_arr, zmom_src_arr,
637 rho_u_forecast_state, rho_v_forecast_state, rho_w_forecast_state,
638 cons_forecast_state);
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
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;
659 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
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;
670 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
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) );
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;
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
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;
702 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
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;
714 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
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) );
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;
731 if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
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);
738 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
739 ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
750 if (w_damping && is_slow_step) {
751 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
756 if (k == domain.smallEnd(2)) {
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) {
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) );
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) );
782 Real wmag = std::abs(rho_w(i,j,k) / rho_on_w_face);
784 Real cflw = wmag * dt * dzInv;
785 if (cflw > cflw_lim) {
787 AMREX_DEVICE_PRINTF(
"w-damping applied at (%d,%d,%d) for w-CFL = %f > %f\n",
788 i,j,k, cflw, cflw_lim);
790 printf(
"w-damping applied at (%d,%d,%d) for w-CFL = %f > %f\n",
791 i,j,k, cflw, cflw_lim);
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);
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::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