Skip to content

Commit

Permalink
Merge branch 'development' into metgrid_interpolator
Browse files Browse the repository at this point in the history
  • Loading branch information
@wiersema1
wiersema1 committed May 24, 2024
2 parents 40b1355 + 517e56c commit 8b83b49
Show file tree
Hide file tree
Showing 11 changed files with 210 additions and 79 deletions.
81 changes: 81 additions & 0 deletions Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 2000
stop_time = 3600.0

amrex.fpe_trap_invalid = 1

fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 20000.0 400.0 10000.0
amr.n_cell = 200 4 100
geometry.is_periodic = 0 1 0
xlo.type = "SlipWall"
xhi.type = "SlipWall"
zlo.type = "SlipWall"
zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.fixed_dt = 0.5
erf.fixed_mri_dt_ratio = 4
#erf.no_substepping = 1
#erf.fixed_dt = 0.1

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
amr.v = 1 # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = 100 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 100 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc qi qrain qsnow qgraup

# SOLVER CHOICES
erf.use_gravity = true
erf.use_coriolis = false

erf.dycore_horiz_adv_type = "Upwind_3rd"
erf.dycore_vert_adv_type = "Upwind_3rd"
erf.dryscal_horiz_adv_type = "Upwind_3rd"
erf.dryscal_vert_adv_type = "Upwind_3rd"
erf.moistscal_horiz_adv_type = "Upwind_3rd"
erf.moistscal_vert_adv_type = "Upwind_3rd"

# PHYSICS OPTIONS
erf.les_type = "None"
erf.pbl_type = "None"
erf.moisture_model = "SAM_NoPrecip_NoIce"
erf.buoyancy_type = 1
erf.use_moist_background = true

erf.molec_diff_type = "ConstantAlpha"
erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
erf.dynamicViscosity = 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
erf.alpha_T = 0.0 # [m^2/s]
erf.alpha_C = 0.0

# INITIAL CONDITIONS
#erf.init_type = "input_sounding"
#erf.input_sounding_file = "BF02_moist_sounding"
#erf.init_sounding_ideal = true

# PROBLEM PARAMETERS (optional)
# warm bubble input
prob.x_c = 10000.0
prob.z_c = 2000.0
prob.x_r = 2000.0
prob.z_r = 2000.0
prob.T_0 = 300.0

prob.do_moist_bubble = true
prob.theta_pert = 2.0
prob.qt_init = 0.02
prob.eq_pot_temp = 320.0
15 changes: 14 additions & 1 deletion Exec/RegTests/Bubble/prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,6 +360,14 @@ Problem::init_custom_pert(
const Real x_r = parms.x_r, y_r = parms.y_r, z_r = parms.z_r;
const Real theta_pert = parms.theta_pert;

int moisture_type = 1;

if(sc.moisture_type == MoistureType::SAM) {
moisture_type = 1;
} else if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
moisture_type = 2;
}

amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k)
{
// Geometry (note we must include these here to get the data on device)
Expand Down Expand Up @@ -410,7 +418,12 @@ Problem::init_custom_pert(
// Cold microphysics are present
int nstate = state_pert.ncomp;
if (nstate == NVAR_max) {
Real omn = std::max(0.0,std::min(1.0,(T-tbgmin)*a_bg));
Real omn;
if(moisture_type == 1) {
omn = std::max(0.0,std::min(1.0,(T-tbgmin)*a_bg));
} else if(moisture_type == 2) {
omn = 1.0;
}
Real qn = state_pert(i, j, k, RhoQ2_comp);
state_pert(i, j, k, RhoQ2_comp) = qn * omn;
state_pert(i, j, k, RhoQ3_comp) = qn * (1.0-omn);
Expand Down
7 changes: 5 additions & 2 deletions Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ enum struct MoistureModelType{
};

enum struct MoistureType {
Kessler, SAM, Kessler_NoRain, None
Kessler, SAM, SAM_NoPrecip_NoIce, Kessler_NoRain, None
};

enum struct WindFarmType {
Expand Down Expand Up @@ -95,6 +95,8 @@ struct SolverChoice {
pp.query("moisture_model", moisture_model_string);
if (moisture_model_string == "SAM") {
moisture_type = MoistureType::SAM;
} else if (moisture_model_string == "SAM_NoPrecip_NoIce") {
moisture_type = MoistureType::SAM_NoPrecip_NoIce;
} else if (moisture_model_string == "Kessler") {
moisture_type = MoistureType::Kessler;
}else if (moisture_model_string == "Kessler_NoRain") {
Expand All @@ -107,7 +109,8 @@ struct SolverChoice {
// Set a different default for moist vs dry
if (moisture_type != MoistureType::None) {
if (moisture_type == MoistureType::Kessler_NoRain ||
moisture_type == MoistureType::SAM) {
moisture_type == MoistureType::SAM ||
moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
buoyancy_type = 1; // asserted in make buoyancy
} else {
buoyancy_type = 2; // uses Tprime
Expand Down
11 changes: 6 additions & 5 deletions Source/IO/ERF_WriteBndryPlanes.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,15 +64,16 @@ WriteBndryPlanes::WriteBndryPlanes (Vector<BoxArray>& grids,
// If the target area is contained at a finer level, use the finest data possible
for (int ilev = 0; ilev < grids.size(); ilev++) {

auto const dxi = geom[ilev].InvCellSizeArray();
const Real* xLo = m_geom[ilev].ProbLo();
auto const dxi = geom[ilev].InvCellSizeArray();
const Box& domain = m_geom[ilev].Domain();

// We create the smallest box that contains all of the cell centers
// in the physical region specified
int ilo = static_cast<int>(Math::floor(box_lo[0] * dxi[0])+.5);
int jlo = static_cast<int>(Math::floor(box_lo[1] * dxi[1])+.5);
int ihi = static_cast<int>(Math::floor(box_hi[0] * dxi[0])+.5)-1;
int jhi = static_cast<int>(Math::floor(box_hi[1] * dxi[1])+.5)-1;
int ilo = static_cast<int>(Math::floor((box_lo[0] - xLo[0]) * dxi[0])+.5);
int jlo = static_cast<int>(Math::floor((box_lo[1] - xLo[1]) * dxi[1])+.5);
int ihi = static_cast<int>(Math::floor((box_hi[0] - xLo[0]) * dxi[0])+.5)-1;
int jhi = static_cast<int>(Math::floor((box_hi[1] - xLo[1]) * dxi[1])+.5)-1;

// Map this to index space -- for now we do no interpolation
target_box.setSmall(IntVect(ilo,jlo,0));
Expand Down
3 changes: 2 additions & 1 deletion Source/Microphysics/EulerianMicrophysics.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,8 @@ public:
{
AMREX_ASSERT( Microphysics::modelType(a_model_type) == MoistureModelType::Eulerian );
m_moist_model.resize(nlev);
if (a_model_type == MoistureType::SAM) {
if (a_model_type == MoistureType::SAM or
a_model_type == MoistureType::SAM_NoPrecip_NoIce) {
SetModel<SAM>();
amrex::Print() << "SAM moisture model!\n";
} else if (a_model_type == MoistureType::Kessler or
Expand Down
1 change: 1 addition & 0 deletions Source/Microphysics/Microphysics.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ public:
static MoistureModelType modelType (const MoistureType a_moisture_type)
{
if ( (a_moisture_type == MoistureType::SAM)
|| (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce)
|| (a_moisture_type == MoistureType::Kessler)
|| (a_moisture_type == MoistureType::Kessler_NoRain)
|| (a_moisture_type == MoistureType::None) ) {
Expand Down
137 changes: 83 additions & 54 deletions Source/Microphysics/SAM/Cloud_SAM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ using namespace amrex;
/**
* Split cloud components according to saturation pressures; source theta from latent heat.
*/
void SAM::Cloud () {
void SAM::Cloud (const SolverChoice& solverChoice) {

constexpr Real an = 1.0/(tbgmax-tbgmin);
constexpr Real bn = tbgmin*an;
Expand All @@ -20,6 +20,16 @@ void SAM::Cloud () {

Real tol = 1.0e-4;

SolverChoice sc = solverChoice;

int moisture_type = 1;

if(solverChoice.moisture_type == MoistureType::SAM) {
moisture_type = 1;
} else if(solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
moisture_type = 2;
}

for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi);
auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi);
Expand Down Expand Up @@ -54,42 +64,57 @@ void SAM::Cloud () {
// This ensures the omn splitting is enforced
// before the Newton iteration, which assumes it is.

// Cloud ice not permitted (melt to form water)
if(tabs_array(i,j,k) >= tbgmax) {
omn = 1.0;
delta_qi = qci_array(i,j,k);
qci_array(i,j,k) = 0.0;
qcl_array(i,j,k) += delta_qi;
tabs_array(i,j,k) -= fac_fus * delta_qi;
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
pres_array(i,j,k) /= 100.0;
omn = 1.0;
if(moisture_type == 1){
// Cloud ice not permitted (melt to form water)
if(tabs_array(i,j,k) >= tbgmax) {
omn = 1.0;
delta_qi = qci_array(i,j,k);
qci_array(i,j,k) = 0.0;
qcl_array(i,j,k) += delta_qi;
tabs_array(i,j,k) -= fac_fus * delta_qi;
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
pres_array(i,j,k) *= 0.01;
}
// Cloud water not permitted (freeze to form ice)
else if(tabs_array(i,j,k) <= tbgmin) {
omn = 0.0;
delta_qc = qcl_array(i,j,k);
qcl_array(i,j,k) = 0.0;
qci_array(i,j,k) += delta_qc;
tabs_array(i,j,k) += fac_fus * delta_qc;
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
pres_array(i,j,k) *= 0.01;
}
// Mixed cloud phase (split according to omn)
else {
omn = an*tabs_array(i,j,k)-bn;
delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn;
delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (1.0 - omn);
qcl_array(i,j,k) = qn_array(i,j,k) * omn;
qci_array(i,j,k) = qn_array(i,j,k) * (1.0 - omn);
tabs_array(i,j,k) += fac_fus * delta_qc;
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
pres_array(i,j,k) *= 0.01;
}
}
// Cloud water not permitted (freeze to form ice)
else if(tabs_array(i,j,k) <= tbgmin) {
omn = 0.0;
delta_qc = qcl_array(i,j,k);
qcl_array(i,j,k) = 0.0;
qci_array(i,j,k) += delta_qc;
tabs_array(i,j,k) += fac_fus * delta_qc;
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
pres_array(i,j,k) /= 100.0;
}
// Mixed cloud phase (split according to omn)
else {
omn = an*tabs_array(i,j,k)-bn;
delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn;
delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (1.0 - omn);
qcl_array(i,j,k) = qn_array(i,j,k) * omn;
qci_array(i,j,k) = qn_array(i,j,k) * (1.0 - omn);
tabs_array(i,j,k) += fac_fus * delta_qc;
else if(moisture_type == 2){
// No ice. ie omn = 1.0
delta_qc = qcl_array(i,j,k) - qn_array(i,j,k);
delta_qi = 0.0;
qcl_array(i,j,k) = qn_array(i,j,k);
qci_array(i,j,k) = 0.0;
tabs_array(i,j,k) += fac_cond * delta_qc;
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
pres_array(i,j,k) /= 100.0;
pres_array(i,j,k) *= 0.01;
}

// Initial guess for temperature & pressure
Expand Down Expand Up @@ -120,19 +145,24 @@ void SAM::Cloud () {
erf_dtqsatw(tabs, pres, dqsatw);
erf_dtqsati(tabs, pres, dqsati);

// Cloud ice not permitted (condensation & fusion)
if(tabs >= tbgmax) {
omn = 1.0;
}
// Cloud water not permitted (sublimation & fusion)
else if(tabs <= tbgmin) {
omn = 0.0;
lstarw = fac_sub;
}
// Mixed cloud phase (condensation & fusion)
else {
omn = an*tabs-bn;
domn = an;
if(moisture_type == 1) {
// Cloud ice not permitted (condensation & fusion)
if(tabs >= tbgmax) {
omn = 1.0;
}
// Cloud water not permitted (sublimation & fusion)
else if(tabs <= tbgmin) {
omn = 0.0;
lstarw = fac_sub;
}
// Mixed cloud phase (condensation & fusion)
else {
omn = an*tabs-bn;
domn = an;
}
} else if(moisture_type == 2) {
omn = 1.0;
domn = 0.0;
}

// Linear combination of each component
Expand All @@ -155,8 +185,7 @@ void SAM::Cloud () {

// Update the pressure
pres = rho_array(i,j,k) * R_d * tabs
* (1.0 + R_v/R_d * qsat);
pres /= 100.0;
* (1.0 + R_v/R_d * qsat) * 0.01;

// Update iteration
niter = niter+1;
Expand Down Expand Up @@ -188,15 +217,15 @@ void SAM::Cloud () {
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);

// Pressure unit conversion
pres_array(i,j,k) /= 100.0;
pres_array(i,j,k) *= 0.01;

}
// We cannot blindly relax to qsat, but we can convert qc/qi -> qv
else {
// Changes in each component
delta_qv = qcl_array(i,j,k) + qci_array(i,j,k);
delta_qc = -qcl_array(i,j,k);
delta_qi = -qci_array(i,j,k);
delta_qv = qcl_array(i,j,k) + qci_array(i,j,k);
delta_qc = qcl_array(i,j,k);
delta_qi = qci_array(i,j,k);

// Partition the change in non-precipitating q
qv_array(i,j,k) += delta_qv;
Expand All @@ -206,7 +235,7 @@ void SAM::Cloud () {
qt_array(i,j,k) = qv_array(i,j,k);

// Update temperature (endothermic since we evap/sublime)
tabs_array(i,j,k) -= fac_fus * delta_qc + fac_sub * delta_qi;
tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi;

// Update pressure
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
Expand All @@ -216,7 +245,7 @@ void SAM::Cloud () {
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);

// Pressure unit conversion
pres_array(i,j,k) /= 100.0;
pres_array(i,j,k) *= 0.01;
}
});
} // mfi
Expand Down
2 changes: 1 addition & 1 deletion Source/Microphysics/SAM/Init_SAM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ void SAM::Copy_State_to_Micro (const MultiFab& cons_in)
tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
states_array(i,j,k,RhoTheta_comp),
qv_array(i,j,k));
pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k))/100.;
pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * 0.01;
});
}
}
Expand Down
Loading

0 comments on commit 8b83b49

Please sign in to comment.