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.

Sign up for GitHub

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

Moist redo #1347

Merged
merged 3 commits into from
Dec 15, 2023
Merged

Moist redo #1347

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
178d21f
Need to do SAM still.
Dec 13, 2023
2d820bb
Ran and restarted BUBBLE with each micro model.
AMLattanzi Dec 14, 2023
dc03e05
NullMoist doesn't have a mic_fab_var data structure. Prevent touching…
AMLattanzi Dec 15, 2023
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
14 changes: 8 additions & 6 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading