Skip to content

Commit 21528d7

Browse files
charleskawczynskicharleskawczynski
committed
Implement and test vertical mass borrowing limiters
1 parent 68cc581 commit 21528d7

File tree

7 files changed

+370
-1
lines changed

7 files changed

+370
-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
v0.14.20
811
-------
912

‎docs/refs.bib‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -257,3 +257,14 @@ @article{Wiin1967
257257
year = {1967},
258258
publisher = {Wiley Online Library}
259259
}
260+
261+
@article{zhang2018impact,
262+
title={Impact of numerical choices on water conservation in the E3SM Atmosphere Model version 1 (EAMv1)},
263+
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},
264+
journal={Geoscientific Model Development},
265+
volume={11},
266+
number={5},
267+
pages={1971--1988},
268+
year={2018},
269+
publisher={Copernicus Publications G{\"o}ttingen, Germany}
270+
}

‎docs/src/api.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -332,14 +332,15 @@ Hypsography.ref_z_to_physical_z
332332
The limiters supertype is
333333
```@docs
334334
Limiters.AbstractLimiter
335+
Limiters.QuasiMonotoneLimiter
336+
Limiters.VerticalMassBorrowingLimiter
335337
```
336338

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

339341

340342
### Interfaces
341343
```@docs
342-
Limiters.QuasiMonotoneLimiter
343344
Limiters.compute_bounds!
344345
Limiters.apply_limiter!
345346
```

‎src/Limiters/Limiters.jl‎

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

2020
# implementations
2121
include("quasimonotone.jl")
22+
include("vertical_mass_borrowing_limiter.jl")
2223

2324
end # end module

‎src/Limiters/vertical_mass_borrowing_limiter.jl‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
import .DataLayouts as DL
2+
3+
"""
4+
VerticalMassBorrowingLimiter(f::Fields.Field, q_min)
5+
6+
A vertical-only mass borrowing limier.
7+
8+
The mass borrower borrows tracer mass from an adjacent layer.
9+
It conserves the 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, lim::VerticalMassBorrowingLimiter)
37+
38+
Apply the vertical mass borrowing
39+
limiter `lim` to field `q`.
40+
"""
41+
apply_limiter!(q::Fields.Field, lim::VerticalMassBorrowingLimiter) =
42+
apply_limiter!(q, axes(q), lim, ClimaComms.device(axes(q)))
43+
44+
function apply_limiter!(
45+
q::Fields.Field,
46+
space::Spaces.FiniteDifferenceSpace,
47+
lim::VerticalMassBorrowingLimiter,
48+
device::ClimaComms.AbstractCPUDevice,
49+
)
50+
cache = (; bmass = lim.bmass, ic = lim.ic, q_min = lim.q_min)
51+
columnwise_massborrow_cpu(q, cache)
52+
end
53+
54+
function apply_limiter!(
55+
q::Fields.Field,
56+
space::Spaces.ExtrudedFiniteDifferenceSpace,
57+
lim::VerticalMassBorrowingLimiter,
58+
device::ClimaComms.AbstractCPUDevice,
59+
)
60+
Fields.bycolumn(axes(q)) do colidx
61+
cache =
62+
(; bmass = lim.bmass[colidx], ic = lim.ic[colidx], q_min = lim.q_min)
63+
columnwise_massborrow_cpu(q[colidx], cache)
64+
end
65+
end
66+
67+
function columnwise_massborrow_cpu(q::Fields.Field, cache) # column fields
68+
(; bmass, ic, q_min) = cache
69+
70+
pdel = Fields.field_values(Fields.Δz_field(q))
71+
# loop over tracers
72+
nlevels = Spaces.nlevels(axes(q))
73+
@. ic = 0
74+
@. bmass = 0
75+
q_vals = Fields.field_values(q)
76+
# top to bottom
77+
for f in 1:DataLayouts.ncomponents(q_vals)
78+
for v in 1:nlevels
79+
CI = CartesianIndex(1, 1, f, v, 1)
80+
# new mass in the current layer
81+
pdel_v = DL.getindex_field(pdel, CI)
82+
nmass = DL.getindex_field(q_vals, CI) + bmass[] / pdel_v
83+
if nmass > q_min[f]
84+
# if new mass in the current layer is positive, don't borrow mass any more
85+
DL.setindex_field!(q_vals, nmass, CI)
86+
bmass[] = 0
87+
else
88+
# set mass to q_min in the current layer, and save bmass
89+
bmass[] = (nmass - q_min[f]) * pdel_v
90+
DL.setindex_field!(q_vals, q_min[f], CI)
91+
ic[] = ic[] + 1
92+
end
93+
end
94+
95+
# bottom to top
96+
for v in nlevels:-1:1
97+
CI = CartesianIndex(1, 1, f, v, 1)
98+
# if the surface layer still needs to borrow mass
99+
if bmass[] < 0
100+
pdel_v = DL.getindex_field(pdel, CI)
101+
# new mass in the current layer
102+
nmass = DL.getindex_field(q_vals, CI) + bmass[] / pdel_v
103+
if nmass > q_min[f]
104+
# if new mass in the current layer is positive, don't borrow mass any more
105+
DL.setindex_field!(q_vals, nmass, CI)
106+
bmass[] = 0
107+
else
108+
# if new mass in the current layer is negative, continue to borrow mass
109+
bmass[] = (nmass - q_min[f]) * pdel_v
110+
DL.setindex_field!(q_vals, q_min[f], CI)
111+
end
112+
end
113+
end
114+
end
115+
116+
return nothing
117+
end

‎test/Limiters/vertical_mass_borrowing_limiter.jl‎

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
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.Grids
10+
using ClimaCore.CommonGrids
11+
using Test
12+
using Random
13+
14+
@testset "Vertical mass borrowing limiter - column" begin
15+
Random.seed!(1234)
16+
FT = Float64
17+
z_elem = 10
18+
z_min = 0
19+
z_max = 1
20+
device = ClimaComms.device()
21+
grid = ColumnGrid(; z_elem, z_min, z_max, device)
22+
cspace = Spaces.FiniteDifferenceSpace(grid, Grids.CellCenter())
23+
context = ClimaComms.context(device)
24+
ArrayType = ClimaComms.array_type(device)
25+
tol = 0.01
26+
perturb = 0.2
27+
28+
# Initialize fields
29+
coords = Fields.coordinate_field(cspace)
30+
q = map(coord -> 0.1, coords)
31+
(; z) = coords
32+
rand_data = ArrayType(rand(size(parent(q))...)) # [0-1]
33+
rand_data = rand_data .- sum(rand_data) / length(rand_data) # make centered about 0 ([-0.5:0.5])
34+
rand_data = (rand_data ./ maximum(rand_data)) .* perturb # rand_data now in [-0.2:0.2]
35+
parent(q) .= parent(q) .+ rand_data # q in [0.1 ± 0.2]
36+
sum_q_init = sum(q)
37+
# Test that the minimum is below 0
38+
@test length(parent(q)) == z_elem == 10
39+
@test 0.3 ≤ count(x -> x < 0, parent(q)) / 10 ≤ 0.5 # ensure reasonable percentage of points are negative
40+
41+
@test -2 * perturb ≤ minimum(q) ≤ -tol
42+
limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,))
43+
Limiters.apply_limiter!(q, limiter)
44+
@test 0 ≤ minimum(q)
45+
@test isapprox(sum(q), sum_q_init; atol = 0.00000000001)
46+
@test isapprox(sum(q), sum_q_init; rtol = 0.01)
47+
end
48+
49+
@testset "Vertical mass borrowing limiter" begin
50+
FT = Float64
51+
z_elem = 10
52+
z_min = 0
53+
z_max = 1
54+
radius = 10
55+
h_elem = 10
56+
n_quad_points = 4
57+
grid = ExtrudedCubedSphereGrid(;
58+
z_elem,
59+
z_min,
60+
z_max,
61+
radius,
62+
h_elem,
63+
n_quad_points,
64+
)
65+
cspace = Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter())
66+
context = ClimaComms.context(device)
67+
ArrayType = ClimaComms.array_type(device)
68+
tol = 0.01
69+
perturb = 0.2
70+
71+
# Initialize fields
72+
coords = Fields.coordinate_field(cspace)
73+
q = map(coord -> 0.1, coords)
74+
75+
rand_data = ArrayType(rand(size(parent(q))...)) # [0-1]
76+
rand_data = rand_data .- sum(rand_data) / length(rand_data) # make centered about 0 ([-0.5:0.5])
77+
rand_data = (rand_data ./ maximum(rand_data)) .* perturb # rand_data now in [-0.2:0.2]
78+
parent(q) .= parent(q) .+ rand_data # q in [0.1 ± 0.2]
79+
sum_q_init = sum(q)
80+
81+
# Test that the minimum is below 0
82+
@test length(parent(q)) == z_elem * h_elem^2 * 6 * n_quad_points^2 == 96000
83+
@test 0.10 ≤ count(x -> x < 0.0001, parent(q)) / 96000 ≤ 1 # ensure 10% of points are negative
84+
85+
@test -0.11 ≤ minimum(q) ≤ -0.08
86+
limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,))
87+
Limiters.apply_limiter!(q, limiter)
88+
@test 0 ≤ minimum(q)
89+
@test isapprox(sum(q), sum_q_init; atol = 0.013)
90+
@test isapprox(sum(q), sum_q_init; rtol = 0.01)
91+
end

0 commit comments

Comments
 (0)