Skip to content

Commit

Permalink
generalize base state to include theta0 and 3 (instead of 1) ghost cells
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Oct 24, 2024
1 parent dda9720 commit 18bb884
Show file tree
Hide file tree
Showing 20 changed files with 158 additions and 91 deletions.
17 changes: 14 additions & 3 deletions Docs/sphinx_doc/Inputs.rst
View file
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
View file
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
View file
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 @@ -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));
}
Expand All @@ -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);
Expand Down
11 changes: 11 additions & 0 deletions Source/ERF_IndexDefines.H
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit 18bb884

Please sign in to comment.