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

More steps towards making moisture a run-time setting #1335

Merged
merged 3 commits into from
Dec 11, 2023
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
5 changes: 5 additions & 0 deletions CMake/BuildERFExe.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,10 @@ function(build_erf_lib erf_lib_name)
${SRC_DIR}/Microphysics/Kessler/Kessler.cpp
${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp
${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp
${SRC_DIR}/Microphysics/FastEddy/Init_FE.cpp
${SRC_DIR}/Microphysics/FastEddy/FastEddy.cpp
${SRC_DIR}/Microphysics/FastEddy/Diagnose_FE.cpp
${SRC_DIR}/Microphysics/FastEddy/Update_FE.cpp
)

if(NOT "${erf_exe_name}" STREQUAL "erf_unit_tests")
Expand Down Expand Up @@ -203,6 +207,7 @@ function(build_erf_lib erf_lib_name)
target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null)
target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM)
target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler)
target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/FastEddy)

if(ERF_ENABLE_RRTMGP)
target_link_libraries(${erf_lib_name} PUBLIC yakl)
Expand Down
3 changes: 2 additions & 1 deletion Docs/sphinx_doc/Inputs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -893,7 +893,8 @@ List of Parameters
| Parameter | Definition | Acceptable | Default |
| | | Values | |
+=============================+==========================+====================+============+
| **erf.moisture_model** | Name of moisture model | "SAM", "Kessler" | "Null" |
| **erf.moisture_model** | Name of moisture model | "SAM", "Kessler", | "Null" |
| | | "FastEddy" | |
+-----------------------------+--------------------------+--------------------+------------+
| **erf.do_cloud** | use basic moisture model | true / false | true |
+-----------------------------+--------------------------+--------------------+------------+
Expand Down
5 changes: 5 additions & 0 deletions Exec/Make.ERF
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,11 @@ include $(ERF_MOISTURE_NULL_DIR)/Make.package
VPATH_LOCATIONS += $(ERF_MOISTURE_NULL_DIR)
INCLUDE_LOCATIONS += $(ERF_MOISTURE_NULL_DIR)

ERF_MOISTURE_FE_DIR = $(ERF_SOURCE_DIR)/Microphysics/FastEddy
include $(ERF_MOISTURE_FE_DIR)/Make.package
VPATH_LOCATIONS += $(ERF_MOISTURE_FE_DIR)
INCLUDE_LOCATIONS += $(ERF_MOISTURE_FE_DIR)

ERF_MOISTURE_SAM_DIR = $(ERF_SOURCE_DIR)/Microphysics/SAM
include $(ERF_MOISTURE_SAM_DIR)/Make.package
VPATH_LOCATIONS += $(ERF_MOISTURE_SAM_DIR)
Expand Down
4 changes: 3 additions & 1 deletion Source/DataStructs/DataStruct.H
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ enum struct TerrainType {
};

enum struct MoistureType {
Kessler, SAM, None
Kessler, SAM, FastEddy, None
};

enum struct Coord {
Expand Down Expand Up @@ -90,6 +90,8 @@ struct SolverChoice {
moisture_type = MoistureType::SAM;
} else if (moisture_model_string == "Kessler") {
moisture_type = MoistureType::Kessler;
} else if (moisture_model_string == "FastEddy") {
moisture_type = MoistureType::FastEddy;
} else {
moisture_type = MoistureType::None;
}
Expand Down
4 changes: 2 additions & 2 deletions Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -905,8 +905,6 @@ ERF::init_only (int lev, Real time)
lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity());
lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity());
lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity());

print_state(vars_new[0][Vars::cons], IntVect(30,30,2));
}

