Skip to content

Commit

Permalink
Moist redo (#1347)
Browse files Browse the repository at this point in the history 
* Need to do SAM still.

* Ran and restarted BUBBLE with each micro model.

* NullMoist doesn't have a mic_fab_var data structure. Prevent touching qmoist without a moisture model.

---------

Co-authored-by: Aaron Lattanzi <[email protected]>
  • Loading branch information
@AMLattanzi
AMLattanzi and Aaron Lattanzi authored Dec 15, 2023
1 parent 6bbf92f commit d100543
Show file tree
Hide file tree
Showing 34 changed files with 1,069 additions and 356 deletions.
14 changes: 8 additions & 6 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ PhysBCFunctNoOp null_bc;
* @param[in] time time at which the data should be filled
* @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data
*/

void
ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
{
Expand Down Expand Up @@ -100,6 +99,10 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)

// We call this even if init_type == real because this routine will fill the vertical bcs
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);

// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None)
micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
}

/*
Expand All @@ -109,7 +112,6 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
* @param[in] time time at which the data should be filled
* @param[out] mf MultiFab to be filled (qmoist[lev])
*/

void
ERF::FillPatchMoistVars (int lev, MultiFab& mf)
{
Expand Down Expand Up @@ -151,7 +153,6 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf)
* @param[in] eddyDiffs diffusion coefficients for LES turbulence models
* @param[in] allow_most_bcs if true then use MOST bcs at the low boundary
*/

void
ERF::FillIntermediatePatch (int lev, Real time,
const Vector<MultiFab*>& mfs,
Expand Down Expand Up @@ -260,9 +261,11 @@ ERF::FillIntermediatePatch (int lev, Real time,
// MOST boundary conditions
if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs)
m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get());
}

//
// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None)
micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
}

/*
* Fill valid and ghost data.
Expand All @@ -274,7 +277,6 @@ ERF::FillIntermediatePatch (int lev, Real time,
* @param[in] time time at which the data should be filled
* @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data
*/

void
ERF::FillCoarsePatch (int lev, Real time, const Vector<MultiFab*>& mfs)
{
Expand Down
6 changes: 5 additions & 1 deletion Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,11 @@ struct SolverChoice {

// Set a different default for moist vs dry
if (moisture_type != MoistureType::None) {
buoyancy_type = 2; // uses Tprime
if (moisture_type == MoistureType::FastEddy) {
buoyancy_type = 1; // asserted in make buoyancy
} else {
buoyancy_type = 2; // uses Tprime
}
}

// Is the terrain static or moving?
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 @@ -525,7 +525,7 @@ private:
amrex::Vector<amrex::MultiFab> rW_new;

Microphysics micro;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> qmoist; // This has up to 6 components: qv, qc, qi, qr, qs, qg
amrex::Vector<amrex::Vector<amrex::MultiFab*>> qmoist; // (lev,ncomp) This has up to 6 components: qv, qc, qi, qr, qs, qg

#if defined(ERF_USE_RRTMGP)
Radiation rad;
Expand Down
50 changes: 16 additions & 34 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,12 @@ ERF::ERF ()
amrex::Print() << "\n";
}

int nlevs_max = max_level + 1;

// NOTE: size micro before readparams (chooses the model at all levels)
micro.ReSize(nlevs_max);
qmoist.resize(nlevs_max);

ReadParameters();
const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
Expand All @@ -128,8 +134,6 @@ ERF::ERF ()
// No valid BoxArray and DistributionMapping have been defined.
// But the arrays for them have been resized.

int nlevs_max = max_level + 1;

istep.resize(nlevs_max, 0);
nsubsteps.resize(nlevs_max, 1);
for (int lev = 1; lev <= max_level; ++lev) {
Expand Down Expand Up @@ -157,8 +161,6 @@ ERF::ERF ()
vars_old[lev].resize(Vars::NumTypes);
}

qmoist.resize(nlevs_max);

mri_integrator_mem.resize(nlevs_max);
physbcs.resize(nlevs_max);

Expand Down Expand Up @@ -217,8 +219,7 @@ ERF::ERF ()
}
}

