Skip to content

Commit f110b57

Browse files
sajjadazimisajjadazimi
authored
Merge pull request #3813 from CliMA/aj/limiters_prep
Move limiter for rhs out of wrappers
2 parents a4fa008 + 19f113f commit f110b57

File tree

3 files changed

+42
-34
lines changed

3 files changed

+42
-34
lines changed

‎reproducibility_tests/ref_counter.jl‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
240
1+
241
22

33
# **README**
44
#
@@ -20,6 +20,9 @@
2020

2121

2222
#=
23+
241
24+
- Add triangle inequality limiter from Horn (2012) and apply it in microphysics
25+
2326
240
2427
- Fix a bug in viscous sponge
2528

‎src/parameterized_tendencies/microphysics/microphysics_wrappers.jl‎

Lines changed: 31 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ function qₚ(q_rain_snow)
2929
return max(FT(0), q_rain_snow)
3030
end
3131

32-
# Helper function to limit the tendency
32+
# Helper function to compute the limit of the tendency
3333
function limit(q, dt, n::Int)
3434
return q / float(dt) / n
3535
end
@@ -85,11 +85,10 @@ function cloud_sources(
8585

8686
S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cm_params, thp, q, ρ, Tₐ(thp, ts))
8787

88-
# keeping the same limiter for now
8988
return ifelse(
9089
S > FT(0),
91-
min(S, limit(qᵥ(thp, ts), dt, 2)),
92-
-min(abs(S), limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 2)),
90+
triangle_inequality_limiter(S, limit(qᵥ(thp, ts), dt, 2)),
91+
-triangle_inequality_limiter(abs(S), limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 2)),
9392
)
9493
end
9594
function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, qₛ, dt) where {FT}
@@ -99,11 +98,10 @@ function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, qₛ, dt) where {FT
9998

10099
S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cm_params, thp, q, ρ, Tₐ(thp, ts))
101100

102-
# keeping the same limiter for now
103101
return ifelse(
104102
S > FT(0),
105-
min(S, limit(qᵥ(thp, ts), dt, 2)),
106-
-min(abs(S), limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 2)),
103+
triangle_inequality_limiter(S, limit(qᵥ(thp, ts), dt, 2)),
104+
-triangle_inequality_limiter(abs(S), limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 2)),
107105
)
108106
end
109107

