Skip to content

Add UpwindBiasedGradient operator #2345

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jul 9, 2025
Merged

Add UpwindBiasedGradient operator #2345

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 98 additions & 0 deletions src/Operators/finitedifference.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3088,6 +3088,104 @@ Base.@propagate_inbounds function stencil_right_boundary(
)
end

"""
UG = UpwindBiasedGradient()
UG.(v, θ)

Compute the gradient of a field `θ` by upwinding it according to the direction
of a vector field `v` on the same space. The gradient stencil is determined by
the sign of the 3rd contravariant component of `v`:
```math
UG(\\boldsymbol{v}, θ)[i] = \\begin{cases}
G(L(θ))[i] \\textrm{, if } v^3[i] > 0 \\\\
G(R(θ))[i] \\textrm{, if } v^3[i] < 0
\\end{cases}
```
where `G` is a gradient operator and `L`/`R` are left/right-bias operators. When
`θ` and `v` are located on centers, `G = GradientF2C()`, `L = LeftBiasedC2F()`,
and `R = RightBiasedC2F()`. When they are located on faces, `G = GradientC2F()`,
`L = LeftBiasedF2C()`, and `R = RightBiasedF2C()`.

No boundary conditions are currently supported. The default behavior on the left
boundary (with index `i_min`) is
```math
UG(\\boldsymbol{v}, θ)[i_min] = G(R(θ))[i_min]
```
and the default behavior on the right boundary (with index `i_max`) is
```math
UG(\\boldsymbol{v}, θ)[i_max] = G(L(θ))[i_max]
```
"""
struct UpwindBiasedGradient{BCS} <: FiniteDifferenceOperator
bcs::BCS
end
function UpwindBiasedGradient(; kwargs...)
assert_no_bcs("UpwindBiasedGradient", kwargs)
return UpwindBiasedGradient(NamedTuple())
end

return_eltype(::UpwindBiasedGradient, velocity, arg) =
Geometry.gradient_result_type(Val((3,)), eltype(arg))

return_space(
::UpwindBiasedGradient,
velocity_space::AllCenterFiniteDifferenceSpace,
arg_space::AllCenterFiniteDifferenceSpace,
) = arg_space
return_space(
::UpwindBiasedGradient,
velocity_space::AllFaceFiniteDifferenceSpace,
arg_space::AllFaceFiniteDifferenceSpace,
) = arg_space

stencil_interior_width(::UpwindBiasedGradient, velocity, arg) =
((0, 0), (-1, 1))
Base.@propagate_inbounds function stencil_interior(
::UpwindBiasedGradient,
space,
idx,
hidx,
velocity,
arg,
)
a⁺ = getidx(space, arg, idx + 1, hidx)
a = getidx(space, arg, idx, hidx)
a⁻ = getidx(space, arg, idx - 1, hidx)
v = Geometry.contravariant3(
getidx(space, velocity, idx, hidx),
Geometry.LocalGeometry(space, idx, hidx),
)
return Geometry.Covariant3Vector(1) ⊗
((1 - sign(v)) / 2 ⊠ a⁺ + sign(v) ⊠ a - (1 + sign(v)) / 2 ⊠ a⁻)
end

boundary_width(::UpwindBiasedGradient, ::AbstractBoundaryCondition) = 1
Base.@propagate_inbounds function stencil_left_boundary(
::UpwindBiasedGradient,
::NullBoundaryCondition,
space,
idx,
hidx,
arg,
)
@assert idx == left_face_boundary_idx(space)
a⁺ = getidx(space, arg, idx + 1, hidx)
a = getidx(space, arg, idx, hidx)
return Geometry.Covariant3Vector(1) ⊗ (a⁺ ⊟ a)
end
Base.@propagate_inbounds function stencil_right_boundary(
::UpwindBiasedGradient,
::NullBoundaryCondition,
space,
idx,
hidx,
arg,
)
@assert idx == right_face_boundary_idx(space)
a = getidx(space, arg, idx, hidx)
a⁻ = getidx(space, arg, idx - 1, hidx)
return Geometry.Covariant3Vector(1) ⊗ (a ⊟ a⁻)
end

abstract type DivergenceOperator <: FiniteDifferenceOperator end
return_eltype(::DivergenceOperator, arg) =
Expand Down
75 changes: 59 additions & 16 deletions test/Operators/finitedifference/convergence_column.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ convergence_rate(err, Δh) =
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