ERF::~ERF ()
= default;
ERF::~ERF () = default;

// advance solution to final time
void
Expand Down Expand Up @@ -529,7 +530,6 @@ ERF::InitData ()
}
}


if (solverChoice.coupling_type == CouplingType::TwoWay) {
int src_comp_reflux = 0;
int num_comp_reflux = vars_new[0][Vars::cons].nComp();
Expand All @@ -548,7 +548,7 @@ ERF::InitData ()
if (solverChoice.moisture_type != MoistureType::None)
{
for (int lev = 0; lev <= finest_level; lev++) {
FillPatchMoistVars(lev, *(qmoist[lev]));
FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component
}
}
}
Expand Down Expand Up @@ -578,23 +578,6 @@ ERF::InitData ()
}
}

// Initialize microphysics here
micro.define(solverChoice);

// Call Init which will call Diagnose to fill qmoist
if (solverChoice.moisture_type != MoistureType::None)
{
for (int lev = 0; lev <= finest_level; ++lev)
{
// If not restarting we need to fill qmoist given qt and qp.
if (restart_chkfile.empty()) {
micro.Init(vars_new[lev][Vars::cons], *(qmoist[lev]),
grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose
micro.Update(vars_new[lev][Vars::cons], *(qmoist[lev]));
}
}
}

// Configure ABLMost params if used MostWall boundary condition
// NOTE: we must set up the MOST routine before calling WritePlotFile because
// WritePlotFile calls FillPatch in order to compute gradients
Expand Down Expand Up @@ -817,6 +800,7 @@ ERF::InitData ()
void
ERF::restart ()
{
// TODO: This could be deleted since ba/dm are not created yet?
for (int lev = 0; lev <= finest_level; ++lev)
{
auto& lev_new = vars_new[lev];
Expand Down Expand Up @@ -859,8 +843,6 @@ ERF::init_only (int lev, Real time)
lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0);
lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0);

qmoist[lev]->setVal(0.);

// Initialize background flow (optional)
if (init_type == "input_sounding") {
// The base state is initialized by integrating vertically through the
Expand Down Expand Up @@ -1177,13 +1159,11 @@ ERF::MakeHorizontalAverages ()

if (use_moisture)
{
MultiFab qv(*(qmoist[lev]), make_alias, 0, 1);

for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
const Box& bx = mfi.validbox();
auto fab_arr = mf.array(mfi);
auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
auto const qv_arr = qv.const_array(mfi);
auto const qv_arr = qmoist[lev][0]->const_array(mfi);

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
Real dens = cons_arr(i, j, k, Rho_comp);
Expand Down Expand Up @@ -1464,6 +1444,12 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
amrex::Print() << "\n";
}

int nlevs_max = max_level + 1;

// NOTE: size micro before readparams (chooses the model at all levels)
micro.ReSize(nlevs_max);
qmoist.resize(nlevs_max);

ReadParameters();
const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
Expand All @@ -1475,8 +1461,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
// No valid BoxArray and DistributionMapping have been defined.
// But the arrays for them have been resized.

int nlevs_max = max_level + 1;

istep.resize(nlevs_max, 0);
nsubsteps.resize(nlevs_max, 1);
for (int lev = 1; lev <= max_level; ++lev) {
Expand Down Expand Up @@ -1504,8 +1488,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
rV_old.resize(nlevs_max);
rW_old.resize(nlevs_max);

qmoist.resize(nlevs_max);

mri_integrator_mem.resize(nlevs_max);
physbcs.resize(nlevs_max);

Expand Down
34 changes: 27 additions & 7 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,15 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
// Microphysics
// *******************************************************************************************
int q_size = micro.Get_Qmoist_Size();
qmoist[lev] = std::make_unique<MultiFab>(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg
qmoist[lev].resize(q_size);
micro.Define(lev, solverChoice);
if (solverChoice.moisture_type != MoistureType::None)
{
micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
}
for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
}

// ********************************************************************************************
// Build the data structures for calculating diffusive/turbulent terms
Expand Down Expand Up @@ -198,9 +206,16 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
//********************************************************************************************
// Microphysics
// *******************************************************************************************
int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1;
int q_size = micro.Get_Qmoist_Size();
qmoist[lev] = std::make_unique<MultiFab>(ba, dm, q_size, ngrow_state);
qmoist[lev].resize(q_size);
micro.Define(lev, solverChoice);
if (solverChoice.moisture_type != MoistureType::None)
{
micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
}
for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
}

init_stuff(lev, ba, dm);

Expand Down Expand Up @@ -305,7 +320,15 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
// Microphysics
// *******************************************************************************************
int q_size = micro.Get_Qmoist_Size();
qmoist[lev]->define(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg
qmoist[lev].resize(q_size);
micro.Define(lev, solverChoice);
if (solverChoice.moisture_type != MoistureType::None)
{
micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
}
for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
}

init_stuff(lev,ba,dm);

Expand Down Expand Up @@ -556,9 +579,6 @@ ERF::ClearLevel (int lev)
rW_new[lev].clear();
rW_old[lev].clear();

// Clears the qmoist memory
qmoist[lev].reset();

// Clears the integrator memory
mri_integrator_mem[lev].reset();
physbcs[lev].reset();
Expand Down
29 changes: 16 additions & 13 deletions Source/IO/Checkpoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,11 +130,14 @@ ERF::WriteCheckpointFile () const

IntVect ng;
if (solverChoice.moisture_type != MoistureType::None) {
// We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
ng = qmoist[lev]->nGrowVect();
MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng);
MultiFab::Copy(moist_vars,*(qmoist[lev]),0,0,qmoist[lev]->nComp(),ng);
VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
// We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
ng = qmoist[lev][mvar]->nGrowVect();
int nvar = qmoist[lev][mvar]->nComp();
MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng);
MultiFab::Copy(moist_vars,*(qmoist[lev][mvar]),0,0,nvar,ng);
VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
}
}

