Skip to content

Commit f40ca5e

Browse files
committed
minor fixes
1 parent 818d6c6 commit f40ca5e

19 files changed

Lines changed: 581 additions & 270 deletions

.buildkite/pipeline.yml

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -727,6 +727,17 @@ steps:
727727
agents:
728728
slurm_gpus: 1
729729

730+
- label: "Unit: FD operator (shmem)"
731+
key: unit_fd_operator_shmem
732+
retry: *retry_policy
733+
command:
734+
- "julia --color=yes --check-bounds=yes --project=.buildkite test/Operators/finitedifference/unit_fd_ops_shared_memory.jl"
735+
- "julia --color=yes --project=.buildkite test/Operators/finitedifference/benchmark_fd_ops_shared_memory.jl"
736+
env:
737+
CLIMACOMMS_DEVICE: "CUDA"
738+
agents:
739+
slurm_gpus: 1
740+
730741
- label: "Unit: gpu columnwise"
731742
key: unit_gpu_columnwise
732743
retry: *retry_policy
@@ -1798,6 +1809,7 @@ steps:
17981809
- "julia --color=yes --project=.buildkite examples/column/advect.jl"
17991810
artifact_paths:
18001811
- "examples/column/output/advect/*"
1812+
- "examples/column/output/advect_step_function/*"
18011813

18021814
- label: ":computer: Column FCT Advection Eq"
18031815
key: "cpu_fct_column_advect"

NEWS.md

Lines changed: 21 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,30 @@
11
ClimaCore.jl Release Notes
22
========================
33

