From 6be0a3c6b4baf3e199ba8776dd5564348026a2a2 Mon Sep 17 00:00:00 2001 From: Kyle Date: Wed, 9 Sep 2026 16:47:40 -0700 Subject: [PATCH] Support reference pressure in anelastic microphysics --- Docs/sphinx_doc/SuperDroplets.rst | 18 + Docs/sphinx_doc/theory/Microphysics.rst | 58 ++ Source/DataStructs/ERF_DataStruct.H | 20 + .../Microphysics/ERF_EulerianMicrophysics.H | 1 + .../Microphysics/Kessler/ERF_InitKessler.cpp | 34 +- Source/Microphysics/Kessler/ERF_Kessler.H | 52 +- .../Morrison/ERF_InitMorrison.cpp | 59 +- Source/Microphysics/Morrison/ERF_Morrison.H | 60 +- Source/Microphysics/Null/ERF_NullMoist.H | 29 +- Source/Microphysics/SAM/ERF_InitSAM.cpp | 102 ++- Source/Microphysics/SAM/ERF_SAM.H | 17 +- Source/Microphysics/SAM/ERF_SAMUtils.H | 108 ++- Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp | 5 +- Source/Microphysics/SatAdj/ERF_SatAdj.H | 1 + .../ERF_SuperDropletsMoist.H | 6 +- Source/Microphysics/WDM6/ERF_AdvanceWDM6.cpp | 27 +- Source/Microphysics/WDM6/ERF_InitWDM6.cpp | 36 +- Source/Microphysics/WDM6/ERF_WDM6.H | 58 +- Source/Microphysics/WSM6/ERF_InitWSM6.cpp | 38 +- Source/Microphysics/WSM6/ERF_UpdateWSM6.cpp | 15 +- Source/Microphysics/WSM6/ERF_WSM6.H | 90 ++- Source/Utils/ERF_MicrophysicsUtils.H | 36 + Tests/Unit/CMakeLists.txt | 3 + .../ERF_GTestSatAdjTemperatureDiagnostics.cpp | 2 + .../ERF_GTestAnelasticMicrophysicsWiring.cpp | 659 ++++++++++++++++++ 25 files changed, 1388 insertions(+), 146 deletions(-) create mode 100644 Tests/Unit/Utils/Microphysics/ERF_GTestAnelasticMicrophysicsWiring.cpp diff --git a/Docs/sphinx_doc/SuperDroplets.rst b/Docs/sphinx_doc/SuperDroplets.rst index e3dbaa7b3c..7e13f5e18c 100644 --- a/Docs/sphinx_doc/SuperDroplets.rst +++ b/Docs/sphinx_doc/SuperDroplets.rst @@ -24,6 +24,24 @@ the same way here. it is read. The physical guidance -- which kernel or which time integrator to prefer for a given problem -- is not covered here. +.. _superdroplets-anelastic-compatibility: + +Anelastic compatibility +----------------------- + +.. warning:: + + ``SuperDroplets`` cannot currently be used in a simulation in which any AMR + level has ``erf.anelastic = 1``. ERF detects this combination during input + validation and exits with an error. + + Use ``SuperDroplets`` with compressible dynamics, or choose a supported + Eulerian microphysics scheme for an anelastic simulation. + +The current Super-Droplet thermodynamic coupling uses the compressible equation +of state and does not consume the hydrostatic base-state pressure used by +anelastic Eulerian microphysics. + Model configuration ------------------- diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index 3d6e6c3677..92e704292b 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -60,6 +60,64 @@ Model overview and transported quantities in ERF number concentration fields: :math:`n_n` (CCN number concentration), :math:`n_c` (cloud droplet number concentration), and :math:`n_r` (rain drop number concentration). +.. _anelastic-microphysics-thermodynamics: + +Anelastic microphysics thermodynamics +------------------------------------- + +ERF supplies the thermodynamic pressure needed by supported Eulerian +microphysics schemes. No additional pressure field or scheme-specific pressure +unit conversion is required in the input file. + +On each anelastic AMR level (``erf.anelastic = 1``), pressure-dependent +microphysics uses the hydrostatic base-state pressure at each cell. Let +:math:`p_{\mathrm{base}}` denote this pressure and :math:`p_{\mathrm{ref}}` +denote ERF's fixed reference pressure used in the definition of potential +temperature. ERF diagnoses + +.. math:: + + \theta = \frac{\rho\theta}{\rho}, \qquad + p_{\mathrm{micro}} = p_{\mathrm{base}}, + +and + +.. math:: + + T = \theta + \left(\frac{p_{\mathrm{base}}}{p_{\mathrm{ref}}}\right)^{R_d/c_p}. + +The ratio :math:`R_d/c_p` uses ``erf.c_p`` for this anelastic +pressure--temperature conversion. Pressure-dependent thermodynamic +calculations use the same hydrostatic base-state pressure consistently through +the microphysics source step. Schemes that convert between temperature and +potential temperature use that same pressure. + +On a compressible AMR level, the existing behavior is unchanged: ERF diagnoses +local pressure and temperature from the prognostic state with the compressible +equation of state. In a mixed hierarchy, this choice is made independently on +each AMR level. + +The reference-pressure treatment applies to the Kessler family +(``Kessler`` and ``Kessler_NoRain``), the SAM family (``SAM``, ``SAM_NoIce``, +and ``SAM_NoPrecip_NoIce``), the Morrison family (``Morrison`` and +``Morrison_NoIce``), ``WSM6``, and ``WDM6``. ``SatAdj`` already uses the +anelastic base-state pressure. ``MoistNoCondensation`` does not have a +pressure-dependent condensation source that requires this treatment. + +ERF performs scheme-specific pressure-unit conversion internally: Kessler and +SAM use hPa internally, while Morrison, WSM6, and WDM6 use Pa. Users should +specify only the usual model and dynamics options. For example: + +.. code-block:: text + + erf.anelastic = 1 + erf.moisture_model = WSM6 + +``erf.anelastic`` may also be specified per AMR level. ``SuperDroplets`` is not +supported when any AMR level is anelastic; see +:ref:`superdroplets-anelastic-compatibility`. + Surface precipitation accumulations ----------------------------------- diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index f931e05507..8d89132609 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -110,6 +110,19 @@ AMREX_ENUM(MoistureType, SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, Morrison, Morrison_NoIce, WSM6, WDM6, SuperDroplets, MoistNoCondensation, None ); +/** + * @brief Whether an anelastic SuperDroplets configuration is unsupported. + * + * The current SuperDroplets implementation diagnoses thermodynamics through + * the compressible equation of state and has no BaseState pressure input. + */ +inline bool anelastic_superdroplets_configuration_invalid ( + const MoistureType moisture_type, + const bool any_anelastic) noexcept +{ + return any_anelastic && moisture_type == MoistureType::SuperDroplets; +} + /** * @brief Wind-farm model. */ @@ -906,6 +919,13 @@ struct SolverChoice { } } + if (anelastic_superdroplets_configuration_invalid(moisture_type, any_anelastic)) { + amrex::Error("erf.moisture_model = SuperDroplets is not supported with " + "erf.anelastic = 1: the Super-Droplet path currently uses " + "the compressible equation of state. Use compressible " + "dynamics or a supported Eulerian moisture model."); + } + // Want to have different immersed forcing defaults depending on anelastic or fully compressible. // We should have different starting values because the stability of the method is dependent on dt. // dt is coarser for anelastic --> need to loosen stiffness of immersed forcing. diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index 320ed58826..7f62afec13 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -61,6 +61,7 @@ public: void Define (const int& lev, /*!< AMR level */ SolverChoice& sc /*!< Solver choice object */) override { + m_moist_model[lev]->SetCurrentLevel(lev); m_moist_model[lev]->Define(sc); } diff --git a/Source/Microphysics/Kessler/ERF_InitKessler.cpp b/Source/Microphysics/Kessler/ERF_InitKessler.cpp index a7a7f91295..e8d790bc8f 100644 --- a/Source/Microphysics/Kessler/ERF_InitKessler.cpp +++ b/Source/Microphysics/Kessler/ERF_InitKessler.cpp @@ -62,6 +62,16 @@ void Kessler::Init (const MultiFab& cons_in, */ void Kessler::Copy_State_to_Micro (const MultiFab& cons_in) { + Copy_State_to_Micro(cons_in, nullptr); +} + +void Kessler::Copy_State_to_Micro (const MultiFab& cons_in, + const MultiFab* base_state) +{ + assert_base_state_available(base_state); + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; + // Get the temperature, density, theta, qt and qp from input for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { const auto& box3d = mfi.growntilebox(); @@ -78,23 +88,21 @@ void Kessler::Copy_State_to_Micro (const MultiFab& cons_in) auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi); + const auto base_array = base_state ? base_state->const_array(mfi) : Array4{}; + const Real rdOcp = m_rdOcp; - // getPgivenRTh returns Pa. Kessler stores pressure in mbar / hPa for the - // qsat helper path, so convert here after forming temperature and density. ParallelFor( box3d, [=] 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); - qp_array(i,j,k) = states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k); - - 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)) * Real(0.01); + kessler_copy_state_to_micro_cell( + states_array, base_array, rho_array, theta_array, qv_array, + qc_array, qp_array, qt_array, tabs_array, pres_array, rdOcp, + use_anelastic_reference_pressure, i, j, k); }); } } +void Kessler::Update_Micro_Vars (MultiFab& cons_in, + const MultiFab* base_state) +{ + Copy_State_to_Micro(cons_in, base_state); +} diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index 485a9fcaf6..05020e728c 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -34,7 +34,44 @@ namespace MicVar_Kess { // derived vars rain_accum, NumVars - }; + }; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void kessler_copy_state_to_micro_cell ( + const amrex::Array4& states, + const amrex::Array4& base, + const amrex::Array4& rho, + const amrex::Array4& theta, + const amrex::Array4& qv, + const amrex::Array4& qc, + const amrex::Array4& qp, + const amrex::Array4& qt, + const amrex::Array4& tabs, + const amrex::Array4& pres, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const int i, const int j, const int k) noexcept +{ + const amrex::Real rho_value = states(i,j,k,Rho_comp); + const amrex::Real rho_theta = states(i,j,k,RhoTheta_comp); + const amrex::Real qv_value = states(i,j,k,RhoQ1_comp) / rho_value; + + rho(i,j,k) = rho_value; + theta(i,j,k) = rho_theta / rho_value; + qv(i,j,k) = qv_value; + qc(i,j,k) = states(i,j,k,RhoQ2_comp) / rho_value; + qp(i,j,k) = states(i,j,k,RhoQ3_comp) / rho_value; + qt(i,j,k) = qv(i,j,k) + qc(i,j,k); + + const amrex::Real p0 = use_anelastic_reference_pressure + ? base(i,j,k,BaseState::p0_comp) : amrex::Real(0.0); + const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state( + rho_value, rho_theta, qv_value, rdOcp, + use_anelastic_reference_pressure, p0); + tabs(i,j,k) = thermo.temperature; + // Kessler's saturation helpers use mbar / hPa. + pres(i,j,k) = amrex::Real(0.01) * thermo.pressure_pa; } class Kessler : public NullMoist { @@ -59,7 +96,9 @@ public: // Use one latent-over-cp factor for both the saturation solve and theta // update. L_v comes from ERF_Constants.H; c_p comes from SolverChoice. m_fac_cond = L_v / sc.c_p; + m_rdOcp = sc.rdOcp; m_do_cond = (!sc.uses_shoc_family()); + set_anelastic_reference_pressure_mode(sc); } // init @@ -82,12 +121,18 @@ public: void Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state); + + void + Update_Micro_Vars (amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state) override; + // Copy state into micro vars void Copy_Micro_to_State (amrex::MultiFab& cons_in) override; - using NullMoist::Update_Micro_Vars; - // update micro vars void Update_Micro_Vars (amrex::MultiFab& cons_in) override @@ -183,6 +228,7 @@ private: // constants amrex::Real m_fac_cond; + amrex::Real m_rdOcp{RdoCp}; bool m_do_cond; MoistureType m_moisture_type = MoistureType::None; diff --git a/Source/Microphysics/Morrison/ERF_InitMorrison.cpp b/Source/Microphysics/Morrison/ERF_InitMorrison.cpp index b51ca8abd4..0f748c14b7 100644 --- a/Source/Microphysics/Morrison/ERF_InitMorrison.cpp +++ b/Source/Microphysics/Morrison/ERF_InitMorrison.cpp @@ -72,6 +72,17 @@ Morrison::Init (const MultiFab& cons_in, void Morrison::Copy_State_to_Micro (const MultiFab& cons_in) { + Copy_State_to_Micro(cons_in, nullptr); +} + +void +Morrison::Copy_State_to_Micro (const MultiFab& cons_in, + const MultiFab* base_state) +{ + assert_base_state_available(base_state); + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; + // Get the temperature, density, theta, qt and qp from input for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { const auto& box3d = mfi.growntilebox(); @@ -101,38 +112,32 @@ Morrison::Copy_State_to_Micro (const MultiFab& cons_in) auto theta_array = mic_fab_vars[MicVar_Morr::theta]->array(mfi); auto tabs_array = mic_fab_vars[MicVar_Morr::tabs]->array(mfi); auto pres_array = mic_fab_vars[MicVar_Morr::pres]->array(mfi); + const auto base_array = base_state ? base_state->const_array(mfi) : Array4{}; + const Real rdOcp = m_rdOcp; // Get pressure, theta, temperature, density, and qt, qp ParallelFor( box3d, [=] 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) = std::max(Real(0),states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp)); - qc_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp)); - qi_array(i,j,k) = std::max(Real(0),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); - qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); - - qpr_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ4_comp)/states_array(i,j,k,Rho_comp)); - qps_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ5_comp)/states_array(i,j,k,Rho_comp)); - qpg_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ6_comp)/states_array(i,j,k,Rho_comp)); - - qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - - nc_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ7_comp) /states_array(i,j,k,Rho_comp)); - ni_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ8_comp) /states_array(i,j,k,Rho_comp)); - nr_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ9_comp) /states_array(i,j,k,Rho_comp)); - ns_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ10_comp)/states_array(i,j,k,Rho_comp)); - ng_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ11_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)); - - // NOTE: the Morrison Fortran version uses Pa not hPa so we don't divideby 100! - pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); // * Real(0.01); + morrison_copy_state_to_micro_cell( + states_array, base_array, rho_array, theta_array, qv_array, + qc_array, qi_array, qn_array, qt_array, qpr_array, qps_array, + qpg_array, qp_array, tabs_array, pres_array, rdOcp, + use_anelastic_reference_pressure, i, j, k); + + const Real rho = rho_array(i,j,k); + + nc_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ7_comp) / rho); + ni_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ8_comp) / rho); + nr_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ9_comp) / rho); + ns_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ10_comp) / rho); + ng_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ11_comp) / rho); }); } } +void +Morrison::Update_Micro_Vars (MultiFab& cons_in, + const MultiFab* base_state) +{ + Copy_State_to_Micro(cons_in, base_state); +} diff --git a/Source/Microphysics/Morrison/ERF_Morrison.H b/Source/Microphysics/Morrison/ERF_Morrison.H index 43e0a5ae9b..2e3891e369 100644 --- a/Source/Microphysics/Morrison/ERF_Morrison.H +++ b/Source/Microphysics/Morrison/ERF_Morrison.H @@ -53,7 +53,52 @@ namespace MicVar_Morr { graup_accum, omega, NumVars - }; + }; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void morrison_copy_state_to_micro_cell ( + const amrex::Array4& states, + const amrex::Array4& base, + const amrex::Array4& rho, + const amrex::Array4& theta, + const amrex::Array4& qv, + const amrex::Array4& qc, + const amrex::Array4& qi, + const amrex::Array4& qn, + const amrex::Array4& qt, + const amrex::Array4& qpr, + const amrex::Array4& qps, + const amrex::Array4& qpg, + const amrex::Array4& qp, + const amrex::Array4& tabs, + const amrex::Array4& pres, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const int i, const int j, const int k) noexcept +{ + const amrex::Real rho_value = states(i,j,k,Rho_comp); + const amrex::Real rho_theta = states(i,j,k,RhoTheta_comp); + rho(i,j,k) = rho_value; + theta(i,j,k) = rho_theta / rho_value; + qv(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ1_comp) / rho_value); + qc(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ2_comp) / rho_value); + qi(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ3_comp) / rho_value); + qn(i,j,k) = qc(i,j,k) + qi(i,j,k); + qt(i,j,k) = qv(i,j,k) + qn(i,j,k); + qpr(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ4_comp) / rho_value); + qps(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ5_comp) / rho_value); + qpg(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ6_comp) / rho_value); + qp(i,j,k) = qpr(i,j,k) + qps(i,j,k) + qpg(i,j,k); + + const amrex::Real p0 = use_anelastic_reference_pressure + ? base(i,j,k,BaseState::p0_comp) : amrex::Real(0.0); + const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state( + rho_value, rho_theta, qv(i,j,k), rdOcp, + use_anelastic_reference_pressure, p0); + tabs(i,j,k) = thermo.temperature; + // Morrison's native pressure unit is Pa. + pres(i,j,k) = thermo.pressure_pa; } class Morrison : public NullMoist { @@ -74,6 +119,7 @@ public: m_moisture_type = sc.moisture_type; m_rdOcp = sc.rdOcp; m_do_cond = (!sc.uses_shoc_family()); + set_anelastic_reference_pressure_mode(sc); } // init @@ -96,12 +142,18 @@ public: void Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state); + + void + Update_Micro_Vars (amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state) override; + // Copy state into micro vars void Copy_Micro_to_State (amrex::MultiFab& cons_in) override; - using NullMoist::Update_Micro_Vars; - void Update_Micro_Vars (amrex::MultiFab& cons_in) override { @@ -198,7 +250,7 @@ private: amrex::Geometry m_geom; // constants - amrex::Real m_rdOcp; + amrex::Real m_rdOcp{RdoCp}; bool m_do_cond; MoistureType m_moisture_type = MoistureType::None; diff --git a/Source/Microphysics/Null/ERF_NullMoist.H b/Source/Microphysics/Null/ERF_NullMoist.H index c11c8da5e7..d92379bfdf 100644 --- a/Source/Microphysics/Null/ERF_NullMoist.H +++ b/Source/Microphysics/Null/ERF_NullMoist.H @@ -111,7 +111,11 @@ public: virtual void - SetCurrentLevel (const int&) {} + SetCurrentLevel (const int& lev) + { + AMREX_ALWAYS_ASSERT(lev >= 0); + m_level = lev; + } virtual void @@ -146,6 +150,29 @@ public: // chooses to implement width-sensitive behavior. Set_RealWidth (const int /*real_width*/) { } +protected: + void + set_anelastic_reference_pressure_mode (const SolverChoice& sc) + { + if (sc.anelastic.empty()) { + m_use_anelastic_reference_pressure = false; + } else { + AMREX_ALWAYS_ASSERT( + m_level < static_cast(sc.anelastic.size())); + m_use_anelastic_reference_pressure = (sc.anelastic[m_level] == 1); + } + } + + void + assert_base_state_available (const amrex::MultiFab* base_state) const + { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !m_use_anelastic_reference_pressure || base_state != nullptr, + "Anelastic microphysics requires the hydrostatic base state at copy-in"); + } + + int m_level{0}; + bool m_use_anelastic_reference_pressure{false}; private: int m_qmoist_size = 0; diff --git a/Source/Microphysics/SAM/ERF_InitSAM.cpp b/Source/Microphysics/SAM/ERF_InitSAM.cpp index 5859d00c36..1a36febf0d 100644 --- a/Source/Microphysics/SAM/ERF_InitSAM.cpp +++ b/Source/Microphysics/SAM/ERF_InitSAM.cpp @@ -93,6 +93,17 @@ SAM::Init (const MultiFab& cons_in, void SAM::Copy_State_to_Micro (const MultiFab& cons_in) { + Copy_State_to_Micro(cons_in, nullptr); +} + +void +SAM::Copy_State_to_Micro (const MultiFab& cons_in, + const MultiFab* base_state) +{ + assert_base_state_available(base_state); + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; + // Get the temperature, density, theta, qt and qp from input for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) { const auto& box3d = mfi.growntilebox(); @@ -116,38 +127,30 @@ SAM::Copy_State_to_Micro (const MultiFab& cons_in) auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi); + const auto base_array = base_state ? base_state->const_array(mfi) : Array4{}; + const Real rdOcp = m_rdOcp; // Get pressure, theta, temperature, density, and qt, qp ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - const SAMPrimitiveCell primitive = - sam_cons_to_primitive(states_array(i,j,k,Rho_comp), - states_array(i,j,k,RhoTheta_comp), - states_array(i,j,k,RhoQ1_comp), - states_array(i,j,k,RhoQ2_comp), - states_array(i,j,k,RhoQ3_comp), - states_array(i,j,k,RhoQ4_comp), - states_array(i,j,k,RhoQ5_comp), - states_array(i,j,k,RhoQ6_comp)); - rho_array(i,j,k) = primitive.rho; - theta_array(i,j,k) = primitive.theta; - qv_array(i,j,k) = primitive.qv; - qc_array(i,j,k) = primitive.qcl; - qi_array(i,j,k) = primitive.qci; - qn_array(i,j,k) = primitive.qn; - qt_array(i,j,k) = primitive.qt; - qpr_array(i,j,k) = primitive.qpr; - qps_array(i,j,k) = primitive.qps; - qpg_array(i,j,k) = primitive.qpg; - qp_array(i,j,k) = primitive.qp; - tabs_array(i,j,k) = primitive.tabs; - pres_array(i,j,k) = primitive.pres_mbar; + sam_copy_state_to_micro_cell( + states_array, base_array, rho_array, theta_array, qv_array, + qc_array, qi_array, qn_array, qt_array, qpr_array, qps_array, + qpg_array, qp_array, tabs_array, pres_array, rdOcp, + use_anelastic_reference_pressure, i, j, k); }); } } +void +SAM::Update_Micro_Vars (MultiFab& cons_in, + const MultiFab* base_state) +{ + Copy_State_to_Micro(cons_in, base_state); + Compute_Coefficients(m_use_anelastic_reference_pressure); +} -void SAM::Compute_Coefficients () +void SAM::Compute_Coefficients (const bool use_anelastic_reference_pressure) { auto accrrc_t = accrrc.table(); auto accrsi_t = accrsi.table(); @@ -178,40 +181,65 @@ void SAM::Compute_Coefficients () PlaneAverage rho_ave(mic_fab_vars[MicVar::rho].get(), m_geom, m_axis); PlaneAverage theta_ave(mic_fab_vars[MicVar::theta].get(), m_geom, m_axis); PlaneAverage qv_ave(mic_fab_vars[MicVar::qv].get(), m_geom, m_axis); + PlaneAverage pres_ave(mic_fab_vars[MicVar::pres].get(), m_geom, m_axis); + PlaneAverage tabs_ave(mic_fab_vars[MicVar::tabs].get(), m_geom, m_axis); rho_ave.compute_averages(ZDir(), rho_ave.field()); - theta_ave.compute_averages(ZDir(), theta_ave.field()); - qv_ave.compute_averages(ZDir(), qv_ave.field()); + if (use_anelastic_reference_pressure) { + pres_ave.compute_averages(ZDir(), pres_ave.field()); + tabs_ave.compute_averages(ZDir(), tabs_ave.field()); + } else { + theta_ave.compute_averages(ZDir(), theta_ave.field()); + qv_ave.compute_averages(ZDir(), qv_ave.field()); + } // get host variable rho, and rhotheta int ncell = rho_ave.ncell_line(); Gpu::HostVector rho_h(ncell), theta_h(ncell), qv_h(ncell); + Gpu::HostVector pres_h(ncell), tabs_h(ncell); rho_ave.line_average(0, rho_h); - theta_ave.line_average(0, theta_h); - qv_ave.line_average(0, qv_h); + if (use_anelastic_reference_pressure) { + pres_ave.line_average(0, pres_h); + tabs_ave.line_average(0, tabs_h); + } else { + theta_ave.line_average(0, theta_h); + qv_ave.line_average(0, qv_h); + } // copy data to device Gpu::DeviceVector rho_d(ncell), theta_d(ncell), qv_d(ncell); + Gpu::DeviceVector pres_d(ncell), tabs_d(ncell); Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); + if (use_anelastic_reference_pressure) { + Gpu::copyAsync(Gpu::hostToDevice, pres_h.begin(), pres_h.end(), pres_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, tabs_h.begin(), tabs_h.end(), tabs_d.begin()); + } else { + Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); + } Gpu::streamSynchronize(); Real* rho_dptr = rho_d.data(); Real* theta_dptr = theta_d.data(); Real* qv_dptr = qv_d.data(); + Real* pres_dptr = pres_d.data(); + Real* tabs_dptr = tabs_d.data(); ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - Real RhoTheta = rho_dptr[k]*theta_dptr[k]; - Real pressure = getPgivenRTh(RhoTheta, qv_dptr[k]); rho1d_t(k) = rho_dptr[k]; - pres1d_t(k) = sam_pa_to_mbar(pressure); - // NOTE: Limit the temperature to the melting point of ice to avoid a divide by - // 0 condition when computing the cold evaporation coefficients. This should - // not affect results since evaporation requires snow/graupel to be present - // and thus TCopy_State_to_Micro(cons_in); - this->Compute_Coefficients(); + this->Compute_Coefficients(m_use_anelastic_reference_pressure); } void @@ -156,7 +163,7 @@ public: } void - Compute_Coefficients (); + Compute_Coefficients (bool use_anelastic_reference_pressure = false); SAMCoefficientRow CoefficientRowAt (int k) const; @@ -361,7 +368,7 @@ private: amrex::Real m_fac_cond; amrex::Real m_fac_fus; amrex::Real m_fac_sub; - amrex::Real m_rdOcp; + amrex::Real m_rdOcp{RdoCp}; bool m_do_cond; MoistureType m_moisture_type = MoistureType::None; diff --git a/Source/Microphysics/SAM/ERF_SAMUtils.H b/Source/Microphysics/SAM/ERF_SAMUtils.H index 2c6018e2cf..1a8aacd92d 100644 --- a/Source/Microphysics/SAM/ERF_SAMUtils.H +++ b/Source/Microphysics/SAM/ERF_SAMUtils.H @@ -178,6 +178,108 @@ amrex::Real sam_theta_from_stored_mbar_converted_to_pa (const amrex::Real& tabs, return getThgivenTandP(tabs, sam_mbar_to_pa(pres_mbar), rdOcp); } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +SAMPrimitiveCell sam_cons_to_primitive_with_thermo ( + const amrex::Real& rho, + const amrex::Real& rho_theta, + const amrex::Real& rho_qv, + const amrex::Real& rho_qcl, + const amrex::Real& rho_qci, + const amrex::Real& rho_qpr, + const amrex::Real& rho_qps, + const amrex::Real& rho_qpg, + const MicrophysicsThermoState& thermo) noexcept +{ + // The pressure and temperature are supplied by ERF's shared thermodynamic + // contract. SAM converts pressure to its established mbar storage here. + SAMPrimitiveCell result{}; + result.rho = rho; + result.theta = rho_theta / rho; + result.qv = amrex::max(amrex::Real(0.0), rho_qv / rho); + result.qcl = amrex::max(amrex::Real(0.0), rho_qcl / rho); + result.qci = amrex::max(amrex::Real(0.0), rho_qci / rho); + result.qn = result.qcl + result.qci; + result.qt = result.qv + result.qn; + result.qpr = amrex::max(amrex::Real(0.0), rho_qpr / rho); + result.qps = amrex::max(amrex::Real(0.0), rho_qps / rho); + result.qpg = amrex::max(amrex::Real(0.0), rho_qpg / rho); + result.qp = result.qpr + result.qps + result.qpg; + result.tabs = thermo.temperature; + result.pres_mbar = sam_pa_to_mbar(thermo.pressure_pa); + return result; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +SAMPrimitiveCell sam_cons_to_primitive_with_base_state ( + const amrex::Real& rho, + const amrex::Real& rho_theta, + const amrex::Real& rho_qv, + const amrex::Real& rho_qcl, + const amrex::Real& rho_qci, + const amrex::Real& rho_qpr, + const amrex::Real& rho_qps, + const amrex::Real& rho_qpg, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const amrex::Real p0) noexcept +{ + const amrex::Real qv = amrex::max(amrex::Real(0.0), rho_qv / rho); + const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state( + rho, rho_theta, qv, rdOcp, + use_anelastic_reference_pressure, p0); + return sam_cons_to_primitive_with_thermo( + rho, rho_theta, rho_qv, rho_qcl, rho_qci, rho_qpr, rho_qps, rho_qpg, + thermo); +} + +// Array-level copy-in used by production and by the wiring test. This keeps +// state-component selection and SAM's primitive working arrays together. +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +SAMPrimitiveCell sam_copy_state_to_micro_cell ( + const amrex::Array4& states, + const amrex::Array4& base, + const amrex::Array4& rho, + const amrex::Array4& theta, + const amrex::Array4& qv, + const amrex::Array4& qc, + const amrex::Array4& qi, + const amrex::Array4& qn, + const amrex::Array4& qt, + const amrex::Array4& qpr, + const amrex::Array4& qps, + const amrex::Array4& qpg, + const amrex::Array4& qp, + const amrex::Array4& tabs, + const amrex::Array4& pres, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const int i, const int j, const int k) noexcept +{ + const amrex::Real rho_value = states(i,j,k,Rho_comp); + const amrex::Real p0 = use_anelastic_reference_pressure + ? base(i,j,k,BaseState::p0_comp) : amrex::Real(0.0); + const SAMPrimitiveCell primitive = sam_cons_to_primitive_with_base_state( + rho_value, states(i,j,k,RhoTheta_comp), + states(i,j,k,RhoQ1_comp), states(i,j,k,RhoQ2_comp), + states(i,j,k,RhoQ3_comp), states(i,j,k,RhoQ4_comp), + states(i,j,k,RhoQ5_comp), states(i,j,k,RhoQ6_comp), rdOcp, + use_anelastic_reference_pressure, p0); + rho(i,j,k) = primitive.rho; + theta(i,j,k) = primitive.theta; + qv(i,j,k) = primitive.qv; + qc(i,j,k) = primitive.qcl; + qi(i,j,k) = primitive.qci; + qn(i,j,k) = primitive.qn; + qt(i,j,k) = primitive.qt; + qpr(i,j,k) = primitive.qpr; + qps(i,j,k) = primitive.qps; + qpg(i,j,k) = primitive.qpg; + qp(i,j,k) = primitive.qp; + tabs(i,j,k) = primitive.tabs; + pres(i,j,k) = primitive.pres_mbar; + return primitive; +} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SAMPrimitiveCell sam_cons_to_primitive (const amrex::Real& rho, const amrex::Real& rho_theta, @@ -188,10 +290,8 @@ SAMPrimitiveCell sam_cons_to_primitive (const amrex::Real& rho, const amrex::Real& rho_qps, const amrex::Real& rho_qpg) noexcept { - // Copy-in reconstructs the SAM primitive cell state from the conserved ERF - // state, clips negative moisture defensively, and diagnoses tabs and - // pres_mbar from the incoming rho*theta/rho*q state. That diagnosed mbar - // pressure is the held pressure for the following source updates. + // Keep the established compressible copy-in formula as its own entry + // point so the compressible path retains its original evaluation order. SAMPrimitiveCell result{}; result.rho = rho; result.theta = rho_theta / rho; diff --git a/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp index ce7be92d57..f2e76b17a1 100644 --- a/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp +++ b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp @@ -49,6 +49,9 @@ void SatAdj::Copy_State_to_Micro (const MultiFab& cons_in) void SatAdj::Update_Micro_Vars (MultiFab& cons_in, const MultiFab* base_state) { + assert_base_state_available(base_state); + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; const Real rdOcp = m_rdOcp; for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) { const auto& tbx = mfi.tilebox(); @@ -79,7 +82,7 @@ void SatAdj::Update_Micro_Vars (MultiFab& cons_in, theta_array(i,j,k) = theta; qv_array(i,j,k) = qv; qc_array(i,j,k) = qc; - if (base_state != nullptr) { + if (use_anelastic_reference_pressure) { const Real p0 = base_array(i,j,k,BaseState::p0_comp); tabs_array(i,j,k) = getTgivenPandTh(p0, theta, rdOcp); pres_array(i,j,k) = p0 * Real(0.01); diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H index d012ea6c19..3006db7fc8 100644 --- a/Source/Microphysics/SatAdj/ERF_SatAdj.H +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -64,6 +64,7 @@ public: m_fac_cond = lcond / sc.c_p; m_rdOcp = sc.rdOcp; m_do_cond = (!sc.uses_shoc_family()); + set_anelastic_reference_pressure_mode(sc); } // init diff --git a/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H b/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H index 49d531d748..854e23508d 100644 --- a/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H +++ b/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H @@ -176,7 +176,11 @@ class SuperDropletsMoist : public NullMoistLagrangian { } /*! \brief Set the current AMR level being processed */ - void SetCurrentLevel (const int& a_lev) override { m_current_lev = a_lev; } + void SetCurrentLevel (const int& a_lev) override + { + NullMoist::SetCurrentLevel(a_lev); + m_current_lev = a_lev; + } /*! \brief Get the diagnostics interval */ int getDiagnosticsInterval () const override { return m_diagnostics_iter; } diff --git a/Source/Microphysics/WDM6/ERF_AdvanceWDM6.cpp b/Source/Microphysics/WDM6/ERF_AdvanceWDM6.cpp index b893c46637..de5c867a3f 100644 --- a/Source/Microphysics/WDM6/ERF_AdvanceWDM6.cpp +++ b/Source/Microphysics/WDM6/ERF_AdvanceWDM6.cpp @@ -897,15 +897,16 @@ void WDM6::Advance(const Real& dt_advance, // ERF stores theta (potential temperature), so we must convert back: theta = T / exner // This matches WRF's conversion: th(i,k,j) = t(i,k) / pii(i,k,j) auto const& theta_arr = mic_fab_vars[MicVar_WDM6::theta]->array(mfi); - constexpr Real p0 = 1.e5; // Reference pressure (Pa) - constexpr Real rdOcp = R_d / Cp_d; // R/cp = 0.286 - + const Real configured_rdOcp = m_rdOcp; + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - // Recompute theta from updated temperature - // exner = (p/p0)^(R/cp) - // theta = T / exner = T * (p0/p)^(R/cp) - Real exner = std::pow(p_arr(i,j,k) / p0, rdOcp); - theta_arr(i,j,k) = t_arr(i,j,k) / exner; + // Use the same held pressure that was passed to WDM6. In + // anelastic mode this is the BaseState p0; in compressible mode + // it is the pressure diagnosed from the incoming state. + theta_arr(i,j,k) = wdm6_theta_from_temperature_and_pressure( + t_arr(i,j,k), p_arr(i,j,k), configured_rdOcp, + use_anelastic_reference_pressure); }); // (Tile-based precipitation diagnostics removed - using global diagnostics instead) @@ -3261,11 +3262,13 @@ void WDM6::Advance(const Real& dt_advance, // POST theta at 100/100 levels, and the resulting // divergence 6.625510053 at k=90 matched the step-1 theta error exactly. { - constexpr Real p0_nat = 1.e5; // Reference pressure (Pa) - constexpr Real rdOcp_nat = R_d / Cp_d; // R/cp = 0.286 + const Real configured_rdOcp = m_rdOcp; + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real exner = std::pow(p_arr(i,j,k) / p0_nat, rdOcp_nat); - w1_theta(i,j,k) = t_arr(i,j,k) / exner; + w1_theta(i,j,k) = wdm6_theta_from_temperature_and_pressure( + t_arr(i,j,k), p_arr(i,j,k), configured_rdOcp, + use_anelastic_reference_pressure); }); } diff --git a/Source/Microphysics/WDM6/ERF_InitWDM6.cpp b/Source/Microphysics/WDM6/ERF_InitWDM6.cpp index ec62d37670..02093fb1db 100644 --- a/Source/Microphysics/WDM6/ERF_InitWDM6.cpp +++ b/Source/Microphysics/WDM6/ERF_InitWDM6.cpp @@ -82,6 +82,17 @@ WDM6::Init(const MultiFab& cons_in, void WDM6::Copy_State_to_Micro(const MultiFab& cons_in) { + Copy_State_to_Micro(cons_in, nullptr); +} + +void +WDM6::Copy_State_to_Micro(const MultiFab& cons_in, + const MultiFab* base_state) +{ + assert_base_state_available(base_state); + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; + for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) { // Match Morrison behavior: refresh microphysics ghost zones from state. const auto& box3d = mfi.growntilebox(); @@ -91,7 +102,7 @@ WDM6::Copy_State_to_Micro(const MultiFab& cons_in) auto theta = mic_fab_vars[MicVar_WDM6::theta]->array(mfi); auto tabs = mic_fab_vars[MicVar_WDM6::tabs]->array(mfi); auto pres = mic_fab_vars[MicVar_WDM6::pres]->array(mfi); - + const auto base_array = base_state ? base_state->const_array(mfi) : Array4{}; auto qv = mic_fab_vars[MicVar_WDM6::qv]->array(mfi); auto qc = mic_fab_vars[MicVar_WDM6::qc]->array(mfi); auto qi = mic_fab_vars[MicVar_WDM6::qi]->array(mfi); @@ -104,17 +115,13 @@ WDM6::Copy_State_to_Micro(const MultiFab& cons_in) auto nr = mic_fab_vars[MicVar_WDM6::nr]->array(mfi); const Real ccn0_local = m_ccn0; // CCN concentration in #/m³ + const Real rdOcp = m_rdOcp; ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) { rho(i,j,k) = states(i,j,k,Rho_comp); - theta(i,j,k) = states(i,j,k,RhoTheta_comp) / states(i,j,k,Rho_comp); - - qv(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ1_comp) / states(i,j,k,Rho_comp)); - qc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ2_comp) / states(i,j,k,Rho_comp)); - qi(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ3_comp) / states(i,j,k,Rho_comp)); - qr(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ4_comp) / states(i,j,k,Rho_comp)); - qs(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ5_comp) / states(i,j,k,Rho_comp)); - qg(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ6_comp) / states(i,j,k,Rho_comp)); + wdm6_copy_state_to_micro_cell( + states, base_array, rho, theta, tabs, pres, qv, qc, qi, qr, + qs, qg, rdOcp, use_anelastic_reference_pressure, i, j, k); // Number concentrations nc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ7_comp) / states(i,j,k,Rho_comp)); @@ -143,10 +150,6 @@ WDM6::Copy_State_to_Micro(const MultiFab& cons_in) // activation from nn. Minimums are enforced in Advance() right before physics, // not during state copying. - tabs(i,j,k) = getTgivenRandRTh(states(i,j,k,Rho_comp), - states(i,j,k,RhoTheta_comp), - qv(i,j,k)); - pres(i,j,k) = getPgivenRTh(states(i,j,k,RhoTheta_comp), qv(i,j,k)); }); } @@ -154,6 +157,13 @@ WDM6::Copy_State_to_Micro(const MultiFab& cons_in) // DON'T clear the flag yet - wait until after Copy_Micro_to_State writes nn to state! } +void +WDM6::Update_Micro_Vars(MultiFab& cons_in, + const MultiFab* base_state) +{ + Copy_State_to_Micro(cons_in, base_state); +} + void WDM6::initialize_coeffs() { diff --git a/Source/Microphysics/WDM6/ERF_WDM6.H b/Source/Microphysics/WDM6/ERF_WDM6.H index 9f0843ec2d..870bdfa1d8 100644 --- a/Source/Microphysics/WDM6/ERF_WDM6.H +++ b/Source/Microphysics/WDM6/ERF_WDM6.H @@ -38,6 +38,56 @@ namespace MicVar_WDM6 { }; } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void wdm6_copy_state_to_micro_cell ( + const amrex::Array4& states, + const amrex::Array4& base, + const amrex::Array4& rho, + const amrex::Array4& theta, + const amrex::Array4& tabs, + const amrex::Array4& pres, + const amrex::Array4& qv, + const amrex::Array4& qc, + const amrex::Array4& qi, + const amrex::Array4& qr, + const amrex::Array4& qs, + const amrex::Array4& qg, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const int i, const int j, const int k) noexcept +{ + const amrex::Real rho_value = states(i,j,k,Rho_comp); + const amrex::Real rho_theta = states(i,j,k,RhoTheta_comp); + rho(i,j,k) = rho_value; + theta(i,j,k) = rho_theta / rho_value; + qv(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ1_comp) / rho_value); + qc(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ2_comp) / rho_value); + qi(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ3_comp) / rho_value); + qr(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ4_comp) / rho_value); + qs(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ5_comp) / rho_value); + qg(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ6_comp) / rho_value); + + const amrex::Real p0 = use_anelastic_reference_pressure + ? base(i,j,k,BaseState::p0_comp) : amrex::Real(0.0); + const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state( + rho_value, rho_theta, qv(i,j,k), rdOcp, + use_anelastic_reference_pressure, p0); + tabs(i,j,k) = thermo.temperature; + pres(i,j,k) = thermo.pressure_pa; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real wdm6_theta_from_temperature_and_pressure ( + const amrex::Real temperature, + const amrex::Real pressure, + const amrex::Real configured_rdOcp, + const bool use_anelastic_reference_pressure) noexcept +{ + const amrex::Real theta_rdOcp = use_anelastic_reference_pressure + ? configured_rdOcp : RdoCp; + return getThgivenTandP(temperature, pressure, theta_rdOcp); +} + // Selects how the WDM6 parameter literals are interpreted. See the LITERAL // PRECISION CONTRACT comment on the constant block below for the full // rationale. Namespace scope, not a class member, because a constexpr member @@ -64,6 +114,8 @@ public: void Define (SolverChoice& sc) override { m_moisture_type = sc.moisture_type; + m_rdOcp = sc.rdOcp; + set_anelastic_reference_pressure_mode(sc); m_axis = sc.ave_plane; m_do_cond = (!sc.uses_shoc_family()); } @@ -166,9 +218,12 @@ public: void Set_Lmask (amrex::iMultiFab* lmask) override { m_lmask = lmask; } void Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + void Copy_State_to_Micro (const amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state); void Copy_Micro_to_State (amrex::MultiFab& cons_in) override; - using NullMoist::Update_Micro_Vars; + void Update_Micro_Vars (amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state) override; void Update_Micro_Vars (amrex::MultiFab& cons_in) override { @@ -226,6 +281,7 @@ private: int nlev{0}, zlo{0}, zhi{0}; int m_axis{2}; bool m_do_cond{true}; + amrex::Real m_rdOcp{RdoCp}; MoistureType m_moisture_type{MoistureType::None}; amrex::MultiFab* m_z_phys_nd{nullptr}; diff --git a/Source/Microphysics/WSM6/ERF_InitWSM6.cpp b/Source/Microphysics/WSM6/ERF_InitWSM6.cpp index 19908cf890..01b0ce1613 100644 --- a/Source/Microphysics/WSM6/ERF_InitWSM6.cpp +++ b/Source/Microphysics/WSM6/ERF_InitWSM6.cpp @@ -46,6 +46,17 @@ WSM6::Init(const MultiFab& cons_in, void WSM6::Copy_State_to_Micro(const MultiFab& cons_in) { + Copy_State_to_Micro(cons_in, nullptr); +} + +void +WSM6::Copy_State_to_Micro(const MultiFab& cons_in, + const MultiFab* base_state) +{ + assert_base_state_available(base_state); + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; + for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) { // Match Morrison behavior: refresh microphysics ghost zones from state. // WSM6 Fortran reads the full (ims:ime, jms:jme, kms:kme) slab. @@ -56,33 +67,30 @@ WSM6::Copy_State_to_Micro(const MultiFab& cons_in) auto theta = mic_fab_vars[MicVar_WSM6::theta]->array(mfi); auto tabs = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi); auto pres = mic_fab_vars[MicVar_WSM6::pres]->array(mfi); - + const auto base_array = base_state ? base_state->const_array(mfi) : Array4{}; auto qv = mic_fab_vars[MicVar_WSM6::qv]->array(mfi); auto qc = mic_fab_vars[MicVar_WSM6::qc]->array(mfi); auto qi = mic_fab_vars[MicVar_WSM6::qi]->array(mfi); auto qr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi); auto qs = mic_fab_vars[MicVar_WSM6::qs]->array(mfi); auto qg = mic_fab_vars[MicVar_WSM6::qg]->array(mfi); + const Real rdOcp = m_rdOcp; ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - rho(i,j,k) = states(i,j,k,Rho_comp); - theta(i,j,k) = states(i,j,k,RhoTheta_comp) / states(i,j,k,Rho_comp); - - qv(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ1_comp) / states(i,j,k,Rho_comp)); - qc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ2_comp) / states(i,j,k,Rho_comp)); - qi(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ3_comp) / states(i,j,k,Rho_comp)); - qr(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ4_comp) / states(i,j,k,Rho_comp)); - qs(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ5_comp) / states(i,j,k,Rho_comp)); - qg(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ6_comp) / states(i,j,k,Rho_comp)); - - tabs(i,j,k) = getTgivenRandRTh(states(i,j,k,Rho_comp), - states(i,j,k,RhoTheta_comp), - qv(i,j,k)); - pres(i,j,k) = getPgivenRTh(states(i,j,k,RhoTheta_comp), qv(i,j,k)); + wsm6_copy_state_to_micro_cell( + states, base_array, rho, theta, tabs, pres, qv, qc, qi, qr, + qs, qg, rdOcp, use_anelastic_reference_pressure, i, j, k); }); } } +void +WSM6::Update_Micro_Vars(MultiFab& cons_in, + const MultiFab* base_state) +{ + Copy_State_to_Micro(cons_in, base_state); +} + void WSM6::initialize_coeffs() { diff --git a/Source/Microphysics/WSM6/ERF_UpdateWSM6.cpp b/Source/Microphysics/WSM6/ERF_UpdateWSM6.cpp index 0e2d7957d8..d3ed49f5f8 100644 --- a/Source/Microphysics/WSM6/ERF_UpdateWSM6.cpp +++ b/Source/Microphysics/WSM6/ERF_UpdateWSM6.cpp @@ -13,6 +13,7 @@ WSM6::Copy_Micro_to_State(MultiFab& cons) auto rho = mic_fab_vars[MicVar_WSM6::rho]->array(mfi); auto theta = mic_fab_vars[MicVar_WSM6::theta]->array(mfi); auto tabs = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi); + auto pres = mic_fab_vars[MicVar_WSM6::pres]->const_array(mfi); auto qv = mic_fab_vars[MicVar_WSM6::qv]->array(mfi); auto qc = mic_fab_vars[MicVar_WSM6::qc]->array(mfi); @@ -20,17 +21,15 @@ WSM6::Copy_Micro_to_State(MultiFab& cons) auto qr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi); auto qs = mic_fab_vars[MicVar_WSM6::qs]->array(mfi); auto qg = mic_fab_vars[MicVar_WSM6::qg]->array(mfi); + const Real rdOcp = m_rdOcp; + const bool use_anelastic_reference_pressure = + m_use_anelastic_reference_pressure; ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - theta(i,j,k) = getThgivenRandT(rho(i,j,k), tabs(i,j,k), RdoCp, qv(i,j,k)); - states(i,j,k,RhoTheta_comp) = rho(i,j,k) * theta(i,j,k); - states(i,j,k,RhoQ1_comp) = rho(i,j,k) * amrex::max(Real(0), qv(i,j,k)); - states(i,j,k,RhoQ2_comp) = rho(i,j,k) * amrex::max(Real(0), qc(i,j,k)); - states(i,j,k,RhoQ3_comp) = rho(i,j,k) * amrex::max(Real(0), qi(i,j,k)); - states(i,j,k,RhoQ4_comp) = rho(i,j,k) * amrex::max(Real(0), qr(i,j,k)); - states(i,j,k,RhoQ5_comp) = rho(i,j,k) * amrex::max(Real(0), qs(i,j,k)); - states(i,j,k,RhoQ6_comp) = rho(i,j,k) * amrex::max(Real(0), qg(i,j,k)); + wsm6_copy_micro_to_state_cell( + states, theta, rho, tabs, pres, qv, qc, qi, qr, qs, qg, + use_anelastic_reference_pressure, rdOcp, i, j, k); }); } diff --git a/Source/Microphysics/WSM6/ERF_WSM6.H b/Source/Microphysics/WSM6/ERF_WSM6.H index df50bed9ff..1d763c3f64 100644 --- a/Source/Microphysics/WSM6/ERF_WSM6.H +++ b/Source/Microphysics/WSM6/ERF_WSM6.H @@ -35,6 +35,88 @@ namespace MicVar_WSM6 { }; } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void wsm6_copy_state_to_micro_cell ( + const amrex::Array4& states, + const amrex::Array4& base, + const amrex::Array4& rho, + const amrex::Array4& theta, + const amrex::Array4& tabs, + const amrex::Array4& pres, + const amrex::Array4& qv, + const amrex::Array4& qc, + const amrex::Array4& qi, + const amrex::Array4& qr, + const amrex::Array4& qs, + const amrex::Array4& qg, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const int i, const int j, const int k) noexcept +{ + const amrex::Real rho_value = states(i,j,k,Rho_comp); + const amrex::Real rho_theta = states(i,j,k,RhoTheta_comp); + rho(i,j,k) = rho_value; + theta(i,j,k) = rho_theta / rho_value; + qv(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ1_comp) / rho_value); + qc(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ2_comp) / rho_value); + qi(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ3_comp) / rho_value); + qr(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ4_comp) / rho_value); + qs(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ5_comp) / rho_value); + qg(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ6_comp) / rho_value); + + const amrex::Real p0 = use_anelastic_reference_pressure + ? base(i,j,k,BaseState::p0_comp) : amrex::Real(0.0); + const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state( + rho_value, rho_theta, qv(i,j,k), rdOcp, + use_anelastic_reference_pressure, p0); + tabs(i,j,k) = thermo.temperature; + pres(i,j,k) = thermo.pressure_pa; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real wsm6_theta_from_temperature_and_pressure ( + const amrex::Real rho, + const amrex::Real temperature, + const amrex::Real pressure, + const amrex::Real qv, + const amrex::Real configured_rdOcp, + const bool use_anelastic_reference_pressure) noexcept +{ + return use_anelastic_reference_pressure + ? getThgivenTandP(temperature, pressure, configured_rdOcp) + : getThgivenRandT(rho, temperature, RdoCp, qv); +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void wsm6_copy_micro_to_state_cell ( + const amrex::Array4& states, + const amrex::Array4& theta, + const amrex::Array4& rho, + const amrex::Array4& tabs, + const amrex::Array4& pres, + const amrex::Array4& qv, + const amrex::Array4& qc, + const amrex::Array4& qi, + const amrex::Array4& qr, + const amrex::Array4& qs, + const amrex::Array4& qg, + const bool use_anelastic_reference_pressure, + const amrex::Real rdOcp, + const int i, const int j, const int k) noexcept +{ + const amrex::Real theta_value = wsm6_theta_from_temperature_and_pressure( + rho(i,j,k), tabs(i,j,k), pres(i,j,k), qv(i,j,k), rdOcp, + use_anelastic_reference_pressure); + theta(i,j,k) = theta_value; + states(i,j,k,RhoTheta_comp) = rho(i,j,k) * theta_value; + states(i,j,k,RhoQ1_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qv(i,j,k)); + states(i,j,k,RhoQ2_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qc(i,j,k)); + states(i,j,k,RhoQ3_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qi(i,j,k)); + states(i,j,k,RhoQ4_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qr(i,j,k)); + states(i,j,k,RhoQ5_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qs(i,j,k)); + states(i,j,k,RhoQ6_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qg(i,j,k)); +} + class WSM6 : public NullMoist { using FabPtr = std::shared_ptr; @@ -45,6 +127,8 @@ public: void Define(SolverChoice& sc) override { m_moisture_type = sc.moisture_type; + m_rdOcp = sc.rdOcp; + set_anelastic_reference_pressure_mode(sc); m_axis = sc.ave_plane; m_do_cond = (!sc.uses_shoc_family()); } @@ -84,9 +168,12 @@ public: void Set_dzmin(const amrex::Real dz_min) override { m_dzmin = dz_min; } void Copy_State_to_Micro(const amrex::MultiFab& cons_in) override; + void Copy_State_to_Micro(const amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state); void Copy_Micro_to_State(amrex::MultiFab& cons_in) override; - using NullMoist::Update_Micro_Vars; + void Update_Micro_Vars(amrex::MultiFab& cons_in, + const amrex::MultiFab* base_state) override; void Update_Micro_Vars(amrex::MultiFab& cons_in) override { @@ -143,6 +230,7 @@ private: int nlev{0}, zlo{0}, zhi{0}; int m_axis{2}; bool m_do_cond{true}; + amrex::Real m_rdOcp{RdoCp}; MoistureType m_moisture_type{MoistureType::None}; amrex::MultiFab* m_z_phys_nd{nullptr}; diff --git a/Source/Utils/ERF_MicrophysicsUtils.H b/Source/Utils/ERF_MicrophysicsUtils.H index 8f57fe2627..2c02be3b39 100644 --- a/Source/Utils/ERF_MicrophysicsUtils.H +++ b/Source/Utils/ERF_MicrophysicsUtils.H @@ -12,6 +12,42 @@ #include #include #include +#include "ERF_EOS.H" + +/** + * Thermodynamic state passed from ERF into an Eulerian microphysics scheme. + * Pressure is always in Pa here; schemes with legacy mbar storage convert it + * once at their copy-in boundary. + */ +struct MicrophysicsThermoState { + amrex::Real pressure_pa; + amrex::Real temperature; +}; + +/** + * Diagnose the state used by one Eulerian microphysics cell. + * + * In anelastic mode, the conserved state supplies theta while the hydrostatic + * reference pressure supplied by BaseState supplies the thermodynamic pressure. + * In compressible mode, retain ERF's existing local-EOS diagnosis. + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +MicrophysicsThermoState diagnose_microphysics_thermo_state ( + const amrex::Real rho, + const amrex::Real rho_theta, + const amrex::Real qv, + const amrex::Real rdOcp, + const bool use_anelastic_reference_pressure, + const amrex::Real p0) noexcept +{ + if (use_anelastic_reference_pressure) { + const amrex::Real theta = rho_theta / rho; + return {p0, getTgivenPandTh(p0, theta, rdOcp)}; + } + + return {getPgivenRTh(rho_theta, qv), + getTgivenRandRTh(rho, rho_theta, qv)}; +} // Positive-argument gamma wrapper. std::lgamma returns log(abs(Gamma(x))) for // negative non-integers, so ERF keeps a positive-only contract here. diff --git a/Tests/Unit/CMakeLists.txt b/Tests/Unit/CMakeLists.txt index a66f1e485a..c7408f93c4 100644 --- a/Tests/Unit/CMakeLists.txt +++ b/Tests/Unit/CMakeLists.txt @@ -144,6 +144,9 @@ target_sources(${erf_exe_name} ${CMAKE_CURRENT_SOURCE_DIR}/Utils/Interpolation/ERF_GTestInterpolationWENOKernel.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Utils/Microphysics/ERF_GTestMicrophysicsUtilsScalar.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Utils/Microphysics/ERF_GTestMicrophysicsUtilsKernel.cpp + # Motivation: supported Eulerian microphysics must consume anelastic p0 + # explicitly while preserving the compressible EOS and SAM plane contract. + ${CMAKE_CURRENT_SOURCE_DIR}/Utils/Microphysics/ERF_GTestAnelasticMicrophysicsWiring.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Microphysics/SatAdj/ERF_GTestSatAdjThermo.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Microphysics/SatAdj/ERF_GTestSatAdjCell.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Microphysics/SatAdj/ERF_GTestSatAdjMultiFab.cpp diff --git a/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjTemperatureDiagnostics.cpp b/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjTemperatureDiagnostics.cpp index 5c8e0a4210..5f4229f327 100644 --- a/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjTemperatureDiagnostics.cpp +++ b/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjTemperatureDiagnostics.cpp @@ -268,6 +268,8 @@ TEST(SatAdjTemperatureDiagnostics, AnelasticPressureContextMatchesReference) SatAdj satadj; SolverChoice sc = make_solver_choice(false); + sc.anelastic = {1}; + satadj.SetCurrentLevel(0); satadj.Define(sc); run_and_sync([&]() { std::unique_ptr z_phys_nd; diff --git a/Tests/Unit/Utils/Microphysics/ERF_GTestAnelasticMicrophysicsWiring.cpp b/Tests/Unit/Utils/Microphysics/ERF_GTestAnelasticMicrophysicsWiring.cpp new file mode 100644 index 0000000000..2ce43ba566 --- /dev/null +++ b/Tests/Unit/Utils/Microphysics/ERF_GTestAnelasticMicrophysicsWiring.cpp @@ -0,0 +1,659 @@ +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "ERF_GTestMicrophysicsCommon.H" + +namespace { + +using amrex::Array4; +using amrex::Box; +using amrex::BoxArray; +using amrex::DistributionMapping; +using amrex::Geometry; +using amrex::IntVect; +using amrex::MFIter; +using amrex::MultiFab; +using amrex::Real; +using namespace microphysics_test; + +constexpr Real kRdOcp = RdoCp; + +Geometry make_test_geometry () +{ + const Box domain(IntVect(0, 0, 0), IntVect(1, 0, 0)); + const amrex::RealBox real_box({AMREX_D_DECL(0.0, 0.0, 0.0)}, + {AMREX_D_DECL(2.0, 1.0, 1.0)}); + amrex::Array periodicity{AMREX_D_DECL(0, 0, 0)}; + return Geometry(domain, &real_box, amrex::CoordSys::cartesian, periodicity.data()); +} + +void initialize_test_state (MultiFab& states, MultiFab& base_state) +{ + states.setVal(Real(0.0)); + base_state.setVal(Real(-12345.0)); + + for (MFIter mfi(states); mfi.isValid(); ++mfi) { + const Box box = mfi.growntilebox(); + const auto state = states.array(mfi); + const auto base = base_state.array(mfi); + + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + const int cell = amrex::max(0, amrex::min(1, i)); + const Real rho = cell == 0 ? Real(1.10) : Real(0.90); + const Real theta = cell == 0 ? Real(303.0) : Real(287.0); + const Real p0 = cell == 0 ? Real(90000.0) : Real(82000.0); + const Real qv = cell == 0 ? Real(0.0030) : Real(0.0070); + const Real qc = cell == 0 ? Real(0.0010) : Real(0.0020); + const Real qi = cell == 0 ? Real(0.0004) : Real(0.0008); + const Real qr = cell == 0 ? Real(0.0007) : Real(0.0011); + const Real qs = cell == 0 ? Real(0.0002) : Real(0.0005); + const Real qg = cell == 0 ? Real(0.0003) : Real(0.0006); + + state(i,j,k,Rho_comp) = rho; + state(i,j,k,RhoTheta_comp) = rho * theta; + state(i,j,k,RhoQ1_comp) = rho * qv; + state(i,j,k,RhoQ2_comp) = rho * qc; + state(i,j,k,RhoQ3_comp) = rho * qi; + state(i,j,k,RhoQ4_comp) = rho * qr; + state(i,j,k,RhoQ5_comp) = rho * qs; + state(i,j,k,RhoQ6_comp) = rho * qg; + state(i,j,k,RhoQ7_comp) = rho * Real(2.0); + state(i,j,k,RhoQ8_comp) = Real(0.0); + state(i,j,k,RhoQ9_comp) = rho * Real(3.0); + state(i,j,k,RhoQ10_comp) = rho * Real(4.0); + state(i,j,k,RhoQ11_comp) = rho * Real(5.0); + + // Only p0 is a physical input to the anelastic microphysics path. + // The other base-state fields deliberately remain sentinels. + base(i,j,k,BaseState::p0_comp) = p0; + }); + } + + amrex::Gpu::streamSynchronize(); +} + +struct TestState { + BoxArray boxes; + DistributionMapping dm; + MultiFab states; + MultiFab base_state; + + TestState () + : boxes(Box(IntVect(0, 0, 0), IntVect(1, 0, 0))), + dm(boxes), + states(boxes, dm, RhoQ11_comp + 1, 1), + base_state(boxes, dm, BaseState::num_comps, 1) + { + initialize_test_state(states, base_state); + } +}; + +struct WorkingArrays { + MultiFab rho; + MultiFab theta; + MultiFab qv; + MultiFab qc; + MultiFab qi; + MultiFab qn; + MultiFab qt; + MultiFab qpr; + MultiFab qps; + MultiFab qpg; + MultiFab qp; + MultiFab tabs; + MultiFab pres; + + explicit WorkingArrays (const TestState& input) + : rho(input.boxes, input.dm, 1, 1), + theta(input.boxes, input.dm, 1, 1), + qv(input.boxes, input.dm, 1, 1), + qc(input.boxes, input.dm, 1, 1), + qi(input.boxes, input.dm, 1, 1), + qn(input.boxes, input.dm, 1, 1), + qt(input.boxes, input.dm, 1, 1), + qpr(input.boxes, input.dm, 1, 1), + qps(input.boxes, input.dm, 1, 1), + qpg(input.boxes, input.dm, 1, 1), + qp(input.boxes, input.dm, 1, 1), + tabs(input.boxes, input.dm, 1, 1), + pres(input.boxes, input.dm, 1, 1) + { + set_all(Real(-999.0)); + } + + void set_all (const Real value) + { + rho.setVal(value); + theta.setVal(value); + qv.setVal(value); + qc.setVal(value); + qi.setVal(value); + qn.setVal(value); + qt.setVal(value); + qpr.setVal(value); + qps.setVal(value); + qpg.setVal(value); + qp.setVal(value); + tabs.setVal(value); + pres.setVal(value); + } +}; + +void copy_kessler (const TestState& input, WorkingArrays& work) +{ + for (MFIter mfi(input.states); mfi.isValid(); ++mfi) { + const Box box = mfi.growntilebox(); + const auto states = input.states.const_array(mfi); + const auto base = input.base_state.const_array(mfi); + const auto rho = work.rho.array(mfi); + const auto theta = work.theta.array(mfi); + const auto qv = work.qv.array(mfi); + const auto qc = work.qc.array(mfi); + const auto qp = work.qp.array(mfi); + const auto qt = work.qt.array(mfi); + const auto tabs = work.tabs.array(mfi); + const auto pres = work.pres.array(mfi); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + kessler_copy_state_to_micro_cell( + states, base, rho, theta, qv, qc, qp, qt, tabs, pres, + kRdOcp, true, i, j, k); + }); + } + amrex::Gpu::streamSynchronize(); +} + +void copy_sam (const TestState& input, WorkingArrays& work) +{ + for (MFIter mfi(input.states); mfi.isValid(); ++mfi) { + const Box box = mfi.growntilebox(); + const auto states = input.states.const_array(mfi); + const auto base = input.base_state.const_array(mfi); + const auto rho = work.rho.array(mfi); + const auto theta = work.theta.array(mfi); + const auto qv = work.qv.array(mfi); + const auto qc = work.qc.array(mfi); + const auto qi = work.qi.array(mfi); + const auto qn = work.qn.array(mfi); + const auto qt = work.qt.array(mfi); + const auto qpr = work.qpr.array(mfi); + const auto qps = work.qps.array(mfi); + const auto qpg = work.qpg.array(mfi); + const auto qp = work.qp.array(mfi); + const auto tabs = work.tabs.array(mfi); + const auto pres = work.pres.array(mfi); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + sam_copy_state_to_micro_cell( + states, base, rho, theta, qv, qc, qi, qn, qt, qpr, qps, + qpg, qp, tabs, pres, kRdOcp, true, i, j, k); + }); + } + amrex::Gpu::streamSynchronize(); +} + +void copy_morrison (const TestState& input, WorkingArrays& work) +{ + for (MFIter mfi(input.states); mfi.isValid(); ++mfi) { + const Box box = mfi.growntilebox(); + const auto states = input.states.const_array(mfi); + const auto base = input.base_state.const_array(mfi); + const auto rho = work.rho.array(mfi); + const auto theta = work.theta.array(mfi); + const auto qv = work.qv.array(mfi); + const auto qc = work.qc.array(mfi); + const auto qi = work.qi.array(mfi); + const auto qn = work.qn.array(mfi); + const auto qt = work.qt.array(mfi); + const auto qpr = work.qpr.array(mfi); + const auto qps = work.qps.array(mfi); + const auto qpg = work.qpg.array(mfi); + const auto qp = work.qp.array(mfi); + const auto tabs = work.tabs.array(mfi); + const auto pres = work.pres.array(mfi); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + morrison_copy_state_to_micro_cell( + states, base, rho, theta, qv, qc, qi, qn, qt, qpr, qps, + qpg, qp, tabs, pres, kRdOcp, true, i, j, k); + }); + } + amrex::Gpu::streamSynchronize(); +} + +void copy_wsm6 (const TestState& input, WorkingArrays& work) +{ + for (MFIter mfi(input.states); mfi.isValid(); ++mfi) { + const Box box = mfi.growntilebox(); + const auto states = input.states.const_array(mfi); + const auto base = input.base_state.const_array(mfi); + const auto rho = work.rho.array(mfi); + const auto theta = work.theta.array(mfi); + const auto tabs = work.tabs.array(mfi); + const auto pres = work.pres.array(mfi); + const auto qv = work.qv.array(mfi); + const auto qc = work.qc.array(mfi); + const auto qi = work.qi.array(mfi); + const auto qr = work.qpr.array(mfi); + const auto qs = work.qps.array(mfi); + const auto qg = work.qpg.array(mfi); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + wsm6_copy_state_to_micro_cell( + states, base, rho, theta, tabs, pres, qv, qc, qi, qr, qs, + qg, kRdOcp, true, i, j, k); + }); + } + amrex::Gpu::streamSynchronize(); +} + +void copy_wdm6 (const TestState& input, WorkingArrays& work) +{ + for (MFIter mfi(input.states); mfi.isValid(); ++mfi) { + const Box box = mfi.growntilebox(); + const auto states = input.states.const_array(mfi); + const auto base = input.base_state.const_array(mfi); + const auto rho = work.rho.array(mfi); + const auto theta = work.theta.array(mfi); + const auto tabs = work.tabs.array(mfi); + const auto pres = work.pres.array(mfi); + const auto qv = work.qv.array(mfi); + const auto qc = work.qc.array(mfi); + const auto qi = work.qi.array(mfi); + const auto qr = work.qpr.array(mfi); + const auto qs = work.qps.array(mfi); + const auto qg = work.qpg.array(mfi); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + wdm6_copy_state_to_micro_cell( + states, base, rho, theta, tabs, pres, qv, qc, qi, qr, qs, + qg, kRdOcp, true, i, j, k); + }); + } + amrex::Gpu::streamSynchronize(); +} + +void expect_reference_pressure_diagnosis (const TestState& input, + const WorkingArrays& work, + const Real pressure_scale) +{ + MultiFab errors(input.boxes, input.dm, 4, 0); + errors.setVal(Real(0.0)); + + for (MFIter mfi(input.states); mfi.isValid(); ++mfi) { + const Box box = mfi.validbox(); + const auto states = input.states.const_array(mfi); + const auto base = input.base_state.const_array(mfi); + const auto rho = work.rho.const_array(mfi); + const auto theta = work.theta.const_array(mfi); + const auto tabs = work.tabs.const_array(mfi); + const auto pres = work.pres.const_array(mfi); + const auto error = errors.array(mfi); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + const Real expected_rho = states(i,j,k,Rho_comp); + const Real expected_theta = states(i,j,k,RhoTheta_comp) / expected_rho; + const Real expected_p0 = base(i,j,k,BaseState::p0_comp); + const Real expected_tabs = getTgivenPandTh(expected_p0, expected_theta, kRdOcp); + error(i,j,k,0) = normalized_error(rho(i,j,k), expected_rho, kValueRelTol); + error(i,j,k,1) = normalized_error(theta(i,j,k), expected_theta, kValueRelTol); + error(i,j,k,2) = normalized_error(tabs(i,j,k), expected_tabs, kValueRelTol); + error(i,j,k,3) = normalized_error(pres(i,j,k), expected_p0 * pressure_scale, + kValueRelTol); + }); + } + + amrex::Gpu::streamSynchronize(); + for (int comp = 0; comp < errors.nComp(); ++comp) { + EXPECT_LE(errors.max(comp), Real(20.0)); + } + EXPECT_LT(work.tabs.min(0), work.tabs.max(0)); + EXPECT_LT(work.pres.min(0), work.pres.max(0)); +} + +SolverChoice make_scheme_choice (const MoistureType moisture_type, + const bool anelastic = false) +{ + SolverChoice sc{}; + sc.c_p = Cp_d; + sc.rdOcp = kRdOcp; + sc.moisture_type = moisture_type; + sc.anelastic = {anelastic ? 1 : 0}; + sc.use_eamxx_shoc = false; + sc.use_native_shoc = false; + sc.ave_plane = 2; + return sc; +} + +void expect_coefficients_near (const SAMCoefficientRow& actual, + const SAMCoefficientRow& expected) +{ + constexpr Real coefficient_tol = Real(1.0e-10); + EXPECT_NEAR(actual.accrrc, expected.accrrc, + scaled_tol(actual.accrrc, expected.accrrc, coefficient_tol)); + EXPECT_NEAR(actual.accrsi, expected.accrsi, + scaled_tol(actual.accrsi, expected.accrsi, coefficient_tol)); + EXPECT_NEAR(actual.accrsc, expected.accrsc, + scaled_tol(actual.accrsc, expected.accrsc, coefficient_tol)); + EXPECT_NEAR(actual.coefice, expected.coefice, + scaled_tol(actual.coefice, expected.coefice, coefficient_tol)); + EXPECT_NEAR(actual.evaps1, expected.evaps1, + scaled_tol(actual.evaps1, expected.evaps1, coefficient_tol)); + EXPECT_NEAR(actual.evaps2, expected.evaps2, + scaled_tol(actual.evaps2, expected.evaps2, coefficient_tol)); + EXPECT_NEAR(actual.accrgi, expected.accrgi, + scaled_tol(actual.accrgi, expected.accrgi, coefficient_tol)); + EXPECT_NEAR(actual.accrgc, expected.accrgc, + scaled_tol(actual.accrgc, expected.accrgc, coefficient_tol)); + EXPECT_NEAR(actual.evapg1, expected.evapg1, + scaled_tol(actual.evapg1, expected.evapg1, coefficient_tol)); + EXPECT_NEAR(actual.evapg2, expected.evapg2, + scaled_tol(actual.evapg2, expected.evapg2, coefficient_tol)); + EXPECT_NEAR(actual.evapr1, expected.evapr1, + scaled_tol(actual.evapr1, expected.evapr1, coefficient_tol)); + EXPECT_NEAR(actual.evapr2, expected.evapr2, + scaled_tol(actual.evapr2, expected.evapr2, coefficient_tol)); +} + +} // namespace + +// Motivation: Each supported Eulerian scheme must receive the anelastic +// reference pressure instead of rebuilding pressure from rho and rho*theta. +TEST(AnelasticMicrophysicsWiring, KesslerCopyInUsesReferencePressure) +{ + TestState input; + WorkingArrays work(input); + copy_kessler(input, work); + expect_reference_pressure_diagnosis(input, work, Real(0.01)); +} + +// Motivation: SAM stores pressure in mbar, but its anelastic temperature must +// still be diagnosed from the BaseState pressure in Pa. +TEST(AnelasticMicrophysicsWiring, SAMCopyInUsesReferencePressure) +{ + TestState input; + WorkingArrays work(input); + copy_sam(input, work); + expect_reference_pressure_diagnosis(input, work, Real(0.01)); +} + +// Motivation: Morrison's source kernels consume pressure in Pa and must see +// the same BaseState pressure used to diagnose their absolute temperature. +TEST(AnelasticMicrophysicsWiring, MorrisonCopyInUsesReferencePressure) +{ + TestState input; + WorkingArrays work(input); + copy_morrison(input, work); + expect_reference_pressure_diagnosis(input, work, Real(1.0)); +} + +// Motivation: WSM6 copy-in and copy-out must use the same held pressure so an +// anelastic source update cannot project theta through a rho-based EOS. +TEST(AnelasticMicrophysicsWiring, WSM6CopyInUsesReferencePressure) +{ + TestState input; + WorkingArrays work(input); + copy_wsm6(input, work); + expect_reference_pressure_diagnosis(input, work, Real(1.0)); + + SolverChoice sc = make_scheme_choice(MoistureType::WSM6, true); + WSM6 wsm6; + wsm6.SetCurrentLevel(0); + wsm6.Define(sc); + std::unique_ptr z_phys_nd; + std::unique_ptr detJ_cc; + const Geometry geom = make_test_geometry(); + wsm6.Init(input.states, input.boxes, geom, Real(1.0), z_phys_nd, detJ_cc); + wsm6.Copy_State_to_Micro(input.states, &input.base_state); + wsm6.Copy_Micro_to_State(input.states); + amrex::Gpu::streamSynchronize(); + + EXPECT_NEAR(input.states.max(RhoTheta_comp), Real(1.10) * Real(303.0), + scaled_tol(Real(1.10) * Real(303.0), Real(1.0), Real(1.0e-11))); +} + +// Motivation: WDM6 must receive p0 during copy-in while retaining the scheme's +// existing pressure units and moisture-component mapping. +TEST(AnelasticMicrophysicsWiring, WDM6CopyInUsesReferencePressure) +{ + TestState input; + WorkingArrays work(input); + copy_wdm6(input, work); + expect_reference_pressure_diagnosis(input, work, Real(1.0)); + + const Real rho_theta_before = input.states.max(RhoTheta_comp); + SolverChoice sc = make_scheme_choice(MoistureType::WDM6, true); + WDM6 wdm6; + wdm6.SetCurrentLevel(0); + wdm6.Define(sc); + std::unique_ptr z_phys_nd; + std::unique_ptr detJ_cc; + const Geometry geom = make_test_geometry(); + wdm6.Init(input.states, input.boxes, geom, Real(1.0), z_phys_nd, detJ_cc); + wdm6.Copy_State_to_Micro(input.states, &input.base_state); + wdm6.Copy_Micro_to_State(input.states); + amrex::Gpu::streamSynchronize(); + + // With no Advance call, copy-in/copy-out must not create a thermodynamic + // tendency merely because the two cells use different reference pressures. + EXPECT_NEAR(input.states.max(RhoTheta_comp), rho_theta_before, + scaled_tol(rho_theta_before, rho_theta_before, Real(1.0e-11))); +} + +// Motivation: WSM6's compressible copy-out must retain the historical RdoCp +// exponent, while anelastic copy-out must use the configured SolverChoice +// exponent with the held reference pressure. +TEST(AnelasticMicrophysicsWiring, WSM6CustomCpKeepsModeSpecificExponent) +{ + constexpr Real custom_cp = Real(900.0); + + { + TestState input; + SolverChoice sc = make_scheme_choice(MoistureType::WSM6, false); + sc.c_p = custom_cp; + sc.rdOcp = R_d / custom_cp; + + WSM6 wsm6; + wsm6.SetCurrentLevel(0); + wsm6.Define(sc); + std::unique_ptr z_phys_nd; + std::unique_ptr detJ_cc; + const Geometry geom = make_test_geometry(); + wsm6.Init(input.states, input.boxes, geom, Real(1.0), z_phys_nd, detJ_cc); + wsm6.Copy_State_to_Micro(input.states); + wsm6.Copy_Micro_to_State(input.states); + amrex::Gpu::streamSynchronize(); + + EXPECT_NEAR(input.states.min(RhoTheta_comp), Real(0.90) * Real(287.0), + scaled_tol(Real(0.90) * Real(287.0), Real(1.0), Real(1.0e-11))); + EXPECT_NEAR(input.states.max(RhoTheta_comp), Real(1.10) * Real(303.0), + scaled_tol(Real(1.10) * Real(303.0), Real(1.0), Real(1.0e-11))); + + const Real temperature = Real(279.0); + const Real pressure = Real(87000.0); + const Real expected = getThgivenRandT( + Real(1.0), temperature, RdoCp, Real(0.004)); + const Real configured = getThgivenRandT( + Real(1.0), temperature, sc.rdOcp, Real(0.004)); + const Real actual = wsm6_theta_from_temperature_and_pressure( + Real(1.0), temperature, pressure, Real(0.004), sc.rdOcp, false); + EXPECT_NEAR(actual, expected, scaled_tol(actual, expected, Real(1.0e-12))); + EXPECT_GT(std::abs(expected - configured), Real(1.0e-2)); + } + + { + TestState input; + SolverChoice sc = make_scheme_choice(MoistureType::WSM6, true); + sc.c_p = custom_cp; + sc.rdOcp = R_d / custom_cp; + + WSM6 wsm6; + wsm6.SetCurrentLevel(0); + wsm6.Define(sc); + std::unique_ptr z_phys_nd; + std::unique_ptr detJ_cc; + const Geometry geom = make_test_geometry(); + wsm6.Init(input.states, input.boxes, geom, Real(1.0), z_phys_nd, detJ_cc); + wsm6.Copy_State_to_Micro(input.states, &input.base_state); + wsm6.Copy_Micro_to_State(input.states); + amrex::Gpu::streamSynchronize(); + + EXPECT_NEAR(input.states.min(RhoTheta_comp), Real(0.90) * Real(287.0), + scaled_tol(Real(0.90) * Real(287.0), Real(1.0), Real(1.0e-11))); + EXPECT_NEAR(input.states.max(RhoTheta_comp), Real(1.10) * Real(303.0), + scaled_tol(Real(1.10) * Real(303.0), Real(1.0), Real(1.0e-11))); + + const Real temperature = Real(279.0); + const Real pressure = Real(87000.0); + const Real expected = getThgivenTandP(temperature, pressure, sc.rdOcp); + const Real actual = wsm6_theta_from_temperature_and_pressure( + Real(1.0), temperature, pressure, Real(0.004), sc.rdOcp, true); + EXPECT_NEAR(actual, expected, scaled_tol(actual, expected, Real(1.0e-12))); + } +} + +// Motivation: Compressible calls must keep the established rho/rho*theta EOS +// diagnosis and must not accidentally consume a caller's base-state pressure. +TEST(AnelasticMicrophysicsWiring, CompressibleDiagnosisIsPreserved) +{ + const Real rho = Real(1.2); + const Real theta = Real(300.0); + const Real rho_theta = rho * theta; + const Real qv = Real(0.004); + const Real expected_pressure = getPgivenRTh(rho_theta, qv); + const Real expected_temperature = getTgivenRandRTh(rho, rho_theta, qv); + + const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state( + rho, rho_theta, qv, kRdOcp, false, Real(82000.0)); + EXPECT_NEAR(thermo.pressure_pa, expected_pressure, + scaled_tol(expected_pressure, expected_pressure, Real(1.0e-12))); + EXPECT_NEAR(thermo.temperature, expected_temperature, + scaled_tol(expected_temperature, expected_temperature, Real(1.0e-12))); + + const SAMPrimitiveCell sam_state = sam_cons_to_primitive( + rho, rho_theta, rho * qv, Real(0.0), Real(0.0), Real(0.0), Real(0.0), Real(0.0)); + EXPECT_NEAR(sam_state.pres_mbar, Real(0.01) * expected_pressure, + scaled_tol(expected_pressure, expected_pressure, Real(1.0e-12))); + EXPECT_NEAR(sam_state.tabs, expected_temperature, + scaled_tol(expected_temperature, expected_temperature, Real(1.0e-12))); +} + +// Motivation: SAM's compressible coefficient path must average rho, theta, and +// qv before applying the nonlinear EOS, rather than averaging diagnosed T. A +// non-null base-state pointer must not override the configured compressible mode. +TEST(AnelasticMicrophysicsWiring, SAMHeterogeneousCompressiblePlaneIgnoresBasePointer) +{ + TestState input; + SolverChoice sc = make_scheme_choice(MoistureType::SAM); + SAM sam; + sam.SetCurrentLevel(0); + sam.Define(sc); + std::unique_ptr z_phys_nd; + std::unique_ptr detJ_cc; + const Geometry geom = make_test_geometry(); + sam.Init(input.states, input.boxes, geom, Real(1.0), z_phys_nd, detJ_cc); + sam.Update_Micro_Vars(input.states, &input.base_state); + + const Real rho_bar = Real(1.0); + const Real theta_bar = Real(295.0); + const Real qv_bar = Real(0.005); + const Real rho_theta_bar = rho_bar * theta_bar; + const Real averaged_state_temperature = getTgivenRandRTh( + rho_bar, rho_theta_bar, qv_bar); + const SAMCoefficientRow expected = sam_compute_coefficient_row( + rho_bar, amrex::min(averaged_state_temperature, Real(273.16)), + erf_gammafff(three + b_rain), + erf_gammafff((Real(5.0) + b_rain) / two), + erf_gammafff(three + b_snow), + erf_gammafff((Real(5.0) + b_snow) / two), + erf_gammafff(three + b_grau), + erf_gammafff((Real(5.0) + b_grau) / two)); + + expect_coefficients_near(sam.CoefficientRowAt(0), expected); + + const Real cell_temperature_0 = getTgivenRandRTh( + Real(1.10), Real(1.10) * Real(303.0), Real(0.0030)); + const Real cell_temperature_1 = getTgivenRandRTh( + Real(0.90), Real(0.90) * Real(287.0), Real(0.0070)); + EXPECT_GT(std::abs(averaged_state_temperature - + Real(0.5) * (cell_temperature_0 + cell_temperature_1)), + Real(1.0e-2)); +} + +// Motivation: WDM6 uses the same held pressure in both conversion paths but +// retains RdoCp for compressible runs and uses configured rdOcp for anelastic. +TEST(AnelasticMicrophysicsWiring, WDM6CustomCpConversionIsModeSpecific) +{ + constexpr Real custom_cp = Real(900.0); + const Real configured_rdOcp = R_d / custom_cp; + const Real temperature = Real(279.0); + const Real pressure = Real(87000.0); + const Real expected_compressible = getThgivenTandP( + temperature, pressure, RdoCp); + const Real expected_anelastic = getThgivenTandP( + temperature, pressure, configured_rdOcp); + + const Real actual_compressible = wdm6_theta_from_temperature_and_pressure( + temperature, pressure, configured_rdOcp, false); + const Real actual_anelastic = wdm6_theta_from_temperature_and_pressure( + temperature, pressure, configured_rdOcp, true); + + EXPECT_NEAR(actual_compressible, expected_compressible, + scaled_tol(actual_compressible, expected_compressible, Real(1.0e-12))); + EXPECT_NEAR(actual_anelastic, expected_anelastic, + scaled_tol(actual_anelastic, expected_anelastic, Real(1.0e-12))); + EXPECT_GT(std::abs(expected_compressible - expected_anelastic), Real(1.0e-2)); +} + +// Motivation: SAM's heterogeneous anelastic plane must average the diagnosed +// pressure/temperature fields, not average theta and then rediagnose a local +// compressible EOS. +TEST(AnelasticMicrophysicsWiring, SAMHeterogeneousPlaneUsesDiagnosedTemperature) +{ + TestState input; + SolverChoice sc = make_scheme_choice(MoistureType::SAM, true); + SAM sam; + sam.SetCurrentLevel(0); + sam.Define(sc); + std::unique_ptr z_phys_nd; + std::unique_ptr detJ_cc; + const Geometry geom = make_test_geometry(); + sam.Init(input.states, input.boxes, geom, Real(1.0), z_phys_nd, detJ_cc); + sam.Update_Micro_Vars(input.states, &input.base_state); + + const Real temp0 = getTgivenPandTh(Real(90000.0), Real(303.0), kRdOcp); + const Real temp1 = getTgivenPandTh(Real(82000.0), Real(287.0), kRdOcp); + const Real tabs = amrex::min(Real(0.5) * (temp0 + temp1), Real(273.16)); + const SAMCoefficientRow expected = sam_compute_coefficient_row( + Real(1.0), tabs, + erf_gammafff(three + b_rain), + erf_gammafff((Real(5.0) + b_rain) / two), + erf_gammafff(three + b_snow), + erf_gammafff((Real(5.0) + b_snow) / two), + erf_gammafff(three + b_grau), + erf_gammafff((Real(5.0) + b_grau) / two)); + + expect_coefficients_near(sam.CoefficientRowAt(0), expected); +} + +// Motivation: The invalid SuperDroplets/anelastic combination must be +// recognized before the simulation reaches a compressible-only source path. +TEST(AnelasticMicrophysicsWiring, SuperDropletsAnelasticConfigurationIsInvalid) +{ + EXPECT_TRUE(anelastic_superdroplets_configuration_invalid( + MoistureType::SuperDroplets, true)); + EXPECT_FALSE(anelastic_superdroplets_configuration_invalid( + MoistureType::SuperDroplets, false)); + EXPECT_FALSE(anelastic_superdroplets_configuration_invalid( + MoistureType::SAM, true)); +}