Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

generalize base state to include theta0 and 3 (instead of 1) ghost cells #1908

Merged
merged 2 commits into from
Oct 24, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 14 additions & 3 deletions Docs/sphinx_doc/Inputs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
================
Expand Down Expand Up @@ -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
=================
Expand Down
2 changes: 1 addition & 1 deletion Source/ERF.H
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::MultiFab>& z_cc);

#ifdef ERF_USE_NETCDF
Expand Down
8 changes: 4 additions & 4 deletions Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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]);
}
Expand Down Expand Up @@ -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());
}

Expand Down Expand Up @@ -1680,7 +1680,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));
}
Expand All @@ -1700,7 +1700,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);
Expand Down
11 changes: 11 additions & 0 deletions Source/ERF_IndexDefines.H
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions Source/ERF_make_new_arrays.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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.);
}

Expand Down
2 changes: 1 addition & 1 deletion Source/ERF_make_new_level.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
10 changes: 4 additions & 6 deletions Source/IO/ERF_Plotfile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -349,17 +349,15 @@ ERF::WritePlotFile (int which, Vector<std::string> 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;
}
Expand Down Expand Up @@ -531,7 +529,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
const Array4<Real const>& 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 {
Expand Down Expand Up @@ -625,7 +623,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
const Array4<Real const>& 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 {
Expand Down
2 changes: 1 addition & 1 deletion Source/IO/ERF_Write1DProfiles.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down
2 changes: 1 addition & 1 deletion Source/IO/ERF_Write1DProfiles_stag.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down
2 changes: 1 addition & 1 deletion Source/IO/ERF_WriteScalarProfiles.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
1 change: 1 addition & 0 deletions Source/Initialization/ERF_Metgrid_utils.H
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::FArrayBox>& NC_ght_fab,
const amrex::Vector<amrex::FArrayBox>& NC_psfc_fab,
Expand Down
67 changes: 43 additions & 24 deletions Source/Initialization/ERF_init1d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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.;

Expand Down Expand Up @@ -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 {

Expand All @@ -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<MultiFab> new_z_phys_cc;
std::unique_ptr<MultiFab> new_z_phys_nd;
Expand All @@ -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());
}
Expand All @@ -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<MultiFab>& z_cc)
{
Real l_gravity = solverChoice.gravity;
Expand Down Expand Up @@ -181,6 +185,7 @@ ERF::erf_enforce_hse (int lev,
Array4<Real> rho_arr = dens.array(mfi);
Array4<Real> pres_arr = pres.array(mfi);
Array4<Real> pi_arr = pi.array(mfi);
Array4<Real> th_arr = theta.array(mfi);
Array4<Real> zcc_arr;
if (l_use_terrain) {
zcc_arr = z_cc->array(mfi);
Expand All @@ -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
Expand All @@ -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;
Expand All @@ -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);
}
}
});
Expand All @@ -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);
});
}

Expand All @@ -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);
});
}

Expand All @@ -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);
});
}

Expand All @@ -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());
}
4 changes: 2 additions & 2 deletions Source/Initialization/ERF_init_custom.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down
Loading
Loading