Skip to content

Implementing double-sided limiters #4037

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
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Implementing double-sided limiters #4037

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
30 changes: 21 additions & 9 deletions src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ end
# The limit is defined as the available amont / n times model time step.
function limit(q, dt, n::Int)
FT = eltype(q)
return max(FT(0), q) / float(dt) / n
return q / float(dt) / n
end

"""
Expand Down Expand Up @@ -85,8 +85,8 @@ function cloud_sources(

return ifelse(
S > FT(0),
triangle_inequality_limiter(S, limit(qᵥ - qₛₗ, dt, 2)),
-triangle_inequality_limiter(abs(S), limit(qₗ, dt, 2)),
triangle_inequality_limiter(S, limit(qᵥ - qₛₗ, dt, 2), limit(qₗ, dt, 2)),
-triangle_inequality_limiter(abs(S), limit(qₗ, dt, 2), limit(qᵥ - qₛₗ, dt, 2)),
)
end
function cloud_sources(
Expand Down Expand Up @@ -135,8 +135,8 @@ function cloud_sources(

return ifelse(
S > FT(0),
triangle_inequality_limiter(S, limit(qᵥ - qₛᵢ, dt, 2)),
-triangle_inequality_limiter(abs(S), limit(qᵢ, dt, 2)),
triangle_inequality_limiter(S, limit(qᵥ - qₛᵢ, dt, 2), limit(qᵢ, dt, 2)),
-triangle_inequality_limiter(abs(S), limit(qᵢ, dt, 2), limit(qᵥ - qₛᵢ, dt, 2)),
)
end

Expand All @@ -154,7 +154,7 @@ defined as Δm_tot / (m_dry + m_tot) for the 0-moment scheme
function q_tot_0M_precipitation_sources(thp, cmp::CMP.Parameters0M, dt, qₜ, ts)
return -triangle_inequality_limiter(
-CM0.remove_precipitation(cmp, PP(thp, ts)),
max(qₜ, 0) / float(dt),
qₜ / float(dt),
)
end

Expand Down Expand Up @@ -224,14 +224,15 @@ function compute_precipitation_sources!(
CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ, true),
CM2.conv_q_liq_to_q_rai(mp.var, qₗ, ρ, mp.Ndp),
)
@. Sᵖ = triangle_inequality_limiter(Sᵖ, limit(qₗ, dt, 5))
@. Sᵖ = triangle_inequality_limiter(Sᵖ, limit(qₗ, dt, 5), limit(qᵣ, dt, 5))
@. Sqₗᵖ -= Sᵖ
@. Sqᵣᵖ += Sᵖ

# snow autoconversion assuming no supersaturation: q_ice -> q_snow
@. Sᵖ = triangle_inequality_limiter(
CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ, true),
limit(qᵢ, dt, 5),
limit(qₛ, dt, 5),
)
@. Sqᵢᵖ -= Sᵖ
@. Sqₛᵖ += Sᵖ
Expand All @@ -240,6 +241,7 @@ function compute_precipitation_sources!(
@. Sᵖ = triangle_inequality_limiter(
CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ, qᵣ, ρ),
limit(qₗ, dt, 5),
limit(qᵣ, dt, 5),
)
@. Sqₗᵖ -= Sᵖ
@. Sqᵣᵖ += Sᵖ
Expand All @@ -248,6 +250,7 @@ function compute_precipitation_sources!(
@. Sᵖ = triangle_inequality_limiter(
CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ, qₛ, ρ),
limit(qᵢ, dt, 5),
limit(qₛ, dt, 5),
)
@. Sqᵢᵖ -= Sᵖ
@. Sqₛᵖ += Sᵖ
Expand All @@ -274,13 +277,15 @@ function compute_precipitation_sources!(
@. Sᵖ = triangle_inequality_limiter(
CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ, qᵣ, ρ),
limit(qᵢ, dt, 5),
limit(qₛ, dt, 5),
)
@. Sqᵢᵖ -= Sᵖ
@. Sqₛᵖ += Sᵖ
# sink of rain via accretion cloud ice - rain
@. Sᵖ = triangle_inequality_limiter(
CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ, qᵣ, ρ),
limit(qᵣ, dt, 5),
limit(qₛ, dt, 5),
)
@. Sqᵣᵖ -= Sᵖ
@. Sqₛᵖ += Sᵖ
Expand All @@ -291,10 +296,12 @@ function compute_precipitation_sources!(
triangle_inequality_limiter(
CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₛ, qᵣ, ρ),
limit(qᵣ, dt, 5),
limit(qₛ, dt, 5),
),
-triangle_inequality_limiter(
CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qᵣ, qₛ, ρ),
limit(qₛ, dt, 5),
limit(qᵣ, dt, 5),
),
)
@. Sqₛᵖ += Sᵖ
Expand Down Expand Up @@ -342,13 +349,15 @@ function compute_precipitation_sinks!(
@. Sᵖ = -triangle_inequality_limiter(
-CM1.evaporation_sublimation(rps..., qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ(thp, ts)),
limit(qᵣ, dt, 5),
limit(qᵥ(thp, ts), dt, 5),
)
@. Sqᵣᵖ += Sᵖ

# melting: q_sno -> q_rai
@. Sᵖ = triangle_inequality_limiter(
CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)),
limit(qₛ, dt, 5),
limit(qᵣ, dt, 5),
)
@. Sqᵣᵖ += Sᵖ
@. Sqₛᵖ -= Sᵖ
Expand All @@ -357,8 +366,8 @@ function compute_precipitation_sinks!(
@. Sᵖ = CM1.evaporation_sublimation(sps..., qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ(thp, ts))
@. Sᵖ = ifelse(
Sᵖ > FT(0),
triangle_inequality_limiter(Sᵖ, limit(qᵥ(thp, ts), dt, 5)),
-triangle_inequality_limiter(FT(-1) * Sᵖ, limit(qₛ, dt, 5)),
triangle_inequality_limiter(Sᵖ, limit(qᵥ(thp, ts), dt, 5), limit(qₛ, dt, 5)),
-triangle_inequality_limiter(FT(-1) * Sᵖ, limit(qₛ, dt, 5), limit(qᵥ(thp, ts), dt, 5)),
)
@. Sqₛᵖ += Sᵖ
#! format: on
Expand Down Expand Up @@ -608,6 +617,7 @@ function compute_warm_precipitation_sources_2M!(
ρ * nₗ,
).dq_rai_dt,
limit(qₗ, dt, 5), # cap rate to at most 20% of qₗ per timestep to ensure stability
limit(qᵣ, dt, 5),
)
@. Sqₗᵖ -= Sᵖ
@. Sqᵣᵖ += Sᵖ
Expand Down Expand Up @@ -654,6 +664,7 @@ function compute_warm_precipitation_sources_2M!(
@. Sᵖ = triangle_inequality_limiter(
CM2.accretion(mp.sb, qₗ, qᵣ, ρ, ρ * nₗ).dq_rai_dt,
limit(qₗ, dt, 5),
limit(qᵣ, dt, 5),
)
@. Sqₗᵖ -= Sᵖ
@. Sqᵣᵖ += Sᵖ
Expand Down Expand Up @@ -683,6 +694,7 @@ function compute_warm_precipitation_sources_2M!(
Tₐ(thp, ts),
).evap_rate_1,
limit(qᵣ, dt, 5),
limit(qᵥ(thp, ts), dt, 5),
)
@. Sqᵣᵖ += Sᵖ

Expand Down
49 changes: 40 additions & 9 deletions src/prognostic_equations/limited_tendencies.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,35 +31,66 @@ NVTX.@annotate function limiters_func!(Y, p, t, ref_Y)
end

"""
triangle_inequality_limiter(force, limit)
triangle_inequality_limiter(force, allowed_source_amount, limit_neg=0)