4-
- Removed unused finite difference operators and boundary condition [2521](https://github.com/CliMA/ClimaCore.jl/pull/2521)
5-
- Removed `SetValue` from `GradientC2F`, `DivergenceC2F`, `CurlC2F`, `UpwindBiasedProductC2F`
6-
- Removed `SetGradient` from `InterpolateC2F` and `WeightedInterpolateC2F`
7-
- Removed `AdvectionC2C`, `AdvectionF2F`, `FluxCorrectionF2F`, and `FluxCorrectionC2C` operators
8-
94
main
105
-------
116

7+
- ![][badge-💥breaking] Removed unused finite difference operators and boundary
8+
conditions [2521](https://github.com/CliMA/ClimaCore.jl/pull/2521)
9+
- Removed `SetValue` from `GradientC2F`, `DivergenceC2F`, `CurlC2F` and
10+
`UpwindBiasedProductC2F`
11+
- Removed `SetGradient` from `InterpolateC2F` and `WeightedInterpolateC2F`
12+
- Removed the `AdvectionC2C`, `AdvectionF2F`, `FluxCorrectionC2C` and
13+
`FluxCorrectionF2F` operators
14+
15+
Each of these can be written in terms of the remaining operators and boundary
16+
conditions; `test/Operators/finitedifference/unit_column.jl` contains a
17+
testset ("Replacements for removed operators and boundary conditions") that
18+
pins the replacements against the stencils they replace. For example, a
19+
`SetValue(x₀)` boundary on `GradientC2F` is the same as
20+
`SetGradient(Covariant3Vector(2 * (x[1] - x₀)))`, and `AdvectionC2C(v, θ)` is
21+
`InterpolateF2C()(dot(Contravariant3Vector(v), GradientC2F()(θ)))`. On GPUs,
22+
the shared memory path for `GradientC2F` now applies to `SetGradient` (rather
23+
than `SetValue`) boundaries, and `InterpolateC2F` no longer lists
24+
`SetGradient` as supported. `MatrixFields.operator_matrix` now reports
25+
`LinVanLeerC2F` as a nonlinear operator instead of failing with a
26+
`MethodError`.
27+
1228
- ![][badge-🔥behavioralΔ] Unified strong/weak spectral element operator variants
1329
via a `FormType` parameter. `Divergence`, `Gradient`, and `Curl` now carry a
1430
second type parameter (`StrongForm` or `WeakForm`), and `WeakDivergence`,

docs/src/examples.md

Lines changed: 13 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -79,25 +79,29 @@ The example code solves the equation for 4 different tendencies with the followi
7979
where ``\partial`` is the [`face-to-center divergence`](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.DivergenceF2C) and $UB$ is the [`center-to-face upwind biased product`](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.UpwindBiasedProductC2F) operator.
8080
- Tendency 2:
8181

82-
$$D = \partial(UB) + \textrm{fcc}(v, \theta),$$
82+
$$D = \partial(UB) - \textrm{fcc}(v, \theta),$$
8383

84-
where $\textrm{fcc}(v, \theta)$ is the [`center-to-center flux correction`](https://github.com/CliMA/ClimaCore.jl/blob/main/src/Operators/finitedifference.jl#L2617) operator.
84+
where $\textrm{fcc}(v, \theta) = \bar{\partial}(|v| G(\theta))$ is a flux correction term, built from the [face-to-center gradient](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.GradientF2C) $\bar{\partial}$ and the [center-to-face gradient](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.GradientC2F) $G$. The gradient of $\theta$ is set to zero on both boundary faces, so that no correction flux passes through them.
8585
- Tendency 3:
8686

87-
$$D = v \cdot G(T)$$
87+
$$D = \bar{I}(v \cdot G(\theta))$$
8888

89-
where``G`` is the [center-to-face gradient](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.GradientC2F) operator, called `gradc2f` in the example code
89+
where ``G`` is the [center-to-face gradient](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.GradientC2F) operator (called `gradc2f` in the example code) and $\bar{I}$ is the [face-to-center interpolation](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.InterpolateF2C) operator.
9090
- Tendency 4:
9191

92-
$$D = A + \textrm{fcc}(v, \theta),$$
92+
$$D = \bar{I}(v \cdot G(\theta)) - \textrm{fcc}(v, \theta),$$
9393

94-
where $\textrm{fcc}(v, \theta)$ is the [`center-to-center flux correction`](https://github.com/CliMA/ClimaCore.jl/blob/main/src/Operators/finitedifference.jl#L2617) operator.
94+
combining the discretizations of tendencies 3 and 2.
9595

9696
#### Set Up
9797

98-
This test case is set up in a 1D column domain ``z \in [0, 4\pi]``, discretized into a mesh of 128 elements. The velocity field is defined as a sinusoidal wave. The boundary conditions are operator dependent, so they depend on the tendency.
99-
* For tendencies 1 and 2 where the upwind biased operator ``UB`` is used, the left boundary is defined as ``sin(a - t)``. The right boundary is ``sin(b - t)``. Here ``a`` and ``b`` are the left and right bounds of the domain.
100-
* For tendencies 3 and 4, where the advection operator ``A`` is used, the left boundary is defined as ``sin(-t)``. The right boundary is extrapolated, meaning its value is set to the closest interior point.
98+
The example runs each of the four tendencies on two initial conditions:
99+
* a smooth one, ``\theta(z, 0) = sin(z)``, on a 1D column domain ``z \in [0, 4\pi]`` discretized into a mesh of 128 elements, and
100+
* a discontinuous one, a step function, on ``z \in [-20, 20]`` discretized into a mesh of 64 elements.
101+
102+
The velocity field is constant and upward. The boundary conditions are operator dependent, so they depend on the tendency.
103+
* For tendencies 1 and 2, where the upwind biased operator ``UB`` is used, the value of ``\theta`` outside of the left boundary is ``sin(a - t)`` and outside of the right boundary is ``sin(b - t)``. Here ``a`` and ``b`` are the left and right bounds of the domain. Since `UpwindBiasedProductC2F` no longer takes a `SetValue` boundary condition, the example evaluates the upwind stencil at the boundary faces itself and imposes the result with a [`SetBoundaryOperator`](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.SetBoundaryOperator).
104+
* For tendencies 3 and 4, the gradient ``G(\theta)`` on the left boundary face is set to ``2 (\theta[1] - sin(-t))``, which is the value it takes when ``\theta = sin(-t)`` outside of the boundary. On the right boundary face it is set to the gradient of the closest interior faces, an extrapolation.
101105

102106
## 2D Cartesian examples
103107

examples/column/advect.jl

Lines changed: 112 additions & 128 deletions
Original file line numberDiff line numberDiff line change
@@ -21,30 +21,93 @@ import TerminalLoggers
2121
Logging.global_logger(TerminalLoggers.TerminalLogger())
2222
const FT = Float64
2323

24-
a = FT(0.0)
25-
b = FT(4pi)
26-
n = 128
2724
α = FT(0.1)
2825

26+
a_sin = FT(0.0)
27+
b_sin = FT(4pi)
2928
domain = Domains.IntervalDomain(
30-
Geometry.ZPoint{FT}(a),
31-
Geometry.ZPoint{FT}(b),
29+
Geometry.ZPoint{FT}(a_sin),
30+
Geometry.ZPoint{FT}(b_sin),
3231
boundary_names = (:left, :right),
3332
)
34-
mesh_sin = Meshes.IntervalMesh(domain, nelems = n)
33+
mesh_sin = Meshes.IntervalMesh(domain, nelems = 128)
3534

36-
a = FT(-20.0)
37-
b = FT(20.0)
38-
n = 64
39-
α = FT(0.1)
35+
a_step = FT(-20.0)
36+
b_step = FT(20.0)
4037
domain = Domains.IntervalDomain(
41-
Geometry.ZPoint(a),
42-
Geometry.ZPoint(b),
38+
Geometry.ZPoint(a_step),
39+
Geometry.ZPoint(b_step),
4340
boundary_names = (:left, :right),
4441
)
45-
mesh_step = Meshes.IntervalMesh(domain, nelems = n)
42+
mesh_step = Meshes.IntervalMesh(domain, nelems = 64)
43+
44+
# `UpwindBiasedProductC2F` no longer takes a `SetValue` boundary condition, so
45+
# evaluate its stencil on the boundary faces here, using the value of θ outside
46+
# of the domain, and impose the result with a `SetBoundaryOperator`.
47+
function upwind_boundary_operator(V, θ, a, b, t)
48+
lg_field = Fields.local_geometry_field(axes(V))
49+
face_bottom = Utilities.PlusHalf(0)
50+
face_top = Fields.nlevels(V) - Utilities.PlusHalf(0)
51+
v_left = Fields.field_values(
52+
Geometry.contravariant3.(
53+
Fields.level(V, face_bottom),
54+
Fields.level(lg_field, face_bottom),
55+
),
56+
)[]
57+
v_right = Fields.field_values(
58+
Geometry.contravariant3.(
59+
Fields.level(V, face_top),
60+
Fields.level(lg_field, face_top),
61+
),
62+
)[]
63+
θ_left = Fields.field_values(Fields.level(θ, 1))[]
64+
θ_right = Fields.field_values(Fields.level(θ, Fields.nlevels(θ)))[]
65+
return Operators.SetBoundaryOperator(;
66+
left = Operators.SetValue(
67+
Geometry.Contravariant3Vector(
68+
Operators.upwind_biased_product(v_left, sin(a - t), θ_left),
69+
),
70+
),
71+
right = Operators.SetValue(
72+
Geometry.Contravariant3Vector(
73+
Operators.upwind_biased_product(v_right, θ_right, sin(b - t)),
74+
),
75+
),
76+
)
77+
end
78+
79+
# The gradient on the left boundary face is the one implied by θ = sin(-t)
80+
# outside of the domain; on the right boundary face it is extrapolated from the
81+
# closest interior faces.
82+
function advection_gradient(θ, t)
83+
θ_1 = Fields.level(θ, 1)
84+
θ_n = Fields.level(θ, Fields.nlevels(θ))
85+
θ_nm1 = Fields.level(θ, Fields.nlevels(θ) - 1)
86+
return Operators.GradientC2F(
87+
left = Operators.SetGradient(
88+
@. lazy(Geometry.Covariant3Vector(2 * (θ_1 - sin(-t))))
89+
),
90+
right = Operators.SetGradient(
91+
@. lazy(Geometry.Covariant3Vector(θ_n - θ_nm1))
92+
),
93+
)
94+
end
95+
96+
# The `Extrapolate` boundary condition of the removed `FluxCorrectionC2C`
97+
# operator drops the term outside of the boundary, i.e. it sets the flux through
98+
# the boundary face -- and hence the inner gradient there -- to zero.
99+
flux_correction_gradient(::Type{FT}) where {FT} = Operators.GradientC2F(
100+
left = Operators.SetGradient(Geometry.Covariant3Vector(FT(0))),
101+
right = Operators.SetGradient(Geometry.Covariant3Vector(FT(0))),
102+
)
103+
46104
device = ClimaComms.device()
47-
for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step))
105+
for (fn, mesh, a, b) in zip(
106+
("sin", "step"),
107+
(mesh_sin, mesh_step),
108+
(a_sin, a_step),
109+
(b_sin, b_step),
110+
)
48111

49112
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
50113
fs = Spaces.FaceFiniteDifferenceSpace(cs)
@@ -61,131 +124,52 @@ for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step))
61124