// read in some parameters from inputs file
Expand Down Expand Up @@ -965,6 +963,8 @@ ERF::ReadParameters ()
micro.SetModel<SAM>();
} else if (solverChoice.moisture_type == MoistureType::Kessler) {
micro.SetModel<Kessler>();
} else if (solverChoice.moisture_type == MoistureType::FastEddy) {
micro.SetModel<FastEddy>();
} else if (solverChoice.moisture_type == MoistureType::None) {
micro.SetModel<NullMoist>();
amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n";
Expand Down
90 changes: 54 additions & 36 deletions Source/IO/Checkpoint.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -106,34 +106,34 @@ ERF::WriteCheckpointFile () const
}
}

// write the MultiFab data to, e.g., chk00010/Level_0/
// Here we make copies of the MultiFab with no ghost cells
for (int lev = 0; lev <= finest_level; ++lev)
{
MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0);
MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0);
VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell"));

MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0);
MultiFab::Copy(xvel,vars_new[lev][Vars::xvel],0,0,1,0);
VisMF::Write(xvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "XFace"));

MultiFab yvel(convert(grids[lev],IntVect(0,1,0)),dmap[lev],1,0);
MultiFab::Copy(yvel,vars_new[lev][Vars::yvel],0,0,1,0);
VisMF::Write(yvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "YFace"));

MultiFab zvel(convert(grids[lev],IntVect(0,0,1)),dmap[lev],1,0);
MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0);
VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace"));

IntVect ng;
#ifdef ERF_USE_MOISTURE
// 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"));
#endif
// write the MultiFab data to, e.g., chk00010/Level_0/
// Here we make copies of the MultiFab with no ghost cells
for (int lev = 0; lev <= finest_level; ++lev)
{
MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0);
MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0);
VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell"));

MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0);
MultiFab::Copy(xvel,vars_new[lev][Vars::xvel],0,0,1,0);
VisMF::Write(xvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "XFace"));

MultiFab yvel(convert(grids[lev],IntVect(0,1,0)),dmap[lev],1,0);
MultiFab::Copy(yvel,vars_new[lev][Vars::yvel],0,0,1,0);
VisMF::Write(yvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "YFace"));

MultiFab zvel(convert(grids[lev],IntVect(0,0,1)),dmap[lev],1,0);
MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0);
VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace"));

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"));
}

// Note that we write the ghost cells of the base state (unlike above)
ng = base_state[lev].nGrowVect();
Expand Down Expand Up @@ -249,7 +249,14 @@ ERF::ReadCheckpointFile ()
// conservative, cell-centered vars
is >> chk_ncomp;
GotoNextLine(is);
AMREX_ASSERT(chk_ncomp == Cons::NumVars);
#if 0
if (solverChoice.moisture_type != MoistureType::None) {
AMREX_ASSERT(chk_ncomp == Cons::NumVars);
} else {
// This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them
AMREX_ASSERT(chk_ncomp == Cons::NumVars-2);
}
#endif

// x-velocity on faces
is >> chk_ncomp;
Expand Down Expand Up @@ -314,7 +321,16 @@ ERF::ReadCheckpointFile ()
{
MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0);
VisMF::Read(cons, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell"));
MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0);

if (solverChoice.moisture_type != MoistureType::None) {
MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0);
} else {
// This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them
MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars-2,0);

// We zero them here so they won't cause problems while we're still debugging
vars_new[lev][Vars::cons].setVal(0.0, RhoQ1_comp, 2);
}

MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0);
VisMF::Read(xvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "XFace"));
Expand All @@ -329,12 +345,14 @@ ERF::ReadCheckpointFile ()
MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0);

IntVect ng;
#ifdef ERF_USE_MOISTURE
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);
#endif
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);
}

// Note that we read the ghost cells of the base state (unlike above)
ng = base_state[lev].nGrowVect();
Expand Down
56 changes: 39 additions & 17 deletions Source/IO/ERF_Write1DProfiles.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -194,9 +194,7 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u , Gpu::HostVector<Rea

MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component

#if defined(ERF_USE_MOISTURE)
MultiFab qv(qmoist[lev], make_alias, 0, 1);
#endif
bool use_moisture = (solverChoice.moisture_type != MoistureType::None);

