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 --project
using Revise; include("test/Operators/finitedifference/convergence_column.jl")
=#
using Test
using ClimaCorePlots
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()

"""
    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)]

@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with FirstOrderOneSided BCs" begin
    FT = Float64
    n_elems_seq = 2 .^ (4, 6, 8, 10)
    stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
    device = ClimaComms.device()

    # for (i, stretch_fn) in enumerate(stretch_fns)
        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)

            centers = getproperty(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))
            # c = sin(z), scalar field defined at the centers
            Δz = FT(2pi / n)
            c = (cos.(centers .- Δz / 2) .- cos.(centers .+ Δz / 2)) ./ Δz
            s = sin.(centers)

            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(
                    Geometry.Contravariant3Vector(FT(0.0)),
                ),
                top = Operators.SetValue(
                    Geometry.Contravariant3Vector(FT(0.0)),
                ),
            )

            corrected_antidiff_flux = @. divf2c(
                C * (third_order_fluxᶠ(w, c) - first_order_fluxᶠ(w, c)),
            )
            adv_wc =
                @. divf2c.(first_order_fluxᶠ(w, c)) + corrected_antidiff_flux

            Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)]

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

            # Error
            err_adv_wc[k] = norm(adv_wc .- cos.(centers))
        end
        Plots.plot!(; legend = :topright)
        Plots.png("convergence_fct.png")

        # Check convergence rate
        conv_adv_wc = convergence_rate(err_adv_wc, Δh)
        # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z)
        @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 0.2006
        @test conv_adv_wc[1] ≈ 0.5 atol = 0.2
        @test conv_adv_wc[2] ≈ 0.5 atol = 0.3
        @test conv_adv_wc[3] ≈ 1.0 atol = 0.55
    end
# end