189 const int ncomp_mf = varnames.size();
193 if (ncomp_mf == 0)
return;
204 for (
int lev = 1; lev <= finest_level; ++lev) {
205 bool fillset =
false;
215 for (
int lev = 0; lev <= finest_level; ++lev) {
216 for (
int mvar(0); mvar<
qmoist[lev].size(); ++mvar) {
217 qmoist[lev][mvar] =
micro->Get_Qmoist_Ptr(lev,mvar);
222 Vector<MultiFab> mf(finest_level+1);
223 for (
int lev = 0; lev <= finest_level; ++lev) {
224 mf[lev].define(grids[lev], dmap[lev], ncomp_mf, 0);
228 Vector<MultiFab> mf_nd(finest_level+1);
230 for (
int lev = 0; lev <= finest_level; ++lev) {
231 BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes();
232 mf_nd[lev].define(nodal_grids, dmap[lev], 3, 0);
233 mf_nd[lev].setVal(0.);
238 Vector<MultiFab> mf_u(finest_level+1);
239 Vector<MultiFab> mf_v(finest_level+1);
240 Vector<MultiFab> mf_w(finest_level+1);
242 for (
int lev = 0; lev <= finest_level; ++lev) {
243 BoxArray grid_stag_u(grids[lev]); grid_stag_u.surroundingNodes(0);
244 BoxArray grid_stag_v(grids[lev]); grid_stag_v.surroundingNodes(1);
245 BoxArray grid_stag_w(grids[lev]); grid_stag_w.surroundingNodes(2);
246 mf_u[lev].define(grid_stag_u, dmap[lev], 1, 0);
247 mf_v[lev].define(grid_stag_v, dmap[lev], 1, 0);
248 mf_w[lev].define(grid_stag_w, dmap[lev], 1, 0);
256 Vector<MultiFab> mf_cc_vel(finest_level+1);
266 for (
int lev = 0; lev <= finest_level; ++lev) {
267 mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,1));
268 average_face_to_cellcenter(mf_cc_vel[lev],0,
280 amrex::Interpolater* mapper = &cell_cons_interp;
281 for (
int lev = 1; lev <= finest_level; ++lev)
283 Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
284 Vector<Real> ftime = {
t_new[lev],
t_new[lev]};
285 Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
286 Vector<Real> ctime = {
t_new[lev],
t_new[lev]};
288 amrex::FillPatchTwoLevels(mf_cc_vel[lev],
t_new[lev], cmf, ctime, fmf, ftime,
289 0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
298 for (
int lev = 0; lev <= finest_level; ++lev)
312 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 0, mf_comp, 1, 0);
316 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 1, mf_comp, 1, 0);
320 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 2, mf_comp, 1, 0);
326 auto calculate_derived = [&](
const std::string& der_name,
331 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
333 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
335 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
337 const Box& bx = mfi.tilebox();
338 auto& dfab = dmf[mfi];
339 auto& sfab = src_mf[mfi];
340 der_function(bx, dfab, 0, 1, sfab, Geom(lev),
t_new[0],
nullptr, lev);
366 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
367 Array<MultiFab const*, AMREX_SPACEDIM> u;
380 MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0);
385 MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0);
390 MultiFab::Copy(mf[lev],th_hse,0,mf_comp,1,0);
397 MultiFab::Copy(mf[lev], p_hse, 0, mf_comp, 1, 0);
400 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
403 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
405 const Box& bx = mfi.tilebox();
406 const Array4<Real >& derdat = mf[lev].array(mfi);
410 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
414 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta,qv_for_p);
424 MultiFab::Copy(mf[lev],
pp_inc[lev], 0, mf_comp, 1, 0);
427 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
430 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
432 const Box& bx = mfi.tilebox();
433 const Array4<Real>& derdat = mf[lev].array(mfi);
434 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
438 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
442 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k);
452 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
454 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
456 const Box& bx = mfi.tilebox();
457 const Array4<Real>& derdat = mf[lev].array(mfi);
459 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
460 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
461 derdat(i, j, k, mf_comp) = S_arr(i,j,k,
Rho_comp) - r0_arr(i,j,k);
470 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
472 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
474 const Box& bx = mfi.tilebox();
475 const Array4<Real>& derdat = mf[lev].array(mfi);
478 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
486 derdat(i, j, k, mf_comp) =
T*std::pow((pressure - pv)/
p_0, -
R_d/fac)*std::exp(
L_v*
qv/(fac*
T)) ;
494 MultiFab* terrain_blank =
m_terrain_drag[lev]->get_terrain_blank_field();
496 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
498 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
500 const Box& bx = mfi.tilebox();
501 const Array4<Real>& derdat = mf[lev].array(mfi);
502 const Array4<Real const>& src = terrain_blank->const_array(mfi);
503 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
504 derdat(i, j, k, mf_comp) = src(i,j,k);
510 #ifdef ERF_USE_WINDFARM
516 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
518 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
520 const Box& bx = mfi.tilebox();
521 const Array4<Real>& derdat = mf[lev].array(mfi);
522 const Array4<Real const>& Nturb_array = Nturb[lev].const_array(mfi);
523 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
524 derdat(i, j, k, mf_comp) = Nturb_array(i,j,k,0);
533 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
535 const Box& bx = mfi.tilebox();
536 const Array4<Real>& derdat = mf[lev].array(mfi);
537 const Array4<Real const>& SMark_array = SMark[lev].const_array(mfi);
538 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
539 derdat(i, j, k, mf_comp) = SMark_array(i,j,k,0);
547 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
549 const Box& bx = mfi.tilebox();
550 const Array4<Real>& derdat = mf[lev].array(mfi);
551 const Array4<Real const>& SMark_array = SMark[lev].const_array(mfi);
552 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
553 derdat(i, j, k, mf_comp) = SMark_array(i,j,k,1);
561 int klo = geom[lev].Domain().smallEnd(2);
562 int khi = geom[lev].Domain().bigEnd(2);
566 auto dxInv = geom[lev].InvCellSizeArray();
569 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
571 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
574 const Box& gbx = mfi.growntilebox(1);
575 const Array4<Real > & p_arr =
pres.array(mfi);
576 const Array4<Real const> & hse_arr =
base_state[lev].const_array(mfi);
579 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
580 p_arr(i,j,k) = hse_arr(i,j,k,1);
583 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
588 pres.FillBoundary(geom[lev].periodicity());
591 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
593 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
596 const Box& bx = mfi.tilebox();
597 const Array4<Real>& derdat = mf[lev].array(mfi);
598 const Array4<Real> & p_arr =
pres.array(mfi);
601 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
603 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
608 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
609 Real gp_zeta_on_iface_lo;
611 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
612 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
613 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
614 }
else if (k == khi) {
615 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
616 p_arr(i-1,j,k ) + p_arr(i,j,k )
617 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
619 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
620 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
621 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
623 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
628 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
629 Real gp_zeta_on_iface_hi;
631 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
632 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
633 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
634 }
else if (k == khi) {
635 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
636 p_arr(i+1,j,k ) + p_arr(i,j,k )
637 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
639 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
640 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
641 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
643 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
646 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
649 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
650 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i+1,j,k) - p_arr(i-1,j,k)) * dxInv[0];
659 auto dxInv = geom[lev].InvCellSizeArray();
663 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
665 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
668 const Box& gbx = mfi.growntilebox(1);
669 const Array4<Real > & p_arr =
pres.array(mfi);
670 const Array4<Real const> & hse_arr =
base_state[lev].const_array(mfi);
673 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
674 p_arr(i,j,k) = hse_arr(i,j,k,1);
677 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
682 pres.FillBoundary(geom[lev].periodicity());
685 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
687 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
690 const Box& bx = mfi.tilebox();
691 const Array4<Real>& derdat = mf[lev].array(mfi);
692 const Array4<Real> & p_arr =
pres.array(mfi);
695 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
697 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
701 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
702 Real gp_zeta_on_jface_lo;
704 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
705 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
706 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
707 }
else if (k == khi) {
708 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
709 p_arr(i,j,k ) + p_arr(i,j-1,k )
710 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
712 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
713 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
714 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
716 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
720 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
721 Real gp_zeta_on_jface_hi;
723 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
724 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
725 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
726 }
else if (k == khi) {
727 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
728 p_arr(i,j+1,k ) + p_arr(i,j,k )
729 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
731 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
732 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
733 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
735 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
737 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
740 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
741 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i,j+1,k) - p_arr(i,j-1,k)) * dxInv[1];
750 auto dxInv = geom[lev].InvCellSizeArray();
752 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
754 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
756 const Box& bx = mfi.tilebox();
757 const Array4<Real >& derdat = mf[lev].array(mfi);
758 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
761 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
763 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
766 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
767 Real gp_zeta_on_iface_lo;
769 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
770 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
771 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
772 }
else if (k == khi) {
773 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
774 p_arr(i-1,j,k ) + p_arr(i,j,k )
775 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
777 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
778 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
779 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
781 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
785 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
786 Real gp_zeta_on_iface_hi;
788 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
789 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
790 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
791 }
else if (k == khi) {
792 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
793 p_arr(i+1,j,k ) + p_arr(i,j,k )
794 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
796 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
797 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
798 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
800 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
802 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
810 auto dxInv = geom[lev].InvCellSizeArray();
812 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
814 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
816 const Box& bx = mfi.tilebox();
817 const Array4<Real >& derdat = mf[lev].array(mfi);
818 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
819 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
820 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
823 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
824 Real gp_zeta_on_jface_lo;
826 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
827 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
828 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
829 }
else if (k == khi) {
830 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
831 p_arr(i,j,k ) + p_arr(i,j-1,k )
832 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
834 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
835 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
836 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
838 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
842 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
843 Real gp_zeta_on_jface_hi;
845 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
846 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
847 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
848 }
else if (k == khi) {
849 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
850 p_arr(i,j+1,k ) + p_arr(i,j,k )
851 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
853 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
854 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
855 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
857 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
859 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
868 MultiFab::Copy(mf[lev],*
z_phys_cc[lev],0,mf_comp,1,0);
874 MultiFab::Copy(mf[lev],*
detJ_cc[lev],0,mf_comp,1,0);
881 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
883 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
885 const Box& bx = mfi.tilebox();
886 const Array4<Real>& derdat = mf[lev].array(mfi);
887 const Array4<Real>& mf_m =
mapfac_m[lev]->array(mfi);
888 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
889 derdat(i ,j ,k, mf_comp) = mf_m(i,j,0);
895 #ifdef ERF_USE_NETCDF
899 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
901 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
903 const Box& bx = mfi.tilebox();
904 const Array4<Real>& derdat = mf[lev].array(mfi);
905 const Array4<Real>& data = lat_m[lev]->array(mfi);
906 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
907 derdat(i, j, k, mf_comp) = data(i,j,0);
914 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
916 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
918 const Box& bx = mfi.tilebox();
919 const Array4<Real>& derdat = mf[lev].array(mfi);
920 const Array4<Real>& data = lon_m[lev]->array(mfi);
921 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
922 derdat(i, j, k, mf_comp) = data(i,j,0);
934 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
936 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
938 const Box& bx = mfi.tilebox();
939 const Array4<Real>& derdat = mf[lev].array(mfi);
940 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
942 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
944 derdat(i ,j ,k, mf_comp) = data(i,j,k,0) / norm;
952 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
954 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
956 const Box& bx = mfi.tilebox();
957 const Array4<Real>& derdat = mf[lev].array(mfi);
958 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
960 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
962 derdat(i ,j ,k, mf_comp) = data(i,j,k,1) / norm;
970 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
972 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
974 const Box& bx = mfi.tilebox();
975 const Array4<Real>& derdat = mf[lev].array(mfi);
976 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
978 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
980 derdat(i ,j ,k, mf_comp) = data(i,j,k,2) / norm;
988 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
990 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
992 const Box& bx = mfi.tilebox();
993 const Array4<Real>& derdat = mf[lev].array(mfi);
994 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
996 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
998 derdat(i ,j ,k, mf_comp) = data(i,j,k,3) / norm;
1006 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
1009 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1011 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
1013 const Box& bx = mfi.tilebox();
1014 auto prim = dmf[mfi].array();
1015 auto const cons = cmf[mfi].const_array();
1017 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1021 prim(i,j,k) = Kmv /
rho;
1049 MultiFab::Copy(mf[lev],*
walldist[lev],0,mf_comp,1,0);
1062 int n_qstate =
micro->Get_Qstate_Size();
1072 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
1105 int n_end = ncomp_cons - 1;
1107 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
1139 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1141 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1143 const Box& bx = mfi.tilebox();
1144 const Array4<Real>& derdat = mf[lev].array(mfi);
1146 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1151 erf_qsatw(
T, pressure, derdat(i,j,k,mf_comp));
1160 MultiFab::Copy(mf[lev],*(
qmoist[lev][4]),0,mf_comp,1,0);
1168 MultiFab::Copy(mf[lev],*(
qmoist[lev][8]),0,mf_comp,1,0);
1173 MultiFab::Copy(mf[lev],*(
qmoist[lev][9]),0,mf_comp,1,0);
1178 MultiFab::Copy(mf[lev],*(
qmoist[lev][10]),0,mf_comp,1,0);
1184 #ifdef ERF_USE_PARTICLES
1185 const auto& particles_namelist( particleData.getNames() );
1186 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
1188 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
1190 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
1191 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1196 Vector<std::string> particle_mesh_plot_names(0);
1197 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
1198 for (
int i = 0; i < particle_mesh_plot_names.size(); i++) {
1199 std::string plot_var_name(particle_mesh_plot_names[i]);
1201 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
1203 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
1204 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1212 MultiFab::Copy(mf[lev], EBFactory(lev).getVolFrac(), 0, mf_comp, 1, 0);
1217 #ifdef ERF_COMPUTE_ERROR
1226 Real H = geom[lev].ProbHi()[2];
1228 Real wavelength = 100.;
1229 Real kp = 2. *
PI / wavelength;
1231 Real
omega = std::sqrt(g * kp);
1234 const auto dx = geom[lev].CellSizeArray();
1237 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1239 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1241 const Box& bx = mfi.validbox();
1242 Box xbx(bx); xbx.surroundingNodes(0);
1246 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1248 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1251 Real
z = 0.25 * (z_nd(i,j,k) + z_nd(i,j+1,k) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1));
1253 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1256 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1258 xvel_arr(i,j,k) -= -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1261 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1263 Real
x = (i + 0.5) * dx[0];
1264 Real
z = 0.25 * ( z_nd(i,j,k) + z_nd(i+1,j,k) + z_nd(i,j+1,k) + z_nd(i+1,j+1,k));
1266 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1269 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1271 zvel_arr(i,j,k) -= Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1275 MultiFab temp_mf(mf[lev].boxArray(), mf[lev].DistributionMap(), AMREX_SPACEDIM, 0);
1276 average_face_to_cellcenter(temp_mf,0,
1280 MultiFab::Copy(mf[lev],temp_mf,0,mf_comp,1,0);
1284 MultiFab::Copy(mf[lev],temp_mf,1,mf_comp,1,0);
1288 MultiFab::Copy(mf[lev],temp_mf,2,mf_comp,1,0);
1294 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1296 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1298 const Box& bx = mfi.validbox();
1299 Box xbx(bx); xbx.surroundingNodes(0);
1304 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1306 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1309 Real
z = 0.25 * (z_nd(i,j,k) + z_nd(i,j+1,k) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1));
1310 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1314 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1315 xvel_arr(i,j,k) += -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1317 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1319 Real
x = (i + 0.5) * dx[0];
1320 Real
z = 0.25 * ( z_nd(i,j,k) + z_nd(i+1,j,k) + z_nd(i,j+1,k) + z_nd(i+1,j+1,k));
1321 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1324 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1326 zvel_arr(i,j,k) += Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1335 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1337 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1339 const Box& bx = mfi.tilebox();
1340 const Array4<Real>& derdat = mf[lev].array(mfi);
1341 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
1344 const auto dx = geom[lev].CellSizeArray();
1345 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1346 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
1348 Real H = geom[lev].ProbHi()[2];
1350 Real wavelength = 100.;
1351 Real kp = 2. *
PI / wavelength;
1353 Real
omega = std::sqrt(g * kp);
1356 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1359 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta) - p0_arr(i,j,k);
1361 Real rho_hse = r0_arr(i,j,k);
1363 Real
x = (i + 0.5) * dx[0];
1364 Real
z = 0.125 * ( z_nd(i,j,k ) + z_nd(i+1,j,k ) + z_nd(i,j+1,k ) + z_nd(i+1,j+1,k )
1365 +z_nd(i,j,k+1) + z_nd(i+1,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1) );
1366 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1369 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1370 Real pprime_exact = -(Ampl *
omega *
omega / kp) * fac *
1371 std::sin(kp *
x - omega_t) * r0_arr(i,j,k);
1373 derdat(i,j,k,mf_comp) -= pprime_exact;
1380 #ifdef ERF_USE_RRTMGP
1382 MultiFab::Copy(mf[lev], *(qheating_rates[lev]), 0, mf_comp, 1, 0);
1386 MultiFab::Copy(mf[lev], *(qheating_rates[lev]), 1, mf_comp, 1, 0);
1393 for (
int lev = 0; lev <= finest_level; ++lev) {
1394 EB_set_covered(mf[lev], 0.0);
1400 for (
int lev(0); lev <= finest_level; ++lev) {
1401 MultiFab::Copy(mf_nd[lev],*
z_phys_nd[lev],0,2,1,0);
1402 Real dz = Geom()[lev].CellSizeArray()[2];
1403 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
1404 const Box& bx = mfi.tilebox();
1405 Array4< Real> mf_arr = mf_nd[lev].array(mfi);
1406 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
1407 mf_arr(i,j,k,2) -= k * dz;
1413 std::string plotfilename;
1414 std::string plotfilenameU;
1415 std::string plotfilenameV;
1416 std::string plotfilenameW;
1422 }
else if (which == 2) {
1434 #ifdef ERF_USE_RRTMGP
1437 if (which==1 && plot_rad) {
1443 if (finest_level == 0)
1445 if (plotfile_type == PlotFileType::Amrex)
1447 Print() <<
"Writing native plotfile " << plotfilename <<
"\n";
1450 GetVecOfConstPtrs(mf),
1451 GetVecOfConstPtrs(mf_nd),
1455 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1456 GetVecOfConstPtrs(mf),
1463 Print() <<
"Writing face velocities" << std::endl;
1464 WriteMultiLevelPlotfile(plotfilenameU, finest_level+1,
1465 GetVecOfConstPtrs(mf_u),
1466 {
"x_velocity_stag"},
1468 WriteMultiLevelPlotfile(plotfilenameV, finest_level+1,
1469 GetVecOfConstPtrs(mf_v),
1470 {
"y_velocity_stag"},
1472 WriteMultiLevelPlotfile(plotfilenameW, finest_level+1,
1473 GetVecOfConstPtrs(mf_w),
1474 {
"z_velocity_stag"},
1478 #ifdef ERF_USE_PARTICLES
1479 particleData.writePlotFile(plotfilename);
1481 #ifdef ERF_USE_NETCDF
1482 }
else if (plotfile_type == PlotFileType::Netcdf) {
1485 writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames,
istep,
t_new[0]);
1489 Print() <<
"Writing no plotfile since plotfile_type is none" << std::endl;
1494 if (plotfile_type == PlotFileType::Amrex) {
1497 int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]);
1498 bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 1) || (ref_ratio[lev0][1] == 1) ) ||
1499 (ref_ratio[lev0][2] == 1) );
1500 for (
int lev = 1; lev < finest_level; lev++) {
1501 any_ratio_one = any_ratio_one ||
1502 ( ( (ref_ratio[lev][0] == 1) || (ref_ratio[lev][1] == 1) ) ||
1503 (ref_ratio[lev][2] == 1) );
1508 Vector<IntVect> r2(finest_level);
1509 Vector<Geometry> g2(finest_level+1);
1510 Vector<MultiFab> mf2(finest_level+1);
1512 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
1515 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
1518 Array<int,AMREX_SPACEDIM> periodicity =
1519 {Geom()[lev0].isPeriodic(0),Geom()[lev0].isPeriodic(1),Geom()[lev0].isPeriodic(2)};
1520 g2[lev0].define(Geom()[lev0].Domain(),&(Geom()[lev0].ProbDomain()),0,periodicity.data());
1522 r2[0] = IntVect(desired_ratio/ref_ratio[lev0][0],
1523 desired_ratio/ref_ratio[lev0][1],
1524 desired_ratio/ref_ratio[lev0][2]);
1526 for (
int lev = 1; lev <= finest_level; ++lev) {
1528 r2[lev-1][0] = r2[lev-2][0] * desired_ratio / ref_ratio[lev-1][0];
1529 r2[lev-1][1] = r2[lev-2][1] * desired_ratio / ref_ratio[lev-1][1];
1530 r2[lev-1][2] = r2[lev-2][2] * desired_ratio / ref_ratio[lev-1][2];
1533 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
1536 Box d2(Geom()[lev].Domain());
1537 d2.refine(r2[lev-1]);
1539 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
1543 for (
int lev = 1; lev <= finest_level; ++lev) {
1544 Interpolater* mapper_c = &pc_interp;
1545 InterpFromCoarseLevel(mf2[lev],
t_new[lev], mf[lev],
1553 Vector<IntVect> rr(finest_level);
1554 for (
int lev = 0; lev < finest_level; ++lev) {
1555 rr[lev] = IntVect(desired_ratio);
1558 Print() <<
"Writing plotfile " << plotfilename <<
"\n";
1561 GetVecOfConstPtrs(mf2),
1562 GetVecOfConstPtrs(mf_nd),
1566 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1567 GetVecOfConstPtrs(mf2), varnames,
1574 GetVecOfConstPtrs(mf),
1575 GetVecOfConstPtrs(mf_nd),
1579 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1580 GetVecOfConstPtrs(mf), varnames,
1584 Print() <<
"Writing face velocities" << std::endl;
1585 WriteMultiLevelPlotfile(plotfilenameU, finest_level+1,
1586 GetVecOfConstPtrs(mf_u),
1587 {
"x_velocity_stag"},
1589 WriteMultiLevelPlotfile(plotfilenameV, finest_level+1,
1590 GetVecOfConstPtrs(mf_v),
1591 {
"y_velocity_stag"},
1593 WriteMultiLevelPlotfile(plotfilenameW, finest_level+1,
1594 GetVecOfConstPtrs(mf_w),
1595 {
"z_velocity_stag"},
1602 #ifdef ERF_USE_PARTICLES
1603 particleData.writePlotFile(plotfilename);
1606 #ifdef ERF_USE_NETCDF
1607 }
else if (plotfile_type == PlotFileType::Netcdf) {
1608 for (
int lev = 0; lev <= finest_level; ++lev) {
1610 writeNCPlotFile(lev, which_box, plotfilename, GetVecOfConstPtrs(mf), varnames,
istep,
t_new[0]);
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
constexpr amrex::Real Cp_l
Definition: ERF_Constants.H:14
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
constexpr amrex::Real L_v
Definition: ERF_Constants.H:16
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:46
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:45
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esatw(amrex::Real t)
Definition: ERF_MicrophysicsUtils.H:64
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:158
PhysBCFunctNoOp null_bc_for_fill
Definition: ERF_Plotfile.cpp:9
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:101
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:87
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:130
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:158
void FillBdyCCVels(amrex::Vector< amrex::MultiFab > &mf_cc_vel)
Definition: ERF_FillBdyCCVels.cpp:11
static amrex::Vector< std::string > PlotFileVarNames(amrex::Vector< std::string > plot_var_names)
Definition: ERF_Plotfile.cpp:174
void WriteMultiLevelPlotfileWithTerrain(const std::string &plotfilename, int nlevels, const amrex::Vector< const amrex::MultiFab * > &mf, const amrex::Vector< const amrex::MultiFab * > &mf_nd, const amrex::Vector< std::string > &varnames, const amrex::Vector< amrex::Geometry > &my_geom, amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &my_ref_ratio, const std::string &versionName="HyperCLaw-V1.1", const std::string &levelPrefix="Level_", const std::string &mfPrefix="Cell", const amrex::Vector< std::string > &extra_dirs=amrex::Vector< std::string >()) const
Definition: ERF_Plotfile.cpp:1619
void writeJobInfo(const std::string &dir) const
Definition: ERF_WriteJobInfo.cpp:9
void Plot_Lsm_Data(amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &ref_ratio)
Definition: ERF_LandSurface.H:84
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:161
@ Mom_h
Definition: ERF_IndexDefines.H:151
@ Theta_h
Definition: ERF_IndexDefines.H:152
void erf_dervortx(const amrex::Box &bx, amrex::FArrayBox &derfab, int dcomp, int ncomp, const amrex::FArrayBox &datfab, const amrex::Geometry &geomdata, amrex::Real, const int *, const int)
Definition: ERF_Derive.cpp:200
void erf_dervorty(const amrex::Box &bx, amrex::FArrayBox &derfab, int dcomp, int ncomp, const amrex::FArrayBox &datfab, const amrex::Geometry &geomdata, amrex::Real, const int *, const int)
Definition: ERF_Derive.cpp:228
void erf_dermagvel(const amrex::Box &bx, amrex::FArrayBox &derfab, int dcomp, int ncomp, const amrex::FArrayBox &datfab, const amrex::Geometry &, amrex::Real, const int *, const int)
Definition: ERF_Derive.cpp:284
void erf_dernull(const Box &, FArrayBox &, int, int, const FArrayBox &, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:39
void erf_dervortz(const amrex::Box &bx, amrex::FArrayBox &derfab, int dcomp, int ncomp, const amrex::FArrayBox &datfab, const amrex::Geometry &geomdata, amrex::Real, const int *, const int)
Definition: ERF_Derive.cpp:256
void erf_dertemp(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:91
void erf_derKE(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:186
void erf_dermoisttemp(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:113
void erf_dersoundspeed(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: ERF_Derive.cpp:58