188 if (plotfile_type == PlotFileType::Amrex) {
189 amrex::Print() <<
"WRITE AMREX " << std::endl;
190 }
else if (plotfile_type == PlotFileType::Netcdf) {
191 amrex::Print() <<
"WRITE NETCDF " << std::endl;
195 const int ncomp_mf = varnames.size();
199 if (ncomp_mf == 0)
return;
205 for (
int lev = 0; lev <= finest_level; ++lev) {
206 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_cc_vel(finest_level+1);
248 for (
int lev = 0; lev <= finest_level; ++lev) {
249 mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,1));
250 average_face_to_cellcenter(mf_cc_vel[lev],0,
262 amrex::Interpolater* mapper = &cell_cons_interp;
263 for (
int lev = 1; lev <= finest_level; ++lev)
265 Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
266 Vector<Real> ftime = {
t_new[lev],
t_new[lev]};
267 Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
268 Vector<Real> ctime = {
t_new[lev],
t_new[lev]};
270 amrex::FillPatchTwoLevels(mf_cc_vel[lev],
t_new[lev], cmf, ctime, fmf, ftime,
271 0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
280 for (
int lev = 0; lev <= finest_level; ++lev)
294 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 0, mf_comp, 1, 0);
298 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 1, mf_comp, 1, 0);
302 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 2, mf_comp, 1, 0);
308 auto calculate_derived = [&](
const std::string& der_name,
313 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
315 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
317 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
319 const Box& bx = mfi.tilebox();
320 auto& dfab = dmf[mfi];
321 auto& sfab = src_mf[mfi];
322 der_function(bx, dfab, 0, 1, sfab, Geom(lev),
t_new[0],
nullptr, lev);
348 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
349 Array<MultiFab const*, AMREX_SPACEDIM> u;
353 computeDivergence(dmf, u, geom[lev]);
357 MultiFab r_hse(
base_state[lev], make_alias, 0, 1);
358 MultiFab p_hse(
base_state[lev], make_alias, 1, 1);
362 MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0);
368 MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0);
375 MultiFab::Copy(mf[lev], p_hse, 0, mf_comp, 1, 0);
378 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
381 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
383 const Box& bx = mfi.tilebox();
384 const Array4<Real >& derdat = mf[lev].array(mfi);
388 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
392 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta,qv_for_p);
402 MultiFab::Copy(mf[lev],
pp_inc[lev], 0, mf_comp, 1, 0);
405 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
408 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
410 const Box& bx = mfi.tilebox();
411 const Array4<Real>& derdat = mf[lev].array(mfi);
412 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
416 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
420 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k);
430 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
432 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
434 const Box& bx = mfi.tilebox();
435 const Array4<Real>& derdat = mf[lev].array(mfi);
437 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
438 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
439 derdat(i, j, k, mf_comp) = S_arr(i,j,k,
Rho_comp) - r0_arr(i,j,k);
448 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
450 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
452 const Box& bx = mfi.tilebox();
453 const Array4<Real>& derdat = mf[lev].array(mfi);
456 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
464 derdat(i, j, k, mf_comp) =
T*std::pow((pressure - pv)/
p_0, -
R_d/fac)*std::exp(
L_v*
qv/(fac*
T)) ;
470 #ifdef ERF_USE_WINDFARM
476 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
478 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
480 const Box& bx = mfi.tilebox();
481 const Array4<Real>& derdat = mf[lev].array(mfi);
482 const Array4<Real const>& Nturb_array = Nturb[lev].const_array(mfi);
483 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
484 derdat(i, j, k, mf_comp) = Nturb_array(i,j,k,0);
492 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
494 const Box& bx = mfi.tilebox();
495 const Array4<Real>& derdat = mf[lev].array(mfi);
496 const Array4<Real const>& SMark_array = SMark[lev].const_array(mfi);
497 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
498 derdat(i, j, k, mf_comp) = SMark_array(i,j,k,0);
506 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
508 const Box& bx = mfi.tilebox();
509 const Array4<Real>& derdat = mf[lev].array(mfi);
510 const Array4<Real const>& SMark_array = SMark[lev].const_array(mfi);
511 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
512 derdat(i, j, k, mf_comp) = SMark_array(i,j,k,1);
520 int klo = geom[lev].Domain().smallEnd(2);
521 int khi = geom[lev].Domain().bigEnd(2);
525 auto dxInv = geom[lev].InvCellSizeArray();
528 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
530 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
533 const Box& gbx = mfi.growntilebox(1);
534 const Array4<Real > & p_arr =
pres.array(mfi);
535 const Array4<Real const> & hse_arr =
base_state[lev].const_array(mfi);
538 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
539 p_arr(i,j,k) = hse_arr(i,j,k,1);
542 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
547 pres.FillBoundary(geom[lev].periodicity());
550 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
552 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
555 const Box& bx = mfi.tilebox();
556 const Array4<Real>& derdat = mf[lev].array(mfi);
557 const Array4<Real> & p_arr =
pres.array(mfi);
560 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
562 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
567 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
568 Real gp_zeta_on_iface_lo;
570 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
571 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
572 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
573 }
else if (k == khi) {
574 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
575 p_arr(i-1,j,k ) + p_arr(i,j,k )
576 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
578 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
579 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
580 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
582 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
587 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
588 Real gp_zeta_on_iface_hi;
590 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
591 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
592 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
593 }
else if (k == khi) {
594 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
595 p_arr(i+1,j,k ) + p_arr(i,j,k )
596 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
598 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
599 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
600 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
602 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
605 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
608 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
609 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i+1,j,k) - p_arr(i-1,j,k)) * dxInv[0];
618 auto dxInv = geom[lev].InvCellSizeArray();
622 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
624 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
627 const Box& gbx = mfi.growntilebox(1);
628 const Array4<Real > & p_arr =
pres.array(mfi);
629 const Array4<Real const> & hse_arr =
base_state[lev].const_array(mfi);
632 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
633 p_arr(i,j,k) = hse_arr(i,j,k,1);
636 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
641 pres.FillBoundary(geom[lev].periodicity());
644 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
646 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
649 const Box& bx = mfi.tilebox();
650 const Array4<Real>& derdat = mf[lev].array(mfi);
651 const Array4<Real> & p_arr =
pres.array(mfi);
654 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
656 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
660 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
661 Real gp_zeta_on_jface_lo;
663 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
664 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
665 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
666 }
else if (k == khi) {
667 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
668 p_arr(i,j,k ) + p_arr(i,j-1,k )
669 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
671 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
672 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
673 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
675 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
679 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
680 Real gp_zeta_on_jface_hi;
682 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
683 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
684 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
685 }
else if (k == khi) {
686 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
687 p_arr(i,j+1,k ) + p_arr(i,j,k )
688 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
690 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
691 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
692 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
694 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
696 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
699 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
700 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i,j+1,k) - p_arr(i,j-1,k)) * dxInv[1];
709 auto dxInv = geom[lev].InvCellSizeArray();
711 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
713 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
715 const Box& bx = mfi.tilebox();
716 const Array4<Real >& derdat = mf[lev].array(mfi);
717 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
720 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
722 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
725 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
726 Real gp_zeta_on_iface_lo;
728 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
729 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
730 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
731 }
else if (k == khi) {
732 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
733 p_arr(i-1,j,k ) + p_arr(i,j,k )
734 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
736 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
737 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
738 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
740 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
744 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
745 Real gp_zeta_on_iface_hi;
747 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
748 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
749 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
750 }
else if (k == khi) {
751 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
752 p_arr(i+1,j,k ) + p_arr(i,j,k )
753 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
755 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
756 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
757 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
759 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
761 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
769 auto dxInv = geom[lev].InvCellSizeArray();
771 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
773 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
775 const Box& bx = mfi.tilebox();
776 const Array4<Real >& derdat = mf[lev].array(mfi);
777 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
778 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
779 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
782 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
783 Real gp_zeta_on_jface_lo;
785 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
786 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
787 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
788 }
else if (k == khi) {
789 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
790 p_arr(i,j,k ) + p_arr(i,j-1,k )
791 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
793 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
794 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
795 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
797 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
801 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
802 Real gp_zeta_on_jface_hi;
804 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
805 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
806 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
807 }
else if (k == khi) {
808 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
809 p_arr(i,j+1,k ) + p_arr(i,j,k )
810 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
812 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
813 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
814 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
816 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
818 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
827 MultiFab::Copy(mf[lev],*
z_phys_cc[lev],0,mf_comp,1,0);
833 MultiFab::Copy(mf[lev],*
detJ_cc[lev],0,mf_comp,1,0);
840 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
842 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
844 const Box& bx = mfi.tilebox();
845 const Array4<Real>& derdat = mf[lev].array(mfi);
846 const Array4<Real>& mf_m =
mapfac_m[lev]->array(mfi);
847 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
848 derdat(i ,j ,k, mf_comp) = mf_m(i,j,0);
854 #ifdef ERF_USE_NETCDF
858 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
860 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
862 const Box& bx = mfi.tilebox();
863 const Array4<Real>& derdat = mf[lev].array(mfi);
864 const Array4<Real>& data = lat_m[lev]->array(mfi);
865 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
866 derdat(i, j, k, mf_comp) = data(i,j,0);
873 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
875 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
877 const Box& bx = mfi.tilebox();
878 const Array4<Real>& derdat = mf[lev].array(mfi);
879 const Array4<Real>& data = lon_m[lev]->array(mfi);
880 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
881 derdat(i, j, k, mf_comp) = data(i,j,0);
893 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
895 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
897 const Box& bx = mfi.tilebox();
898 const Array4<Real>& derdat = mf[lev].array(mfi);
899 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
901 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
903 derdat(i ,j ,k, mf_comp) = data(i,j,k,0) / norm;
911 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
913 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
915 const Box& bx = mfi.tilebox();
916 const Array4<Real>& derdat = mf[lev].array(mfi);
917 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
919 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
921 derdat(i ,j ,k, mf_comp) = data(i,j,k,1) / norm;
929 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
931 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
933 const Box& bx = mfi.tilebox();
934 const Array4<Real>& derdat = mf[lev].array(mfi);
935 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
937 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
939 derdat(i ,j ,k, mf_comp) = data(i,j,k,2) / norm;
947 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
949 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
951 const Box& bx = mfi.tilebox();
952 const Array4<Real>& derdat = mf[lev].array(mfi);
953 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
955 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
957 derdat(i ,j ,k, mf_comp) = data(i,j,k,3) / norm;
994 int n_qstate =
micro->Get_Qstate_Size();
1004 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
1037 int n_end = ncomp_cons - 1;
1039 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
1071 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1073 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1075 const Box& bx = mfi.tilebox();
1076 const Array4<Real>& derdat = mf[lev].array(mfi);
1078 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1083 erf_qsatw(
T, pressure, derdat(i,j,k,mf_comp));
1092 MultiFab::Copy(mf[lev],*(
qmoist[lev][4]),0,mf_comp,1,0);
1100 MultiFab::Copy(mf[lev],*(
qmoist[lev][8]),0,mf_comp,1,0);
1105 MultiFab::Copy(mf[lev],*(
qmoist[lev][9]),0,mf_comp,1,0);
1110 MultiFab::Copy(mf[lev],*(
qmoist[lev][10]),0,mf_comp,1,0);
1116 #ifdef ERF_USE_PARTICLES
1117 const auto& particles_namelist( particleData.getNames() );
1118 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
1120 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
1122 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
1123 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1128 Vector<std::string> particle_mesh_plot_names(0);
1129 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
1130 for (
int i = 0; i < particle_mesh_plot_names.size(); i++) {
1131 std::string plot_var_name(particle_mesh_plot_names[i]);
1133 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
1135 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
1136 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1144 MultiFab::Copy(mf[lev], EBFactory(lev).getVolFrac(), 0, mf_comp, 1, 0);
1149 #ifdef ERF_COMPUTE_ERROR
1158 Real H = geom[lev].ProbHi()[2];
1160 Real wavelength = 100.;
1161 Real kp = 2. *
PI / wavelength;
1163 Real
omega = std::sqrt(g * kp);
1166 const auto dx = geom[lev].CellSizeArray();
1169 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1171 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1173 const Box& bx = mfi.validbox();
1174 Box xbx(bx); xbx.surroundingNodes(0);
1178 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1180 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1183 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));
1185 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1188 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1190 xvel_arr(i,j,k) -= -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1193 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1195 Real
x = (i + 0.5) * dx[0];
1196 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));
1198 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1201 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1203 zvel_arr(i,j,k) -= Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1207 MultiFab temp_mf(mf[lev].boxArray(), mf[lev].DistributionMap(), AMREX_SPACEDIM, 0);
1208 average_face_to_cellcenter(temp_mf,0,
1212 MultiFab::Copy(mf[lev],temp_mf,0,mf_comp,1,0);
1216 MultiFab::Copy(mf[lev],temp_mf,1,mf_comp,1,0);
1220 MultiFab::Copy(mf[lev],temp_mf,2,mf_comp,1,0);
1226 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1228 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1230 const Box& bx = mfi.validbox();
1231 Box xbx(bx); xbx.surroundingNodes(0);
1236 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1238 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1241 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));
1242 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1246 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1247 xvel_arr(i,j,k) += -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1249 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1251 Real
x = (i + 0.5) * dx[0];
1252 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));
1253 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1256 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1258 zvel_arr(i,j,k) += Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1267 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1269 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1271 const Box& bx = mfi.tilebox();
1272 const Array4<Real>& derdat = mf[lev].array(mfi);
1273 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
1276 const auto dx = geom[lev].CellSizeArray();
1277 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1278 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
1280 Real H = geom[lev].ProbHi()[2];
1282 Real wavelength = 100.;
1283 Real kp = 2. *
PI / wavelength;
1285 Real
omega = std::sqrt(g * kp);
1288 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1291 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta) - p0_arr(i,j,k);
1293 Real rho_hse = r0_arr(i,j,k);
1295 Real
x = (i + 0.5) * dx[0];
1296 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 )
1297 +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) );
1298 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1301 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1302 Real pprime_exact = -(Ampl *
omega *
omega / kp) * fac *
1303 std::sin(kp *
x - omega_t) * r0_arr(i,j,k);
1305 derdat(i,j,k,mf_comp) -= pprime_exact;
1312 #ifdef ERF_USE_RRTMGP
1314 MultiFab::Copy(mf[lev], *(qheating_rates[lev]), 0, mf_comp, 1, 0);
1318 MultiFab::Copy(mf[lev], *(qheating_rates[lev]), 1, mf_comp, 1, 0);
1325 for (
int lev = 0; lev <= finest_level; ++lev) {
1326 EB_set_covered(mf[lev], 0.0);
1332 for (
int lev(0); lev <= finest_level; ++lev) {
1333 MultiFab::Copy(mf_nd[lev],*
z_phys_nd[lev],0,2,1,0);
1334 Real dz = Geom()[lev].CellSizeArray()[2];
1335 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
1336 const Box& bx = mfi.tilebox();
1337 Array4< Real> mf_arr = mf_nd[lev].array(mfi);
1338 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
1339 mf_arr(i,j,k,2) -= k * dz;
1345 std::string plotfilename;
1348 }
else if (which == 2) {
1357 #ifdef ERF_USE_RRTMGP
1360 if (which==1 && plot_rad) {
1366 if (finest_level == 0)
1368 if (plotfile_type == PlotFileType::Amrex)
1370 Print() <<
"Writing native plotfile " << plotfilename <<
"\n";
1373 GetVecOfConstPtrs(mf),
1374 GetVecOfConstPtrs(mf_nd),
1378 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1379 GetVecOfConstPtrs(mf),
1385 #ifdef ERF_USE_PARTICLES
1386 particleData.writePlotFile(plotfilename);
1389 }
else if (plotfile_type == PlotFileType::HDF5) {
1390 Print() <<
"Writing plotfile " << plotfilename+
"d01.h5" <<
"\n";
1391 WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1,
1392 GetVecOfConstPtrs(mf),
1396 #ifdef ERF_USE_NETCDF
1397 }
else if (plotfile_type == PlotFileType::Netcdf) {
1400 writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames,
istep,
t_new[0]);
1404 Print() <<
"Writing no plotfile since plotfile_type is none" << std::endl;
1409 if (plotfile_type == PlotFileType::Amrex) {
1412 int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]);
1413 bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 0) || (ref_ratio[lev0][1] == 0) ) ||
1414 (ref_ratio[lev0][2] == 0) );
1418 Vector<IntVect> r2(finest_level);
1419 Vector<Geometry> g2(finest_level+1);
1420 Vector<MultiFab> mf2(finest_level+1);
1422 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
1425 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
1428 Array<int,AMREX_SPACEDIM> periodicity =
1429 {Geom()[lev0].isPeriodic(0),Geom()[lev0].isPeriodic(1),Geom()[lev0].isPeriodic(2)};
1430 g2[lev0].define(Geom()[lev0].Domain(),&(Geom()[lev0].ProbDomain()),0,periodicity.data());
1432 r2[0] = IntVect(desired_ratio/ref_ratio[lev0][0],
1433 desired_ratio/ref_ratio[lev0][1],
1434 desired_ratio/ref_ratio[lev0][2]);
1436 for (
int lev = 1; lev <= finest_level; ++lev) {
1438 r2[lev-1][0] *= desired_ratio / ref_ratio[lev-1][0];
1439 r2[lev-1][1] *= desired_ratio / ref_ratio[lev-1][1];
1440 r2[lev-1][2] *= desired_ratio / ref_ratio[lev-1][2];
1443 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
1446 Box d2(Geom()[lev].Domain());
1447 d2.refine(r2[lev-1]);
1449 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
1453 for (
int lev = 1; lev <= finest_level; ++lev) {
1454 Interpolater* mapper_c = &pc_interp;
1455 InterpFromCoarseLevel(mf2[lev],
t_new[lev], mf[lev],
1463 Vector<IntVect> rr(finest_level);
1464 for (
int lev = 0; lev < finest_level; ++lev) {
1465 rr[lev] = IntVect(desired_ratio);
1468 Print() <<
"Writing plotfile " << plotfilename <<
"\n";
1471 GetVecOfConstPtrs(mf2),
1472 GetVecOfConstPtrs(mf_nd),
1476 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1477 GetVecOfConstPtrs(mf2), varnames,
1484 GetVecOfConstPtrs(mf),
1485 GetVecOfConstPtrs(mf_nd),
1489 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1490 GetVecOfConstPtrs(mf), varnames,
1497 #ifdef ERF_USE_PARTICLES
1498 particleData.writePlotFile(plotfilename);
1501 #ifdef ERF_USE_NETCDF
1502 }
else if (plotfile_type == PlotFileType::Netcdf) {
1503 for (
int lev = 0; lev <= finest_level; ++lev) {
1505 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:42
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:45
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esatw(amrex::Real t)
Definition: ERF_Microphysics_Utils.H:52
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_Microphysics_Utils.H:142
PhysBCFunctNoOp null_bc_for_fill
Definition: ERF_Plotfile.cpp:10
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:99
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:85
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:128
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:156
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:1514
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
@ omega
Definition: ERF_SAM.H:49
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