diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce new file mode 100644 index 000000000..fc18f8f23 --- /dev/null +++ b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble_SAM_NoPrecip_NoIce @@ -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 diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index 8957ab3f1..b6382081b 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -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) @@ -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); diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index b738ec678..0ab5ea5cc 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -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 { @@ -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") { @@ -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 diff --git a/Source/IO/ERF_WriteBndryPlanes.cpp b/Source/IO/ERF_WriteBndryPlanes.cpp index fc8d493e0..eda3353c9 100644 --- a/Source/IO/ERF_WriteBndryPlanes.cpp +++ b/Source/IO/ERF_WriteBndryPlanes.cpp @@ -64,15 +64,16 @@ WriteBndryPlanes::WriteBndryPlanes (Vector& 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(Math::floor(box_lo[0] * dxi[0])+.5); - int jlo = static_cast(Math::floor(box_lo[1] * dxi[1])+.5); - int ihi = static_cast(Math::floor(box_hi[0] * dxi[0])+.5)-1; - int jhi = static_cast(Math::floor(box_hi[1] * dxi[1])+.5)-1; + int ilo = static_cast(Math::floor((box_lo[0] - xLo[0]) * dxi[0])+.5); + int jlo = static_cast(Math::floor((box_lo[1] - xLo[1]) * dxi[1])+.5); + int ihi = static_cast(Math::floor((box_hi[0] - xLo[0]) * dxi[0])+.5)-1; + int jhi = static_cast(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)); diff --git a/Source/Microphysics/EulerianMicrophysics.H b/Source/Microphysics/EulerianMicrophysics.H index 4df1af178..4e1ac1af2 100644 --- a/Source/Microphysics/EulerianMicrophysics.H +++ b/Source/Microphysics/EulerianMicrophysics.H @@ -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(); amrex::Print() << "SAM moisture model!\n"; } else if (a_model_type == MoistureType::Kessler or diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index efd601460..ea96f0e79 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -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) ) { diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 8c9264f3e..89483d887 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -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; @@ -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); @@ -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 @@ -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 @@ -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; @@ -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; @@ -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) @@ -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 diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index 0fc219432..1d593044f 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -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; }); } } diff --git a/Source/Microphysics/SAM/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp index e28aadfae..296cc776c 100644 --- a/Source/Microphysics/SAM/Precip.cpp +++ b/Source/Microphysics/SAM/Precip.cpp @@ -167,12 +167,12 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta + // Latent heat source for theta (endothermic fusion & exothermic melting) tabs_array(i,j,k) += fac_fus * ( dqca * (1.0 - omp) - dqia * omp ); 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; } //================================================== @@ -188,12 +188,12 @@ void SAM::Precip () { dqpg = evapg1_t(k)*sqrt(qpg) + evapg2_t(k)*pow(qpg,powg2); // NOTE: This is always a sink for precipitating comps - // since qv0. If we are // in a super-saturated state (qv>qsat) the Newton // iterations in Cloud() will have handled condensation. - dqpr *= -dtn * (qv_array(i,j,k)/qsat - 1.0); - dqps *= -dtn * (qv_array(i,j,k)/qsat - 1.0); - dqpg *= -dtn * (qv_array(i,j,k)/qsat - 1.0); + dqpr *= dtn * (1.0 - qv_array(i,j,k)/qsat); + dqps *= dtn * (1.0 - qv_array(i,j,k)/qsat); + dqpg *= dtn * (1.0 - qv_array(i,j,k)/qsat); // Limit to avoid negative moisture fractions dqpr = std::min(qpr_array(i,j,k),dqpr); @@ -211,8 +211,8 @@ void SAM::Precip () { qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - // Latent heat source for theta - tabs_array(i,j,k) += fac_cond * dqpr + fac_sub * (dqps + dqpg); + // Latent heat source for theta (endothermic) + tabs_array(i,j,k) -= fac_cond * dqpr + fac_sub * (dqps + dqpg); 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); diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 998712d89..3b4d9f280 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -63,7 +63,7 @@ public: virtual ~SAM () = default; // cloud physics - void Cloud (); + void Cloud (const SolverChoice& sc); // ice physics void IceFall (); @@ -121,14 +121,16 @@ public: // wrapper to do all the updating void Advance (const amrex::Real& dt_advance, - const SolverChoice& /*solverChoice*/) override + const SolverChoice& sc) override { dt = dt_advance; - this->Cloud(); - this->IceFall(); - this->Precip(); - this->PrecipFall(); + this->Cloud(sc); + if(sc.moisture_type == MoistureType::SAM) { + this->IceFall(); + this->Precip(); + this->PrecipFall(); + } } amrex::MultiFab* diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index 5dabd32dd..65533381e 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -136,7 +136,7 @@ void make_buoyancy (Vector& S_data, AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); } - if (solverChoice.moisture_type == MoistureType::SAM) { + if (solverChoice.moisture_type == MoistureType::SAM or solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) { AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); }