Skip to content

Commit

Permalink
add RhoQv_comp and RhoQr_comp as integers that we can test on to dete… (
Browse files Browse the repository at this point in the history 
#1780)

* add RhoQv_comp and RhoQr_comp as integers that we can test on to determine moisture model; defaults are -1

* remove unused

* remove unused

* fixes for last commit

* fix gpu oops

* fix typo

* fix more oops

* allow NSCALARS > 1
  • Loading branch information
@asalmgren
asalmgren authored Sep 2, 2024
1 parent aee27dd commit 66784a7
Show file tree
Hide file tree
Showing 30 changed files with 458 additions and 204 deletions.
12 changes: 5 additions & 7 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,16 +47,14 @@ public:
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Hwave,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Lwave,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& eddyDiffs,
int n_qstate,
amrex::Real start_bdy_time = 0.0,
amrex::Real bdy_time_interval = 0.0)
: m_exp_most(use_exp_most),
m_rotate(use_rot_most),
m_start_bdy_time(start_bdy_time),
m_bdy_time_interval(bdy_time_interval),
m_geom(geom),
m_ma(geom,vars_old,Theta_prim,Qv_prim,Qr_prim,z_phys_nd),
m_n_qstate(n_qstate)
m_ma(geom,vars_old,Theta_prim,Qv_prim,Qr_prim,z_phys_nd)
{
// We have a moisture model if Qv_prim is a valid pointer
use_moisture = (Qv_prim[0].get());
Expand Down Expand Up @@ -369,14 +367,16 @@ public:
void
update_pblh (const int& lev,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
amrex::MultiFab* z_phys_cc);
amrex::MultiFab* z_phys_cc,
const int RhoQv_comp, const int RhoQr_comp);

template <typename PBLHeightEstimator>
void
compute_pblh (const int& lev,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
amrex::MultiFab* z_phys_cc,
const PBLHeightEstimator& est);
const PBLHeightEstimator& est,
const int RhoQv_comp, const int RhoQr_comp);

void
read_custom_roughness (const int& lev,
Expand Down Expand Up @@ -529,8 +529,6 @@ private:
amrex::Vector<amrex::MultiFab*> m_Hwave_lev;
amrex::Vector<amrex::MultiFab*> m_Lwave_lev;
amrex::Vector<amrex::MultiFab*> m_eddyDiffs_lev;

int m_n_qstate;
};

#endif /* ABLMOST_H */
15 changes: 6 additions & 9 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -560,11 +560,12 @@ ABLMost::get_lsm_tsurf (const int& lev)
void
ABLMost::update_pblh (const int& lev,
Vector<Vector<MultiFab>>& vars,
MultiFab* z_phys_cc)
MultiFab* z_phys_cc,
int RhoQv_comp, int RhoQr_comp)
{
if (pblh_type == PBLHeightCalcType::MYNN25) {
MYNNPBLH estimator;
compute_pblh(lev, vars, z_phys_cc, estimator);
compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQr_comp);
} else if (pblh_type == PBLHeightCalcType::YSU) {
amrex::Error("YSU PBLH calc not implemented yet");
}
Expand All @@ -575,16 +576,12 @@ void
ABLMost::compute_pblh (const int& lev,
Vector<Vector<MultiFab>>& vars,
MultiFab* z_phys_cc,
const PBLHeightEstimator& est)
const PBLHeightEstimator& est,
int RhoQv_comp, int RhoQr_comp)
{
int moist_flag = (m_n_qstate > 0) ? 1 : 0;
if (m_n_qstate > 3) {
moist_flag = RhoQ4_comp;
}

est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
vars[lev][Vars::cons],m_lmask_lev[lev][0],
moist_flag);
RhoQv_comp, RhoQr_comp);
}

