Skip to content

Commit

Permalink
Update Kessler microphysics (#1373)
Browse files Browse the repository at this point in the history 
* Updating Kessler microphysics

* Updating inputs for Kessler

---------

Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Ann Almgren <[email protected]>
Co-authored-by: Aaron M. Lattanzi <[email protected]>
  • Loading branch information
@nataraj2 @asalmgren @AMLattanzi
4 people authored Jan 12, 2024
1 parent ba76913 commit fe70d5f
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 77 deletions.
16 changes: 4 additions & 12 deletions Exec/SquallLine_2D/inputs_moisture
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 9000
max_step = 14400
stop_time = 90000.0

amrex.fpe_trap_invalid = 1
Expand All @@ -9,23 +9,14 @@ fabarray.mfiter_tile_size = 2048 1024 2048
# PROBLEM SIZE & GEOMETRY
geometry.prob_lo = -25000. 0. 0.
geometry.prob_hi = 25000. 400. 20000.
amr.n_cell = 192 4 81 # dx=dy=dz=100 m
amr.n_cell = 192 4 128 # dx=dy=dz=100 m

# periodic in x to match WRF setup
# - as an alternative, could use symmetry at x=0 and outflow at x=25600
geometry.is_periodic = 1 1 0
zlo.type = "SlipWall"
zhi.type = "SlipWall"

erf.sponge_strength = 2.0
#erf.use_zhi_sponge_damping = true
erf.zhi_sponge_start = 12000.0

erf.sponge_density = 1.2
erf.sponge_x_velocity = 0.0
erf.sponge_y_velocity = 0.0
erf.sponge_z_velocity = 0.0

# TIME STEP CONTROL
erf.use_native_mri = 1
erf.fixed_dt = 1.0 # fixed time step [s] -- Straka et al 1993
Expand All @@ -46,7 +37,7 @@ amr.check_int = 1000 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # root name of plotfile
erf.plot_int_1 = 100 # number of timesteps between plotfiles
erf.plot_int_1 = 60 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens

# SOLVER CHOICE
Expand Down Expand Up @@ -74,6 +65,7 @@ erf.dryscal_vert_adv_type = "Upwind_3rd"
erf.moistscal_horiz_adv_type = "Upwind_3rd"
erf.moistscal_vert_adv_type = "Upwind_3rd"


# PROBLEM PARAMETERS (optional)
prob.z_tr = 12000.0
prob.height = 1200.0
Expand Down
33 changes: 0 additions & 33 deletions Source/Microphysics/Kessler/Diagnose_Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,38 +7,5 @@
*/
void Kessler::Diagnose ()
{
auto qt = mic_fab_vars[MicVar_Kess::qt];
auto qp = mic_fab_vars[MicVar_Kess::qp];
auto qv = mic_fab_vars[MicVar_Kess::qv];
auto qn = mic_fab_vars[MicVar_Kess::qn];
auto qcl = mic_fab_vars[MicVar_Kess::qcl];
auto qci = mic_fab_vars[MicVar_Kess::qci];
auto qpl = mic_fab_vars[MicVar_Kess::qpl];
auto qpi = mic_fab_vars[MicVar_Kess::qpi];
auto tabs = mic_fab_vars[MicVar_Kess::tabs];

for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qt_array = qt->array(mfi);
auto qp_array = qp->array(mfi);
auto qv_array = qv->array(mfi);
auto qn_array = qn->array(mfi);
auto qcl_array = qcl->array(mfi);
auto qci_array = qci->array(mfi);
auto qpl_array = qpl->array(mfi);
auto qpi_array = qpi->array(mfi);

const auto& box3d = mfi.tilebox();

amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
amrex::Real omn = 1.0;
qcl_array(i,j,k) = qn_array(i,j,k)*omn;
qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn);
amrex::Real omp = 1.0;
qpl_array(i,j,k) = qp_array(i,j,k)*omp;
qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp);
});
}
}

3 changes: 0 additions & 3 deletions Source/Microphysics/Kessler/Init_Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,6 @@ void Kessler::Copy_State_to_Micro (const MultiFab& cons_in)

auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
auto qi_array = mic_fab_vars[MicVar_Kess::qci]->array(mfi);
auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);

auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
Expand All @@ -84,7 +82,6 @@ void Kessler::Copy_State_to_Micro (const MultiFab& cons_in)
qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
qp_array(i,j,k) = states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp);
qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);

tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
states_array(i,j,k,RhoTheta_comp),
Expand Down
47 changes: 18 additions & 29 deletions Source/Microphysics/Kessler/Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,15 +13,12 @@ using namespace amrex;
*/
void Kessler::AdvanceKessler ()
{
auto qt = mic_fab_vars[MicVar_Kess::qt];
auto qv = mic_fab_vars[MicVar_Kess::qv];
auto qc = mic_fab_vars[MicVar_Kess::qcl];
auto qp = mic_fab_vars[MicVar_Kess::qp];
auto qn = mic_fab_vars[MicVar_Kess::qn];
auto tabs = mic_fab_vars[MicVar_Kess::tabs];
auto pres = mic_fab_vars[MicVar_Kess::pres];

auto qcl = mic_fab_vars[MicVar_Kess::qcl];
auto theta = mic_fab_vars[MicVar_Kess::theta];
auto qv = mic_fab_vars[MicVar_Kess::qv];
auto rho = mic_fab_vars[MicVar_Kess::rho];

auto dz = m_geom.CellSize(2);
Expand Down Expand Up @@ -69,16 +66,13 @@ void Kessler::AdvanceKessler ()

Real dtn = dt;

// get the temperature, dentisy, theta, qt and qp from input
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi);
auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi);
auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi);
auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);

auto theta_array = theta->array(mfi);
auto qv_array = qv->array(mfi);
auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);

const auto& box3d = mfi.tilebox();
Expand All @@ -90,14 +84,9 @@ void Kessler::AdvanceKessler ()

ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));

if(qt_array(i,j,k) == 0.0){
qv_array(i,j,k) = 0.0;
qn_array(i,j,k) = 0.0;
}

//------- Autoconversion/accretion
Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat;
Expand All @@ -115,8 +104,8 @@ void Kessler::AdvanceKessler ()
dq_clwater_to_vapor = 0.0;


//Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
//Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));

// If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
if(qv_array(i,j,k) > qsat){
Expand All @@ -126,8 +115,8 @@ void Kessler::AdvanceKessler ()

// If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and
// reducing temperature
if(qv_array(i,j,k) < qsat and qn_array(i,j,k) > 0.0){
dq_clwater_to_vapor = std::min(qn_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
if(qv_array(i,j,k) < qsat and qc_array(i,j,k) > 0.0){
dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
}

if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) {
Expand All @@ -143,8 +132,8 @@ void Kessler::AdvanceKessler ()

// If there is cloud water present then do accretion and autoconversion to rain

if (qn_array(i,j,k) > 0.0) {
qcc = qn_array(i,j,k);
if (qc_array(i,j,k) > 0.0) {
qcc = qc_array(i,j,k);

autor = 0.0;
if (qcc > qcw0) {
Expand All @@ -156,7 +145,7 @@ void Kessler::AdvanceKessler ()
dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001));

// If the amount of change is more than the amount of qc present, then dq = qc
dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k));
dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k));
}

if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0;
Expand All @@ -165,15 +154,15 @@ void Kessler::AdvanceKessler ()
if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0;
//dq_sed = 0.0;

qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain;
qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor;
qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor;
qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;

theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);

qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));
});
}
}

0 comments on commit fe70d5f

Please sign in to comment.