cent_field_exact = zeros(FT, cs)
face_field = zeros(FT, fs)
Expand Down Expand Up @@ -98,7 +98,7 @@ end
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

face_field_exact = zeros(FT, fs)
cent_field = zeros(FT, cs)
Expand Down Expand Up @@ -150,7 +150,7 @@ end
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

face_field_exact = Geometry.Covariant3Vector.(zeros(FT, fs))
cent_field = zeros(FT, cs)
Expand Down Expand Up @@ -202,7 +202,7 @@ end
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

cent_field_exact = Geometry.Covariant3Vector.(zeros(FT, cs))
cent_field = Geometry.Covariant3Vector.(zeros(FT, cs))
Expand Down Expand Up @@ -258,7 +258,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
faces = getproperty(Fields.coordinate_field(fs), :z)
Expand Down Expand Up @@ -409,6 +409,49 @@ end
@test conv_curl_sin_f[1] ≤ conv_curl_sin_f[2] ≤ conv_curl_sin_f[3]
end

@testset "UpwindBiasedGradient on (uniform) periodic mesh, random w" begin
FT = Float64
device = ClimaComms.device()

n_elems_seq = 2 .^ (5, 6, 7, 8)
center_errors = zeros(FT, length(n_elems_seq))
face_errors = zeros(FT, length(n_elems_seq))
Δh = zeros(FT, length(n_elems_seq))

for (k, n) in enumerate(n_elems_seq)
domain = Domains.IntervalDomain(
Geometry.ZPoint{FT}(-pi),
Geometry.ZPoint{FT}(pi);
periodic = true,
)
mesh = Meshes.IntervalMesh(domain; nelems = n)

center_space = Spaces.CenterFiniteDifferenceSpace(device, mesh)
face_space = Spaces.face_space(center_space)

ᶜw = Geometry.WVector.(map(_ -> 2 * rand() - 1, ones(center_space)))
ᶠw = Geometry.WVector.(map(_ -> 2 * rand() - 1, ones(face_space)))

ᶜz = Fields.coordinate_field(center_space).z
ᶠz = Fields.coordinate_field(face_space).z

upwind_biased_grad = Operators.UpwindBiasedGradient()
ᶜ∇sinz = Geometry.WVector.(upwind_biased_grad.(ᶜw, sin.(ᶜz)))
ᶠ∇sinz = Geometry.WVector.(upwind_biased_grad.(ᶠw, sin.(ᶠz)))

center_errors[k] = norm(ᶜ∇sinz .- Geometry.WVector.(cos.(ᶜz)))
face_errors[k] = norm(ᶠ∇sinz .- Geometry.WVector.(cos.(ᶠz)))
Δh[k] = Spaces.local_geometry_data(face_space).J[vindex(1)]
end

@test all(error -> error < 0.1, center_errors)
@test all(error -> error < 0.1, face_errors)

center_convergence_rates = convergence_rate(center_errors, Δh)
face_convergence_rates = convergence_rate(face_errors, Δh)
@test all(rate -> isapprox(rate, 1; atol = 0.01), center_convergence_rates)
@test all(rate -> isapprox(rate, 1; atol = 0.01), face_convergence_rates)
end

@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform) periodic mesh, constant w" begin
FT = Float64
Expand All @@ -428,7 +471,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)

Expand Down Expand Up @@ -479,7 +522,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
faces = getproperty(Fields.coordinate_field(fs), :z)
Expand Down Expand Up @@ -531,7 +574,7 @@ end
mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)

Expand Down Expand Up @@ -593,7 +636,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
faces = getproperty(Fields.coordinate_field(fs), :z)
Expand Down Expand Up @@ -650,7 +693,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
Expand Down Expand Up @@ -706,7 +749,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)

Expand Down Expand Up @@ -761,7 +804,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
Expand Down Expand Up @@ -817,7 +860,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
Expand Down Expand Up @@ -873,7 +916,7 @@ end
mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
Expand Down Expand Up @@ -944,7 +987,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

centers = getproperty(Fields.coordinate_field(cs), :z)
C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
Expand Down Expand Up @@ -1016,7 +1059,7 @@ end
mesh = Meshes.IntervalMesh(domain; nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
fs = Spaces.face_space(cs)

# advective velocity
c = Geometry.WVector.(ones(Float64, fs),)
Expand Down
Loading