void
Expand Down
15 changes: 15 additions & 0 deletions Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,14 +104,22 @@ struct SolverChoice {
pp.query("moisture_model", moisture_model_string);
if (moisture_model_string == "SAM") {
moisture_type = MoistureType::SAM;
RhoQv_comp = RhoQ1_comp;
RhoQr_comp = RhoQ4_comp;
} else if (moisture_model_string == "SAM_NoIce") {
moisture_type = MoistureType::SAM_NoIce;
RhoQv_comp = RhoQ1_comp;
RhoQr_comp = RhoQ4_comp;
} else if (moisture_model_string == "SAM_NoPrecip_NoIce") {
moisture_type = MoistureType::SAM_NoPrecip_NoIce;
RhoQv_comp = RhoQ1_comp;
} else if (moisture_model_string == "Kessler") {
moisture_type = MoistureType::Kessler;
RhoQv_comp = RhoQ1_comp;
RhoQr_comp = RhoQ3_comp;
}else if (moisture_model_string == "Kessler_NoRain") {
moisture_type = MoistureType::Kessler_NoRain;
RhoQv_comp = RhoQ1_comp;
} else {
moisture_type = MoistureType::None;
}
Expand Down Expand Up @@ -601,6 +609,13 @@ struct SolverChoice {
bool do_cloud {true};
bool do_precip {true};
bool use_moist_background {false};
int RhoQv_comp {-1};

// This component will be model-dependent:
// if a model with no rain, this will stay -1
// if Kessler, then it will be set to RhoQ3
// if SAM, then it will be set to RhoQ4
int RhoQr_comp {-1};

amrex::Real latitude_lo=-1e10, longitude_lo=-1e10;
std::string windfarm_loc_table, windfarm_spec_table;
Expand Down
9 changes: 6 additions & 3 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -465,7 +465,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
const Geometry& geom,
const MultiFab& mapfac_u, const MultiFab& mapfac_v,
const std::unique_ptr<MultiFab>& z_phys_nd,
const TurbChoice& turbChoice, const Real const_grav,
const SolverChoice& solverChoice,
std::unique_ptr<ABLMost>& most,
const bool& exp_most,
const bool& use_moisture,
Expand All @@ -484,6 +484,9 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
// ComputeTurbulentViscosityPBL computes the PBL viscosity just for the vertical component.
//

TurbChoice turbChoice = solverChoice.turbChoice[level];
const Real const_grav = solverChoice.gravity;

if (most) {
bool l_use_turb = ( turbChoice.les_type == LESType::Smagorinsky ||
turbChoice.les_type == LESType::Deardorff ||
Expand All @@ -508,11 +511,11 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
if (turbChoice.pbl_type == PBLType::MYNN25) {
ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity,
geom, turbChoice, most, use_moisture,
level, bc_ptr, vert_only, z_phys_nd);
level, bc_ptr, vert_only, z_phys_nd,
solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
} else if (turbChoice.pbl_type == PBLType::YSU) {
ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity,
geom, turbChoice, most, use_moisture,
level, bc_ptr, vert_only, z_phys_nd);
}

}
16 changes: 6 additions & 10 deletions Source/Diffusion/Diffusion.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,13 +61,12 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain,
amrex::Array4< amrex::Real>& qfx2_z,
amrex::Array4< amrex::Real>& diss,
const amrex::Array4<const amrex::Real>& mu_turb,
const DiffChoice& diffChoice,
const TurbChoice& turbChoice,
const SolverChoice& solverChoice,
const int level,
const amrex::Array4<const amrex::Real>& tm_arr,
const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> grav_gpu,
const amrex::BCRec* bc_ptr,
const bool use_most,
const int n_qstate);
const bool use_most);

void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
int start_comp, int num_comp,
Expand Down Expand Up @@ -100,15 +99,12 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
amrex::Array4< amrex::Real>& qfx2_z,
amrex::Array4< amrex::Real>& diss,
const amrex::Array4<const amrex::Real>& mu_turb,
const DiffChoice& diffChoice,
const TurbChoice& turbChoice,
const SolverChoice& solverChoice,
const int level,
const amrex::Array4<const amrex::Real>& tm_arr,
const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> grav_gpu,
const amrex::BCRec* bc_ptr,
const bool use_most,
const int n_qstate);


const bool use_most);

