Add UpwindBiasedGradient operator

Author: dennisYatunin

src/Operators/finitedifference.jl

@@ -3088,6 +3088,107 @@ Base.@propagate_inbounds function stencil_right_boundary(
     )
 end
 
+"""
+    UG = UpwindBiasedGradient()
+    UG.(v, θ)
+
+Compute the vertical gradient of the field `θ` by upwinding it according to the
+direction of a vector field `v` on the same space.
+
+More precisely, the gradient is computed based on the sign of the 3rd
+contravariant component of `v`:
+```math
+UG(\\boldsymbol{v}, θ)[i] = \\begin{cases}
+    G(LB(θ))[i] \\textrm{, if } v^3[i] > 0 \\\\
+    G(RB(θ))[i] \\textrm{, if } v^3[i] < 0
+\\end{cases}
+```
+where `G` is a gradient operator and `LB`/`RB` are left/right-bias operators.
+When `θ` and `v` are located on centers, `G = GradientF2C()`,
+`LB = LeftBiasedC2F()`, and `RB = RightBiasedC2F()`. When they are located on
+faces, `G = GradientC2F()`, `LB = LeftBiasedF2C()`, and `RB = RightBiasedF2C()`.
+
+No boundary conditions are currently supported. The default behavior on the left
+boundary is
+```math
+UG(\\boldsymbol{v}, θ)[i_min] = G(RB(θ))[i_min]
+```
+and the default behavior on the right boundary is
+```math
+UG(\\boldsymbol{v}, θ)[i_max] = G(LB(θ))[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,
+)
+    FT = Spaces.undertype(space)
+    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),
+    )
+    ∂a∂ξ₃_times_2 = (1 - sign(v)) ⊠ a⁺ + 2 * sign(v) ⊠ a - (1 + sign(v)) ⊠ a⁻
+    return Geometry.Covariant3Vector(FT(1) / 2) ⊗ ∂a∂ξ₃_times_2
+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) =

test/Operators/finitedifference/convergence_column.jl

@@ -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.FaceFiniteDifferenceSpace(center_space)
+
+        ᶜw = Geometry.WVector.(map(_ -> rand(), ones(center_space)))
+        ᶠw = Geometry.WVector.(map(_ -> rand(), 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