Skip to content

Source/Utils: microphysics saturation and HSE numerics #3497

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

Open
ffmares wants to merge 2 commits into
base: development
Choose a base branch
from
Open

Source/Utils: microphysics saturation and HSE numerics #3497

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a788880
Source/Utils: microphysics saturation and HSE numerics
Jul 30, 2026
b4f85ac
Merge branch 'development' into pr/float-bubble-eos-tuning
asalmgren Aug 1, 2026
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
193 changes: 181 additions & 12 deletions Source/Utils/ERF_HSEUtils.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,24 @@ namespace HSEutils
const int MAX_ITER = 10;
#ifdef AMREX_USE_FLOAT
const amrex::Real TOL = Real(1.e-4);
const amrex::Real TEMP_STEP_REL_TOL = Real(1.e-4);
const amrex::Real TEMP_STEP_ABS_TOL = Real(1.e-3);
const amrex::Real TEMP_RES_REL_TOL = Real(1.e-6);
const amrex::Real TEMP_RES_ABS_TOL = Real(1.e-3);
const amrex::Real PRES_STEP_REL_TOL = Real(1.e-4);
const amrex::Real PRES_STEP_ABS_TOL = Real(1.e-2);
const amrex::Real PRES_RES_REL_TOL = Real(1.e-6);
const amrex::Real PRES_RES_ABS_TOL = Real(1.e-1);
#else
const amrex::Real TOL = Real(1.e-8);
const amrex::Real TEMP_STEP_REL_TOL = Real(1.e-8);
const amrex::Real TEMP_STEP_ABS_TOL = Real(1.e-8);
const amrex::Real TEMP_RES_REL_TOL = Real(1.e-10);
const amrex::Real TEMP_RES_ABS_TOL = Real(1.e-8);
const amrex::Real PRES_STEP_REL_TOL = Real(1.e-8);
const amrex::Real PRES_STEP_ABS_TOL = Real(1.e-6);
const amrex::Real PRES_RES_REL_TOL = Real(1.e-10);
const amrex::Real PRES_RES_ABS_TOL = Real(1.e-6);
#endif