// Note that we write the ghost cells of the base state (unlike above)
Expand Down Expand Up @@ -300,7 +303,6 @@ ERF::ReadCheckpointFile ()
}

for (int lev = 0; lev <= finest_level; ++lev) {

// read in level 'lev' BoxArray from Header
BoxArray ba;
ba.readFrom(is);
Expand Down Expand Up @@ -338,13 +340,14 @@ ERF::ReadCheckpointFile ()
MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0);

IntVect ng;
if (solverChoice.moisture_type == MoistureType::None) {
qmoist[lev]->setVal(0.0);
} else {
ng = qmoist[lev]->nGrowVect();
MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng);
VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
MultiFab::Copy(*(qmoist[lev]),moist_vars,0,0,qmoist[lev]->nComp(),ng);
if (solverChoice.moisture_type != MoistureType::None) {
for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
ng = qmoist[lev][mvar]->nGrowVect();
int nvar = qmoist[lev][mvar]->nComp();
MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng);
VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
MultiFab::Copy(*(qmoist[lev][mvar]),moist_vars,0,0,nvar,ng);
}
}

// Note that we read the ghost cells of the base state (unlike above)
Expand Down
5 changes: 2 additions & 3 deletions Source/IO/ERF_Write1DProfiles.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,18 +246,17 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u , Gpu::HostVector<Rea

if (use_moisture)
{
MultiFab qv(*(qmoist[lev]), make_alias, 0, 1);

for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Array4<Real>& fab_arr = mf_out.array(mfi);
const Array4<Real>& fab_arr = mf_out.array(mfi);
const Array4<Real>& cons_arr = mf_cons.array(mfi);
const Array4<Real>& u_cc_arr = u_cc.array(mfi);
const Array4<Real>& v_cc_arr = v_cc.array(mfi);
const Array4<Real>& w_cc_arr = w_cc.array(mfi);
const Array4<Real>& p0_arr = p_hse.array(mfi);
const Array4<Real>& qv_arr = qv.array(mfi);
const Array4<Real>& qv_arr = qmoist[0][0]->array(mfi); // TODO: Is this written only on lev 0?

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Expand Down
Loading

0 comments on commit d100543

Please sign in to comment.