Incorrect stencil in Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F convergence test

[charleskawczynski]

I spoke with @akshaysridhar today about this FCT test, and plotted the computed vs exact solutions:

Here's the full modified script in case it's helpful:

#=
julia --check-bounds=yes --project=.buildkite
using Revise; include("../cc_fd_biased_prod.jl")
=#
using Test
using ClimaCorePlots
using LazyBroadcast: @lazy
using Plots
using StaticArrays, IntervalSets, LinearAlgebra

using ClimaComms
ClimaComms.@import_required_backends
import ClimaCore: slab, Domains, Meshes, Topologies, Spaces, Fields, Operators
import ClimaCore.Domains: Geometry
import ClimaCore.DataLayouts: vindex

device = ClimaComms.device()

function stencil_location(space, bc, idx)
    (li, lw, rw, ri) = Operators.window_bounds(space, bc)
    if !Topologies.isperiodic(Spaces.vertical_topology(space))
        if idx in li:(lw - 1)
            return Operators.LeftBoundaryWindow{Spaces.left_boundary_name(space)}()
        end
    end
    if idx in lw:rw
        return Interior()
    end
    if !Topologies.isperiodic(Spaces.vertical_topology(space))
        if idx in (rw + 1):ri
            return Operators.RightBoundaryWindow{Spaces.right_boundary_name(space)}()
        end
    end
    error("Uncaught case")
end

"""
    convergence_rate(err, Δh)

Estimate convergence rate given vectors `err` and `Δh`

    err = C Δh^p+ H.O.T
    err_k ≈ C Δh_k^p
    err_k/err_m ≈ Δh_k^p/Δh_m^p
    log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
    log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
    log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p

"""
convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

FT = Float64
n_elems_seq = 2 .^ (4, 6, 8, 10)
stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
device = ClimaComms.device()

i = 1
stretch_fn = stretch_fns[1]
err_adv_wc = zeros(FT, length(n_elems_seq))
Δh = zeros(FT, length(n_elems_seq))
Plots.plot()
for (k, n) in enumerate(n_elems_seq)
    domain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(-pi),
        Geometry.ZPoint{FT}(pi);
        boundary_names = (:bottom, :top),
    )
    mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)

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

    zc = Fields.coordinate_field(cs).z
    C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)

    # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator
    # Unitary, constant advective velocity
    w = Geometry.WVector.(ones(fs))
    # parent(w) .= parent(w) * -1
    Δz = FT(2pi / n)
    c = sin.(zc)

    zcontra3 = Geometry.Contravariant3Vector(FT(0.0))

    first_order_fluxᶠ = Operators.UpwindBiasedProductC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
        bottom = Operators.FirstOrderOneSided(),
        top = Operators.FirstOrderOneSided(),
    )

    divf2c = Operators.DivergenceF2C(
        bottom = Operators.SetValue(
            zcontra3,
        ),
        top = Operators.SetValue(
            zcontra3,
        ),
    )

    adv_wc = @lazy @. divf2c.(third_order_fluxᶠ(w, c))
    zero_bc = Operators.SetBoundaryOperator(;
        bottom = Operators.SetValue(zcontra3), top = Operators.SetValue(zcontra3)
    )
    zbc = zero_bc.(third_order_fluxᶠ.(w, c)).components.data.:1
    zbc .= zbc ./ maximum(zbc) .* maximum(c)
    Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)]

    exact = cos.(zc)
    err = parent(Fields.field_values(abs.(Base.materialize(adv_wc) .- exact)))

    if k==1
        Plots.plot!(cos.(zc); label = "exact")
    end
    Plots.plot!(Base.materialize(adv_wc); label = "computed, dh=$(round(Δh[k]; digits=4))")
    Plots.plot!(c; label = "c, dh=$(round(Δh[k]; digits=4))")
    Plots.plot!(zbc; label = "zbc, dh=$(round(Δh[k]; digits=4))")

    # Error
    err_adv_wc[k] = norm(Base.materialize(adv_wc) .- cos.(zc))
end
Plots.plot!(; legend = :topleft)
Plots.png("convergence_ubp.png")

[charleskawczynski]

If we only look at @. divf2c.(first_order_fluxᶠ(w, c)), which is UpwindBiasedProductC2F with or without extrapolate BCS (they're not used since the divergence operator is composed with it), then we have:

which also seems wrong on the boundaries, unless we somehow account for that in ClimaAtmos.

cc @szy21, @trontrytel (since this operator appears to be used in ClimaAtmos with edmf)

Actually, I think this one may be fine, only the boundary points look incorrect, but I don't think those are used when composed ( it would still be good to add unit tests to confirm this ). But the 3rd order one appears incorrect for the boundary and first interior points (i.e., first and second cell centers).

[szy21]

@dennisYatunin?

[dennisYatunin]

Could you please provide plots of c, w, and zero_bc.(first_order_fluxᶠ.(w, c)), where zero_bc = SetBoundaryOperator(; bottom = SetValue(0), top = SetValue(0))? That should help us narrow down the issue.

[charleskawczynski]

I updated the reproducer to include c. w is just a ones field. So, velocity is upward. For zero_bc, I normalized it so that it would fit nicely on the plot (via zbc .= zbc ./ maximum(zbc) .* maximum(c)):

I'm still not 100% clear on how our recursive Operators.getidx modifies stencils as it starts reaching boundaries. It seems like we only handle two cases: interior, and one point inside the boundary window, but there can be many points inside the boundary window, and a different number of points will be available depending on the idx.

[charleskawczynski]

I just met with @dennisYatunin, and we found that the BC used is reduced order, so this may actually be okay. Also, the Upwind3rdOrderBiasedProductC2F is required to be composed, so the boundary values are (known to be) invalid.

[charleskawczynski]

I'll leave this issue open until we update the convergence test to get the expected order.