Skip to content

Support reference pressure in anelastic microphysics #3564

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
pressel wants to merge 1 commit into
base: development
Choose a base branch
from
Draft

Support reference pressure in anelastic microphysics #3564

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 18 additions & 0 deletions Docs/sphinx_doc/SuperDroplets.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
-------------------

Expand Down
58 changes: 58 additions & 0 deletions Docs/sphinx_doc/theory/Microphysics.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
-----------------------------------

Expand Down
20 changes: 20 additions & 0 deletions Source/DataStructs/ERF_DataStruct.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*/
Expand Down Expand Up @@ -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.
Expand Down
1 change: 1 addition & 0 deletions Source/Microphysics/ERF_EulerianMicrophysics.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down
34 changes: 21 additions & 13 deletions Source/Microphysics/Kessler/ERF_InitKessler.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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<Real const>{};
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);
}
52 changes: 49 additions & 3 deletions Source/Microphysics/Kessler/ERF_Kessler.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<const amrex::Real>& states,
const amrex::Array4<const amrex::Real>& base,
const amrex::Array4<amrex::Real>& rho,
const amrex::Array4<amrex::Real>& theta,
const amrex::Array4<amrex::Real>& qv,
const amrex::Array4<amrex::Real>& qc,
const amrex::Array4<amrex::Real>& qp,
const amrex::Array4<amrex::Real>& qt,
const amrex::Array4<amrex::Real>& tabs,
const amrex::Array4<amrex::Real>& 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 {
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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;

Expand Down
59 changes: 32 additions & 27 deletions Source/Microphysics/Morrison/ERF_InitMorrison.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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<Real const>{};
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);
}
Loading
Loading