void ComputeStressConsVisc_N (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz, amrex::Box tbxyz, amrex::Real mu_eff,
const amrex::Array4<const amrex::Real>& cell_data,
Expand Down
18 changes: 9 additions & 9 deletions Source/Diffusion/DiffusionSrcForState_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,16 +58,18 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
Array4< Real>& qfx2_z,
Array4< Real>& diss,
const Array4<const Real>& mu_turb,
const DiffChoice &diffChoice,
const TurbChoice &turbChoice,
const SolverChoice &solverChoice,
const int level,
const Array4<const Real>& tm_arr,
const GpuArray<Real,AMREX_SPACEDIM> grav_gpu,
const BCRec* bc_ptr,
const bool use_most,
const int n_qstate)
const bool use_most)
{
BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N);

DiffChoice diffChoice = solverChoice.diffChoice;
TurbChoice turbChoice = solverChoice.turbChoice[level];

amrex::ignore_unused(use_most);

const Real dx_inv = cellSizeInv[0];
Expand Down Expand Up @@ -667,10 +669,8 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
int qty_index = RhoQKE_comp;
auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;

int moist_comp = (n_qstate > 0) ? 1 : 0;
if (n_qstate > 3) {
moist_comp = RhoQ4_comp;
}
const int rhoqv_comp = solverChoice.RhoQv_comp;
const int rhoqr_comp = solverChoice.RhoQr_comp;

ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Expand All @@ -685,7 +685,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
mu_turb,cellSizeInv,domain,
pbl_mynn_B1_l,tm_arr(i,j,0),
moist_comp,
rhoqv_comp, rhoqr_comp,
c_ext_dir_on_zlo, c_ext_dir_on_zhi,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi);
Expand Down
18 changes: 9 additions & 9 deletions Source/Diffusion/DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,16 +71,18 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
Array4< Real>& /*qfx2_z*/,
Array4< Real>& diss,
const Array4<const Real>& mu_turb,
const DiffChoice &diffChoice,
const TurbChoice &turbChoice,
const SolverChoice &solverChoice,
const int level,
const Array4<const Real>& tm_arr,
const GpuArray<Real,AMREX_SPACEDIM> grav_gpu,
const BCRec* bc_ptr,
const bool use_most,
const int n_qstate)
const bool use_most)
{
BL_PROFILE_VAR("DiffusionSrcForState_T()",DiffusionSrcForState_T);

DiffChoice diffChoice = solverChoice.diffChoice;
TurbChoice turbChoice = solverChoice.turbChoice[level];

amrex::ignore_unused(use_most);

const Real dx_inv = cellSizeInv[0];
Expand Down Expand Up @@ -733,10 +735,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
int qty_index = RhoQKE_comp;
auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;

int moist_comp = (n_qstate > 0) ? 1 : 0;
if (n_qstate > 3) {
moist_comp = RhoQ4_comp;
}
const int rhoqv_comp = solverChoice.RhoQv_comp;
const int rhoqr_comp = solverChoice.RhoQr_comp;

ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Expand All @@ -752,7 +752,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
mu_turb,cellSizeInv,domain,
pbl_mynn_B1_l,tm_arr(i,j,0),
moist_comp,
rhoqv_comp, rhoqr_comp,
c_ext_dir_on_zlo, c_ext_dir_on_zhi,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi,
Expand Down
2 changes: 1 addition & 1 deletion Source/Diffusion/EddyViscosity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
const amrex::Geometry& geom,
const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v,
const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
const TurbChoice& turbChoice, const amrex::Real const_grav,
const SolverChoice& solverChoice,
std::unique_ptr<ABLMost>& most,
const bool& exp_most,
const bool& use_moisture,
Expand Down
31 changes: 18 additions & 13 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -980,12 +980,10 @@ ERF::InitData ()
}
}

const int n_qstate = micro->Get_Qstate_Size();

