From e9cec7237bbbf2ddad76621e8cdd5fde01b25860 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 12 Dec 2024 14:37:00 -0800 Subject: [PATCH] Saturation Adjustment Micro Model (#2020) * SatAdj micro model tested on moist bubble rise. * Add docs. --------- Co-authored-by: Ann Almgren --- CMake/BuildERFExe.cmake | 4 + Docs/sphinx_doc/Inputs.rst | 2 +- Docs/sphinx_doc/theory/Microphysics.rst | 5 + Exec/Make.ERF | 5 + Exec/MoistRegTests/Bubble/ERF_Prob.cpp | 3 - Source/DataStructs/ERF_DataStruct.H | 14 +- Source/IO/ERF_Plotfile.cpp | 2 +- .../Microphysics/ERF_EulerianMicrophysics.H | 13 +- Source/Microphysics/ERF_Microphysics.H | 6 +- Source/Microphysics/Kessler/ERF_Kessler.H | 6 +- Source/Microphysics/SAM/ERF_SAM.H | 8 +- Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp | 76 +++++++ Source/Microphysics/SatAdj/ERF_SatAdj.H | 214 ++++++++++++++++++ Source/Microphysics/SatAdj/ERF_SatAdj.cpp | 83 +++++++ .../Microphysics/SatAdj/ERF_UpdateSatAdj.cpp | 37 +++ Source/Microphysics/SatAdj/Make.package | 4 + 16 files changed, 459 insertions(+), 23 deletions(-) create mode 100644 Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp create mode 100644 Source/Microphysics/SatAdj/ERF_SatAdj.H create mode 100644 Source/Microphysics/SatAdj/ERF_SatAdj.cpp create mode 100644 Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp create mode 100644 Source/Microphysics/SatAdj/Make.package diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 47e1e8cab..dc4888ff7 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -207,6 +207,9 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/Kessler/ERF_InitKessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_UpdateKessler.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_InitSatAdj.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_SatAdj.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -255,6 +258,7 @@ endif() target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) + target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 46ec4148a..f654250d3 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -1443,7 +1443,7 @@ List of Parameters | | | Values | | +=============================+==========================+====================+============+ | **erf.moisture_model** | Name of moisture model | "SAM", "Kessler", | "Null" | -| | | "FastEddy" | | +| | | "SatAdj" | | +-----------------------------+--------------------------+--------------------+------------+ | **erf.do_cloud** | use basic moisture model | true / false | true | +-----------------------------+--------------------------+--------------------+------------+ diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index b62d6704a..afe563eff 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -170,3 +170,8 @@ The evaporation rate of rain is .. math:: Q_{revp} = 2\pi(S-1)n_{0R}[0.78\lambda_{R}^{-2}+0.31S_{c}^{1/3}\Gamma[(b+5)/2]a^{1/2}\mu^{-1/2}(\frac{\rho_{0}}{\rho})^{1/4}\lambda_{R}^{(b+5)/2}](\frac{1}{\rho})(\frac{L_{v}^{2}}{K_{0}R_{w}T^{2}}+\frac{1}{\rho r_{s}\psi})^{-1} + + Saturation Adjustment Microphysics Model +---------------------------------- +The saturation adjustment microphysics model is the simplest possible moisture model and only transports the +water vapor mixing ratio, :math:`q_v`, and the cloud water mixing ration, :math:`q_c`. Evaporation, :math:`q_v \longrightarrow q_c`, and condensation, :math:`q_c \longrightarrow q_v`, are the only relevant mechanisms. The final saturation state, :math:`q_v = q_{vs}(T)` is obtained from Newton-Raphson iterations on the thermal temperature :math:`T`. diff --git a/Exec/Make.ERF b/Exec/Make.ERF index d90ccc172..d25135f68 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -149,6 +149,11 @@ include $(ERF_MOISTURE_KESSLER_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) +ERF_MOISTURE_SATADJ_DIR = $(ERF_SOURCE_DIR)/Microphysics/SatAdj +include $(ERF_MOISTURE_SATADJ_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_SATADJ_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_SATADJ_DIR) + # If using windfarm parametrization, then compile all models and choose # at runtime from the inputs ifeq ($(USE_WINDFARM), TRUE) diff --git a/Exec/MoistRegTests/Bubble/ERF_Prob.cpp b/Exec/MoistRegTests/Bubble/ERF_Prob.cpp index a28cdfdaa..f079ba5d9 100644 --- a/Exec/MoistRegTests/Bubble/ERF_Prob.cpp +++ b/Exec/MoistRegTests/Bubble/ERF_Prob.cpp @@ -303,7 +303,6 @@ Problem::init_custom_pert( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - state_pert(i, j, k, RhoQ3_comp) = 0.0; } }); } // do_moist_bubble @@ -389,7 +388,6 @@ Problem::init_custom_pert( // mean states state_pert(i, j, k, RhoQ1_comp) = rho*q_v_hot; state_pert(i, j, k, RhoQ2_comp) = rho*(parms_d.qt_init - q_v_hot); - state_pert(i, j, k, RhoQ3_comp) = 0.0; // Cold microphysics are present int nstate = state_pert.ncomp; @@ -429,7 +427,6 @@ Problem::init_custom_pert( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - state_pert(i, j, k, RhoQ3_comp) = 0.0; } }); } // do_moist_bubble diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 0be9f91e5..c7e9e6b32 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -42,7 +42,7 @@ AMREX_ENUM(MoistureModelType, ); AMREX_ENUM(MoistureType, - Kessler, SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler_NoRain, None + SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, None ); AMREX_ENUM(WindFarmType, @@ -130,15 +130,19 @@ struct SolverChoice { } else if (moisture_type == MoistureType::Kessler_NoRain) { RhoQv_comp = RhoQ1_comp; RhoQc_comp = RhoQ2_comp; + } else if (moisture_type == MoistureType::SatAdj) { + RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; } // TODO: should we set default for dry?? // 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_NoIce || - moisture_type == MoistureType::SAM_NoPrecip_NoIce) + if (moisture_type == MoistureType::Kessler_NoRain || + moisture_type == MoistureType::SAM || + moisture_type == MoistureType::SAM_NoIce || + moisture_type == MoistureType::SAM_NoPrecip_NoIce || + moisture_type == MoistureType::SatAdj) { buoyancy_type = 1; // uses Rhoprime } else { diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index fb1d331e6..26ed76e16 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -1072,7 +1072,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p // Precipitating components //-------------------------------------------------------------------------- - if(containerHasElement(plot_var_names, "qp")) + if(containerHasElement(plot_var_names, "qp") && (n_qstate >= 3)) { int n_start = (n_qstate > 3) ? RhoQ4_comp : RhoQ3_comp; int n_end = ncomp_cons - 1; diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index cb91854b4..edf9c34cd 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -8,6 +8,7 @@ #include "ERF_NullMoist.H" #include "ERF_SAM.H" #include "ERF_Kessler.H" +#include "ERF_SatAdj.H" #include "ERF_Microphysics.H" /*! \brief Eulerian microphysics interface */ @@ -34,6 +35,8 @@ public: } else if (a_model_type == MoistureType::Kessler || a_model_type == MoistureType::Kessler_NoRain) { SetModel(); + } else if (a_model_type == MoistureType::SatAdj) { + SetModel(); } else if (a_model_type == MoistureType::None) { SetModel(); } else { @@ -108,12 +111,12 @@ public: /*! \brief get the indices and names of moisture model variables for restart at a given level */ - void Get_Qmoist_Restart_Vars ( const int a_lev, /*!< level */ - const SolverChoice& a_sc, /*!< Solver choice object */ - std::vector& a_idx, /*!< indices */ - std::vector& a_names /*!< names */ ) const override + void Get_Qmoist_Restart_Vars (const int a_lev, /*!< level */ + const SolverChoice& a_sc, /*!< Solver choice object */ + std::vector& a_idx, /*!< indices */ + std::vector& a_names /*!< names */ ) const override { - m_moist_model[a_lev]->Qmoist_Restart_Vars( a_sc, a_idx, a_names ); + m_moist_model[a_lev]->Qmoist_Restart_Vars(a_sc, a_idx, a_names); } protected: diff --git a/Source/Microphysics/ERF_Microphysics.H b/Source/Microphysics/ERF_Microphysics.H index 29eab027d..084572e0f 100644 --- a/Source/Microphysics/ERF_Microphysics.H +++ b/Source/Microphysics/ERF_Microphysics.H @@ -58,7 +58,10 @@ public: /*! \brief get the indices and names of moisture model variables for restart at a given level */ - virtual void Get_Qmoist_Restart_Vars ( int, const SolverChoice&, std::vector&, std::vector& ) const = 0; + virtual void Get_Qmoist_Restart_Vars (int, + const SolverChoice&, + std::vector&, + std::vector& ) const = 0; /*! \brief query if a specified moisture model is Eulerian or Lagrangian */ static MoistureModelType modelType (const MoistureType a_moisture_type) @@ -68,6 +71,7 @@ public: || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce) || (a_moisture_type == MoistureType::Kessler) || (a_moisture_type == MoistureType::Kessler_NoRain) + || (a_moisture_type == MoistureType::SatAdj) || (a_moisture_type == MoistureType::None) ) { return MoistureModelType::Eulerian; } else { diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index f361f013a..f890ef35b 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -125,9 +125,9 @@ public: Qstate_Size () override { return Kessler::m_qstate_size; } void - Qmoist_Restart_Vars ( const SolverChoice& a_sc, - std::vector& a_idx, - std::vector& a_names) const override + Qmoist_Restart_Vars (const SolverChoice& a_sc, + std::vector& a_idx, + std::vector& a_names) const override { a_idx.clear(); a_names.clear(); diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index e514b509a..7652b3093 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -264,9 +264,9 @@ public: } void - Qmoist_Restart_Vars ( const SolverChoice& a_sc, - std::vector& a_idx, - std::vector& a_names) const override + Qmoist_Restart_Vars (const SolverChoice& a_sc, + std::vector& a_idx, + std::vector& a_names) const override { a_idx.clear(); a_names.clear(); @@ -283,7 +283,7 @@ private: // Number of qmoist variables (qt, qv, qcl, qci, qp, qpr, qps, qpg) int m_qmoist_size = 11; - // Number of qmoist variables + // Number of qstate variables int m_qstate_size = 6; // CFL MAX for vertical advection diff --git a/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp new file mode 100644 index 000000000..5ce0b67c9 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp @@ -0,0 +1,76 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void SatAdj::Init (const MultiFab& cons_in, + const BoxArray& /*grids*/, + const Geometry& geom, + const Real& dt_advance, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*detJ_cc*/) +{ + dt = dt_advance; + m_geom = geom; + + MicVarMap.resize(m_qmoist_size); + MicVarMap = {MicVar_SatAdj::qv, MicVar_SatAdj::qc}; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_SatAdj::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), + 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } +} + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + */ +void SatAdj::Copy_State_to_Micro (const MultiFab& cons_in) +{ + // Get the temperature, density, theta, qt and qp from input + for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) { + const auto& tbx = mfi.tilebox(); + + auto states_array = cons_in.array(mfi); + + auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + + auto rho_array = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi); + + // Get pressure, theta, temperature, density, and qt, qp + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + rho_array(i,j,k) = states_array(i,j,k,Rho_comp); + theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); + 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); + + 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)); + + // Pressure in [mbar] for qsat evaluation + 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/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H new file mode 100644 index 000000000..8b42743b5 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -0,0 +1,214 @@ +#ifndef ERF_SATADJ_H +#define ERF_SATADJ_H + +/* + * Implementation saturation adjustment microphysics model + * The model transports qv and qc and does Newton iterations + * to complete condensation/evaporation. + */ + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_EOS.H" +#include "ERF_Constants.H" +#include "ERF_MicrophysicsUtils.H" +#include "ERF_IndexDefines.H" +#include "ERF_DataStruct.H" +#include "ERF_NullMoist.H" +#include "ERF_TileNoZ.H" + +namespace MicVar_SatAdj { + enum { + // independent variables + rho=0, // density + theta, // liquid/ice water potential temperature + tabs, // temperature + pres, // pressure + // non-precipitating vars + qv, // cloud vapor + qc, // cloud water + NumVars + }; +} + +class SatAdj : public NullMoist { + + using FabPtr = std::shared_ptr; + +public: + // constructor + SatAdj () {} + + // destructor + virtual ~SatAdj () = default; + + // cloud physics + void AdvanceSatAdj (const SolverChoice& /*solverChoice*/); + + // Set up for first time + void + Define (SolverChoice& sc) override + { + m_fac_cond = lcond / sc.c_p; + m_rdOcp = sc.rdOcp; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + const amrex::BoxArray& /*grids*/, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*detJ_cc*/) override; + + // Copy state into micro vars + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + + // Copy state into micro vars + void + Copy_Micro_to_State (amrex::MultiFab& cons_in) override; + + // update micro vars + void + Update_Micro_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_State_to_Micro(cons_in); + } + + // update state vars + void + Update_State_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_Micro_to_State(cons_in); + } + + // wrapper to advance micro vars + void + Advance (const amrex::Real& dt_advance, + const SolverChoice& solverChoice) override + { + dt = dt_advance; + + this->AdvanceSatAdj(solverChoice); + } + + amrex::MultiFab* + Qmoist_Ptr (const int& varIdx) override + { + AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size); + return mic_fab_vars[MicVarMap[varIdx]].get(); + } + + int + Qmoist_Size () override { return SatAdj::m_qmoist_size; } + + int + Qstate_Size () override { return SatAdj::m_qstate_size; } + + void + Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/, + std::vector& a_idx, + std::vector& a_names) const override + { + a_idx.clear(); + a_names.clear(); + } + + AMREX_GPU_HOST_DEVICE + AMREX_FORCE_INLINE + static amrex::Real + NewtonIterSat (int& i, int& j, int& k, + const amrex::Real& fac_cond, + const amrex::Array4& tabs_array, + const amrex::Array4& pres_array, + const amrex::Array4& qv_array, + const amrex::Array4& qc_array) + { + // Solution tolerance + amrex::Real tol = 1.0e-8; + + // Saturation moisture fractions + amrex::Real qsat, dqsat; + + // Newton iteration vars + int niter; + amrex::Real fff, dfff; + amrex::Real dtabs, delta_qv; + + // Initial guess for temperature & pressure + amrex::Real tabs = tabs_array(i,j,k); + amrex::Real pres = pres_array(i,j,k); + + niter = 0; + dtabs = 1; + + //================================================== + // Newton iteration to qv=qsat (cloud phase only) + //================================================== + do { + // Saturation moisture fractions + erf_qsatw(tabs, pres, qsat); + erf_dtqsatw(tabs, pres, dqsat); + + // Function for root finding: + // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) + fff = -tabs + tabs_array(i,j,k) + fac_cond*(qv_array(i,j,k) - qsat); + + // Derivative of function (T_new iterated on) + dfff = -1.0 - fac_cond*dqsat; + + // Update the temperature + dtabs = -fff/dfff; + tabs += dtabs; + + // Update iteration + niter = niter+1; + } while(std::abs(dtabs) > tol && niter < 20); + + // Update qsat from last iteration (dq = dq/dt * dt) + qsat += dqsat*dtabs; + + // Changes in each component + delta_qv = qv_array(i,j,k) - qsat; + + // Partition the change in non-precipitating q + qv_array(i,j,k) = qsat; + qc_array(i,j,k) += delta_qv; + + // Return to temperature + return tabs; + } + +private: + // Number of qmoist variables (qv, qc) + int m_qmoist_size = 2; + + // Number of qstate variables + int m_qstate_size = 2; + + // MicVar map (Qmoist indices -> MicVar enum) + amrex::Vector MicVarMap; + + // geometry + amrex::Geometry m_geom; + + // timestep + amrex::Real dt; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_rdOcp ; + + // independent variables + amrex::Array mic_fab_vars; +}; +#endif diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp new file mode 100644 index 000000000..0956e4156 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp @@ -0,0 +1,83 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + +/** + * Compute Precipitation-related Microphysics quantities. + */ +void SatAdj::AdvanceSatAdj (const SolverChoice& /*solverChoice*/) +{ + auto tabs = mic_fab_vars[MicVar_SatAdj::tabs]; + + // Expose for GPU + Real d_fac_cond = m_fac_cond; + Real rdOcp = m_rdOcp; + + // get the temperature, dentisy, theta, qt and qc from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const auto& tbx = mfi.tilebox(); + + auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); + + //------- Evaporation/condensation + Real qsat; + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat); + + // There is enough moisutre to drive to equilibrium + if ((qv_array(i,j,k)+qc_array(i,j,k)) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , + d_fac_cond, tabs_array, pres_array, + qv_array , qc_array ); + + // Update theta (constant pressure) + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // + // We cannot blindly relax to qsat, but we can convert qc/qi -> qv. + // The concept here is that if we put all the moisture into qv and modify + // the temperature, we can then check if qv > qsat occurs (for final T/P/qv). + // If the reduction in T/qsat and increase in qv does trigger the + // aforementioned condition, we can do Newton iteration to drive qv = qsat. + // + } else { + // Changes in each component + Real delta_qc = qc_array(i,j,k); + + // Partition the change in non-precipitating q + qv_array(i,j,k) += qc_array(i,j,k); + qc_array(i,j,k) = 0.0; + + // Update temperature (endothermic since we evap/sublime) + tabs_array(i,j,k) -= d_fac_cond * delta_qc; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // Verify assumption that qv > qsat does not occur + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat); + if (qv_array(i,j,k) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , + d_fac_cond , tabs_array, pres_array, + qv_array , qc_array ); + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + } + } + }); + } +} diff --git a/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp new file mode 100644 index 000000000..b45d619ef --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp @@ -0,0 +1,37 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc + */ +void SatAdj::Copy_Micro_to_State (MultiFab& cons) +{ + // Get the temperature, density, theta, qt and qp from input + for (MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& tbx = mfi.tilebox(); + + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + + // get potential total density, temperature, qt, qp + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); +} + diff --git a/Source/Microphysics/SatAdj/Make.package b/Source/Microphysics/SatAdj/Make.package new file mode 100644 index 000000000..5b9ff43e4 --- /dev/null +++ b/Source/Microphysics/SatAdj/Make.package @@ -0,0 +1,4 @@ +CEXE_sources += ERF_InitSatAdj.cpp +CEXE_sources += ERF_UpdateSatAdj.cpp +CEXE_sources += ERF_SatAdj.cpp +CEXE_headers += ERF_SatAdj.H