From fe70d5f9fce04491eda005c497def180edcfbbc2 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 11 Jan 2024 16:12:18 -0800 Subject: [PATCH] Update Kessler microphysics (#1373) * Updating Kessler microphysics * Updating inputs for Kessler --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Ann Almgren Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- Exec/SquallLine_2D/inputs_moisture | 16 ++----- .../Microphysics/Kessler/Diagnose_Kessler.cpp | 33 ------------- Source/Microphysics/Kessler/Init_Kessler.cpp | 3 -- Source/Microphysics/Kessler/Kessler.cpp | 47 +++++++------------ 4 files changed, 22 insertions(+), 77 deletions(-) diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 175352278..18cbe950c 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 9000 +max_step = 14400 stop_time = 90000.0 amrex.fpe_trap_invalid = 1 @@ -9,7 +9,7 @@ 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 @@ -17,15 +17,6 @@ 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 @@ -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 @@ -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 diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index e63a4a84d..c25762e80 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -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); - }); - } } diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp index 6f626550f..08ce5a893 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -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); @@ -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), diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index c7077fdd1..9f01ed269 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -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); @@ -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(); @@ -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; @@ -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){ @@ -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) { @@ -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) { @@ -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; @@ -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)); }); } }