/**
Expand Down Expand Up @@ -473,18 +489,107 @@ Real compute_F_for_temp_in_zone(const Real T_b, const Real p_b, const Real q_t,
return eq_pot_temp - T_b*std::pow((p_b - p_v)/p_0, -R_d/fac)*std::exp(L_v*q_v/(fac*T_b));
}

struct TempNewtonEval {
Real F;
Real dF;
Real dFdp;
};

AMREX_FORCE_INLINE
AMREX_GPU_HOST_DEVICE
Real compute_temperature (const Real p_b, const Real q_t, const Real eq_pot_temp, const bool use_empirical,
const int which_zone, const Real scaled_height)
TempNewtonEval compute_F_and_dF_for_temp_in_zone(const Real T_b, const Real p_b, const Real q_t,
const Real eq_pot_temp, const bool use_empirical,
const int which_zone, const Real scaled_height)
{
Real T_b = Real(200.0), delta_T; // Initial guess
Real fac = Cp_d + Cp_l*q_t;
Real RH = compute_relative_humidity(p_b, T_b, use_empirical, which_zone, scaled_height);
Real q_v = vapor_mixing_ratio(p_b, T_b, RH, use_empirical, which_zone);
Real p_s = compute_saturation_pressure(T_b, use_empirical);
Real p_v = compute_vapor_pressure(p_s, RH);
Real G = std::pow((p_b - p_v)/p_0, -R_d/fac) * std::exp(L_v*q_v/(fac*T_b));
Real F = eq_pot_temp - T_b*G;

Real dp_s_dT = Real(100.0) * erf_dtesatw(T_b, use_empirical);
Real dp_v_dT = Real(0.0);
Real dq_v_dT = Real(0.0);
Real dp_v_dp = Real(0.0);
Real dq_v_dp = Real(0.0);

if (which_zone == 1) {
Real q_s;
Real dq_s_dT;
erf_qsatw(T_b, p_b, q_s);
erf_dtqsatw(T_b, p_b, dq_s_dT);

Real q_ref = Real(0.014);
Real RH_ref = q_ref / q_s;
Real dRH_dT = -q_ref * dq_s_dT / (q_s * q_s);
dp_v_dT = dp_s_dT * RH_ref + p_s * dRH_dT;
dp_v_dp = q_ref / RdoRv;
} else {
dp_v_dT = RH * dp_s_dT;
dq_v_dT = RdoRv * p_b * dp_v_dT / ((p_b - p_v) * (p_b - p_v));
dq_v_dp = -RdoRv * p_v / ((p_b - p_v) * (p_b - p_v));
}

Real dlnA_dT = (R_d / fac) * (dp_v_dT / (p_b - p_v));
Real dB_dT = (L_v / fac) * (dq_v_dT / T_b - q_v / (T_b * T_b));
Real dF = -G * (Real(1.0) + T_b * (dlnA_dT + dB_dT));

Real dlnA_dp = - (R_d / fac) * (Real(1.0) - dp_v_dp) / (p_b - p_v);
Real dB_dp = (L_v / fac) * (dq_v_dp / T_b);
Real dFdp = -T_b * G * (dlnA_dp + dB_dp);

return {F, dF, dFdp};
}

AMREX_FORCE_INLINE
AMREX_GPU_HOST_DEVICE
Real compute_temperature (const Real p_b, const Real q_t, const Real eq_pot_temp, const bool use_empirical,
const int which_zone, const Real scaled_height, const Real T_guess = Real(200.0))
{
#ifdef AMREX_USE_FLOAT
Real eps = Real(1.e-4) * T_b;
Real T_b = (T_guess > Real(0.0)) ? T_guess : Real(200.0); // Initial guess
Real delta_T = Real(0.0);
Real residual_T = Real(0.0);
Real max_rel_step = Real(0.1);
bool converged = false;

for (int iter=0; iter<20; iter++)
{
TempNewtonEval newton_eval = compute_F_and_dF_for_temp_in_zone(T_b, p_b, q_t, eq_pot_temp,
use_empirical, which_zone, scaled_height);
residual_T = newton_eval.F;
delta_T = -residual_T / newton_eval.dF;

Real step_tol = TEMP_STEP_ABS_TOL + TEMP_STEP_REL_TOL * std::abs(T_b);
Real resid_tol = TEMP_RES_ABS_TOL + TEMP_RES_REL_TOL * std::abs(eq_pot_temp);
if (std::abs(delta_T) <= step_tol && std::abs(residual_T) <= resid_tol) {
converged = true;
break;
}

Real max_step = max_rel_step * std::abs(T_b);
if (std::abs(delta_T) > max_step) {
delta_T = std::copysign(max_step, delta_T);
}

T_b += delta_T;
}

if (!converged) {
amrex::Print() << "Newton Raphson for temperature could not converge; last residual = "
<< residual_T << ", last delta_T = " << delta_T << std::endl;
amrex::Abort("Newton Raphson for temperature could not converge");
}

return T_b;
#else
ignore_unused(T_guess);
Real T_b = Real(200.0), delta_T; // Initial guess

Real eps = Real(1.e-10) * T_b;
#endif

for (int iter=0; iter<20; iter++)
{
Real F = compute_F_for_temp_in_zone(T_b , p_b, q_t, eq_pot_temp, use_empirical, which_zone, scaled_height);
Expand All @@ -495,10 +600,11 @@ Real compute_temperature (const Real p_b, const Real q_t, const Real eq_pot_temp
}

if (std::fabs(delta_T) > TOL * T_b) {
amrex::Abort("Newton Raphson for temperature could not converge");
amrex::Abort("Newton Raphson for temperature could not converge");
}

return T_b;
#endif
}

AMREX_FORCE_INLINE
Expand Down Expand Up @@ -540,7 +646,10 @@ void compute_rho (const Real& pressure, Real& theta, Real& rho, Real& q_v, Real&
theta = compute_theta(scaled_height, theta_0, theta_tr, z_tr, T_tr);
T_b = getTgivenPandTh(pressure, theta, RdoCp);
} else {
T_b = compute_temperature(pressure, q_t, eq_pot_temp, use_empirical, which_zone, scaled_height);
if (!(T_b > Real(0.0))) {
T_b = eq_pot_temp;
}
T_b = compute_temperature(pressure, q_t, eq_pot_temp, use_empirical, which_zone, scaled_height, T_b);
theta = getThgivenTandP(T_b, pressure, RdoCp);
}

Expand Down Expand Up @@ -592,13 +701,72 @@ Real compute_p_k_in_zone (Real &p_k, const Real p_k_minus_1, Real &theta_k, Real
const int which_zone, const Real scaled_height, const bool T_from_theta,
const Real theta_0, const Real theta_tr, const Real z_tr, const Real T_tr)
{
Real delta_p_k;

#ifdef AMREX_USE_FLOAT
Real eps = Real(1e-4);
Real delta_p_k = Real(0.0);
Real residual_p = Real(0.0);
Real max_rel_step = Real(0.1);
Real rel_eps = Real(1e-6);
bool converged = false;

for(int iter=0; iter<20; iter++)
{
residual_p = compute_F(p_k, p_k_minus_1, theta_k, rho_k, q_v_k, T_dp, T_b, dz, rho_k_minus_1,
q_t, eq_pot_temp, use_empirical, which_zone, scaled_height, T_from_theta,
theta_0, theta_tr, z_tr, T_tr);

Real dp = rel_eps * std::max(Real(1.0), std::abs(p_k));
Real p_plus = p_k + dp;
Real p_minus = p_k - dp;

Real theta_tmp = theta_k;
Real rho_tmp = rho_k;
Real qv_tmp = q_v_k;
Real Tdp_tmp = T_dp;
Real Tb_tmp = T_b;
Real F_plus = compute_F(p_plus, p_k_minus_1, theta_tmp, rho_tmp, qv_tmp, Tdp_tmp, Tb_tmp, dz,
rho_k_minus_1, q_t, eq_pot_temp, use_empirical, which_zone,
scaled_height, T_from_theta, theta_0, theta_tr, z_tr, T_tr);

theta_tmp = theta_k;
rho_tmp = rho_k;
qv_tmp = q_v_k;
Tdp_tmp = T_dp;
Tb_tmp = T_b;
Real F_minus = compute_F(p_minus, p_k_minus_1, theta_tmp, rho_tmp, qv_tmp, Tdp_tmp, Tb_tmp, dz,
rho_k_minus_1, q_t, eq_pot_temp, use_empirical, which_zone,
scaled_height, T_from_theta, theta_0, theta_tr, z_tr, T_tr);

Real F_prime = (F_plus - F_minus) / (Real(2.0) * dp);
if (std::abs(F_prime) < Real(1e-14)) {
F_prime = std::copysign(Real(1e-14), F_prime == Real(0.0) ? Real(1.0) : F_prime);
}

delta_p_k = -residual_p / F_prime;
Real step_tol = PRES_STEP_ABS_TOL + PRES_STEP_REL_TOL * std::abs(p_k);
Real resid_scale = std::max(std::abs(p_k), std::max(std::abs(p_k_minus_1), Real(1.0)));
Real resid_tol = PRES_RES_ABS_TOL + PRES_RES_REL_TOL * resid_scale;
if (std::abs(delta_p_k) <= step_tol && std::abs(residual_p) <= resid_tol) {
converged = true;
break;
}

Real max_step = max_rel_step * std::abs(p_k);
if (std::abs(delta_p_k) > max_step) {
delta_p_k = std::copysign(max_step, delta_p_k);
}
p_k += delta_p_k;
}

if (!converged) {
amrex::Print() << "Newton Raphson for pressure could not converge; last residual = "
<< residual_p << ", last delta_p_k = " << delta_p_k << std::endl;
amrex::Abort("Newton Raphson for pressure could not converge");
}
return p_k;
#else
Real delta_p_k;

Real eps = Real(1e-10);
#endif

for(int iter=0; iter<20; iter++)
{
Expand All @@ -618,6 +786,7 @@ Real compute_p_k_in_zone (Real &p_k, const Real p_k_minus_1, Real &theta_k, Real
}

return p_k;
#endif
}

AMREX_FORCE_INLINE
Expand All @@ -630,7 +799,7 @@ init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v,
const Real z_tr_1 = -one, const Real z_tr_2 = -one,
const Real theta_0 = amrex::Real(0), const Real theta_tr = amrex::Real(0), const Real T_tr = amrex::Real(0))
{
Real T_b, T_dp;
Real T_b = eq_pot_temp, T_dp;

int which_zone = -1;
Real scaled_height = amrex::Real(0);
Expand Down
53 changes: 33 additions & 20 deletions Source/Utils/ERF_MicrophysicsUtils.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,33 +71,32 @@ amrex::Real erf_esatw_cc (amrex::Real t) {
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real erf_esatw_flatau_poly (amrex::Real dtt)
{
amrex::Real const a0 = amrex::Real(6.11239921);
amrex::Real const a1 = amrex::Real(0.443987641);
amrex::Real const a2 = amrex::Real(0.142986287e-1);
amrex::Real const a3 = amrex::Real(0.264847430e-3);
amrex::Real const a4 = amrex::Real(0.302950461e-5);
amrex::Real const a5 = amrex::Real(0.206739458e-7);
amrex::Real const a6 = amrex::Real(0.640689451e-10);
amrex::Real const a7 = -amrex::Real(0.952447341e-13);
amrex::Real const a8 = -amrex::Real(0.976195544e-15);
constexpr amrex::Real a0 = amrex::Real(6.11239921);
constexpr amrex::Real a1 = amrex::Real(0.443987641);
constexpr amrex::Real a2 = amrex::Real(0.142986287e-1);
constexpr amrex::Real a3 = amrex::Real(0.264847430e-3);
constexpr amrex::Real a4 = amrex::Real(0.302950461e-5);
constexpr amrex::Real a5 = amrex::Real(0.206739458e-7);
constexpr amrex::Real a6 = amrex::Real(0.640689451e-10);
constexpr amrex::Real a7 = -amrex::Real(0.952447341e-13);
constexpr amrex::Real a8 = -amrex::Real(0.976195544e-15);

return a0 + dtt*(a1+dtt*(a2+dtt*(a3+dtt*(a4+dtt*(a5+dtt*(a6+dtt*(a7+a8*dtt)))))));
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real erf_dtesatw_flatau_poly (amrex::Real dtt)
{
amrex::Real const a0 = amrex::Real(0.443956472);
amrex::Real const a1 = amrex::Real(0.285976452e-1);
amrex::Real const a2 = amrex::Real(0.794747212e-3);
amrex::Real const a3 = amrex::Real(0.121167162e-4);
amrex::Real const a4 = amrex::Real(0.103167413e-6);
amrex::Real const a5 = amrex::Real(0.385208005e-9);
amrex::Real const a6 = -amrex::Real(0.604119582e-12);
amrex::Real const a7 = -amrex::Real(0.792933209e-14);
amrex::Real const a8 = -amrex::Real(0.599634321e-17);
constexpr amrex::Real a1 = amrex::Real(0.443987641);
constexpr amrex::Real a2 = amrex::Real(0.142986287e-1);
constexpr amrex::Real a3 = amrex::Real(0.264847430e-3);
constexpr amrex::Real a4 = amrex::Real(0.302950461e-5);
constexpr amrex::Real a5 = amrex::Real(0.206739458e-7);
constexpr amrex::Real a6 = amrex::Real(0.640689451e-10);
constexpr amrex::Real a7 = -amrex::Real(0.952447341e-13);
constexpr amrex::Real a8 = -amrex::Real(0.976195544e-15);

return a0 + dtt*(a1+dtt*(a2+dtt*(a3+dtt*(a4+dtt*(a5+dtt*(a6+dtt*(a7+a8*dtt)))))));
return a1 + dtt*(amrex::Real(2.0)*a2 + dtt*(amrex::Real(3.0)*a3 + dtt*(amrex::Real(4.0)*a4 + dtt*(amrex::Real(5.0)*a5 + dtt*(amrex::Real(6.0)*a6 + dtt*(amrex::Real(7.0)*a7 + amrex::Real(8.0)*a8*dtt))))));
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand Down Expand Up @@ -183,7 +182,21 @@ amrex::Real erf_dtesatw_cc (amrex::Real t) {
// at the switch because the Flatau polynomial and Magnus-style exponential
// approximation are independent fits rather than a matched composite model.
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real erf_dtesatw (amrex::Real t) {
amrex::Real erf_dtesatw (amrex::Real t, bool use_empirical=false) {
if (use_empirical) {
constexpr amrex::Real t0 = amrex::Real(273.15);
constexpr amrex::Real b = amrex::Real(34.494);
constexpr amrex::Real c1 = amrex::Real(4924.99);
constexpr amrex::Real c2 = amrex::Real(237.1);
constexpr amrex::Real d = amrex::Real(105.0);
constexpr amrex::Real n = amrex::Real(1.57);

amrex::Real x = t - t0;
amrex::Real esatw = amrex::Real(0.01) * std::exp(b - c1/(x + c2)) / std::pow(x + d, n);
amrex::Real dln = c1 / ((x + c2) * (x + c2)) - n / (x + d);
return esatw * dln;
}

if (erf_use_positive_esatw_poly(t)) {
return erf_dtesatw_flatau_poly(t - amrex::Real(273.16));
}
Expand Down
Loading