Skip to content

Commit 2d4e94d

Browse files
charleskawczynskicharleskawczynski
committed
Implement and test vertical mass borrowing limiters
Update src/Limiters/vertical_mass_borrowing_limiter.jl Co-authored-by: Tapio Schneider <[email protected]> Update src/Limiters/vertical_mass_borrowing_limiter.jl Co-authored-by: Tapio Schneider <[email protected]> Update src/Limiters/vertical_mass_borrowing_limiter.jl Co-authored-by: Tapio Schneider <[email protected]> Use density-dz for pressure thickness
1 parent 8aa9af9 commit 2d4e94d

File tree

7 files changed

+446
-1
lines changed

7 files changed

+446
-1
lines changed

‎NEWS.md

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,9 @@ ClimaCore.jl Release Notes
44
main
55
-------
66

7+
- A new limiter, `VerticalMassBorrowingLimiter`, was added. Which redistributes all negative values from a given field while preserving mass.
8+
PR [2084](https://github.com/CliMA/ClimaCore.jl/pull/2084)
9+
710
### ![][badge-✨feature/enhancement] Various improvements to `Remapper` [2060](https://github.com/CliMA/ClimaCore.jl/pull/2060)
811

912
The `ClimaCore.Remapping` module received two improvements. First, `Remapper` is

‎docs/refs.bib

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -271,3 +271,14 @@ @article{Wiin1967
271271
year = {1967},
272272
publisher = {Wiley Online Library}
273273
}
274+
275+
@article{zhang2018impact,
276+
title={Impact of numerical choices on water conservation in the E3SM Atmosphere Model version 1 (EAMv1)},
277+
author={Zhang, Kai and Rasch, Philip J and Taylor, Mark A and Wan, Hui and Leung, Ruby and Ma, Po-Lun and Golaz, Jean-Christophe and Wolfe, Jon and Lin, Wuyin and Singh, Balwinder and others},
278+
journal={Geoscientific Model Development},
279+
volume={11},
280+
number={5},
281+
pages={1971--1988},
282+
year={2018},
283+
publisher={Copernicus Publications G{\"o}ttingen, Germany}
284+
}

‎docs/src/api.md

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -356,14 +356,15 @@ Hypsography.ref_z_to_physical_z
356356
The limiters supertype is
357357
```@docs
358358
Limiters.AbstractLimiter
359+
Limiters.QuasiMonotoneLimiter
360+
Limiters.VerticalMassBorrowingLimiter
359361
```
360362

361363
This class of flux-limiters is applied only in the horizontal direction (on spectral advection operators).
362364

363365

364366
### Interfaces
365367
```@docs
366-
Limiters.QuasiMonotoneLimiter
367368
Limiters.compute_bounds!
368369
Limiters.apply_limiter!
369370
```

‎src/Limiters/Limiters.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,5 +20,6 @@ abstract type AbstractLimiter end
2020

2121
# implementations
2222
include("quasimonotone.jl")
23+
include("vertical_mass_borrowing_limiter.jl")
2324

2425
end # end module

‎src/Limiters/vertical_mass_borrowing_limiter.jl

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
import .DataLayouts as DL
2+
3+
"""
4+
VerticalMassBorrowingLimiter(f::Fields.Field, q_min)
5+
6+
A vertical-only mass borrowing limiter.
7+
8+
The mass borrower borrows tracer mass from an adjacent, lower layer.
9+
It conserves the total tracer mass and can avoid negative tracers.
10+
11+
At level k, it will first borrow the mass from the layer k+1 (lower level).
12+
If the mass is not sufficient in layer k+1, it will borrow mass from
13+
layer k+2. The borrower will proceed this process until the bottom layer.
14+
If the tracer mass in the bottom layer goes negative, it will repeat the
15+
process from the bottom to the top. In this way, the borrower works for
16+
any shape of mass profiles.
17+
18+
This code was adapted from [E3SM](https://github.com/E3SM-Project/E3SM/blob/2c377c5ec9a5585170524b366ad85074ab1b1a5c/components/eam/src/physics/cam/massborrow.F90)
19+
20+
References:
21+
- [zhang2018impact](@cite)
22+
"""
23+
struct VerticalMassBorrowingLimiter{F, T}
24+
bmass::F
25+
ic::F
26+
q_min::T
27+
end
28+
function VerticalMassBorrowingLimiter(f::Fields.Field, q_min)
29+
bmass = similar(Spaces.level(f, 1))
30+
ic = similar(Spaces.level(f, 1))
31+
return VerticalMassBorrowingLimiter(bmass, ic, q_min)
32+
end
33+
34+
35+
"""
36+
apply_limiter!(q::Fields.Field, ρ::Fields.Field, lim::VerticalMassBorrowingLimiter)
37+
38+
Apply the vertical mass borrowing
39+
limiter `lim` to field `q`, given
40+
density `ρ`.
41+
"""
42+
apply_limiter!(
43+
q::Fields.Field,
44+
ρ::Fields.Field,
45+
lim::VerticalMassBorrowingLimiter,
46+
) = apply_limiter!(q, ρ, axes(q), lim, ClimaComms.device(axes(q)))
47+
48+
function apply_limiter!(
49+
q::Fields.Field,
50+
ρ::Fields.Field,
51+
space::Spaces.FiniteDifferenceSpace,
52+
lim::VerticalMassBorrowingLimiter,
53+
device::ClimaComms.AbstractCPUDevice,
54+
)
55+
cache = (; bmass = lim.bmass, ic = lim.ic, q_min = lim.q_min)
56+
columnwise_massborrow_cpu(q, ρ, cache)
57+
end
58+
59+
function apply_limiter!(
60+
q::Fields.Field,
61+
ρ::Fields.Field,
62+
space::Spaces.ExtrudedFiniteDifferenceSpace,
63+
lim::VerticalMassBorrowingLimiter,
64+
device::ClimaComms.AbstractCPUDevice,
65+
)
66+
Fields.bycolumn(axes(q)) do colidx
67+
cache = (;
68+
bmass = lim.bmass[colidx],
69+
ic = lim.ic[colidx],
70+
q_min = lim.q_min,
71+
)
72+
columnwise_massborrow_cpu(q[colidx], ρ[colidx], cache)
73+
end
74+
end
75+
76+
# TODO: can we improve the performance?
77+
# `bycolumn` on the CPU may be better here since we could multithread it.
78+
function columnwise_massborrow_cpu(q::Fields.Field, ρ::Fields.Field, cache) # column fields
79+
(; bmass, ic, q_min) = cache
80+
81+
Δz = Fields.Δz_field(q)
82+
Δz_vals = Fields.field_values(Δz)
83+
ρ_vals = Fields.field_values(ρ)
84+
# loop over tracers
85+
nlevels = Spaces.nlevels(axes(q))
86+
@. ic = 0
87+
@. bmass = 0
88+
q_vals = Fields.field_values(q)
89+
# top to bottom
90+
for f in 1:DataLayouts.ncomponents(q_vals)
91+
for v in 1:nlevels
92+
CI = CartesianIndex(1, 1, f, v, 1)
93+
# new mass in the current layer
94+
ρΔz_v =
95+
DL.getindex_field(Δz_vals, CI) * DL.getindex_field(ρ_vals, CI)
96+
nmass = DL.getindex_field(q_vals, CI) + bmass[] / ρΔz_v
97+
if nmass > q_min[f]
98+
# if new mass in the current layer is positive, don't borrow mass any more
99+
DL.setindex_field!(q_vals, nmass, CI)
100+
bmass[] = 0
101+
else
102+
# set mass to q_min in the current layer, and save bmass
103+
bmass[] = (nmass - q_min[f]) * ρΔz_v
104+
DL.setindex_field!(q_vals, q_min[f], CI)
105+
ic[] = ic[] + 1
106+
end
107+
end
108+
109+
# bottom to top
110+
for v in nlevels:-1:1
111+
CI = CartesianIndex(1, 1, f, v, 1)
112+
# if the surface layer still needs to borrow mass
113+
if bmass[] < 0
114+
ρΔz_v =
115+
DL.getindex_field(Δz_vals, CI) *
116+
DL.getindex_field(ρ_vals, CI)
117+
# new mass in the current layer
118+
nmass = DL.getindex_field(q_vals, CI) + bmass[] / ρΔz_v
119+
if nmass > q_min[f]
120+
# if new mass in the current layer is positive, don't borrow mass any more
121+
DL.setindex_field!(q_vals, nmass, CI)
122+
bmass[] = 0
123+
else
124+
# if new mass in the current layer is negative, continue to borrow mass
125+
bmass[] = (nmass - q_min[f]) * ρΔz_v
126+
DL.setindex_field!(q_vals, q_min[f], CI)
127+
end
128+
end
129+
end
130+
end
131+
132+
return nothing
133+
end

‎test/Limiters/vertical_mass_borrowing_limiter.jl

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
#=
2+
julia --project
3+
using Revise; include(joinpath("test", "Limiters", "vertical_mass_borrowing_limiter.jl"))
4+
=#
5+
using ClimaComms
6+
ClimaComms.@import_required_backends
7+
using ClimaCore: Fields, Spaces, Limiters
8+
using ClimaCore.RecursiveApply
9+
using ClimaCore.Geometry
10+
using ClimaCore.Grids
11+
using ClimaCore.CommonGrids
12+
using Test
13+
using Random
14+
15+
function perturb_field!(f::Fields.Field; perturb_radius)
16+
device = ClimaComms.device(f)
17+
ArrayType = ClimaComms.array_type(device)
18+
rand_data = ArrayType(rand(size(parent(f))...)) # [0-1]
19+
rand_data = rand_data .- sum(rand_data) / length(rand_data) # make centered about 0 ([-0.5:0.5])
20+
rand_data = (rand_data ./ maximum(rand_data)) .* perturb_radius # rand_data now in [-perturb_radius:perturb_radius]
21+
parent(f) .= parent(f) .+ rand_data # f in [f ± perturb_radius]
22+
return nothing
23+
end
24+
25+
@testset "Vertical mass borrowing limiter - column" begin
26+
Random.seed!(1234)
27+
z_elem = 10
28+
z_min = 0
29+
z_max = 1
30+
device = ClimaComms.device()
31+
grid = ColumnGrid(; z_elem, z_min, z_max, device)
32+
cspace = Spaces.FiniteDifferenceSpace(grid, Grids.CellCenter())
33+
tol = 0.01
34+
perturb_q = 0.3
35+
perturb_ρ = 0.2
36+
37+
# Initialize fields
38+
coords = Fields.coordinate_field(cspace)
39+
ρ = map(coord -> 1.0, coords)
40+
q = map(coord -> 0.1, coords)
41+
(; z) = coords
42+
perturb_field!(q; perturb_radius = perturb_q)
43+
perturb_field!(ρ; perturb_radius = perturb_ρ)
44+
ρq_init = ρ .⊠ q
45+
sum_ρq_init = sum(ρq_init)
46+
47+
# Test that the minimum is below 0
48+
@test length(parent(q)) == z_elem == 10
49+
@test 0.3 ≤ count(x -> x < 0, parent(q)) / 10 ≤ 0.5 # ensure reasonable percentage of points are negative
50+
51+
@test -2 * perturb_q ≤ minimum(q) ≤ -tol
52+
limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,))
53+
Limiters.apply_limiter!(q, ρ, limiter)
54+
@test 0 ≤ minimum(q)
55+
ρq = ρ .⊠ q
56+
@test isapprox(sum(ρq), sum_ρq_init; atol = 1e-15)
57+
@test isapprox(sum(ρq), sum_ρq_init; rtol = 1e-10)
58+
# @show sum(ρq) # 0.07388931313511024
59+
# @show sum_ρq_init # 0.07388931313511025
60+
end
61+
62+
@testset "Vertical mass borrowing limiter - sphere" begin
63+
z_elem = 10
64+
z_min = 0
65+
z_max = 1
66+
radius = 10
67+
h_elem = 10
68+
n_quad_points = 4
69+
grid = ExtrudedCubedSphereGrid(;
70+
z_elem,
71+
z_min,
72+
z_max,
73+
radius,
74+
h_elem,
75+
n_quad_points,
76+
)
77+
cspace = Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter())
78+
tol = 0.01
79+
perturb_q = 0.2
80+
perturb_ρ = 0.1
81+
82+
# Initialize fields
83+
coords = Fields.coordinate_field(cspace)
84+
ρ = map(coord -> 1.0, coords)
85+
q = map(coord -> 0.1, coords)
86+
87+
perturb_field!(q; perturb_radius = perturb_q)
88+
perturb_field!(ρ; perturb_radius = perturb_ρ)
89+
ρq_init = ρ .⊠ q
90+
sum_ρq_init = sum(ρq_init)
91+
92+
# Test that the minimum is below 0
93+
@test length(parent(q)) == z_elem * h_elem^2 * 6 * n_quad_points^2 == 96000
94+
@test 0.10 ≤ count(x -> x < 0.0001, parent(q)) / 96000 ≤ 1 # ensure 10% of points are negative
95+
96+
@test -2 * perturb_q ≤ minimum(q) ≤ -tol
97+
limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,))
98+
Limiters.apply_limiter!(q, ρ, limiter)
99+
@test 0 ≤ minimum(q)
100+
ρq = ρ .⊠ q
101+
@test isapprox(sum(ρq), sum_ρq_init; atol = 0.029)
102+
@test isapprox(sum(ρq), sum_ρq_init; rtol = 0.00023)
103+
# @show sum(ρq) # 125.5483524237572
104+
# @show sum_ρq_init # 125.52052986152977
105+
end
106+
107+
@testset "Vertical mass borrowing limiter - deep atmosphere" begin
108+
z_elem = 10
109+
z_min = 0
110+
z_max = 1
111+
radius = 10
112+
h_elem = 10
113+
n_quad_points = 4
114+
grid = ExtrudedCubedSphereGrid(;
115+
z_elem,
116+
z_min,
117+
z_max,
118+
radius,
119+
global_geometry = Geometry.DeepSphericalGlobalGeometry(radius),
120+
h_elem,
121+
n_quad_points,
122+
)
123+
cspace = Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter())
124+
tol = 0.01
125+
perturb_q = 0.2
126+
perturb_ρ = 0.1
127+
128+
# Initialize fields
129+
coords = Fields.coordinate_field(cspace)
130+
ρ = map(coord -> 1.0, coords)
131+
q = map(coord -> 0.1, coords)
132+
133+
perturb_field!(q; perturb_radius = perturb_q)
134+
perturb_field!(ρ; perturb_radius = perturb_ρ)
135+
ρq_init = ρ .⊠ q
136+
sum_ρq_init = sum(ρq_init)
137+
138+
# Test that the minimum is below 0
139+
@test length(parent(q)) == z_elem * h_elem^2 * 6 * n_quad_points^2 == 96000
140+
@test 0.10 ≤ count(x -> x < 0.0001, parent(q)) / 96000 ≤ 1 # ensure 10% of points are negative
141+
142+
@test -2 * perturb_q ≤ minimum(q) ≤ -tol
143+
limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,))
144+
Limiters.apply_limiter!(q, ρ, limiter)
145+
@test 0 ≤ minimum(q)
146+
ρq = ρ .⊠ q
147+
@test isapprox(sum(ρq), sum_ρq_init; atol = 0.269)
148+
@test isapprox(sum(ρq), sum_ρq_init; rtol = 0.00199)
149+
# @show sum(ρq) # 138.90494931641584
150+
# @show sum_ρq_init # 139.1735731377394
151+
end

0 commit comments

Comments
 (0)