157 const int ncomp_mf = varnames.size();
161 if (ncomp_mf == 0)
return;
166 for (
int lev = 0; lev <= finest_level; ++lev) {
167 bool fillset =
false;
176 for (
int lev = 0; lev <= finest_level; ++lev) {
177 for (
int mvar(0); mvar<
qmoist[lev].size(); ++mvar) {
178 qmoist[lev][mvar] =
micro->Get_Qmoist_Ptr(lev,mvar);
183 Vector<MultiFab> mf(finest_level+1);
184 for (
int lev = 0; lev <= finest_level; ++lev) {
185 mf[lev].define(grids[lev], dmap[lev], ncomp_mf, 0);
189 Vector<MultiFab> mf_nd(finest_level+1);
191 for (
int lev = 0; lev <= finest_level; ++lev) {
192 BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes();
193 mf_nd[lev].define(nodal_grids, dmap[lev], 3, 0);
194 mf_nd[lev].setVal(0.);
199 Vector<MultiFab> mf_cc_vel(finest_level+1);
209 for (
int lev = 0; lev <= finest_level; ++lev) {
210 mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,1));
211 average_face_to_cellcenter(mf_cc_vel[lev],0,
223 amrex::Interpolater* mapper = &cell_cons_interp;
224 for (
int lev = 1; lev <= finest_level; ++lev)
226 Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
227 Vector<Real> ftime = {
t_new[lev],
t_new[lev]};
228 Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
229 Vector<Real> ctime = {
t_new[lev],
t_new[lev]};
231 amrex::FillPatchTwoLevels(mf_cc_vel[lev],
t_new[lev], cmf, ctime, fmf, ftime,
232 0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
241 for (
int lev = 0; lev <= finest_level; ++lev)
246 AMREX_ALWAYS_ASSERT(
cons_names.size() >= ncomp_cons);
247 for (
int i = 0; i < ncomp_cons; ++i) {
256 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 0, mf_comp, 1, 0);
260 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 1, mf_comp, 1, 0);
264 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 2, mf_comp, 1, 0);
270 auto calculate_derived = [&](
const std::string& der_name,
275 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
277 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
279 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
281 const Box& bx = mfi.tilebox();
282 auto& dfab = dmf[mfi];
283 auto& sfab = src_mf[mfi];
284 der_function(bx, dfab, 0, 1, sfab, Geom(lev),
t_new[0],
nullptr, lev);
305 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
306 Array<MultiFab const*, AMREX_SPACEDIM> u;
310 computeDivergence(dmf, u, geom[lev]);
314 MultiFab r_hse(
base_state[lev], make_alias, 0, 1);
315 MultiFab p_hse(
base_state[lev], make_alias, 1, 1);
319 MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0);
325 MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0);
332 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
334 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
336 const Box& bx = mfi.tilebox();
337 const Array4<Real >& derdat = mf[lev].array(mfi);
341 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
345 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta,qv_for_p);
353 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
355 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
357 const Box& bx = mfi.tilebox();
358 const Array4<Real>& derdat = mf[lev].array(mfi);
359 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
363 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
367 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k);
375 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
377 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
379 const Box& bx = mfi.tilebox();
380 const Array4<Real>& derdat = mf[lev].array(mfi);
382 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
383 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
384 derdat(i, j, k, mf_comp) = S_arr(i,j,k,
Rho_comp) - r0_arr(i,j,k);
393 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
395 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
397 const Box& bx = mfi.tilebox();
398 const Array4<Real>& derdat = mf[lev].array(mfi);
401 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
409 derdat(i, j, k, mf_comp) =
T*std::pow((pressure - pv)/
p_0, -
R_d/fac)*std::exp(
L_v*
qv/(fac*
T)) ;
415 #ifdef ERF_USE_WINDFARM
418 std::cout <<
"Plotting num_turb" <<
"\n";
420 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
422 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
424 const Box& bx = mfi.tilebox();
425 const Array4<Real>& derdat = mf[lev].array(mfi);
426 const Array4<Real const>& Nturb_array =
Nturb[lev].const_array(mfi);
427 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
428 derdat(i, j, k, mf_comp) = Nturb_array(i,j,k,0);
435 int klo = geom[lev].Domain().smallEnd(2);
436 int khi = geom[lev].Domain().bigEnd(2);
440 auto dxInv = geom[lev].InvCellSizeArray();
443 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
445 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
448 const Box& gbx = mfi.growntilebox(1);
449 const Array4<Real> & p_arr =
pres.array(mfi);
451 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
455 pres.FillBoundary(geom[lev].periodicity());
458 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
460 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
463 const Box& bx = mfi.tilebox();
464 const Array4<Real>& derdat = mf[lev].array(mfi);
465 const Array4<Real> & p_arr =
pres.array(mfi);
468 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
470 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
475 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
476 Real gp_zeta_on_iface_lo;
478 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
479 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
480 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
481 }
else if (k == khi) {
482 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
483 p_arr(i-1,j,k ) + p_arr(i,j,k )
484 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
486 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
487 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
488 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
490 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
495 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
496 Real gp_zeta_on_iface_hi;
498 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
499 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
500 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
501 }
else if (k == khi) {
502 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
503 p_arr(i+1,j,k ) + p_arr(i,j,k )
504 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
506 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
507 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
508 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
510 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
513 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
516 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
517 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i+1,j,k) - p_arr(i-1,j,k)) * dxInv[0];
526 auto dxInv = geom[lev].InvCellSizeArray();
530 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
532 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
535 const Box& gbx = mfi.growntilebox(1);
536 const Array4<Real> & p_arr =
pres.array(mfi);
538 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept {
542 pres.FillBoundary(geom[lev].periodicity());
545 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
547 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
550 const Box& bx = mfi.tilebox();
551 const Array4<Real>& derdat = mf[lev].array(mfi);
552 const Array4<Real> & p_arr =
pres.array(mfi);
555 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
557 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
561 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
562 Real gp_zeta_on_jface_lo;
564 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
565 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
566 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
567 }
else if (k == khi) {
568 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
569 p_arr(i,j,k ) + p_arr(i,j-1,k )
570 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
572 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
573 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
574 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
576 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
580 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
581 Real gp_zeta_on_jface_hi;
583 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
584 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
585 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
586 }
else if (k == khi) {
587 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
588 p_arr(i,j+1,k ) + p_arr(i,j,k )
589 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
591 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
592 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
593 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
595 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
597 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
600 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
601 derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i,j+1,k) - p_arr(i,j-1,k)) * dxInv[1];
610 auto dxInv = geom[lev].InvCellSizeArray();
612 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
614 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
616 const Box& bx = mfi.tilebox();
617 const Array4<Real >& derdat = mf[lev].array(mfi);
618 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
621 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
623 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
626 Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
627 Real gp_zeta_on_iface_lo;
629 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
630 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
631 - p_arr(i-1,j,k ) - p_arr(i,j,k ) );
632 }
else if (k == khi) {
633 gp_zeta_on_iface_lo = 0.5 * dxInv[2] * (
634 p_arr(i-1,j,k ) + p_arr(i,j,k )
635 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
637 gp_zeta_on_iface_lo = 0.25 * dxInv[2] * (
638 p_arr(i-1,j,k+1) + p_arr(i,j,k+1)
639 - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) );
641 Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo;
645 Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
646 Real gp_zeta_on_iface_hi;
648 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
649 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
650 - p_arr(i+1,j,k ) - p_arr(i,j,k ) );
651 }
else if (k == khi) {
652 gp_zeta_on_iface_hi = 0.5 * dxInv[2] * (
653 p_arr(i+1,j,k ) + p_arr(i,j,k )
654 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
656 gp_zeta_on_iface_hi = 0.25 * dxInv[2] * (
657 p_arr(i+1,j,k+1) + p_arr(i,j,k+1)
658 - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) );
660 Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi;
662 derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi);
670 auto dxInv = geom[lev].InvCellSizeArray();
672 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
674 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
676 const Box& bx = mfi.tilebox();
677 const Array4<Real >& derdat = mf[lev].array(mfi);
678 const Array4<Real const>& p_arr = p_hse.const_array(mfi);
679 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
680 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
683 Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
684 Real gp_zeta_on_jface_lo;
686 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
687 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
688 - p_arr(i,j,k ) - p_arr(i,j-1,k ) );
689 }
else if (k == khi) {
690 gp_zeta_on_jface_lo = 0.5 * dxInv[2] * (
691 p_arr(i,j,k ) + p_arr(i,j-1,k )
692 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
694 gp_zeta_on_jface_lo = 0.25 * dxInv[2] * (
695 p_arr(i,j,k+1) + p_arr(i,j-1,k+1)
696 - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) );
698 Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo;
702 Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
703 Real gp_zeta_on_jface_hi;
705 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
706 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
707 - p_arr(i,j+1,k ) - p_arr(i,j,k ) );
708 }
else if (k == khi) {
709 gp_zeta_on_jface_hi = 0.5 * dxInv[2] * (
710 p_arr(i,j+1,k ) + p_arr(i,j,k )
711 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
713 gp_zeta_on_jface_hi = 0.25 * dxInv[2] * (
714 p_arr(i,j+1,k+1) + p_arr(i,j,k+1)
715 - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) );
717 Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi;
719 derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi);
728 MultiFab::Copy(mf[lev],*
z_phys_cc[lev],0,mf_comp,1,0);
734 MultiFab::Copy(mf[lev],*
detJ_cc[lev],0,mf_comp,1,0);
741 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
743 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
745 const Box& bx = mfi.tilebox();
746 const Array4<Real>& derdat = mf[lev].array(mfi);
747 const Array4<Real>& mf_m =
mapfac_m[lev]->array(mfi);
748 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
749 derdat(i ,j ,k, mf_comp) = mf_m(i,j,0);
755 #ifdef ERF_USE_NETCDF
759 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
761 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
763 const Box& bx = mfi.tilebox();
764 const Array4<Real>& derdat = mf[lev].array(mfi);
765 const Array4<Real>& data = lat_m[lev]->array(mfi);
766 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
767 derdat(i, j, k, mf_comp) = data(i,j,0);
774 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
776 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
778 const Box& bx = mfi.tilebox();
779 const Array4<Real>& derdat = mf[lev].array(mfi);
780 const Array4<Real>& data = lon_m[lev]->array(mfi);
781 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
782 derdat(i, j, k, mf_comp) = data(i,j,0);
794 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
796 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
798 const Box& bx = mfi.tilebox();
799 const Array4<Real>& derdat = mf[lev].array(mfi);
800 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
802 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
804 derdat(i ,j ,k, mf_comp) = data(i,j,k,0) / norm;
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>& data =
vel_t_avg[lev]->array(mfi);
820 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
822 derdat(i ,j ,k, mf_comp) = data(i,j,k,1) / norm;
830 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
832 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
834 const Box& bx = mfi.tilebox();
835 const Array4<Real>& derdat = mf[lev].array(mfi);
836 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
838 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
840 derdat(i ,j ,k, mf_comp) = data(i,j,k,2) / norm;
848 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
850 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
852 const Box& bx = mfi.tilebox();
853 const Array4<Real>& derdat = mf[lev].array(mfi);
854 const Array4<Real>& data =
vel_t_avg[lev]->array(mfi);
856 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
858 derdat(i ,j ,k, mf_comp) = data(i,j,k,3) / norm;
895 int n_qstate =
micro->Get_Qstate_Size();
905 MultiFab::Copy( mf[lev], Sm, n_start, mf_comp, 1, 0);
906 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
907 MultiFab::Add( mf[lev], Sm, n_comp, mf_comp, 1, 0);
909 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
916 MultiFab::Copy( mf[lev], Sm,
RhoQ1_comp, mf_comp, 1, 0);
917 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
924 MultiFab::Copy( mf[lev], Sm,
RhoQ2_comp, mf_comp, 1, 0);
925 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
932 MultiFab::Copy( mf[lev], Sm,
RhoQ3_comp, mf_comp, 1, 0);
933 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
942 int n_end = ncomp_cons - 1;
945 MultiFab::Copy( mf[lev], Sm, n_start, mf_comp, 1, 0);
946 for (
int n_comp(n_start+1); n_comp <= n_end; ++n_comp) {
947 MultiFab::Add( mf[lev], Sm, n_comp, mf_comp, 1, 0);
949 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
958 MultiFab::Copy( mf[lev], Sm, n_start , mf_comp, 1, 0);
959 MultiFab::Divide(mf[lev], Sm,
Rho_comp, mf_comp, 1, 0);
966 MultiFab::Copy( mf[lev], Sm,
RhoQ5_comp, mf_comp, 1, 0);
967 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
974 MultiFab::Copy( mf[lev], Sm,
RhoQ6_comp, mf_comp, 1, 0);
975 MultiFab::Divide(mf[lev], Sm,
Rho_comp , mf_comp, 1, 0);
982 MultiFab rain_accum_mf(*(
qmoist[lev][4]), make_alias, 0, 1);
983 MultiFab::Copy(mf[lev],rain_accum_mf,0,mf_comp,1,0);
991 MultiFab rain_accum_mf(*(
qmoist[lev][8]), make_alias, 0, 1);
992 MultiFab::Copy(mf[lev],rain_accum_mf,0,mf_comp,1,0);
997 MultiFab snow_accum_mf(*(
qmoist[lev][9]), make_alias, 0, 1);
998 MultiFab::Copy(mf[lev],snow_accum_mf,0,mf_comp,1,0);
1003 MultiFab graup_accum_mf(*(
qmoist[lev][10]), make_alias, 0, 1);
1004 MultiFab::Copy(mf[lev],graup_accum_mf,0,mf_comp,1,0);
1010 #ifdef ERF_USE_PARTICLES
1011 const auto& particles_namelist( particleData.getNames() );
1012 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
1014 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
1016 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
1017 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1022 Vector<std::string> particle_mesh_plot_names(0);
1023 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
1024 for (
int i = 0; i < particle_mesh_plot_names.size(); i++) {
1025 std::string plot_var_name(particle_mesh_plot_names[i]);
1027 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
1029 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
1030 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
1038 MultiFab::Copy(mf[lev], EBFactory(lev).getVolFrac(), 0, mf_comp, 1, 0);
1043 #ifdef ERF_COMPUTE_ERROR
1052 Real H = geom[lev].ProbHi()[2];
1054 Real wavelength = 100.;
1055 Real kp = 2. *
PI / wavelength;
1057 Real
omega = std::sqrt(g * kp);
1060 const auto dx = geom[lev].CellSizeArray();
1063 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1065 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1067 const Box& bx = mfi.validbox();
1068 Box xbx(bx); xbx.surroundingNodes(0);
1072 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1074 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1077 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));
1079 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1082 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1084 xvel_arr(i,j,k) -= -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1087 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1089 Real
x = (i + 0.5) * dx[0];
1090 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));
1092 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1095 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1097 zvel_arr(i,j,k) -= Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1101 MultiFab temp_mf(mf[lev].boxArray(), mf[lev].DistributionMap(), AMREX_SPACEDIM, 0);
1102 average_face_to_cellcenter(temp_mf,0,
1106 MultiFab::Copy(mf[lev],temp_mf,0,mf_comp,1,0);
1110 MultiFab::Copy(mf[lev],temp_mf,1,mf_comp,1,0);
1114 MultiFab::Copy(mf[lev],temp_mf,2,mf_comp,1,0);
1120 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1122 for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
1124 const Box& bx = mfi.validbox();
1125 Box xbx(bx); xbx.surroundingNodes(0);
1130 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1132 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1135 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));
1136 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1140 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1141 xvel_arr(i,j,k) += -Ampl *
omega * fac * std::sin(kp *
x - omega_t);
1143 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
1145 Real
x = (i + 0.5) * dx[0];
1146 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));
1147 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1150 Real fac = std::sinh( kp * (
z - H) ) / std::sinh(kp * H);
1152 zvel_arr(i,j,k) += Ampl *
omega * fac * std::cos(kp *
x - omega_t);
1161 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
1163 for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
1165 const Box& bx = mfi.tilebox();
1166 const Array4<Real>& derdat = mf[lev].array(mfi);
1167 const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
1170 const auto dx = geom[lev].CellSizeArray();
1171 const Array4<Real const>& z_nd =
z_phys_nd[lev]->const_array(mfi);
1172 const Array4<Real const>& r0_arr = r_hse.const_array(mfi);
1174 Real H = geom[lev].ProbHi()[2];
1176 Real wavelength = 100.;
1177 Real kp = 2. *
PI / wavelength;
1179 Real
omega = std::sqrt(g * kp);
1182 ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
1185 derdat(i, j, k, mf_comp) =
getPgivenRTh(rhotheta) - p0_arr(i,j,k);
1187 Real rho_hse = r0_arr(i,j,k);
1189 Real
x = (i + 0.5) * dx[0];
1190 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 )
1191 +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) );
1192 Real z_base = Ampl * std::sin(kp *
x - omega_t);
1195 Real fac = std::cosh( kp * (
z - H) ) / std::sinh(kp * H);
1196 Real pprime_exact = -(Ampl *
omega *
omega / kp) * fac *
1197 std::sin(kp *
x - omega_t) * r0_arr(i,j,k);
1199 derdat(i,j,k,mf_comp) -= pprime_exact;
1208 for (
int lev = 0; lev <= finest_level; ++lev) {
1209 EB_set_covered(mf[lev], 0.0);
1215 for (
int lev(0); lev <= finest_level; ++lev) {
1216 MultiFab::Copy(mf_nd[lev],*
z_phys_nd[lev],0,2,1,0);
1217 Real dz = Geom()[lev].CellSizeArray()[2];
1218 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
1219 const Box& bx = mfi.tilebox();
1220 Array4< Real> mf_arr = mf_nd[lev].array(mfi);
1221 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
1222 mf_arr(i,j,k,2) -= k * dz;
1228 std::string plotfilename;
1231 else if (which == 2)
1239 if (finest_level == 0)
1242 Print() <<
"Writing native plotfile " << plotfilename <<
"\n";
1245 GetVecOfConstPtrs(mf),
1246 GetVecOfConstPtrs(mf_nd),
1250 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1251 GetVecOfConstPtrs(mf),
1257 #ifdef ERF_USE_PARTICLES
1258 particleData.Checkpoint(plotfilename);
1262 Print() <<
"Writing plotfile " << plotfilename+
"d01.h5" <<
"\n";
1263 WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1,
1264 GetVecOfConstPtrs(mf),
1268 #ifdef ERF_USE_NETCDF
1272 writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames,
istep,
t_new[0]);
1275 Print() <<
"User specified plot_filetype = " <<
plotfile_type << std::endl;
1276 Abort(
"Dont know this plot_filetype");
1281 Vector<IntVect> r2(finest_level);
1282 Vector<Geometry> g2(finest_level+1);
1283 Vector<MultiFab> mf2(finest_level+1);
1285 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
1288 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
1291 Array<int,AMREX_SPACEDIM> periodicity =
1292 {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)};
1293 g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data());
1296 r2[0] = IntVect(1,1,ref_ratio[0][0]);
1297 for (
int lev = 1; lev <= finest_level; ++lev) {
1301 r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0];
1304 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
1307 Box d2(Geom()[lev].Domain());
1308 d2.refine(r2[lev-1]);
1310 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
1314 for (
int lev = 1; lev <= finest_level; ++lev) {
1315 Interpolater* mapper_c = &pc_interp;
1316 InterpFromCoarseLevel(mf2[lev],
t_new[lev], mf[lev],
1324 Vector<IntVect> rr(finest_level);
1325 for (
int lev = 0; lev < finest_level; ++lev) {
1326 rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]);
1329 Print() <<
"Writing plotfile " << plotfilename <<
"\n";
1332 GetVecOfConstPtrs(mf),
1333 GetVecOfConstPtrs(mf_nd),
1337 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
1338 GetVecOfConstPtrs(mf2), varnames,
1344 #ifdef ERF_USE_PARTICLES
1345 particleData.Checkpoint(plotfilename);
1347 #ifdef ERF_USE_NETCDF
1349 for (
int lev = 0; lev <= finest_level; ++lev) {
1351 writeNCPlotFile(lev, which_box, plotfilename, GetVecOfConstPtrs(mf), varnames,
istep,
t_new[0]);
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: EOS.H:29
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
#define RhoQ6_comp
Definition: IndexDefines.H:22
#define RhoQ5_comp
Definition: IndexDefines.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esatw(amrex::Real t)
Definition: Microphysics_Utils.H:44
PhysBCFunctNoOp null_bc_for_fill
Definition: 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: TerrainMetrics.H:93
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: TerrainMetrics.H:79
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: TerrainMetrics.H:122
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: TerrainMetrics.H:150
void FillBdyCCVels(amrex::Vector< amrex::MultiFab > &mf_cc_vel)
Definition: ERF_FillPatch.cpp:548
static amrex::Vector< std::string > PlotFileVarNames(amrex::Vector< std::string > plot_var_names)
Definition: Plotfile.cpp:142
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, amrex::Real time, const amrex::Vector< int > &level_steps, 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: Plotfile.cpp:1360
void writeJobInfo(const std::string &dir) const
Definition: writeJobInfo.cpp:9
void Plot_Lsm_Data(amrex::Real time, const amrex::Vector< int > &level_steps, const amrex::Vector< amrex::IntVect > &ref_ratio)
Definition: LandSurface.H:80
@ Theta_v
Definition: IndexDefines.H:128
@ Mom_h
Definition: IndexDefines.H:121
@ Mom_v
Definition: IndexDefines.H:127
@ PBL_lengthscale
Definition: IndexDefines.H:133
@ Theta_h
Definition: IndexDefines.H:122
@ omega
Definition: 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: Derive.cpp:198
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: Derive.cpp:226
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: Derive.cpp:282
void erf_dernull(const Box &, FArrayBox &, int, int, const FArrayBox &, const Geometry &, Real, const int *, const int)
Definition: 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: Derive.cpp:254
void erf_dertemp(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: Derive.cpp:91
void erf_derKE(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: Derive.cpp:163
void erf_derQKE(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: Derive.cpp:184
void erf_dersoundspeed(const Box &bx, FArrayBox &derfab, int, int, const FArrayBox &datfab, const Geometry &, Real, const int *, const int)
Definition: Derive.cpp:58