for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Expand All @@ -206,9 +204,6 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u , Gpu::HostVector<Rea
const Array4<Real>& v_cc_arr = v_cc.array(mfi);
const Array4<Real>& w_cc_arr = w_cc.array(mfi);
const Array4<Real>& cons_arr = mf_cons.array(mfi);
#if defined(ERF_USE_MOISTURE)
const Array4<Real>& qv_arr = qv.array(mfi);
#endif
const Array4<Real>& p0_arr = p_hse.array(mfi);

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
Expand Down Expand Up @@ -237,18 +232,45 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u , Gpu::HostVector<Rea
fab_arr(i, j, k,14) = tke * u_cc_arr(i,j,k); // ku
fab_arr(i, j, k,15) = tke * v_cc_arr(i,j,k); // kv
fab_arr(i, j, k,16) = tke * w_cc_arr(i,j,k); // kw
#if defined(ERF_USE_MOISTURE)
Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
#else
Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp));
#endif
p -= p0_arr(i,j,k);
fab_arr(i, j, k,17) = p; // p'
fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k); // p'u
fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k); // p'v
fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k); // p'w

if (!use_moisture) {
Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp));
p -= p0_arr(i,j,k);
fab_arr(i, j, k,17) = p; // p'
fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k); // p'u
fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k); // p'v
fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k); // p'w
}
});
}
} // mfi

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>& 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);

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));

p -= p0_arr(i,j,k);
fab_arr(i, j, k,17) = p; // p'
fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k); // p'u
fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k); // p'v
fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k); // p'w
});
} // mfi
} // use_moisture

h_avg_rho = sumToLine(mf_out, 0,1,domain,zdir);
h_avg_th = sumToLine(mf_out, 1,1,domain,zdir);
Expand Down
5 changes: 3 additions & 2 deletions Source/IO/Plotfile.cpp.convergence
Original file line number Diff line number Diff line change
Expand Up @@ -544,7 +544,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
mf_comp ++;
}

#if defined(ERF_USE_MOISTURE)
if (solverChoice.moisture_type != MoistureType::None) {
calculate_derived("qt", derived::erf_derQt);
calculate_derived("qp", derived::erf_derQp);

Expand Down Expand Up @@ -590,7 +590,8 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0);
mf_comp += 1;
}
#elif defined(ERF_USE_WARM_NO_PRECIP)
}
#if defined(ERF_USE_WARM_NO_PRECIP)
calculate_derived("qv", derived::erf_derQv);
calculate_derived("qc", derived::erf_derQc);
#endif
Expand Down
19 changes: 9 additions & 10 deletions Source/IO/ReadFromWRFInput.cpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "NCWpsFile.H"
#include "AMReX_FArrayBox.H"
#include "DataStruct.H"

using namespace amrex;

Expand All @@ -16,14 +17,12 @@ read_from_wrfinput (int lev,
FArrayBox& NC_MSFM_fab, FArrayBox& NC_SST_fab,
FArrayBox& NC_C1H_fab , FArrayBox& NC_C2H_fab,
FArrayBox& NC_RDNW_fab,
#if defined(ERF_USE_MOISTURE)
FArrayBox& NC_QVAPOR_fab,
FArrayBox& NC_QCLOUD_fab,
FArrayBox& NC_QRAIN_fab,
#elif defined(ERF_USE_WARM_NO_PRECIP)
#endif
FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab,
FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab)
FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab,
MoistureType moisture_type)
{
amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;

Expand Down Expand Up @@ -51,12 +50,12 @@ read_from_wrfinput (int lev,
NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 15
NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 16
NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17
#if defined(ERF_USE_MOISTURE)
NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 18
NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 19
NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20
# elif defined(ERF_USE_WARM_NO_PRECIP)
#endif

if (moisture_type != MoistureType::None) {
NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 18
NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 19
NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20
}

// Read the netcdf file and fill these FABs
amrex::Print() << "Building initial FABS from file " << fname << std::endl;
Expand Down
Loading
Loading