Skip to content

Commit 13f55b2

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

File tree

4 files changed

+336
-0
lines changed

4 files changed

+336
-0
lines changed

‎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: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
#=
2+
3+
From E3SM:
4+
https://github.com/E3SM-Project/E3SM/blob/2c377c5ec9a5585170524b366ad85074ab1b1a5c/components/eam/src/physics/cam/massborrow.F90
5+
Described in Appendix B of
6+
https://gmd.copernicus.org/articles/11/1971/2018/gmd-11-1971-2018.pdf
7+
8+
=#
9+
10+
import .DataLayouts as DL
11+
12+
struct VerticalMassBorrowingLimiter{F, T}
13+
bmass::F
14+
ic::F
15+
qmin::T
16+
end
17+
function VerticalMassBorrowingLimiter(f::Fields.Field, qmin)
18+
bmass = similar(Spaces.level(f, 1))
19+
ic = similar(Spaces.level(f, 1))
20+
return VerticalMassBorrowingLimiter(bmass, ic, qmin)
21+
end
22+
23+
apply_limiter!(q::Fields.Field, lim::VerticalMassBorrowingLimiter) =
24+
apply_limiter!(q, axes(q), lim, ClimaComms.device(axes(q)))
25+
26+
function apply_limiter!(
27+
q::Fields.Field,
28+
space::Spaces.FiniteDifferenceSpace,
29+
lim::VerticalMassBorrowingLimiter,
30+
device::ClimaComms.AbstractCPUDevice,
31+
)
32+
cache = (; bmass = lim.bmass, ic = lim.ic, qmin = lim.qmin)
33+
columnwise_massborrow_cpu(q, cache)
34+
end
35+
36+
function apply_limiter!(
37+
q::Fields.Field,
38+
space::Spaces.ExtrudedFiniteDifferenceSpace,
39+
lim::VerticalMassBorrowingLimiter,
40+
device::ClimaComms.AbstractCPUDevice,
41+
)
42+
Fields.bycolumn(axes(q)) do colidx
43+
cache =
44+
(; bmass = lim.bmass[colidx], ic = lim.ic[colidx], qmin = lim.qmin)
45+
columnwise_massborrow_cpu(q[colidx], cache)
46+
end
47+
end
48+
49+
function columnwise_massborrow_cpu(q::Fields.Field, cache) # column fields
50+
(; bmass, ic, qmin) = cache
51+
52+
pdel = Fields.field_values(Fields.Δz_field(q))
53+
# loop over tracers
54+
nlevels = Spaces.nlevels(axes(q))
55+
@. ic = 0
56+
@. bmass = 0
57+
q_vals = Fields.field_values(q)
58+
# top to bottom
59+
for f in 1:DataLayouts.ncomponents(q_vals)
60+
for v in 1:nlevels
61+
CI = CartesianIndex(1, 1, f, v, 1)
62+
# new mass in the current layer
63+
pdel_v = DL.getindex_field(pdel, CI)
64+
nmass = DL.getindex_field(q_vals, CI) + bmass[] / pdel_v
65+
if nmass > qmin[f]
66+
# if new mass in the current layer is positive, don't borrow mass any more
67+
DL.setindex_field!(q_vals, nmass, CI)
68+
bmass[] = 0
69+
else
70+
# set mass to qmin in the current layer, and save bmass
71+
bmass[] = (nmass - qmin[f]) * pdel_v
72+
DL.setindex_field!(q_vals, qmin[f], CI)
73+
ic[] = ic[] + 1
74+
end
75+
end
76+
77+
# bottom to top
78+
for v in nlevels:-1:1
79+
CI = CartesianIndex(1, 1, f, v, 1)
80+
# if the surface layer still needs to borrow mass
81+
if bmass[] < 0
82+
pdel_v = DL.getindex_field(pdel, CI)
83+
# new mass in the current layer
84+
nmass = DL.getindex_field(q_vals, CI) + bmass[] / pdel_v
85+
if nmass > qmin[f]
86+
# if new mass in the current layer is positive, don't borrow mass any more
87+
DL.setindex_field!(q_vals, nmass, CI)
88+
bmass[] = 0
89+
else
90+
# if new mass in the current layer is negative, continue to borrow mass
91+
bmass[] = (nmass - qmin[f]) * pdel_v
92+
DL.setindex_field!(q_vals, qmin[f], CI)
93+
end
94+
end
95+
end
96+
end
97+
98+
return nothing
99+
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