@@ -119,9 +117,9 @@ Returns the qₜ source term due to precipitation formation
119117
defined as Δm_tot / (m_dry + m_tot) for the 0-moment scheme
120118
"""
121119
function q_tot_0M_precipitation_sources(thp, cmp::CMP.Parameters0M, dt, qₜ, ts)
122-
return -min(
123-
max(qₜ, 0) / float(dt),
120+
return -triangle_inequality_limiter(
124121
-CM0.remove_precipitation(cmp, PP(thp, ts)),
122+
max(qₜ, 0) / float(dt),
125123
)
126124
end
127125

@@ -190,77 +188,77 @@ function compute_precipitation_sources!(
190188
CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts, qₚ(qᵣ)), true),
191189
CM2.conv_q_liq_to_q_rai(mp.var, qₗ(thp, ts, qₚ(qᵣ)), ρ, mp.Ndp),
192190
)
193-
@. Sᵖ = min(limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5), Sᵖ)
191+
@. Sᵖ = triangle_inequality_limiter(Sᵖ, limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5))
194192
@. Sqₗᵖ -= Sᵖ
195193
@. Sqᵣᵖ += Sᵖ
196194

197195
# snow autoconversion assuming no supersaturation: q_ice -> q_snow
198-
@. Sᵖ = min(
199-
limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
196+
@. Sᵖ = triangle_inequality_limiter(
200197
CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ(thp, ts, qₚ(qₛ)), true),
198+
limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
201199
)
202200
@. Sqᵢᵖ -= Sᵖ
203201
@. Sqₛᵖ += Sᵖ
204202

205203
# accretion: q_liq + q_rain -> q_rain
206-
@. Sᵖ = min(
207-
limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
204+
@. Sᵖ = triangle_inequality_limiter(
208205
CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ(thp, ts, qₚ(qᵣ)), qₚ(qᵣ), ρ),
206+
limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
209207
)
210208
@. Sqₗᵖ -= Sᵖ
211209
@. Sqᵣᵖ += Sᵖ
212210

213211
# accretion: q_ice + q_snow -> q_snow
214-
@. Sᵖ = min(
215-
limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
212+
@. Sᵖ = triangle_inequality_limiter(
216213
CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ(thp, ts, qₚ(qₛ)), qₚ(qₛ), ρ),
214+
limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
217215
)
218216
@. Sqᵢᵖ -= Sᵖ
219217
@. Sqₛᵖ += Sᵖ
220218

221219
# accretion: q_liq + q_sno -> q_sno or q_rai
222220
# sink of cloud water via accretion cloud water + snow
223-
@. Sᵖ = min(
224-
limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
221+
@. Sᵖ = triangle_inequality_limiter(
225222
CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts, qₚ(qᵣ)), qₚ(qₛ), ρ),
223+
limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
226224
)
227225
# if T < T_freeze cloud droplets freeze to become snow
228226
# else the snow melts and both cloud water and snow become rain
229227
α(thp, ts) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze)
230228
@. Sᵖ_snow = ifelse(
231229
Tₐ(thp, ts) < mp.ps.T_freeze,
232230
Sᵖ,
233-
FT(-1) * min(Sᵖ * α(thp, ts), limit(qₚ(qₛ), dt, 5)),
231+
FT(-1) * triangle_inequality_limiter(Sᵖ * α(thp, ts), limit(qₚ(qₛ), dt, 5)),
234232
)
235233
@. Sqₛᵖ += Sᵖ_snow
236234
@. Sqₗᵖ -= Sᵖ
237235
@. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow)
238236

239237
# accretion: q_ice + q_rai -> q_sno
240-
@. Sᵖ = min(
241-
limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
238+
@. Sᵖ = triangle_inequality_limiter(
242239
CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₚ(qₛ)), qₚ(qᵣ), ρ),
240+
limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
243241
)
244242
@. Sqᵢᵖ -= Sᵖ
245243
@. Sqₛᵖ += Sᵖ
246244
# sink of rain via accretion cloud ice - rain
247-
@. Sᵖ = min(
248-
limit(qₚ(qᵣ), dt, 5),
245+
@. Sᵖ = triangle_inequality_limiter(
249246
CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₚ(qₛ)), qₚ(qᵣ), ρ),
247+
limit(qₚ(qᵣ), dt, 5),
250248
)
251249
@. Sqᵣᵖ -= Sᵖ
252250
@. Sqₛᵖ += Sᵖ
253251

254252
# accretion: q_rai + q_sno -> q_rai or q_sno
255253
@. Sᵖ = ifelse(
256254
Tₐ(thp, ts) < mp.ps.T_freeze,
257-
min(
258-
limit(qₚ(qᵣ), dt, 5),
255+
triangle_inequality_limiter(
259256
CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₚ(qₛ), qₚ(qᵣ), ρ),
257+
limit(qₚ(qᵣ), dt, 5),
260258
),
261-
-min(
262-
limit(qₚ(qₛ), dt, 5),
259+
-triangle_inequality_limiter(
263260
CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qₚ(qᵣ), qₚ(qₛ), ρ),
261+
limit(qₚ(qₛ), dt, 5),
264262
),
265263
)
266264
@. Sqₛᵖ += Sᵖ
@@ -302,16 +300,16 @@ function compute_precipitation_sinks!(
302300

303301
#! format: off
304302
# evaporation: q_rai -> q_vap
305-
@. Sᵖ = -min(
306-
limit(qₚ(qᵣ), dt, 5),
303+
@. Sᵖ = -triangle_inequality_limiter(
307304
-CM1.evaporation_sublimation(rps..., PP(thp, ts), qₚ(qᵣ), ρ, Tₐ(thp, ts)),
305+
limit(qₚ(qᵣ), dt, 5),
308306
)
309307
@. Sqᵣᵖ += Sᵖ
310308

311309
# melting: q_sno -> q_rai
312-
@. Sᵖ = min(
313-
limit(qₚ(qₛ), dt, 5),
310+
@. Sᵖ = triangle_inequality_limiter(
314311
CM1.snow_melt(sps..., qₚ(qₛ), ρ, Tₐ(thp, ts)),
312+
limit(qₚ(qₛ), dt, 5),
315313
)
316314
@. Sqᵣᵖ += Sᵖ
317315
@. Sqₛᵖ -= Sᵖ
@@ -320,8 +318,8 @@ function compute_precipitation_sinks!(
320318
@. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₚ(qₛ), ρ, Tₐ(thp, ts))
321319
@. Sᵖ = ifelse(
322320
Sᵖ > FT(0),
323-
min(limit(qᵥ(thp, ts), dt, 5), Sᵖ),
324-
-min(limit(qₚ(qₛ), dt, 5), FT(-1) * Sᵖ),
321+
triangle_inequality_limiter(Sᵖ, limit(qᵥ(thp, ts), dt, 5)),
322+
-triangle_inequality_limiter(FT(-1) * Sᵖ, limit(qₚ(qₛ), dt, 5)),
325323
)
326324
@. Sqₛᵖ += Sᵖ
327325
#! format: on

‎src/prognostic_equations/limited_tendencies.jl‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11

2+
# Interacting with ClimaCore quasimonotone limiter for advection
23
NVTX.@annotate function limiters_func!(Y, p, t, ref_Y)
34
(; limiter) = p.numerics
45
if !isnothing(limiter)
@@ -9,3 +10,9 @@ NVTX.@annotate function limiters_func!(Y, p, t, ref_Y)
910
end
1011
return nothing
1112
end
13+
14+
# Helper function to limit source term tendencies the trianlge inequality-based
15+
# limiter from Horn (2012, doi:10.5194/gmd-5-345-2012)
16+
function triangle_inequality_limiter(force, limit)
17+
return force + limit - sqrt(force^2 + limit^2)
18+
end

0 commit comments

Comments
 (0)