Skip to content

Commit 406bfeb

Browse files
Merge pull request #1662 from CliMA/as/limiter-template
TVD slope limiters
2 parents 6539b89 + 084f693 commit 406bfeb

File tree

8 files changed

+1110
-29
lines changed

8 files changed

+1110
-29
lines changed

.buildkite/pipeline.yml

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1454,6 +1454,20 @@ steps:
14541454
- "julia --color=yes --project=.buildkite examples/column/fct_advection.jl"
14551455
artifact_paths:
14561456
- "examples/column/output/fct_advection/*"
1457+
1458+
- label: ":computer: Column TVD Slope-limited Advection Eq"
1459+
key: "cpu_tvd_column_advect"
1460+
command:
1461+
- "julia --color=yes --project=.buildkite examples/column/tvd_advection.jl"
1462+
artifact_paths:
1463+
- "examples/column/output/tvd_advection/*"
1464+
1465+
- label: ":computer: Column Lin vanLeer Limiter Advection Eq"
1466+
key: "cpu_lvl_column_advect"
1467+
command:
1468+
- "julia --color=yes --project=.buildkite examples/column/vanleer_advection.jl"
1469+
artifact_paths:
1470+
- "examples/column/output/vanleer_advection/*"
14571471

14581472
- label: ":computer: Column BB FCT Advection Eq"
14591473
key: "cpu_bb_fct_column_advect"

docs/refs.bib

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -124,6 +124,20 @@ @article{GubaOpt2014
124124
doi = {10.1016/j.jcp.2014.02.029}
125125
}
126126

