Skip to content

Commit 122ee93

Browse files
dennisYatunindennisYatunin
authored
Merge pull request #2345 from CliMA/dy/upwind_gradient
Add UpwindBiasedGradient operator
2 parents 5a4af69 + e34abc1 commit 122ee93

File tree

2 files changed

+157
-16
lines changed

2 files changed

+157
-16
lines changed

‎src/Operators/finitedifference.jl

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3088,6 +3088,104 @@ Base.@propagate_inbounds function stencil_right_boundary(
30883088
)
30893089
end
30903090

3091+
"""
3092+
UG = UpwindBiasedGradient()
3093+
UG.(v, θ)
3094+
3095+
Compute the gradient of a field `θ` by upwinding it according to the direction
3096+
of a vector field `v` on the same space. The gradient stencil is determined by
3097+
the sign of the 3rd contravariant component of `v`:
3098+
```math
3099+
UG(\\boldsymbol{v}, θ)[i] = \\begin{cases}
3100+
G(L(θ))[i] \\textrm{, if } v^3[i] > 0 \\\\
3101+
G(R(θ))[i] \\textrm{, if } v^3[i] < 0
3102+
\\end{cases}
3103+
```
3104+
where `G` is a gradient operator and `L`/`R` are left/right-bias operators. When
3105+
`θ` and `v` are located on centers, `G = GradientF2C()`, `L = LeftBiasedC2F()`,
3106+
and `R = RightBiasedC2F()`. When they are located on faces, `G = GradientC2F()`,
3107+
`L = LeftBiasedF2C()`, and `R = RightBiasedF2C()`.
3108+
3109+
No boundary conditions are currently supported. The default behavior on the left
3110+
boundary (with index `i_min`) is
3111+
```math
3112+
UG(\\boldsymbol{v}, θ)[i_min] = G(R(θ))[i_min]
3113+
```
3114+
and the default behavior on the right boundary (with index `i_max`) is
3115+
```math
3116+
UG(\\boldsymbol{v}, θ)[i_max] = G(L(θ))[i_max]
3117+
```
3118+
"""
3119+
struct UpwindBiasedGradient{BCS} <: FiniteDifferenceOperator
3120+
bcs::BCS
3121+
end
3122+
function UpwindBiasedGradient(; kwargs...)
3123+
assert_no_bcs("UpwindBiasedGradient", kwargs)
3124+
return UpwindBiasedGradient(NamedTuple())
3125+
end
3126+
3127+
return_eltype(::UpwindBiasedGradient, velocity, arg) =
3128+
Geometry.gradient_result_type(Val((3,)), eltype(arg))
3129+
3130+
return_space(
3131+
::UpwindBiasedGradient,
3132+
velocity_space::AllCenterFiniteDifferenceSpace,
3133+
arg_space::AllCenterFiniteDifferenceSpace,
3134+
) = arg_space
3135+
return_space(
3136+
::UpwindBiasedGradient,
3137+
velocity_space::AllFaceFiniteDifferenceSpace,
3138+
arg_space::AllFaceFiniteDifferenceSpace,
3139+
) = arg_space
3140+
3141+
stencil_interior_width(::UpwindBiasedGradient, velocity, arg) =
3142+
((0, 0), (-1, 1))
3143+
Base.@propagate_inbounds function stencil_interior(
3144+
::UpwindBiasedGradient,
3145+
space,
3146+
idx,
3147+
hidx,
3148+
velocity,
3149+
arg,
3150+
)
3151+
a⁺ = getidx(space, arg, idx + 1, hidx)
3152+
a = getidx(space, arg, idx, hidx)
3153+
a⁻ = getidx(space, arg, idx - 1, hidx)
3154+
v = Geometry.contravariant3(
3155+
getidx(space, velocity, idx, hidx),
3156+
Geometry.LocalGeometry(space, idx, hidx),
3157+
)
3158+
return Geometry.Covariant3Vector(1) ⊗
3159+
((1 - sign(v)) / 2 ⊠ a⁺ + sign(v) ⊠ a - (1 + sign(v)) / 2 ⊠ a⁻)
3160+
end
3161+
3162+
boundary_width(::UpwindBiasedGradient, ::AbstractBoundaryCondition) = 1
3163+
Base.@propagate_inbounds function stencil_left_boundary(
3164+
::UpwindBiasedGradient,
3165+
::NullBoundaryCondition,
3166+
space,
3167+
idx,
3168+
hidx,
3169+
arg,
3170+
)
3171+
@assert idx == left_face_boundary_idx(space)
3172+
a⁺ = getidx(space, arg, idx + 1, hidx)
3173+
a = getidx(space, arg, idx, hidx)
3174+
return Geometry.Covariant3Vector(1) ⊗ (a⁺ ⊟ a)
3175+
end
3176+
Base.@propagate_inbounds function stencil_right_boundary(
3177+
::UpwindBiasedGradient,
3178+
::NullBoundaryCondition,
3179+
space,
3180+
idx,
3181+
hidx,
3182+
arg,
3183+
)
3184+
@assert idx == right_face_boundary_idx(space)
3185+
a = getidx(space, arg, idx, hidx)
3186+
a⁻ = getidx(space, arg, idx - 1, hidx)
3187+
return Geometry.Covariant3Vector(1) ⊗ (a ⊟ a⁻)
3188+
end
30913189

30923190
abstract type DivergenceOperator <: FiniteDifferenceOperator end
30933191
return_eltype(::DivergenceOperator, arg) =

‎test/Operators/finitedifference/convergence_column.jl

Lines changed: 59 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ convergence_rate(err, Δh) =
5050
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)
5151

