diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index b3559ff68..3ce8902f8 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -23,11 +23,21 @@ Governing Equations | | equations (instead of | | | | | the compressible equations) | | | +--------------------------+-----------------------------+---------------+-------------+ +| **erf.use_fft** | use FFT rather than | true / false | false | +| | multigrid to solve the | | | +| | the Poisson equations | | | ++--------------------------+-----------------------------+---------------+-------------+ +| **erf.mg_v** | verbosity of the multigrid | Integer >= 0 | 0 | +| | solver if used | | | +| | the Poisson equations | | | ++--------------------------+-----------------------------+---------------+-------------+ .. note:: - To solve the anelastic equations, you must set ERF_USE_POISSON_SOLVE = TRUE if using - gmake or ERF_ENABLE_POISSON_SOLVE if using cmake. + To solve the anelastic equations, you must set USE_POISSON_SOLVE = TRUE if using + gmake or ERF_ENABLE_POISSON_SOLVE if using cmake. This will enable use of the + AMReX-based Poisson solver. To optionally use the FFT solver, you must additionally + set USE_FFT = TRUE if using gmake. Problem Geometry ================ @@ -1221,7 +1231,8 @@ integrating the hydrostatic equation from the surface. If **erf.init_type = custom** or **erf.init_type = input_sounding**, ``erf.nc_init_file`` and ``erf.nc_bdy_file`` do not need to be set. -Setting **erf.project_initial_velocity = 1** will have no effect if the code is not built with **ERF_USE_POISSON_SOLVE** defined. +Setting **erf.project_initial_velocity = 1** will have no effect if the code is not built with **ERF_USE_POISSON_SOLVE** defined +if using cmake or with **USE_POISSON_SOLVE = TRUE** if using gmake. Map Scale Factors ================= diff --git a/Source/ERF.H b/Source/ERF.H index a977fa35a..6c0532a9d 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -348,7 +348,7 @@ public: const std::string &mfPrefix) const; void erf_enforce_hse (int lev, - amrex::MultiFab& dens, amrex::MultiFab& pres, amrex::MultiFab& pi, + amrex::MultiFab& dens, amrex::MultiFab& pres, amrex::MultiFab& pi, amrex::MultiFab& th, std::unique_ptr& z_cc); #ifdef ERF_USE_NETCDF diff --git a/Source/ERF.cpp b/Source/ERF.cpp index ef6ec5bfc..176d8e07a 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -552,7 +552,7 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) // Copy z_phs_nd and detJ_cc at end of timestep MultiFab::Copy(*z_phys_nd[lev], *z_phys_nd_new[lev], 0, 0, 1, z_phys_nd[lev]->nGrowVect()); MultiFab::Copy( *detJ_cc[lev], *detJ_cc_new[lev], 0, 0, 1, detJ_cc[lev]->nGrowVect()); - MultiFab::Copy(base_state[lev],base_state_new[lev],0,0,3,1); + MultiFab::Copy(base_state[lev],base_state_new[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect()); make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]); } @@ -944,7 +944,7 @@ ERF::InitData_post () // For moving terrain only if (solverChoice.terrain_type != TerrainType::Static) { - MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,3,1); + MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect()); base_state_new[lev].FillBoundary(geom[lev].periodicity()); } @@ -1675,7 +1675,7 @@ ERF::MakeHorizontalAverages () fab_arr(i, j, k, 1) = cons_arr(i, j, k, RhoTheta_comp) / dens; if (!use_moisture) { if (is_anelastic) { - fab_arr(i,j,k,2) = hse_arr(i,j,k,1); + fab_arr(i,j,k,2) = hse_arr(i,j,k,BaseState::p0_comp); } else { fab_arr(i,j,k,2) = getPgivenRTh(cons_arr(i,j,k,RhoTheta_comp)); } @@ -1695,7 +1695,7 @@ ERF::MakeHorizontalAverages () ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real dens = cons_arr(i, j, k, Rho_comp); if (is_anelastic) { - fab_arr(i,j,k,2) = hse_arr(i,j,k,1); + fab_arr(i,j,k,2) = hse_arr(i,j,k,BaseState::p0_comp); } else { Real qv = cons_arr(i, j, k, RhoQ1_comp) / dens; fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv); diff --git a/Source/ERF_IndexDefines.H b/Source/ERF_IndexDefines.H index 811946464..f7728d1bc 100644 --- a/Source/ERF_IndexDefines.H +++ b/Source/ERF_IndexDefines.H @@ -59,6 +59,17 @@ #define PrimQ5_comp (RhoQ5_comp-1) #define PrimQ6_comp (RhoQ6_comp-1) +// Base state variables +namespace BaseState { + enum { + r0_comp = 0, + p0_comp, + pi0_comp, + th0_comp, + num_comps + }; +} + // NOTE: we still use this indexing even if no moisture // NOTE: We assume a single boundary condition for all passive scalars // NOTE: and a single boundary condition for all components of the base state diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 78a249141..2c839db59 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -28,11 +28,11 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // ******************************************************************************************** // Base state holds r_0, pres_0, pi_0 (in that order) // ******************************************************************************************** - tmp_base_state.define(ba,dm,3,1); + tmp_base_state.define(ba,dm,BaseState::num_comps,3); tmp_base_state.setVal(0.); if (solverChoice.use_terrain && solverChoice.terrain_type != TerrainType::Static) { - base_state_new[lev].define(ba,dm,3,1); + base_state_new[lev].define(ba,dm,BaseState::num_comps,3); base_state_new[lev].setVal(0.); } diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 64e458e80..3ed7ad42e 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -232,7 +232,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // InterpFromCoarseLevel(base_state[lev], base_state[lev].nGrowVect(), IntVect(0,0,0), // do not fill ghost cells outside the domain - base_state[lev-1], 0, 0, 3, + base_state[lev-1], 0, 0, base_state[lev].nComp(), geom[lev-1], geom[lev], refRatio(lev-1), &cell_cons_interp, domain_bcs_type, BCVars::cons_bc); diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 75557acd4..7b6635b4d 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -349,17 +349,15 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp += 1; } - MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component + MultiFab r_hse(base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse(base_state[lev], make_alias, BaseState::p0_comp, 1); if (containerHasElement(plot_var_names, "pres_hse")) { - // p_0 is second component of base_state MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0); mf_comp += 1; } if (containerHasElement(plot_var_names, "dens_hse")) { - // r_0 is first component of base_state MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0); mf_comp += 1; } @@ -531,7 +529,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); if (solverChoice.anelastic[lev] == 1) { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - p_arr(i,j,k) = hse_arr(i,j,k,1); + p_arr(i,j,k) = hse_arr(i,j,k,BaseState::p0_comp); }); } else { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -625,7 +623,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); if (solverChoice.anelastic[lev] == 1) { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - p_arr(i,j,k) = hse_arr(i,j,k,1); + p_arr(i,j,k) = hse_arr(i,j,k,BaseState::p0_comp); }); } else { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index 450515a7a..94d560404 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -263,7 +263,7 @@ ERF::derive_diag_profiles(Real /*time*/, int nvars = vars_new[lev][Vars::cons].nComp(); MultiFab mf_cons(vars_new[lev][Vars::cons], make_alias, 0, nvars); - MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component + MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); bool use_moisture = (solverChoice.moisture_type != MoistureType::None); diff --git a/Source/IO/ERF_Write1DProfiles_stag.cpp b/Source/IO/ERF_Write1DProfiles_stag.cpp index 81fa117f5..1e6d1c37f 100644 --- a/Source/IO/ERF_Write1DProfiles_stag.cpp +++ b/Source/IO/ERF_Write1DProfiles_stag.cpp @@ -345,7 +345,7 @@ ERF::derive_diag_profiles_stag (Real /*time*/, int nvars = vars_new[lev][Vars::cons].nComp(); MultiFab mf_cons(vars_new[lev][Vars::cons], make_alias, 0, nvars); - MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component + MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); bool use_moisture = (solverChoice.moisture_type != MoistureType::None); diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index b531afc85..854b0136b 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -39,7 +39,7 @@ ERF::sum_integrated_quantities (Real time) MultiFab pert_dens(vars_new[lev][Vars::cons].boxArray(), vars_new[lev][Vars::cons].DistributionMap(), 1,0); - MultiFab r_hse (base_state[lev], make_alias, 0, 1); // r_0 is first component + MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); for ( MFIter mfi(pert_dens,TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); diff --git a/Source/Initialization/ERF_Metgrid_utils.H b/Source/Initialization/ERF_Metgrid_utils.H index b0ed1f4b2..4ebfb0a18 100644 --- a/Source/Initialization/ERF_Metgrid_utils.H +++ b/Source/Initialization/ERF_Metgrid_utils.H @@ -94,6 +94,7 @@ init_base_state_from_metgrid (const bool use_moisture, amrex::FArrayBox& r_hse_fab, amrex::FArrayBox& p_hse_fab, amrex::FArrayBox& pi_hse_fab, + amrex::FArrayBox& th_hse_fab, amrex::FArrayBox& z_phys_cc_fab, const amrex::Vector& NC_ght_fab, const amrex::Vector& NC_psfc_fab, diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 40612b852..745d8b469 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -19,17 +19,18 @@ using namespace amrex; void ERF::initHSE (int lev) { - // This integrates up through column to update p_hse, pi_hse; + // This integrates up through column to update p_hse, pi_hse, th_hse; // r_hse is not const b/c FillBoundary is called at the end for r_hse and p_hse - MultiFab r_hse (base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component - MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component + MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); + MultiFab pi_hse(base_state[lev], make_alias, BaseState::pi0_comp, 1); + MultiFab th_hse(base_state[lev], make_alias, BaseState::th0_comp, 1); bool all_boxes_touch_bottom = true; Box domain(geom[lev].Domain()); - int icomp = 0; int bccomp = BCVars::base_bc; int ncomp = 3; + int icomp = 0; int bccomp = BCVars::base_bc; int ncomp = BaseState::num_comps; Real time = 0.; @@ -75,7 +76,7 @@ ERF::initHSE (int lev) prob->erf_init_dens_hse(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); } - erf_enforce_hse(lev, r_hse, p_hse, pi_hse, z_phys_cc[lev]); + erf_enforce_hse(lev, r_hse, p_hse, pi_hse, th_hse, z_phys_cc[lev]); } else { @@ -85,12 +86,14 @@ ERF::initHSE (int lev) DistributionMapping dm_new(ba_new); - MultiFab new_base_state(ba_new, dm_new, 3, 1); - new_base_state.ParallelCopy(base_state[lev],0,0,3,1,1); + MultiFab new_base_state(ba_new, dm_new, BaseState::num_comps, base_state[lev].nGrowVect()); + new_base_state.ParallelCopy(base_state[lev],0,0,base_state[lev].nComp(), + base_state[lev].nGrowVect(),base_state[lev].nGrowVect()); - MultiFab new_r_hse (new_base_state, make_alias, 0, 1); // r_0 is first component - MultiFab new_p_hse (new_base_state, make_alias, 1, 1); // p_0 is second component - MultiFab new_pi_hse(new_base_state, make_alias, 2, 1); // pi_0 is third component + MultiFab new_r_hse (new_base_state, make_alias, BaseState::r0_comp, 1); + MultiFab new_p_hse (new_base_state, make_alias, BaseState::p0_comp, 1); + MultiFab new_pi_hse(new_base_state, make_alias, BaseState::pi0_comp, 1); + MultiFab new_th_hse(new_base_state, make_alias, BaseState::th0_comp, 1); std::unique_ptr new_z_phys_cc; std::unique_ptr new_z_phys_nd; @@ -111,17 +114,18 @@ ERF::initHSE (int lev) prob->erf_init_dens_hse(new_r_hse, new_z_phys_nd, new_z_phys_cc, geom[lev]); } - erf_enforce_hse(lev, new_r_hse, new_p_hse, new_pi_hse, new_z_phys_cc); + erf_enforce_hse(lev, new_r_hse, new_p_hse, new_pi_hse, new_th_hse, new_z_phys_cc); // Now copy back into the original arrays - base_state[lev].ParallelCopy(new_base_state,0,0,3,1,1); + base_state[lev].ParallelCopy(new_base_state,0,0,base_state[lev].nComp(), + base_state[lev].nGrowVect(),base_state[lev].nGrowVect()); } // // Impose physical bc's on the base state -- the values outside the fine region // but inside the domain have already been filled in the call above to InterpFromCoarseLevel // - (*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect(),time,bccomp); + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect(),time,bccomp); base_state[lev].FillBoundary(geom[lev].periodicity()); } @@ -147,7 +151,7 @@ ERF::initHSE () */ void ERF::erf_enforce_hse (int lev, - MultiFab& dens, MultiFab& pres, MultiFab& pi, + MultiFab& dens, MultiFab& pres, MultiFab& pi, MultiFab& theta, std::unique_ptr& z_cc) { Real l_gravity = solverChoice.gravity; @@ -181,6 +185,7 @@ ERF::erf_enforce_hse (int lev, Array4 rho_arr = dens.array(mfi); Array4 pres_arr = pres.array(mfi); Array4 pi_arr = pi.array(mfi); + Array4 th_arr = theta.array(mfi); Array4 zcc_arr; if (l_use_terrain) { zcc_arr = z_cc->array(mfi); @@ -200,12 +205,14 @@ ERF::erf_enforce_hse (int lev, hz = 0.5*dz; } - pres_arr(i,j,klo ) = p_0 - hz * rho_arr(i,j,klo) * l_gravity; - pi_arr(i,j,klo ) = getExnergivenP(pres_arr(i,j,klo ), rdOcp); + pres_arr(i,j,klo) = p_0 - hz * rho_arr(i,j,klo) * l_gravity; + pi_arr(i,j,klo) = getExnergivenP(pres_arr(i,j,klo), rdOcp); + th_arr(i,j,klo) =getRhoThetagivenP(pres_arr(i,j,klo)) / rho_arr(i,j,klo); // Set ghost cell with dz and rho at boundary pres_arr(i,j,klo-1) = p_0 + hz * rho_arr(i,j,klo) * l_gravity; - pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); + pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); + th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1); } else { // If level > 0 and klo > 0, we need to use the value of pres_arr(i,j,klo-1) which was @@ -222,6 +229,10 @@ ERF::erf_enforce_hse (int lev, pi_arr(i,j,klo ) = getExnergivenP(pres_arr(i,j,klo ), rdOcp); pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); + + th_arr(i,j,klo ) = getRhoThetagivenP(pres_arr(i,j,klo )) / rho_arr(i,j,klo ); + th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1); + } Real dens_interp; @@ -231,12 +242,14 @@ ERF::erf_enforce_hse (int lev, dens_interp = 0.5*(rho_arr(i,j,k) + rho_arr(i,j,k-1)); pres_arr(i,j,k) = pres_arr(i,j,k-1) - dz_loc * dens_interp * l_gravity; pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + th_arr(i,j,k) = getRhoThetagivenP(pres_arr(i,j,k)) / rho_arr(i,j,k); } } else { for (int k = klo+1; k <= khi; k++) { dens_interp = 0.5*(rho_arr(i,j,k) + rho_arr(i,j,k-1)); pres_arr(i,j,k) = pres_arr(i,j,k-1) - dz * dens_interp * l_gravity; pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + th_arr(i,j,k) = getRhoThetagivenP(pres_arr(i,j,k)) / rho_arr(i,j,k); } } }); @@ -251,7 +264,8 @@ ERF::erf_enforce_hse (int lev, bx.setBig(0,domlo_x-1); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(domlo_x,j,k); - pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + pi_arr(i,j,k) = pi_arr(domlo_x,j,k); + th_arr(i,j,k) = th_arr(domlo_x,j,k); }); } @@ -262,7 +276,8 @@ ERF::erf_enforce_hse (int lev, bx.setBig(0,domhi_x+1); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(domhi_x,j,k); - pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + pi_arr(i,j,k) = pi_arr(domhi_x,j,k); + th_arr(i,j,k) = th_arr(domhi_x,j,k); }); } @@ -273,7 +288,8 @@ ERF::erf_enforce_hse (int lev, bx.setBig(1,domlo_y-1); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(i,domlo_y,k); - pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + pi_arr(i,j,k) = pi_arr(i,domlo_y,k); + th_arr(i,j,k) = th_arr(i,domlo_y,k); }); } @@ -284,11 +300,14 @@ ERF::erf_enforce_hse (int lev, bx.setBig(1,domhi_y+1); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(i,domhi_y,k); - pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + pi_arr(i,j,k) = pi_arr(i,domhi_y,k); + th_arr(i,j,k) = th_arr(i,domhi_y,k); }); } } - dens.FillBoundary(geom[lev].periodicity()); - pres.FillBoundary(geom[lev].periodicity()); + dens.FillBoundary(geom[lev].periodicity()); + pres.FillBoundary(geom[lev].periodicity()); + pi.FillBoundary(geom[lev].periodicity()); + theta.FillBoundary(geom[lev].periodicity()); } diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 984c57004..4fcbc9104 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -27,8 +27,8 @@ ERF::init_custom (int lev) { auto& lev_new = vars_new[lev]; - MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component + MultiFab r_hse(base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse(base_state[lev], make_alias, BaseState::p0_comp, 1); MultiFab cons_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), lev_new[Vars::cons].nComp() , lev_new[Vars::cons].nGrow()); diff --git a/Source/Initialization/ERF_init_from_hse.cpp b/Source/Initialization/ERF_init_from_hse.cpp index 9623cb8e7..6dbc3c768 100644 --- a/Source/Initialization/ERF_init_from_hse.cpp +++ b/Source/Initialization/ERF_init_from_hse.cpp @@ -6,8 +6,6 @@ #include #include -#include "AMReX_Print.H" - using namespace amrex; /** @@ -23,8 +21,8 @@ using namespace amrex; * - save r_hse * - call ERF::enforce_hse(...), calculates p_hse from saved r_hse (redundant, * but needed because p_hse is not necessarily calculated by the Problem - * implementation) and pi_hse -- note: this pressure does not exactly match - * the p_hse from before because what is calculated by init_isentropic_hse + * implementation) and pi_hse and th_hse -- note: this pressure does not exactly + * match the p_hse from before because what is calculated by init_isentropic_hse * comes from the EOS whereas what is calculated here comes from the hydro- * static equation * @@ -35,8 +33,8 @@ ERF::init_from_hse (int lev) { auto& lev_new = vars_new[lev]; - MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component + MultiFab r_hse(base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse(base_state[lev], make_alias, BaseState::p0_comp, 1); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) diff --git a/Source/Initialization/ERF_init_from_input_sounding.cpp b/Source/Initialization/ERF_init_from_input_sounding.cpp index abbf40565..ec31635e5 100644 --- a/Source/Initialization/ERF_init_from_input_sounding.cpp +++ b/Source/Initialization/ERF_init_from_input_sounding.cpp @@ -23,6 +23,7 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx, Array4 const &r_hse_arr, Array4 const &p_hse_arr, Array4 const &pi_hse_arr, + Array4 const &th_hse_arr, GeometryData const &geomdata, Array4 const &z_cc_arr, const Real& l_gravity, @@ -79,7 +80,7 @@ ERF::init_from_input_sounding (int lev) // InterpFromCoarseLevel(base_state[lev], base_state[lev].nGrowVect(), IntVect(0,0,0), // do not fill ghost cells outside the domain - base_state[lev-1], 0, 0, 3, + base_state[lev-1], 0, 0, base_state[lev].nComp(), geom[lev-1], geom[lev], refRatio(lev-1), &cell_cons_interp, domain_bcs_type, BCVars::base_bc); @@ -88,16 +89,17 @@ ERF::init_from_input_sounding (int lev) // when the corners need to be filled by, for example, reflection of the fine ghost // cell outside the fine region but inide the domain. int bccomp = BCVars::base_bc; - int icomp = 0; int ncomp = 3; Real dummy_time = 0.; - (*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect(),dummy_time,bccomp); + Real dummy_time = 0.; + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect(),dummy_time,bccomp); } auto& lev_new = vars_new[lev]; // update if init_sounding_ideal == true - MultiFab r_hse (base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component - MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component + MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); + MultiFab pi_hse(base_state[lev], make_alias, BaseState::pi0_comp, 1); + MultiFab th_hse(base_state[lev], make_alias, BaseState::th0_comp, 1); const Real l_gravity = solverChoice.gravity; const Real l_rdOcp = solverChoice.rdOcp; @@ -112,9 +114,10 @@ ERF::init_from_input_sounding (int lev) const auto &xvel_arr = lev_new[Vars::xvel].array(mfi); const auto &yvel_arr = lev_new[Vars::yvel].array(mfi); const auto &zvel_arr = lev_new[Vars::zvel].array(mfi); - Array4 r_hse_arr = r_hse.array(mfi); - Array4 p_hse_arr = p_hse.array(mfi); + Array4 r_hse_arr = r_hse.array(mfi); + Array4 p_hse_arr = p_hse.array(mfi); Array4 pi_hse_arr = pi_hse.array(mfi); + Array4 th_hse_arr = th_hse.array(mfi); Array4 z_cc_arr = (solverChoice.use_terrain) ? z_phys_cc[lev]->const_array(mfi) : Array4{}; Array4 z_nd_arr = (solverChoice.use_terrain) ? z_phys_nd[lev]->const_array(mfi) : Array4{}; @@ -125,7 +128,7 @@ ERF::init_from_input_sounding (int lev) // calculated by calc_rho_p() init_bx_scalars_from_input_sounding_hse( bx, cons_arr, - r_hse_arr, p_hse_arr, pi_hse_arr, + r_hse_arr, p_hse_arr, pi_hse_arr, th_hse_arr, geom[lev].data(), z_cc_arr, l_gravity, l_rdOcp, l_moist, input_sounding_data); } @@ -211,6 +214,7 @@ init_bx_scalars_from_input_sounding (const Box &bx, * @param r_hse_arr Array4 specifying the density HSE base state data we are to initialize * @param p_hse_arr Array4 specifying the pressure HSE base state data we are to initialize * @param pi_hse_arr Array4 specifying the Exner pressure HSE base state data we are to initialize + * @param th_hse_arr Array4 specifying the base state potential temperature we are to initialize * @param geomdata GeometryData object specifying the domain geometry * @param l_gravity Real number specifying the gravitational acceleration constant * @param l_rdOcp Real number specifying the Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$) @@ -222,6 +226,7 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx, Array4 const &r_hse_arr, Array4 const &p_hse_arr, Array4 const &pi_hse_arr, + Array4 const &th_hse_arr, GeometryData const &geomdata, Array4 const &z_cc_arr, const Real& /*l_gravity*/, @@ -269,9 +274,10 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx, // Update hse quantities with values calculated from InputSoundingData.calc_rho_p() qv_k = (l_moist) ? interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size) : 0.0; - r_hse_arr (i, j, k) = rho_k * (1.0 + qv_k); - p_hse_arr (i, j, k) = getPgivenRTh(rhoTh_k, qv_k); - pi_hse_arr(i, j, k) = getExnergivenRTh(rhoTh_k, l_rdOcp); + r_hse_arr (i,j,k) = rho_k * (1.0 + qv_k); + p_hse_arr (i,j,k) = getPgivenRTh(rhoTh_k, qv_k); + pi_hse_arr(i,j,k) = getExnergivenRTh(rhoTh_k, l_rdOcp); + th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k)) / r_hse_arr(i,j,k); // FOEXTRAP hse arrays if (k==kbot) @@ -279,12 +285,14 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx, r_hse_arr (i, j, k-1) = r_hse_arr (i,j,k); p_hse_arr (i, j, k-1) = p_hse_arr (i,j,k); pi_hse_arr(i, j, k-1) = pi_hse_arr(i,j,k); + th_hse_arr(i, j, k-1) = th_hse_arr(i,j,k); } else if (k==ktop) { r_hse_arr (i, j, k+1) = r_hse_arr (i,j,k); p_hse_arr (i, j, k+1) = p_hse_arr (i,j,k); pi_hse_arr(i, j, k+1) = pi_hse_arr(i,j,k); + th_hse_arr(i, j, k+1) = th_hse_arr(i,j,k); } // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 79f26c3a5..d6e2660e6 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -330,12 +330,14 @@ ERF::init_from_metgrid (int lev) } // mf - MultiFab r_hse (base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component - MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component + MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); + MultiFab pi_hse(base_state[lev], make_alias, BaseState::pi0_comp, 1); + MultiFab th_hse(base_state[lev], make_alias, BaseState::th0_comp, 1); for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { FArrayBox& p_hse_fab = p_hse[mfi]; FArrayBox& pi_hse_fab = pi_hse[mfi]; + FArrayBox& th_hse_fab = th_hse[mfi]; FArrayBox& r_hse_fab = r_hse[mfi]; FArrayBox& cons_fab = lev_new[Vars::cons][mfi]; FArrayBox& z_phys_nd_fab = (*z_phys)[mfi]; @@ -344,11 +346,12 @@ ERF::init_from_metgrid (int lev) // p_hse calculate dry pressure // r_hse calculate dry density // pi_hse calculate Exner term given pressure + // th_hse calculate potential temperature const Box valid_bx = mfi.validbox(); init_base_state_from_metgrid(use_moisture, l_rdOcp, valid_bx, flag_psfc, - cons_fab, r_hse_fab, p_hse_fab, pi_hse_fab, + cons_fab, r_hse_fab, p_hse_fab, pi_hse_fab, th_hse_fab, z_phys_nd_fab, NC_ght_fab, NC_psfc_fab, fabs_for_bcs); } // mf @@ -357,6 +360,7 @@ ERF::init_from_metgrid (int lev) r_hse.FillBoundary(geom[lev].periodicity()); p_hse.FillBoundary(geom[lev].periodicity()); pi_hse.FillBoundary(geom[lev].periodicity()); + th_hse.FillBoundary(geom[lev].periodicity()); // NOTE: fabs_for_bcs is defined over the whole domain on each rank. // However, the operations needed to define the data on the ERF @@ -763,6 +767,7 @@ init_state_from_metgrid (const bool use_moisture, * @param r_hse_fab FArrayBox holding the hydrostatic base state density we are initializing * @param p_hse_fab FArrayBox holding the hydrostatic base state pressure we are initializing * @param pi_hse_fab FArrayBox holding the hydrostatic base Exner pressure we are initializing + * @param th_hse_fab FArrayBox holding the base state potential temperature we are initializing * @param z_phys_nd_fab FArrayBox holding nodal z coordinate data for terrain * @param NC_ght_fab Vector of FArrayBox objects holding metgrid data for height of cell centers * @param NC_psfc_fab Vector of FArrayBox objects holding metgrid data for surface pressure @@ -777,6 +782,7 @@ init_base_state_from_metgrid (const bool use_moisture, FArrayBox& r_hse_fab, FArrayBox& p_hse_fab, FArrayBox& pi_hse_fab, + FArrayBox& th_hse_fab, FArrayBox& z_phys_cc_fab, const Vector& /*NC_ght_fab*/, const Vector& NC_psfc_fab, @@ -819,6 +825,7 @@ init_base_state_from_metgrid (const bool use_moisture, const Array4& r_hse_arr = r_hse_fab.array(); const Array4& p_hse_arr = p_hse_fab.array(); const Array4& pi_hse_arr = pi_hse_fab.array(); + const Array4& th_hse_arr = th_hse_fab.array(); auto psfc_flag = flag_psfc_vec[0]; // ******************************************************** @@ -934,6 +941,7 @@ init_base_state_from_metgrid (const bool use_moisture, r_hse_arr(i,j,k) *= (1.0 + Qv); pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp); + th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k)) / r_hse_arr(i,j,k); }); // FOEXTRAP hse arrays @@ -942,9 +950,10 @@ init_base_state_from_metgrid (const bool use_moisture, { int jj = amrex::max(j ,valid_bx.smallEnd(1)); jj = amrex::min(jj,valid_bx.bigEnd(1)); - r_hse_arr(i,j,k) = r_hse_arr(i+1,jj,k); - p_hse_arr(i,j,k) = p_hse_arr(i+1,jj,k); + r_hse_arr(i,j,k) = r_hse_arr(i+1,jj,k); + p_hse_arr(i,j,k) = p_hse_arr(i+1,jj,k); pi_hse_arr(i,j,k) = pi_hse_arr(i+1,jj,k); + th_hse_arr(i,j,k) = th_hse_arr(i+1,jj,k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -953,39 +962,43 @@ init_base_state_from_metgrid (const bool use_moisture, r_hse_arr(i,j,k) = r_hse_arr(i-1,jj,k); p_hse_arr(i,j,k) = p_hse_arr(i-1,jj,k); pi_hse_arr(i,j,k) = pi_hse_arr(i-1,jj,k); + th_hse_arr(i,j,k) = th_hse_arr(i-1,jj,k); }); ParallelFor(gvbx_ylo, gvbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - r_hse_arr(i,j,k) = r_hse_arr(i,j+1,k); - p_hse_arr(i,j,k) = p_hse_arr(i,j+1,k); + r_hse_arr(i,j,k) = r_hse_arr(i,j+1,k); + p_hse_arr(i,j,k) = p_hse_arr(i,j+1,k); pi_hse_arr(i,j,k) = pi_hse_arr(i,j+1,k); + th_hse_arr(i,j,k) = th_hse_arr(i,j+1,k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - r_hse_arr(i,j,k) = r_hse_arr(i,j-1,k); - p_hse_arr(i,j,k) = p_hse_arr(i,j-1,k); + r_hse_arr(i,j,k) = r_hse_arr(i,j-1,k); + p_hse_arr(i,j,k) = p_hse_arr(i,j-1,k); pi_hse_arr(i,j,k) = pi_hse_arr(i,j-1,k); + th_hse_arr(i,j,k) = th_hse_arr(i,j-1,k); }); ParallelFor(gvbx_zlo, gvbx_zhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - r_hse_arr(i,j,k) = r_hse_arr(i,j,k+1); - p_hse_arr(i,j,k) = p_hse_arr(i,j,k+1); + r_hse_arr(i,j,k) = r_hse_arr(i,j,k+1); + p_hse_arr(i,j,k) = p_hse_arr(i,j,k+1); pi_hse_arr(i,j,k) = pi_hse_arr(i,j,k+1); + th_hse_arr(i,j,k) = th_hse_arr(i,j,k+1); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - r_hse_arr(i,j,k) = r_hse_arr(i,j,k-1); - p_hse_arr(i,j,k) = p_hse_arr(i,j,k-1); + r_hse_arr(i,j,k) = r_hse_arr(i,j,k-1); + p_hse_arr(i,j,k) = p_hse_arr(i,j,k-1); pi_hse_arr(i,j,k) = pi_hse_arr(i,j,k-1); + th_hse_arr(i,j,k) = th_hse_arr(i,j,k-1); }); } int ntimes = NC_psfc_fab.size(); for (int itime=0; itime& NC_ALB_fab, const Vector& NC_PB_fab, @@ -270,9 +271,10 @@ ERF::init_from_wrfinput (int lev) } // use_terrain - MultiFab r_hse (base_state[lev], make_alias, 0, 1); // r_0 is first component - MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component - MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component + MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); + MultiFab pi_hse(base_state[lev], make_alias, BaseState::pi0_comp, 1); + MultiFab th_hse(base_state[lev], make_alias, BaseState::th0_comp, 1); IntVect ng = p_hse.nGrowVect(); const Real l_rdOcp = solverChoice.rdOcp; @@ -283,11 +285,12 @@ ERF::init_from_wrfinput (int lev) FArrayBox& cons_fab = lev_new[Vars::cons][mfi]; FArrayBox& p_hse_fab = p_hse[mfi]; FArrayBox& pi_hse_fab = pi_hse[mfi]; + FArrayBox& th_hse_fab = th_hse[mfi]; FArrayBox& r_hse_fab = r_hse[mfi]; const Box gtbx = mfi.tilebox(IntVect(0), ng); init_base_state_from_wrfinput(lev, gtbx, domain, l_rdOcp, solverChoice.moisture_type, n_qstate, - cons_fab, p_hse_fab, pi_hse_fab, r_hse_fab, + cons_fab, p_hse_fab, pi_hse_fab, th_hse_fab, r_hse_fab, NC_ALB_fab, NC_PB_fab, NC_P_fab); } @@ -295,6 +298,7 @@ ERF::init_from_wrfinput (int lev) r_hse.FillBoundary(geom[lev].periodicity()); p_hse.FillBoundary(geom[lev].periodicity()); pi_hse.FillBoundary(geom[lev].periodicity()); + th_hse.FillBoundary(geom[lev].periodicity()); } if (init_type == InitType::Real && (lev == 0)) { @@ -446,6 +450,7 @@ init_msfs_from_wrfinput (int /*lev*/, FArrayBox& msfu_fab, * @param l_rdOcp Real constant specifying Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$) * @param p_hse FArrayBox specifying the hydrostatic base state pressure we initialize * @param pi_hse FArrayBox specifying the hydrostatic base state Exner pressure we initialize + * @param th_hse FArrayBox specifying the hydrostatic base state potential temperature * @param r_hse FArrayBox specifying the hydrostatic base state density we initialize * @param NC_ALB_fab Vector of FArrayBox objects containing WRF data specifying 1/density * @param NC_PB_fab Vector of FArrayBox objects containing WRF data specifying pressure @@ -460,6 +465,7 @@ init_base_state_from_wrfinput (int /*lev*/, FArrayBox& cons_fab, FArrayBox& p_hse, FArrayBox& pi_hse, + FArrayBox& th_hse, FArrayBox& r_hse, const Vector& NC_ALB_fab, const Vector& NC_PB_fab, @@ -478,6 +484,7 @@ init_base_state_from_wrfinput (int /*lev*/, const Array4& cons_arr = cons_fab.array(); const Array4& p_hse_arr = p_hse.array(); const Array4& pi_hse_arr = pi_hse.array(); + const Array4& th_hse_arr = th_hse.array(); const Array4& r_hse_arr = r_hse.array(); //const Array4& alpha_arr = NC_ALB_fab[idx].const_array(); const Array4& nc_pb_arr = NC_PB_fab[idx].const_array(); @@ -508,9 +515,10 @@ init_base_state_from_wrfinput (int /*lev*/, Real Rhse_Sum = cons_arr(ii,jj,kk,Rho_comp); for (int q_offset(0); q_offset& mom_mf, Mult MLPoisson mlpoisson(geom_tmp, ba_tmp, dm_tmp, info); - MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component + MultiFab r_hse(base_state[lev], make_alias, Basestate::r0_comp, 1); auto bclo = get_projection_bc(Orientation::low); auto bchi = get_projection_bc(Orientation::high);