‎test/Limiters/vertical_mass_borrowing_limiter_advection.jl

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
#=
2+
julia --project=.buildkite
3+
using Revise; include(joinpath("test", "Limiters", "vertical_mass_borrowing_limiter_advection.jl"))
4+
=#
5+
using Test
6+
using LinearAlgebra
7+
import ClimaComms
8+
ClimaComms.@import_required_backends
9+
using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
10+
using ClimaTimeSteppers
11+
12+
import ClimaCore:
13+
Fields,
14+
Domains,
15+
Limiters,
16+
Topologies,
17+
Meshes,
18+
DataLayouts,
19+
Operators,
20+
Geometry,
21+
Spaces
22+
23+
24+
# Advection Equation, with constant advective velocity (so advection form = flux form)
25+
# ∂_t y + w ∂_z y = 0
26+
# the solution translates to the right at speed w,
27+
# so at time t, the solution is y(z - w * t)
28+
29+
# visualization artifacts
30+
ENV["GKSwstype"] = "nul"
31+
using ClimaCorePlots, Plots
32+
Plots.GRBackend()
33+
dir = "vert_mass_borrow_advection"
34+
path = joinpath(@__DIR__, "output", dir)
35+
mkpath(path)
36+
37+
function lim!(y, parameters, t, y_ref)
38+
(; w, Δt, limiter) = parameters
39+
Limiters.apply_limiter!(y.q, limiter)
40+
return nothing
41+
end
42+
43+
function tendency!(yₜ, y, parameters, t)
44+
(; w, Δt, limiter) = parameters
45+
FT = Spaces.undertype(axes(y.q))
46+
bcvel = pulse(-π, t, z₀, zₕ, z₁)
47+
divf2c = Operators.DivergenceF2C(
48+
bottom = Operators.SetValue(Geometry.WVector(FT(bcvel))),
49+
top = Operators.SetValue(Geometry.WVector(FT(0))),
50+
)
51+
upwind1 = Operators.UpwindBiasedProductC2F(
52+
bottom = Operators.Extrapolate(),
53+
top = Operators.Extrapolate(),
54+
)
55+
upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
56+
bottom = Operators.ThirdOrderOneSided(),
57+
top = Operators.ThirdOrderOneSided(),
58+
)
59+
If = Operators.InterpolateC2F()
60+
@. yₜ.q =
61+
# -divf2c(w * If(y.q))
62+
-divf2c(upwind1(w, y.q) * If(y.q))
63+
# -divf2c(upwind3(w, y.q) * If(y.q))
64+
return nothing
65+
end
66+
67+
# Define a pulse wave or square wave
68+
69+
FT = Float64
70+
t₀ = FT(0.0)
71+
Δt = 0.0001 * 25
72+
t₁ = 200Δt
73+
z₀ = FT(0.0)
74+
zₕ = FT(1.0)
75+
z₁ = FT(1.0)
76+
speed = FT(-1.0)
77+
pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀
78+
79+
n = 2 .^ 6
80+
81+
domain = Domains.IntervalDomain(
82+
Geometry.ZPoint{FT}(-π),
83+
Geometry.ZPoint{FT}(π);
84+
boundary_names = (:bottom, :top),
85+
)
86+
87+
# stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(FT(7.0)))
88+
stretch_fns = (Meshes.Uniform(),)
89+
plot_string = ["uniform", "stretched"]
90+
91+
@testset "VerticalMassBorrowingLimiter on advection" begin
92+
for (i, stretch_fn) in enumerate(stretch_fns)
93+
mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)
94+
cent_space = Spaces.CenterFiniteDifferenceSpace(mesh)
95+
face_space = Spaces.FaceFiniteDifferenceSpace(cent_space)
96+
z = Fields.coordinate_field(cent_space).z
97+
O = ones(FT, face_space)
98+
99+
# Initial condition
100+
q_init = pulse.(z, 0.0, z₀, zₕ, z₁)
101+
q = q_init
102+
y = Fields.FieldVector(; q)
103+
limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,))
104+
105+
# Unitary, constant advective velocity
106+
w = Geometry.WVector.(speed .* O)
107+
108+
# Solve the ODE
109+
parameters = (; w, Δt, limiter)
110+
prob = ODEProblem(
111+
ClimaODEFunction(; lim!, T_lim! = tendency!),
112+
y,
113+
(t₀, t₁),
114+
parameters,
115+
)
116+
sol = solve(
117+
prob,
118+
ExplicitAlgorithm(SSP33ShuOsher()),
119+
dt = Δt,
120+
saveat = Δt,
121+
)
122+
123+
q_init = sol.u[1].q
124+
q_final = sol.u[end].q
125+
q_analytic = pulse.(z, t₁, z₀, zₕ, z₁)
126+
err = norm(q_final .- q_analytic)
127+
rel_mass_err = norm((sum(q_final) - sum(q_init)) / sum(q_init))
128+
129+
130+
p = plot()
131+
Plots.plot!(q_init, label = "init")
132+
Plots.plot!(q_final, label = "computed")
133+
Plots.plot!(q_analytic, label = "analytic")
134+
Plots.plot!(; legend = :topright)
135+
Plots.plot!(; xlabel = "q", title = "VerticalMassBorrowingLimiter")
136+
f = joinpath(
137+
path,
138+
"VerticalMassBorrowingLimiter_comparison_$(plot_string[i]).png",
139+
)
140+
Plots.png(p, f)
141+
@test err ≤ 0.25
142+
@test rel_mass_err ≤ 10eps()
143+
@test minimum(q_final) ≥ 0
144+
end
145+
end

0 commit comments

Comments
 (0)