5252
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
53-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
53+
fs = Spaces.face_space(cs)
5454

5555
cent_field_exact = zeros(FT, cs)
5656
face_field = zeros(FT, fs)
@@ -98,7 +98,7 @@ end
9898
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)
9999

100100
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
101-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
101+
fs = Spaces.face_space(cs)
102102

103103
face_field_exact = zeros(FT, fs)
104104
cent_field = zeros(FT, cs)
@@ -150,7 +150,7 @@ end
150150
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)
151151

152152
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
153-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
153+
fs = Spaces.face_space(cs)
154154

155155
face_field_exact = Geometry.Covariant3Vector.(zeros(FT, fs))
156156
cent_field = zeros(FT, cs)
@@ -202,7 +202,7 @@ end
202202
mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)
203203

204204
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
205-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
205+
fs = Spaces.face_space(cs)
206206

207207
cent_field_exact = Geometry.Covariant3Vector.(zeros(FT, cs))
208208
cent_field = Geometry.Covariant3Vector.(zeros(FT, cs))
@@ -258,7 +258,7 @@ end
258258
mesh = Meshes.IntervalMesh(domain; nelems = n)
259259

260260
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
261-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
261+
fs = Spaces.face_space(cs)
262262

263263
centers = getproperty(Fields.coordinate_field(cs), :z)
264264
faces = getproperty(Fields.coordinate_field(fs), :z)
@@ -409,6 +409,49 @@ end
409409
@test conv_curl_sin_f[1] ≤ conv_curl_sin_f[2] ≤ conv_curl_sin_f[3]
410410
end
411411

412+
@testset "UpwindBiasedGradient on (uniform) periodic mesh, random w" begin
413+
FT = Float64
414+
device = ClimaComms.device()
415+
416+
n_elems_seq = 2 .^ (5, 6, 7, 8)
417+
center_errors = zeros(FT, length(n_elems_seq))
418+
face_errors = zeros(FT, length(n_elems_seq))
419+
Δh = zeros(FT, length(n_elems_seq))
420+
421+
for (k, n) in enumerate(n_elems_seq)
422+
domain = Domains.IntervalDomain(
423+
Geometry.ZPoint{FT}(-pi),
424+
Geometry.ZPoint{FT}(pi);
425+
periodic = true,
426+
)
427+
mesh = Meshes.IntervalMesh(domain; nelems = n)
428+
429+
center_space = Spaces.CenterFiniteDifferenceSpace(device, mesh)
430+
face_space = Spaces.face_space(center_space)
431+
432+
ᶜw = Geometry.WVector.(map(_ -> 2 * rand() - 1, ones(center_space)))
433+
ᶠw = Geometry.WVector.(map(_ -> 2 * rand() - 1, ones(face_space)))
434+
435+
ᶜz = Fields.coordinate_field(center_space).z
436+
ᶠz = Fields.coordinate_field(face_space).z
437+
438+
upwind_biased_grad = Operators.UpwindBiasedGradient()
439+
ᶜ∇sinz = Geometry.WVector.(upwind_biased_grad.(ᶜw, sin.(ᶜz)))
440+
ᶠ∇sinz = Geometry.WVector.(upwind_biased_grad.(ᶠw, sin.(ᶠz)))
441+
442+
center_errors[k] = norm(ᶜ∇sinz .- Geometry.WVector.(cos.(ᶜz)))
443+
face_errors[k] = norm(ᶠ∇sinz .- Geometry.WVector.(cos.(ᶠz)))
444+
Δh[k] = Spaces.local_geometry_data(face_space).J[vindex(1)]
445+
end
446+
447+
@test all(error -> error < 0.1, center_errors)
448+
@test all(error -> error < 0.1, face_errors)
449+
450+
center_convergence_rates = convergence_rate(center_errors, Δh)
451+
face_convergence_rates = convergence_rate(face_errors, Δh)
452+
@test all(rate -> isapprox(rate, 1; atol = 0.01), center_convergence_rates)
453+
@test all(rate -> isapprox(rate, 1; atol = 0.01), face_convergence_rates)
454+
end
412455

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

430473
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
431-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
474+
fs = Spaces.face_space(cs)
432475

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

@@ -479,7 +522,7 @@ end
479522
mesh = Meshes.IntervalMesh(domain; nelems = n)
480523

481524
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
482-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
525+
fs = Spaces.face_space(cs)
483526

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

533576
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
534-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
577+
fs = Spaces.face_space(cs)
535578

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

@@ -593,7 +636,7 @@ end
593636
mesh = Meshes.IntervalMesh(domain; nelems = n)
594637

595638
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
596-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
639+
fs = Spaces.face_space(cs)
597640

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

652695
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
653-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
696+
fs = Spaces.face_space(cs)
654697

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

708751
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
709-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
752+
fs = Spaces.face_space(cs)
710753

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

@@ -761,7 +804,7 @@ end
761804
mesh = Meshes.IntervalMesh(domain; nelems = n)
762805

763806
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
764-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
807+
fs = Spaces.face_space(cs)
765808

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

819862
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
820-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
863+
fs = Spaces.face_space(cs)
821864

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

875918
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
876-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
919+
fs = Spaces.face_space(cs)
877920

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

946989
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
947-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
990+
fs = Spaces.face_space(cs)
948991

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

10181061
cs = Spaces.CenterFiniteDifferenceSpace(device, mesh)
1019-
fs = Spaces.FaceFiniteDifferenceSpace(cs)
1062+
fs = Spaces.face_space(cs)
10201063

10211064
# advective velocity
10221065
c = Geometry.WVector.(ones(Float64, fs),)

0 commit comments

Comments
 (0)