Limits a `force` (or source term) based on a maximum allowable `limit` using a
Limits a `force` (or source term) based on maximum allowable limits using a
formula derived from the triangle inequality, as proposed by Horn (2012).

The formula used is: `L = F + M - sqrt(F² + M²)`, where `F` is the `force` and
`M` is the `limit`.
`M` is the `allowed_source_amount`.

This limiter is designed to smoothly reduce the `force` as it approaches or
exceeds the `limit`, ensuring the result `L` satisfies `0 ≤ L ≤ M` if `F ≥ 0`
exceeds the `allowed_source_amount`, ensuring the result `L` satisfies `0 ≤ L ≤ M` if `F ≥ 0`
and `M > 0`. It also preserves `L ≤ F`. It's particularly useful for ensuring
that source terms (e.g., emissions, chemical production rates) do not become
unphysically large or lead to numerical instability, while being continuously
differentiable.


Due to numerical errors, `allowed_source_amount` can be negative.
This spurious behavior leads to returned `force` switching signs.
In this case the source and destination categories are swapped.
We try to limit this new force by the amount available in the old destination category.
If both source and destination are negative, we return zero, as there is nothing the limiter can do.

Arguments:
- `force`: The original force or source term value.
- `limit`: The maximum permissible positive value for the limited force.
- `force`: The original force or source term value (positive or negative).
- `intended_source_amount`: The available mass in the source category (positive normally,
but can be negative due to numerical errors).
- `limit_neg`: The available mass in the destination category for reverse conversions
(defaults to 0). Only used when `intended_source_amount ≤ 0`.

Returns:
- The limited force value.
- The limited force value, guaranteed not to deplete more mass than available in either category.

Reference:
- Horn, M. (2012). "ASAMgpu V1.0 – a moist fully compressible atmospheric model using
graphics processing units (GPUs)". Geoscientific Model Development,
5, 345–353. https://doi.org/10.5194/gmd-5-345-2012
"""
function triangle_inequality_limiter(force, allowed_source_amount, limit_neg = 0)

function triangle_inequality_limiter(force, limit)
FT = eltype(force)
return force == FT(0) ? force : force + limit - sqrt(force^2 + limit^2)

if force == FT(0) || allowed_source_amount == FT(0)
# If the force is zero or the intended source amount is zero, return zero.
# (We have a separate case for it, because of the numerics issues we were sometimes getting
# small values after limiting even when force was zero.)
return FT(0)
elseif allowed_source_amount > FT(0)
# If the intended source amount is positive, limit the force by the intended source amount.
return force + allowed_source_amount - sqrt(force^2 + allowed_source_amount^2)
elseif limit_neg > FT(0)
# Due to numerics we might end up with negative numbers in the intended_source_amount.
# This results in the tendecy switching signs. In that case we need to try to limit by the
# amound available in the destination category.
force1 = force + allowed_source_amount - sqrt(force^2 + allowed_source_amount^2)
force2 = -force1 + limit_neg - sqrt(force1^2 + limit_neg^2)
return -force2
else
# If both the source and destination are negatve, we set the force to zero
# (There is nothing the limiter can do in this case.)
return FT(0)
end
end

Loading