62125
# Solve advection Equation: ∂θ/dt = -∂(vθ)
63126

64-
# upwinding
65-
function tendency1!(dθ, θ, _, t)
66-
lg_field = Fields.local_geometry_field(fs)
67-
lg_left = Fields.level(lg_field, Utilities.PlusHalf(0))
68-
lg_right = Fields.level(lg_field, Fields.nlevels(lg_field) - Utilities.PlusHalf(0))
69-
v_left = Fields.field_values(
70-
Geometry.contravariant3.(Fields.level(V, Utilities.PlusHalf(0)), lg_left),
71-
)[]
72-
aᴸᴮ = sin(a - t)
73-
aᴸ = Fields.field_values(Fields.level(θ, 1))[]
74-
left_bc = Operators.SetValue(
75-
Geometry.Contravariant3Vector(Operators.upwind_biased_product(v_left, aᴸᴮ, aᴸ)),
76-
)
77-
v_right = Fields.field_values(
78-
Geometry.contravariant3.(
79-
Fields.level(V, Fields.nlevels(V) - Utilities.PlusHalf(0)),
80-
lg_right,
81-
),
82-
)[]
83-
aᴿᴮ = sin(b - t)
84-
aᴿ = Fields.field_values(Fields.level(θ, Fields.nlevels(θ)))[]
85-
right_bc = Operators.SetValue(
86-
Geometry.Contravariant3Vector(
87-
Operators.upwind_biased_product(v_right, aᴿ, aᴿᴮ),
127+
lg_field = Fields.local_geometry_field(fs)
128+
= Operators.DivergenceF2C()
129+
UB = Operators.UpwindBiasedProductC2F()
130+
gradf2c = Operators.GradientF2C()
131+
interpf2c = Operators.InterpolateF2C()
132+
gradc2f_fcc = flux_correction_gradient(FT)
133+
134+
# the flux correction term, equivalent to the removed `FluxCorrectionC2C`
135+
# operator with `Extrapolate` boundary conditions
136+
fcc(θ) = @. lazy(
137+
adjoint(
138+
gradf2c(
139+
adjoint(gradc2f_fcc(θ)) * Geometry.Contravariant3Vector(
140+
abs(Geometry.contravariant3(V, lg_field)),
141+
),
88142
),
89-
)
90-
set_bcs = Operators.SetBoundaryOperator(; left = left_bc, right = right_bc)
91-
UB = Operators.UpwindBiasedProductC2F()
92-
= Operators.DivergenceF2C()
143+
) * Geometry.Contravariant3Vector(1),
144+
)
93145

146+
# upwinding
147+
function tendency1!(dθ, θ, _, t)
148+
set_bcs = upwind_boundary_operator(V, θ, a, b, t)
94149
return @.= -(set_bcs(UB(V, θ)))
95150
end
151+
# upwinding, with flux correction
96152
function tendency2!(dθ, θ, _, t)
97-
lg_field = Fields.local_geometry_field(fs)
98-
lg_left = Fields.level(lg_field, Utilities.PlusHalf(0))
99-
lg_right = Fields.level(lg_field, Fields.nlevels(lg_field) - Utilities.PlusHalf(0))
100-
v_left = Fields.field_values(
101-
Geometry.contravariant3.(Fields.level(V, Utilities.PlusHalf(0)), lg_left),
102-
)[]
103-
aᴸᴮ = sin(a - t)
104-
aᴸ = Fields.field_values(Fields.level(θ, 1))[]
105-
left_bc = Operators.SetValue(
106-
Geometry.Contravariant3Vector(Operators.upwind_biased_product(v_left, aᴸᴮ, aᴸ)),
107-
)
108-
v_right = Fields.field_values(
109-
Geometry.contravariant3.(
110-
Fields.level(V, Fields.nlevels(V) - Utilities.PlusHalf(0)),
111-
lg_right,
112-
),
113-
)[]
114-
aᴿᴮ = sin(b - t)
115-
aᴿ = Fields.field_values(Fields.level(θ, Fields.nlevels(θ)))[]
116-
right_bc = Operators.SetValue(
117-
Geometry.Contravariant3Vector(
118-
Operators.upwind_biased_product(v_right, aᴿ, aᴿᴮ),
119-
),
120-
)
121-
set_bcs = Operators.SetBoundaryOperator(; left = left_bc, right = right_bc)
122-
UB = Operators.UpwindBiasedProductC2F()
123-
= Operators.DivergenceF2C()
124-
left_center = Fields.level(θ, 1)
125-
θ_top = Fields.level(θ, Fields.nlevels(θ))
126-
θ_top_m1 = Fields.level(θ, Fields.nlevels(θ) - 1)
127-
right_center_left_biased_grad =
128-
@. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1))
129-
right_gradient_extrapolate = Operators.SetGradient(right_center_left_biased_grad)
130-
θ_2 = Fields.level(θ, 2)
131-
_left_lazy_grad2 = @. lazy(Geometry.Covariant3Vector(θ_2 - left_center))
132-
left_gradient_extrapolate = Operators.SetGradient(_left_lazy_grad2)
133-
gradc2f_fcc = Operators.GradientC2F(
134-
left = left_gradient_extrapolate,
135-
right = right_gradient_extrapolate,
136-
)
137-
gradf2c = Operators.GradientF2C()
138-
return @.=
139-
-(set_bcs(UB(V, θ))) +
140-
parent(
141-
gradf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f_fcc(θ))),
142-
).data.:1
153+
set_bcs = upwind_boundary_operator(V, θ, a, b, t)
154+
correction = fcc(θ)
155+
return @.= -(set_bcs(UB(V, θ))) + correction
143156
end
144-
# use the advection operator
157+
# advection, written as an interpolated center-to-face gradient
145158
function tendency3!(dθ, θ, _, t)
146-
left_center = Fields.level(θ, 1)
147-
θ_top = Fields.level(θ, Fields.nlevels(θ))
148-
θ_top_m1 = Fields.level(θ, Fields.nlevels(θ) - 1)
149-
right_center_left_biased_grad =
150-
@. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1))
151-
left_gradient =
152-
Operators.SetGradient(
153-
@. lazy(Geometry.Covariant3Vector(2 * (left_center - sin(-t))))
154-
)
155-
right_gradient = Operators.SetGradient(right_center_left_biased_grad)
156-
gradc2f = Operators.GradientC2F(left = left_gradient, right = right_gradient)
157-
interpf2c = Operators.InterpolateF2C()
159+
gradc2f = advection_gradient(θ, t)
158160
return @.=
159-
-1 * interpf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f(θ)))
161+
-interpf2c(
162+
Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f(θ)),
163+
)
160164
end
161-
# use the advection operator
165+
# advection, with flux correction
162166
function tendency4!(dθ, θ, _, t)
163-
left_center = Fields.level(θ, 1)
164-
θ_top = Fields.level(θ, Fields.nlevels(θ))
165-
θ_top_m1 = Fields.level(θ, Fields.nlevels(θ) - 1)
166-
right_center_left_biased_grad =
167-
@. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1))
168-
left_gradient =
169-
Operators.SetGradient(
170-
@. lazy(Geometry.Covariant3Vector(2 * (left_center - sin(-t))))
171-
)
172-
right_gradient_extrapolate = Operators.SetGradient(right_center_left_biased_grad)
173-
gradc2f =
174-
Operators.GradientC2F(left = left_gradient, right = right_gradient_extrapolate)
175-
interpf2c = Operators.InterpolateF2C()
176-
θ_2 = Fields.level(θ, 2)
177-
_left_lazy_grad4 = @. lazy(Geometry.Covariant3Vector(θ_2 - left_center))
178-
left_gradient_extrapolate = Operators.SetGradient(_left_lazy_grad4)
179-
gradc2f_fcc = Operators.GradientC2F(
180-
left = left_gradient_extrapolate,
181-
right = right_gradient_extrapolate,
182-
)
183-
gradf2c = Operators.GradientF2C()
167+
gradc2f = advection_gradient(θ, t)
168+
correction = fcc(θ)
184169
return @.=
185-
-1 * interpf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f(θ))) +
186-
parent(
187-
gradf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f_fcc(θ))),
188-
).data.:1
170+
-interpf2c(
171+
Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f(θ)),
172+
) + correction
189173
end
190174

191175
# use the advection operator

examples/column/heat.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,8 @@ T = Fields.zeros(FT, cs)
3737
# Solve Heat Equation: ∂_t T = α ∇²T
3838
function ∑tendencies!(dT, T, _, t)
3939

40-
bcs_bottom = Operators.SetValue(FT(0.0))
40+
# a Dirichlet condition, T = 0, on the bottom boundary face: the gradient
41+
# there is 2 (T[1] - 0) / Δz
4142
bottom_level_T = Fields.level(T, 1)
4243
bottom_grad = @. lazy(2.0 * Geometry.Covariant3Vector(bottom_level_T))
4344
bcs_bottom = Operators.SetGradient(bottom_grad)

0 commit comments

Comments
 (0)