Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum.
60 BL_PROFILE_REGION(
"erf_make_mom_sources()");
62 Box domain(geom.Domain());
63 const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
94 if (solverChoice.
terrain_type == TerrainType::ImmersedForcing) {
96 amrex::Error(
" Currently forest canopy cannot be used with immersed forcing");
106 auto cosphi = solverChoice.
cosphi;
107 auto sinphi = solverChoice.
sinphi;
139 Real rhoUA_target{0};
140 Real rhoVA_target{0};
147 Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
148 TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;
151 enforce_massflux_x || enforce_massflux_y))
154 const int u_offset = 1;
155 const int v_offset = 1;
164 r_ave.compute_averages(
ZDir(), r_ave.field());
166 int ncell = r_ave.ncell_line();
167 Gpu::HostVector< Real> r_plane_h(ncell);
168 Gpu::DeviceVector< Real> r_plane_d(ncell);
170 r_ave.line_average(
Rho_comp, r_plane_h);
172 Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());
174 Real* dptr_r = r_plane_d.data();
176 Box tdomain = domain; tdomain.grow(2,ng_c[2]);
177 r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});
179 dptr_r_plane = r_plane_tab.table();
180 ParallelFor(ncell, [=] AMREX_GPU_DEVICE (
int k) noexcept
182 dptr_r_plane(k-
offset) = dptr_r[k];
186 IntVect ng_u = S_data[
IntVars::xmom].nGrowVect(); ng_u[2] = u_offset;
189 IntVect ng_v = S_data[
IntVars::ymom].nGrowVect(); ng_v[2] = v_offset;
192 u_ave.compute_averages(
ZDir(), u_ave.field());
193 v_ave.compute_averages(
ZDir(), v_ave.field());
195 int u_ncell = u_ave.ncell_line();
196 int v_ncell = v_ave.ncell_line();
197 Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
198 Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);
200 u_ave.line_average(0, u_plane_h);
201 v_ave.line_average(0, v_plane_h);
203 Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
204 Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());
206 Real* dptr_u = u_plane_d.data();
207 Real* dptr_v = v_plane_d.data();
209 Box udomain = domain; udomain.grow(2,ng_u[2]);
210 Box vdomain = domain; vdomain.grow(2,ng_v[2]);
211 u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
212 v_plane_tab.resize({vdomain.smallEnd(2)}, {vdomain.bigEnd(2)});
214 dptr_u_plane = u_plane_tab.table();
215 ParallelFor(u_ncell, [=] AMREX_GPU_DEVICE (
int k) noexcept
217 dptr_u_plane(k-u_offset) = dptr_u[k];
220 dptr_v_plane = v_plane_tab.table();
221 ParallelFor(v_ncell, [=] AMREX_GPU_DEVICE (
int k) noexcept
223 dptr_v_plane(k-v_offset) = dptr_v[k];
227 if (enforce_massflux_x || enforce_massflux_y) {
228 Real Lx = geom.ProbHi(0) - geom.ProbLo(0);
229 Real Ly = geom.ProbHi(1) - geom.ProbLo(1);
231 if (solverChoice.
mesh_type == MeshType::ConstantDz) {
233 rhoUA = std::accumulate(u_plane_h.begin() + u_offset + massflux_klo,
234 u_plane_h.begin() + u_offset + massflux_khi+1, 0.0);
235 rhoVA = std::accumulate(v_plane_h.begin() + v_offset + massflux_klo,
236 v_plane_h.begin() + v_offset + massflux_khi+1, 0.0);
237 rhoUA_target = std::accumulate(r_plane_h.begin() +
offset + massflux_klo,
238 r_plane_h.begin() +
offset + massflux_khi+1, 0.0);
239 rhoVA_target = rhoUA_target;
241 rhoUA *= geom.CellSize(2) * Ly;
242 rhoVA *= geom.CellSize(2) * Lx;
243 rhoUA_target *= geom.CellSize(2) * Ly;
244 rhoVA_target *= geom.CellSize(2) * Lx;
246 }
else if (solverChoice.
mesh_type == MeshType::StretchedDz) {
248 for (
int k=massflux_klo; k < massflux_khi; ++k) {
249 rhoUA += u_plane_h[k + u_offset] * stretched_dz_h[k];
250 rhoVA += v_plane_h[k + v_offset] * stretched_dz_h[k];
251 rhoUA_target += r_plane_h[k +
offset] * stretched_dz_h[k];
253 rhoVA_target = rhoUA_target;
262 rhoUA_target *= U_target;
263 rhoVA_target *= V_target;
265 Print() <<
"Integrated mass flux : " << rhoUA <<
" " << rhoVA
266 <<
" (target: " << rhoUA_target <<
" " << rhoVA_target <<
")"
274 for ( MFIter mfi(S_data[
IntVars::cons]); mfi.isValid(); ++mfi)
276 Box tbx = mfi.nodaltilebox(0);
277 Box tby = mfi.nodaltilebox(1);
278 Box tbz = mfi.nodaltilebox(2);
279 if (tbz.bigEnd(2) == domain.bigEnd(2)+1) tbz.growHi(2,-1);
281 const Array4<const Real>& cell_data = S_data[
IntVars::cons].array(mfi);
282 const Array4<const Real>& rho_u = S_data[
IntVars::xmom].array(mfi);
283 const Array4<const Real>& rho_v = S_data[
IntVars::ymom].array(mfi);
284 const Array4<const Real>& rho_w = S_data[
IntVars::zmom].array(mfi);
286 const Array4<const Real>& u =
xvel.array(mfi);
287 const Array4<const Real>& v =
yvel.array(mfi);
288 const Array4<const Real>& w = wvel.array(mfi);
290 const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
291 const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
292 const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
294 const Array4<const Real>& r0 = r_hse.const_array(mfi);
296 const Array4<const Real>& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) :
297 Array4<const Real>{};
298 const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
299 Array4<const Real>{};
301 const Array4<const Real>& cphi_arr = (cosPhi_mf) ? cosPhi_mf->const_array(mfi) :
302 Array4<const Real>{};
303 const Array4<const Real>& sphi_arr = (sinPhi_mf) ? sinPhi_mf->const_array(mfi) :
304 Array4<const Real>{};
306 const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
307 const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
312 if (use_coriolis && is_slow_step) {
313 if (var_coriolis && has_lat_lon) {
314 ParallelFor(tbx, tby, tbz,
315 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
317 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));
318 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));
319 Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i-1,j,0));
320 Real cphi_loc = 0.5 * (cphi_arr(i,j,0) + cphi_arr(i-1,j,0));
321 xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sphi_loc - rho_w_loc * cphi_loc);
323 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
324 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));
325 Real sphi_loc = 0.5 * (sphi_arr(i,j,0) + sphi_arr(i,j-1,0));
326 ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sphi_loc;
328 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
329 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));
330 zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cphi_arr(i,j,0);
333 ParallelFor(tbx, tby, tbz,
334 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
336 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));
337 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));
338 xmom_src_arr(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi);
340 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
341 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));
342 ymom_src_arr(i, j, k) += -coriolis_factor * rho_u_loc * sinphi;
344 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
345 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));
346 zmom_src_arr(i, j, k) += coriolis_factor * rho_u_loc * cosphi;
354 Real zlo = geom.ProbLo(2);
355 Real dz = geom.CellSize(2);
360 if ((is_slow_step && !use_Rayleigh_fast) || (!is_slow_step && use_Rayleigh_fast)) {
361 if (rayleigh_damp_U) {
362 ParallelFor(tbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
364 Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
365 Real zfrac = 1 - (ztop - zcc) / zdamp;
368 Real uu = rho_u(i,j,k) / rho_on_u_face;
370 xmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (uu -
ubar[k]) * rho_on_u_face;
375 if (rayleigh_damp_V) {
376 ParallelFor(tby, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
378 Real zcc = (z_cc_arr) ? z_cc_arr(i,j,k) : zlo + (k+0.5)*dz;
379 Real zfrac = 1 - (ztop - zcc) / zdamp;
382 Real vv = rho_v(i,j,k) / rho_on_v_face;
384 ymom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (vv -
vbar[k]) * rho_on_v_face;
389 if (rayleigh_damp_W) {
390 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
392 Real zstag = (z_nd_arr) ? z_nd_arr(i,j,k) : zlo + k*dz;
393 Real zfrac = 1 - (ztop - zstag) / zdamp;
396 Real ww = rho_w(i,j,k) / rho_on_w_face;
398 zmom_src_arr(i, j, k) -= dampcoef*sinefac*sinefac * (ww -
wbar[k]) * rho_on_w_face;
408 ParallelFor(tbx, tby, tbz,
409 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
412 xmom_src_arr(i, j, k) += rho_on_u_face * abl_geo_forcing[0];
414 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
417 ymom_src_arr(i, j, k) += rho_on_v_face * abl_geo_forcing[1];
419 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
422 zmom_src_arr(i, j, k) += rho_on_w_face * abl_geo_forcing[2];
429 if (geo_wind_profile && is_slow_step) {
430 ParallelFor(tbx, tby,
431 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
434 xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
436 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
439 ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
449 ParallelFor(tbx, tby,
450 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
452 Real dzInv = 0.5*dxInv[2];
454 Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
455 + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
456 Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
457 + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
458 dzInv = 1.0 / (z_xf_hi - z_xf_lo);
460 Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,
nr) + cell_data(i-1,j,k,
nr) );
461 Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
462 Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
463 Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
464 xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
466 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
468 Real dzInv = 0.5*dxInv[2];
470 Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
471 + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
472 Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
473 + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
474 dzInv = 1.0 / (z_yf_hi - z_yf_lo);
476 Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,
nr) + cell_data(i,j-1,k,
nr) );
477 Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
478 Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
479 Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
480 ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
483 ParallelFor(tbx, tby,
484 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
486 Real dzInv = 0.5*dxInv[2];
488 Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
489 + z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
490 Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
491 + z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
492 dzInv = 1.0 / (z_xf_hi - z_xf_lo);
494 Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
495 Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
496 Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
497 xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
499 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
501 Real dzInv = 0.5*dxInv[2];
503 Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
504 + z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
505 Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
506 + z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
507 dzInv = 1.0 / (z_yf_hi - z_yf_lo);
509 Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
510 Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
511 Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
512 ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
531 for (
int nt = 1; nt < n_sounding_times; nt++) {
534 if (itime_n == n_sounding_times-1) {
537 itime_np1 = itime_n+1;
540 coeff_n =
Real(1.0) - coeff_np1;
545 const Real* u_inp_sound_n = input_sounding_data.
U_inp_sound_d[itime_n].dataPtr();
546 const Real* u_inp_sound_np1 = input_sounding_data.
U_inp_sound_d[itime_np1].dataPtr();
547 const Real* v_inp_sound_n = input_sounding_data.
V_inp_sound_d[itime_n].dataPtr();
548 const Real* v_inp_sound_np1 = input_sounding_data.
V_inp_sound_d[itime_np1].dataPtr();
549 ParallelFor(tbx, tby,
550 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
552 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));
554 xmom_src_arr(i, j, k) += cell_data(i, j, k,
nr) * nudge_u;
556 [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
558 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));
560 ymom_src_arr(i, j, k) += cell_data(i, j, k,
nr) * nudge_v;
569 const Array4<const Real>& mf_ux = mapfac[MapFac::ux]->const_array(mfi);
570 const Array4<const Real>& mf_uy = mapfac[MapFac::uy]->const_array(mfi);
571 const Array4<const Real>& mf_vx = mapfac[MapFac::vx]->const_array(mfi);
572 const Array4<const Real>& mf_vy = mapfac[MapFac::vy]->const_array(mfi);
574 u, cell_data, xmom_src_arr, mf_ux, mf_uy);
576 v, cell_data, ymom_src_arr, mf_vx, mf_vy);
587 z_cc_arr, xmom_src_arr, ymom_src_arr,
588 rho_u, rho_v, d_sponge_ptrs_at_lev);
593 xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w,
594 r0, z_nd_arr, z_cc_arr);
602 ((is_slow_step && !use_canopy_fast) || (!is_slow_step && use_canopy_fast))) {
603 ParallelFor(tbx, tby, tbz,
604 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
606 const Real ux = u(i, j, k);
607 const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
608 + v(i, j+1, k ) + v(i-1, j+1, k ) );
609 const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
610 + w(i, j , k+1) + w(i-1, j , k+1) );
611 const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
612 const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i-1, j, k));
613 xmom_src_arr(i, j, k) -= f_drag * ux * windspeed;
615 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
617 const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
618 + u(i+1, j , k ) + u(i+1, j-1, k ) );
619 const Real uy = v(i, j, k);
620 const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
621 + w(i , j , k+1) + w(i , j-1, k+1) );
622 const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
623 const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j-1, k));
624 ymom_src_arr(i, j, k) -= f_drag * uy * windspeed;
626 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
628 const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k )
629 + u(i , j , k-1) + u(i+1, j , k-1) );
630 const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k )
631 + v(i , j , k-1) + v(i , j+1, k-1) );
633 const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
634 const Real f_drag = 0.5 * (f_drag_arr(i, j, k) + f_drag_arr(i, j, k-1));
635 zmom_src_arr(i, j, k) -= f_drag * uz * windspeed;
641 if (solverChoice.
terrain_type == TerrainType::ImmersedForcing &&
642 ((is_slow_step && !use_ImmersedForcing_fast) || (!is_slow_step && use_ImmersedForcing_fast))) {
643 const Real drag_coefficient=10.0/dz;
645 ParallelFor(tbx, tby, tbz,
646 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
648 const Real ux = u(i, j, k);
649 const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k )
650 + v(i, j+1, k ) + v(i-1, j+1, k ) );
651 const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k )
652 + w(i, j , k+1) + w(i-1, j , k+1) );
653 const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
654 const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
655 const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
656 xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed;
658 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
660 const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k )
661 + u(i+1, j , k ) + u(i+1, j-1, k ) );
662 const Real uy = v(i, j, k);
663 const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k )
664 + w(i , j , k+1) + w(i , j-1, k+1) );
665 const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
666 const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
667 const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
668 ymom_src_arr(i, j, k) -= t_blank * CdM * 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 t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
679 const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0);
680 zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed;
687 if (is_slow_step && (enforce_massflux_x || enforce_massflux_y)) {
690 ParallelFor(tbx, tby,
691 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
692 xmom_src_arr(i, j, k) += tau_inv * (rhoUA_target - rhoUA);
694 [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
695 ymom_src_arr(i, j, k) += tau_inv * (rhoVA_target - rhoVA);
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 PIoTwo
Definition: ERF_Constants.H:7
@ ubar
Definition: ERF_DataStruct.H:91
@ wbar
Definition: ERF_DataStruct.H:91
@ vbar
Definition: ERF_DataStruct.H:91
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:16
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 coriolis_factor
Definition: ERF_DataStruct.H:838
static MeshType mesh_type
Definition: ERF_DataStruct.H:777
bool rayleigh_damp_V
Definition: ERF_DataStruct.H:810
bool rayleigh_damp_substep
Definition: ERF_DataStruct.H:816
amrex::Real rayleigh_zdamp
Definition: ERF_DataStruct.H:814
amrex::Real const_massflux_v
Definition: ERF_DataStruct.H:908
amrex::Real cosphi
Definition: ERF_DataStruct.H:839
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:881
int massflux_klo
Definition: ERF_DataStruct.H:912
bool custom_w_subsidence
Definition: ERF_DataStruct.H:845
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:851
bool rayleigh_damp_U
Definition: ERF_DataStruct.H:809
amrex::Real rayleigh_ztop
Definition: ERF_DataStruct.H:815
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:819
amrex::Real sinphi
Definition: ERF_DataStruct.H:840
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:883
amrex::Real const_massflux_u
Definition: ERF_DataStruct.H:907
bool use_coriolis
Definition: ERF_DataStruct.H:806
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:867
bool variable_coriolis
Definition: ERF_DataStruct.H:886
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:847
static TerrainType terrain_type
Definition: ERF_DataStruct.H:768
bool rayleigh_damp_W
Definition: ERF_DataStruct.H:811
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:787
bool has_lat_lon
Definition: ERF_DataStruct.H:885
bool do_forest_drag
Definition: ERF_DataStruct.H:904
amrex::Real const_massflux_tau
Definition: ERF_DataStruct.H:909
int massflux_khi
Definition: ERF_DataStruct.H:913
bool forest_substep
Definition: ERF_DataStruct.H:820
amrex::Real rayleigh_dampcoef
Definition: ERF_DataStruct.H:813
int ave_plane
Definition: ERF_DataStruct.H:888
std::string sponge_type
Definition: ERF_SpongeStruct.H:58