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)) ;
495 MultiFab::Copy(mf[lev],*terrain_blank,0,mf_comp,1,0);
499 #ifdef ERF_USE_WINDFARM
504 MultiFab::Copy(mf[lev],Nturb[lev],0,mf_comp,1,0);
512 MultiFab::Copy(mf[lev],SMark[lev],0,mf_comp,1,0);
519 MultiFab::Copy(mf[lev],SMark[lev],1,mf_comp,1,0);
525 int klo = geom[lev].Domain().smallEnd(2);
526 int khi = geom[lev].Domain().bigEnd(2);
530 auto dxInv = geom[lev].InvCellSizeArray();
533 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
535 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
538 const Box& gbx = mfi.growntilebox(1);
539 const Array4<Real > & p_arr =
pres.array(mfi);
540 const Array4<Real const> & hse_arr =
base_state[lev].const_array(mfi);
543 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
544 p_arr(i,j,k) = hse_arr(i,j,k,1);
547 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
552 pres.FillBoundary(geom[lev].periodicity());
555 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
557 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
560 const Box& bx = mfi.tilebox();
561 const Array4<Real>& derdat = mf[lev].array(mfi);
562 const Array4<Real> & p_arr =
pres.array(mfi);
565 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
567 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
572 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
573 Real gp_zeta_on_iface_lo;
575 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
576 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
577 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
578 }
else if (k == khi) {
579 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
580 p_arr(i-1,j,k ) + p_arr(i,j,k )
581 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
583 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
584 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
585 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
587 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
592 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
593 Real gp_zeta_on_iface_hi;
595 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
596 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
597 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
598 }
else if (k == khi) {
599 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
600 p_arr(i+1,j,k ) + p_arr(i,j,k )
601 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
603 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
604 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
605 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
607 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
610 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
613 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
614 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i+1,j,k) - p_arr(i-1,j,k)) * dxInv[0];
623 auto dxInv = geom[lev].InvCellSizeArray();
627 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
629 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
632 const Box& gbx = mfi.growntilebox(1);
633 const Array4<Real > & p_arr =
pres.array(mfi);
634 const Array4<Real const> & hse_arr =
base_state[lev].const_array(mfi);
637 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
638 p_arr(i,j,k) = hse_arr(i,j,k,1);
641 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
646 pres.FillBoundary(geom[lev].periodicity());
649 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
651 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
654 const Box& bx = mfi.tilebox();
655 const Array4<Real>& derdat = mf[lev].array(mfi);
656 const Array4<Real> & p_arr =
pres.array(mfi);
659 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
661 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
665 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
666 Real gp_zeta_on_jface_lo;
668 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
669 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
670 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
671 }
else if (k == khi) {
672 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
673 p_arr(i,j,k ) + p_arr(i,j-1,k )
674 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
676 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
677 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
678 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
680 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
684 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
685 Real gp_zeta_on_jface_hi;
687 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
688 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
689 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
690 }
else if (k == khi) {
691 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
692 p_arr(i,j+1,k ) + p_arr(i,j,k )
693 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
695 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
696 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
697 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
699 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
701 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
704 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
705 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i,j+1,k) - p_arr(i,j-1,k)) * dxInv[1];
714 auto dxInv = geom[lev].InvCellSizeArray();
716 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
718 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
720 const Box& bx = mfi.tilebox();
721 const Array4<Real >& derdat = mf[lev].array(mfi);
722 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
725 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
727 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
730 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
731 Real gp_zeta_on_iface_lo;
733 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
734 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
735 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
736 }
else if (k == khi) {
737 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
738 p_arr(i-1,j,k ) + p_arr(i,j,k )
739 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
741 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
742 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
743 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
745 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
749 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
750 Real gp_zeta_on_iface_hi;
752 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
753 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
754 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
755 }
else if (k == khi) {
756 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
757 p_arr(i+1,j,k ) + p_arr(i,j,k )
758 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
760 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
761 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
762 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
764 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
766 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
774 auto dxInv = geom[lev].InvCellSizeArray();
776 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
778 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
780 const Box& bx = mfi.tilebox();
781 const Array4<Real >& derdat = mf[lev].array(mfi);
782 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
783 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
784 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
787 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
788 Real gp_zeta_on_jface_lo;
790 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
791 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
792 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
793 }
else if (k == khi) {
794 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
795 p_arr(i,j,k ) + p_arr(i,j-1,k )
796 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
798 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
799 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
800 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
802 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
806 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
807 Real gp_zeta_on_jface_hi;
809 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
810 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
811 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
812 }
else if (k == khi) {
813 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
814 p_arr(i,j+1,k ) + p_arr(i,j,k )
815 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
817 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
818 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
819 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
821 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
823 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
832 MultiFab::Copy(mf[lev],*
z_phys_cc[lev],0,mf_comp,1,0);
838 MultiFab::Copy(mf[lev],*
detJ_cc[lev],0,mf_comp,1,0);
845 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
847 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
849 const Box& bx = mfi.tilebox();
850 const Array4<Real>& derdat = mf[lev].array(mfi);
851 const Array4<Real>& mf_m =
mapfac_m[lev]->array(mfi);
852 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
853 derdat(i ,j ,k, mf_comp) = mf_m(i,j,0);
859 #ifdef ERF_USE_NETCDF
863 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
865 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
867 const Box& bx = mfi.tilebox();
868 const Array4<Real>& derdat = mf[lev].array(mfi);
869 const Array4<Real>& data = lat_m[lev]->array(mfi);
870 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
871 derdat(i, j, k, mf_comp) = data(i,j,0);
878 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
880 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
882 const Box& bx = mfi.tilebox();
883 const Array4<Real>& derdat = mf[lev].array(mfi);
884 const Array4<Real>& data = lon_m[lev]->array(mfi);
885 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
886 derdat(i, j, k, mf_comp) = data(i,j,0);
898 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
900 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
902 const Box& bx = mfi.tilebox();
903 const Array4<Real>& derdat = mf[lev].array(mfi);
904 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
906 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
908 derdat(i ,j ,k, mf_comp) = data(i,j,k,0) / norm;
916 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
918 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
920 const Box& bx = mfi.tilebox();
921 const Array4<Real>& derdat = mf[lev].array(mfi);
922 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
924 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
926 derdat(i ,j ,k, mf_comp) = data(i,j,k,1) / norm;
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,2) / 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,3) / norm;
970 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
973 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
975 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
977 const Box& bx = mfi.tilebox();
978 auto prim = dmf[mfi].array();
979 auto const cons = cmf[mfi].const_array();
981 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
985 prim(i,j,k) = Kmv /
rho;
1013 MultiFab::Copy(mf[lev],*
walldist[lev],0,mf_comp,1,0);
1017 MultiFab::Copy(mf[lev],*
SFS_diss_lev[lev],0,mf_comp,1,0);
1030 int n_qstate =
micro->Get_Qstate_Size();
1040 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
1073 int n_end = ncomp_cons - 1;
1075 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
1107 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1109 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1111 const Box& bx = mfi.tilebox();
1112 const Array4<Real>& derdat = mf[lev].array(mfi);
1114 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1119 erf_qsatw(
T, pressure, derdat(i,j,k,mf_comp));
1128 MultiFab::Copy(mf[lev],*(
qmoist[lev][4]),0,mf_comp,1,0);
1136 MultiFab::Copy(mf[lev],*(
qmoist[lev][8]),0,mf_comp,1,0);
1141 MultiFab::Copy(mf[lev],*(
qmoist[lev][9]),0,mf_comp,1,0);
1146 MultiFab::Copy(mf[lev],*(
qmoist[lev][10]),0,mf_comp,1,0);
1152 #ifdef ERF_USE_PARTICLES
1153 const auto& particles_namelist( particleData.getNames() );
1154 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
1156 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
1158 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
1159 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1164 Vector<std::string> particle_mesh_plot_names(0);
1165 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
1166 for (
int i = 0; i < particle_mesh_plot_names.size(); i++) {
1167 std::string plot_var_name(particle_mesh_plot_names[i]);
1169 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
1171 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
1172 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1182 MultiFab::Copy(mf[lev],
EBFactory(lev).getVolFrac(), 0, mf_comp, 1, 0);
1184 mf[lev].setVal(1.0, mf_comp, 1, 0);
1189 #ifdef ERF_COMPUTE_ERROR
1198 Real H = geom[lev].ProbHi()[2];
1200 Real wavelength = 100.;
1201 Real kp = 2. *
PI / wavelength;
1203 Real
omega = std::sqrt(g * kp);
1206 const auto dx = geom[lev].CellSizeArray();
1209 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1211 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1213 const Box& bx = mfi.validbox();
1214 Box xbx(bx); xbx.surroundingNodes(0);
1218 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1220 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1223 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));
1225 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1228 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1230 xvel_arr(i,j,k) -= -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1233 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1235 Real
x = (i + 0.5) * dx[0];
1236 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));
1238 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1241 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1243 zvel_arr(i,j,k) -= Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1247 MultiFab temp_mf(mf[lev].boxArray(), mf[lev].DistributionMap(), AMREX_SPACEDIM, 0);
1248 average_face_to_cellcenter(temp_mf,0,
1252 MultiFab::Copy(mf[lev],temp_mf,0,mf_comp,1,0);
1256 MultiFab::Copy(mf[lev],temp_mf,1,mf_comp,1,0);
1260 MultiFab::Copy(mf[lev],temp_mf,2,mf_comp,1,0);
1266 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1268 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1270 const Box& bx = mfi.validbox();
1271 Box xbx(bx); xbx.surroundingNodes(0);
1276 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1278 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1281 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));
1282 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1286 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1287 xvel_arr(i,j,k) += -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1289 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1291 Real
x = (i + 0.5) * dx[0];
1292 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));
1293 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1296 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1298 zvel_arr(i,j,k) += Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1307 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1309 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1311 const Box& bx = mfi.tilebox();
1312 const Array4<Real>& derdat = mf[lev].array(mfi);
1313 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
1316 const auto dx = geom[lev].CellSizeArray();
1317 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1318 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
1320 Real H = geom[lev].ProbHi()[2];
1322 Real wavelength = 100.;
1323 Real kp = 2. *
PI / wavelength;
1325 Real
omega = std::sqrt(g * kp);
1328 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1331 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta) - p0_arr(i,j,k);
1333 Real rho_hse = r0_arr(i,j,k);
1335 Real
x = (i + 0.5) * dx[0];
1336 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 )
1337 +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) );
1338 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1341 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1342 Real pprime_exact = -(Ampl *
omega *
omega / kp) * fac *
1343 std::sin(kp *
x - omega_t) * r0_arr(i,j,k);
1345 derdat(i,j,k,mf_comp) -= pprime_exact;
1352 #ifdef ERF_USE_RRTMGP
1354 MultiFab::Copy(mf[lev], *(qheating_rates[lev]), 0, mf_comp, 1, 0);
1358 MultiFab::Copy(mf[lev], *(qheating_rates[lev]), 1, mf_comp, 1, 0);
1366 for (
int lev = 0; lev <= finest_level; ++lev) {
1367 EB_set_covered(mf[lev], 0.0);
1373 for (
int lev(0); lev <= finest_level; ++lev) {
1374 MultiFab::Copy(mf_nd[lev],*
z_phys_nd[lev],0,2,1,0);
1375 Real dz = Geom()[lev].CellSizeArray()[2];
1376 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
1377 const Box& bx = mfi.tilebox();
1378 Array4<Real> mf_arr = mf_nd[lev].array(mfi);
1379 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
1380 mf_arr(i,j,k,2) -= k * dz;
1386 std::string plotfilename;
1387 std::string plotfilenameU;
1388 std::string plotfilenameV;
1389 std::string plotfilenameW;
1395 }
else if (which == 2) {
1407 #ifdef ERF_USE_RRTMGP
1418 if (finest_level == 0)
1420 if (plotfile_type == PlotFileType::Amrex)
1422 Print() <<
"Writing native plotfile " << plotfilename <<
"\n";
1425 GetVecOfConstPtrs(mf),
1426 GetVecOfConstPtrs(mf_nd),
1430 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1431 GetVecOfConstPtrs(mf),
1438 Print() <<
"Writing face velocities" << std::endl;
1439 WriteMultiLevelPlotfile(plotfilenameU, finest_level+1,
1440 GetVecOfConstPtrs(mf_u),
1441 {
"x_velocity_stag"},
1443 WriteMultiLevelPlotfile(plotfilenameV, finest_level+1,
1444 GetVecOfConstPtrs(mf_v),
1445 {
"y_velocity_stag"},
1447 WriteMultiLevelPlotfile(plotfilenameW, finest_level+1,
1448 GetVecOfConstPtrs(mf_w),
1449 {
"z_velocity_stag"},
1453 #ifdef ERF_USE_PARTICLES
1454 particleData.writePlotFile(plotfilename);
1456 #ifdef ERF_USE_NETCDF
1457 }
else if (plotfile_type == PlotFileType::Netcdf) {
1460 writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames,
istep,
t_new[0]);
1464 Print() <<
"Writing no plotfile since plotfile_type is none" << std::endl;
1469 if (plotfile_type == PlotFileType::Amrex) {
1472 int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]);
1473 bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 1) || (ref_ratio[lev0][1] == 1) ) ||
1474 (ref_ratio[lev0][2] == 1) );
1475 for (
int lev = 1; lev < finest_level; lev++) {
1476 any_ratio_one = any_ratio_one ||
1477 ( ( (ref_ratio[lev][0] == 1) || (ref_ratio[lev][1] == 1) ) ||
1478 (ref_ratio[lev][2] == 1) );
1483 Vector<IntVect> r2(finest_level);
1484 Vector<Geometry> g2(finest_level+1);
1485 Vector<MultiFab> mf2(finest_level+1);
1487 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
1490 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
1493 Array<int,AMREX_SPACEDIM> periodicity =
1494 {Geom()[lev0].isPeriodic(0),Geom()[lev0].isPeriodic(1),Geom()[lev0].isPeriodic(2)};
1495 g2[lev0].define(Geom()[lev0].Domain(),&(Geom()[lev0].ProbDomain()),0,periodicity.data());
1497 r2[0] = IntVect(desired_ratio/ref_ratio[lev0][0],
1498 desired_ratio/ref_ratio[lev0][1],
1499 desired_ratio/ref_ratio[lev0][2]);
1501 for (
int lev = 1; lev <= finest_level; ++lev) {
1503 r2[lev-1][0] = r2[lev-2][0] * desired_ratio / ref_ratio[lev-1][0];
1504 r2[lev-1][1] = r2[lev-2][1] * desired_ratio / ref_ratio[lev-1][1];
1505 r2[lev-1][2] = r2[lev-2][2] * desired_ratio / ref_ratio[lev-1][2];
1508 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
1511 Box d2(Geom()[lev].Domain());
1512 d2.refine(r2[lev-1]);
1514 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
1522 Vector<BCRec> temp_domain_bcs_type;
1523 temp_domain_bcs_type.resize(ncomp_mf);
1528 for (
int lev = 1; lev <= finest_level; ++lev) {
1529 Interpolater* mapper_c = &pc_interp;
1530 InterpFromCoarseLevel(mf2[lev],
t_new[lev], mf[lev],
1534 r2[lev-1], mapper_c, temp_domain_bcs_type, 0);
1538 Vector<IntVect> rr(finest_level);
1539 for (
int lev = 0; lev < finest_level; ++lev) {
1540 rr[lev] = IntVect(desired_ratio);
1543 Print() <<
"Writing plotfile " << plotfilename <<
"\n";
1546 GetVecOfConstPtrs(mf2),
1547 GetVecOfConstPtrs(mf_nd),
1551 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1552 GetVecOfConstPtrs(mf2), varnames,
1559 GetVecOfConstPtrs(mf),
1560 GetVecOfConstPtrs(mf_nd),
1564 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1565 GetVecOfConstPtrs(mf), varnames,
1569 Print() <<
"Writing face velocities" << std::endl;
1570 WriteMultiLevelPlotfile(plotfilenameU, finest_level+1,
1571 GetVecOfConstPtrs(mf_u),
1572 {
"x_velocity_stag"},
1574 WriteMultiLevelPlotfile(plotfilenameV, finest_level+1,
1575 GetVecOfConstPtrs(mf_v),
1576 {
"y_velocity_stag"},
1578 WriteMultiLevelPlotfile(plotfilenameW, finest_level+1,
1579 GetVecOfConstPtrs(mf_w),
1580 {
"z_velocity_stag"},
1587 #ifdef ERF_USE_PARTICLES
1588 particleData.writePlotFile(plotfilename);
1591 #ifdef ERF_USE_NETCDF
1592 }
else if (plotfile_type == PlotFileType::Netcdf) {
1593 for (
int lev = 0; lev <= finest_level; ++lev) {
1595 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:1604
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