m_most = std::make_unique<ABLMost>(geom, use_exp_most, use_rot_most,
vars_old, Theta_prim, Qv_prim, Qr_prim, z_phys_nd,
sst_lev, lmask_lev, lsm_data, lsm_flux,
Hwave, Lwave, eddyDiffs_lev, n_qstate
Hwave, Lwave, eddyDiffs_lev
#ifdef ERF_USE_NETCDF
,start_bdy_time, bdy_time_interval
#endif
Expand All @@ -1003,25 +1001,32 @@ ERF::InitData ()
{
Real time = t_new[lev];
IntVect ng = Theta_prim[lev]->nGrowVect();
MultiFab S(vars_new[lev][Vars::cons],make_alias,0,RhoTheta_comp+1);
MultiFab::Copy( *Theta_prim[lev], S, RhoTheta_comp, 0, 1, ng);
MultiFab::Divide(*Theta_prim[lev], S, Rho_comp , 0, 1, ng);

MultiFab::Copy( *Theta_prim[lev], vars_new[lev][Vars::cons], RhoTheta_comp, 0, 1, ng);
MultiFab::Divide(*Theta_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);

if (solverChoice.moisture_type != MoistureType::None) {
ng = Qv_prim[lev]->nGrowVect();
int RhoQr_comp = (micro->Get_Qstate_Size() > 3) ? RhoQ4_comp : RhoQ3_comp;
MultiFab Sm(vars_new[lev][Vars::cons],make_alias,0,RhoQr_comp+1);
MultiFab::Copy( *Qv_prim[lev], Sm, RhoQ1_comp, 0, 1, ng);
MultiFab::Copy( *Qr_prim[lev], Sm, RhoQr_comp, 0, 1, ng);
MultiFab::Divide(*Qv_prim[lev], Sm, Rho_comp , 0, 1, ng);
MultiFab::Divide(*Qr_prim[lev], Sm, Rho_comp , 0, 1, ng);

MultiFab::Copy( *Qv_prim[lev], vars_new[lev][Vars::cons], RhoQ1_comp, 0, 1, ng);
MultiFab::Divide(*Qv_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);

int rhoqr_comp = solverChoice.RhoQr_comp;
if (rhoqr_comp > -1) {
MultiFab::Copy( *Qr_prim[lev], vars_new[lev][Vars::cons], rhoqr_comp, 0, 1, ng);
MultiFab::Divide(*Qr_prim[lev], vars_new[lev][Vars::cons], Rho_comp, 0, 1, ng);
} else {
Qr_prim[lev]->setVal(0.0);
}
}
m_most->update_mac_ptrs(lev, vars_new, Theta_prim, Qv_prim, Qr_prim);

if (restart_chkfile == "") {
// Only do this if starting from scratch; if restarting, then
// we don't want to call update_fluxes multiple times because
// it will change u* and theta* from their previous values
m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get());
m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get(),
solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
m_most->update_fluxes(lev, time);
}
}
Expand Down
10 changes: 4 additions & 6 deletions Source/ERF_Tagging.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,14 @@ ERF::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)

// This allows dynamic refinement based on the value of qv
} else if ( ref_tags[j].Field() == "qv" ) {
MultiFab qv(vars_new[levc][Vars::cons],make_alias,0,RhoQ1_comp+1);
MultiFab::Copy( *mf, qv, RhoQ1_comp, 0, 1, 0);
MultiFab::Divide(*mf, qv, Rho_comp , 0, 1, 0);
MultiFab::Copy( *mf, vars_new[levc][Vars::cons], RhoQ1_comp, 0, 1, 0);
MultiFab::Divide(*mf, vars_new[levc][Vars::cons], Rho_comp, 0, 1, 0);


// This allows dynamic refinement based on the value of qc
} else if (ref_tags[j].Field() == "qc" ) {
MultiFab qc(vars_new[levc][Vars::cons],make_alias,0,RhoQ2_comp+1);
MultiFab::Copy( *mf, qc, RhoQ2_comp, 0, 1, 0);
MultiFab::Divide(*mf, qc, Rho_comp , 0, 1, 0);
MultiFab::Copy( *mf, vars_new[levc][Vars::cons], RhoQ2_comp, 0, 1, 0);
MultiFab::Divide(*mf, vars_new[levc][Vars::cons], Rho_comp, 0, 1, 0);


// This allows dynamic refinement based on the value of the scalar/pressure/theta
Expand Down
Loading

0 comments on commit 66784a7

Please sign in to comment.