127+
@article{Lin1994,
128+
author = {Shian-Jiann Lin and Winston C. Chao and Y. C. Sud and G. K. Walker},
129+
title = {A Class of the van Leer-type Transport Schemes and Its Application to the Moisture Transport in a General Circulation Model},
130+
journal = {Monthly Weather Review},
131+
year = {1994},
132+
publisher = {American Meteorological Society},
133+
address = {Boston MA, USA},
134+
volume = {122},
135+
number = {7},
136+
doi = {10.1175/1520-0493(1994)122<1575:ACOTVL>2.0.CO;2},
137+
pages= {1575 - 1593},
138+
url = {https://journals.ametsoc.org/view/journals/mwre/122/7/1520-0493_1994_122_1575_acotvl_2_0_co_2.xml}
139+
}
140+
127141
@article{Nair2005,
128142
title = {A Discontinuous Galerkin Transport Scheme on the Cubed Sphere},
129143
author = {Nair, Ramachandran D and Thomas, Stephen J and Loft, Richard D},

docs/src/operators.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@ LeftBiasedC2F
6464
RightBiasedC2F
6565
LeftBiasedF2C
6666
RightBiasedF2C
67+
AbstractTVDSlopeLimiter
6768
```
6869

6970
### Derivative operators

examples/column/tvd_advection.jl

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
using Test
2+
using LinearAlgebra
3+
import ClimaComms
4+
ClimaComms.@import_required_backends
5+
using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
6+
using ClimaTimeSteppers
7+
8+
import ClimaCore:
9+
Fields,
10+
Domains,
11+
Topologies,
12+
Meshes,
13+
DataLayouts,
14+
Operators,
15+
Geometry,
16+
Spaces
17+
18+
19+
# Advection Equation, with constant advective velocity (so advection form = flux form)
20+
# ∂_t y + w ∂_z y = 0
21+
# the solution translates to the right at speed w,
22+
# so at time t, the solution is y(z - w * t)
23+
24+
# visualization artifacts
25+
ENV["GKSwstype"] = "nul"
26+
using ClimaCorePlots, Plots
27+
Plots.GRBackend()
28+
dir = "tvd_advection"
29+
path = joinpath(@__DIR__, "output", dir)
30+
mkpath(path)
31+
32+
33+
function tendency!(yₜ, y, parameters, t)
34+
(; w, Δt, limiter_method) = parameters
35+
FT = Spaces.undertype(axes(y.q))
36+
bcvel = pulse(-π, t, z₀, zₕ, z₁)
37+
divf2c = Operators.DivergenceF2C(
38+
bottom = Operators.SetValue(Geometry.WVector(FT(0))),
39+
top = Operators.SetValue(Geometry.WVector(FT(0))),
40+
)
41+
upwind1 = Operators.UpwindBiasedProductC2F(
42+
bottom = Operators.Extrapolate(),
43+
top = Operators.Extrapolate(),
44+
)
45+
upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
46+
bottom = Operators.ThirdOrderOneSided(),
47+
top = Operators.ThirdOrderOneSided(),
48+
)
49+
FCTZalesak = Operators.FCTZalesak(
50+
bottom = Operators.FirstOrderOneSided(),
51+
top = Operators.FirstOrderOneSided(),
52+
)
53+
TVDSlopeLimited = Operators.TVDLimitedFluxC2F(
54+
bottom = Operators.FirstOrderOneSided(),
55+
top = Operators.FirstOrderOneSided(),
56+
method = limiter_method,
57+
)
58+
59+
If = Operators.InterpolateC2F()
60+
61+
if limiter_method == "Zalesak"
62+
@. yₜ.q =
63+
-divf2c(
64+
upwind1(w, y.q) + FCTZalesak(
65+
upwind3(w, y.q) - upwind1(w, y.q),
66+
y.q / Δt,
67+
y.q / Δt - divf2c(upwind1(w, y.q)),
68+
),
69+
)
70+
else
71+
Δfluxₕ = @. w * If(y.q)
72+
Δfluxₗ = @. upwind1(w, y.q)
73+
@. yₜ.q =
74+
-divf2c(
75+
upwind1(w, y.q) +
76+
TVDSlopeLimited(upwind3(w, y.q) - upwind1(w, y.q), y.q / Δt, w),
77+
)
78+
end
79+
end
80+
81+
# Define a pulse wave or square wave
82+
83+
FT = Float64
84+
t₀ = FT(0.0)
85+
t₁ = FT(6)
86+
z₀ = FT(0.0)
87+
zₕ = FT(2π)
88+
z₁ = FT(1.0)
89+
speed = FT(-1.0)
90+
pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) zₕ ? z₁ : z₀
91+
92+
n = 2 .^ 8
93+
elemlist = 2 .^ [3, 4, 5, 6, 7, 8, 9, 10]
94+
Δt = FT(0.3) * (20π / n)
95+
@info "Timestep Δt[s]: $(Δt)"
96+
97+
domain = Domains.IntervalDomain(
98+
Geometry.ZPoint{FT}(-10π),
99+
Geometry.ZPoint{FT}(10π);
100+
boundary_names = (:bottom, :top),
101+
)
102+
103+
stretch_fns = [Meshes.Uniform()]
104+
plot_string = ["uniform"]
105+
106+
for (i, stretch_fn) in enumerate(stretch_fns)
107+
limiter_methods = (
108+
Operators.RZeroLimiter(),
109+
Operators.RMaxLimiter(),
110+
Operators.KorenLimiter(),
111+
Operators.SuperbeeLimiter(),
112+
Operators.MonotonizedCentralLimiter(),
113+
"Zalesak",
114+
)
115+
for (j, limiter_method) in enumerate(limiter_methods)
116+
@info (limiter_method, stretch_fn)
117+
mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)
118+
cent_space = Spaces.CenterFiniteDifferenceSpace(mesh)
119+
face_space = Spaces.FaceFiniteDifferenceSpace(cent_space)
120+
z = Fields.coordinate_field(cent_space).z
121+
O = ones(FT, face_space)
122+
123+
# Initial condition
124+
q_init = pulse.(z, 0.0, z₀, zₕ, z₁)
125+
y = Fields.FieldVector(q = q_init)
126+
127+
# Unitary, constant advective velocity
128+
w = Geometry.WVector.(speed .* O)
129+
130+
# Solve the ODE
131+
parameters = (; w, Δt, limiter_method)
132+
prob = ODEProblem(
133+
ClimaODEFunction(; T_exp! = tendency!),
134+
y,
135+
(t₀, t₁),
136+
parameters,
137+
)
138+
sol = solve(
139+
prob,
140+
ExplicitAlgorithm(SSP33ShuOsher()),
141+
dt = Δt,
142+
saveat = Δt,
143+
)
144+
145+
q_final = sol.u[end].q
146+
q_analytic = pulse.(z, t₁, z₀, zₕ, z₁)
147+
148+
err = norm(q_final .- q_analytic)
149+
rel_mass_err = norm((sum(q_final) - sum(q_init)) / sum(q_init))
150+
151+
if j == 1
152+
fig = Plots.plot(q_analytic; label = "Exact", color = :red)
153+
end
154+
linstyl = [:dash, :dot, :dashdot, :dashdotdot, :dash, :solid]
155+
clrs = [:orange, :gray, :green, :maroon, :pink, :blue]
156+
if limiter_method == "Zalesak"
157+
fig = plot!(
158+
q_final;
159+
label = "Zalesak",
160+
linestyle = linstyl[j],
161+
color = clrs[j],
162+
dpi = 400,
163+
xlim = (-0.5, 1.1),
164+
ylim = (-15, 10),
165+
)
166+
else
167+
fig = plot!(
168+
q_final;
169+
label = "$(typeof(limiter_method))"[21:end],
170+
linestyle = linstyl[j],
171+
color = clrs[j],
172+
dpi = 400,
173+
xlim = (-0.5, 1.1),
174+
ylim = (-15, 10),
175+
)
176+
end
177+
fig = plot!(legend = :outerbottom, legendcolumns = 2)
178+
if j == length(limiter_methods)
179+
Plots.png(
180+
fig,
181+
joinpath(
182+
path,
183+
"SlopeLimitedFluxSolution_" *
184+
"$(typeof(limiter_method))"[21:end] *
185+
".png",
186+
),
187+
)
188+
end
189+
@test err 0.25
190+
@test rel_mass_err 10eps()
191+
end
192+
end

0 commit comments

Comments
 (0)