Introduce abstractions for tvd limiters

Update numerical flux stencils to use tvd limiters
Update column examples
Update deformation flow example
+ format

.buildkite/pipeline.yml

@@ -1450,6 +1450,13 @@ steps:
           - "julia --color=yes --project=.buildkite examples/column/fct_advection.jl"
         artifact_paths:
           - "examples/column/output/fct_advection/*"
+
+      - label: ":computer: Column TVD Slope-limited Advection Eq"
+        key: "cpu_tvd_column_advect"
+        command:
+          - "julia --color=yes --project=.buildkite examples/column/tvd_advection.jl"
+        artifact_paths:
+          - "examples/column/output/tvd_advection/*"
 
       - label: ":computer: Column BB FCT Advection Eq"
         key: "cpu_bb_fct_column_advect"

docs/src/operators.md

@@ -64,6 +64,7 @@ LeftBiasedC2F
 RightBiasedC2F
 LeftBiasedF2C
 RightBiasedF2C
+AbstractTVDSlopeLimiter
 ```
 
 ### Derivative operators

examples/column/tvd_advection.jl

@@ -0,0 +1,156 @@
+using Test
+using LinearAlgebra
+import ClimaComms
+ClimaComms.@import_required_backends
+using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
+using ClimaTimeSteppers
+
+import ClimaCore:
+    Fields,
+    Domains,
+    Topologies,
+    Meshes,
+    DataLayouts,
+    Operators,
+    Geometry,
+    Spaces
+
+
+# Advection Equation, with constant advective velocity (so advection form = flux form)
+# ∂_t y + w ∂_z y  = 0
+# the solution translates to the right at speed w,
+# so at time t, the solution is y(z - w * t)
+
+# visualization artifacts
+ENV["GKSwstype"] = "nul"
+using ClimaCorePlots, Plots
+Plots.GRBackend()
+dir = "tvd_advection"
+path = joinpath(@__DIR__, "output", dir)
+mkpath(path)
+
+function linkfig(figpath, alt = "")
+    # buildkite-agent upload figpath
+    # link figure in logs if we are running on CI
+    if get(ENV, "BUILDKITE", "") == "true"
+        artifact_url = "artifact://$figpath"
+        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
+    end
+end
+
+function tendency!(yₜ, y, parameters, t)
+    (; w, Δt, limiter_method) = parameters
+    FT = Spaces.undertype(axes(y.q))
+    bcvel = pulse(-π, t, z₀, zₕ, z₁)
+    divf2c = Operators.DivergenceF2C(
+        bottom = Operators.SetValue(Geometry.WVector(FT(bcvel))),
+        top = Operators.SetValue(Geometry.WVector(FT(0))),
+    )
+    upwind1 = Operators.UpwindBiasedProductC2F(
+        bottom = Operators.Extrapolate(),
+        top = Operators.Extrapolate(),
+    )
+    upwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
+        bottom = Operators.ThirdOrderOneSided(),
+        top = Operators.ThirdOrderOneSided(),
+    )
+    LimitedFlux = Operators.TVDSlopeLimitedFlux(
+        bottom = Operators.FirstOrderOneSided(),
+        top = Operators.FirstOrderOneSided(),
+        method = limiter_method,
+    )
+    If = Operators.InterpolateC2F()
+    @. yₜ.q =
+        -divf2c(
+            upwind1(w, y.q) +
+            LimitedFlux(upwind3(w, y.q) - upwind1(w, y.q), y.q / Δt),
+        )
+    # @. yₜ.q =
+    #     -divf2c(w * If(y.q))
+end
+
+# Define a pulse wave or square wave
+
+FT = Float64
+t₀ = FT(0.0)
+Δt = 0.0001 * 25
+t₁ = 200Δt
+z₀ = FT(0.0)
+zₕ = FT(1.0)
+z₁ = FT(1.0)
+speed = FT(-1.0)
+pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀
+
+n = 2 .^ 6
+
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(-π),
+    Geometry.ZPoint{FT}(π);
+    boundary_names = (:bottom, :top),
+)
+
+stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(FT(7.0)))
+plot_string = ["uniform", "stretched"]
+
+for (i, stretch_fn) in enumerate(stretch_fns)
+    limiter_methods = (
+        Operators.RZeroLimiter(),
+        Operators.RHalfLimiter(),
+        Operators.RMaxLimiter(),
+        Operators.MinModLimiter(),
+        Operators.KorenLimiter(),
+        Operators.SuperbeeLimiter(),
+        Operators.MonotonizedCentralLimiter(),
+        Operators.VanLeerLimiter(),
+    )
+    for limiter_method in limiter_methods
+        @info (limiter_method, stretch_fn)
+        mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n)
+        cent_space = Spaces.CenterFiniteDifferenceSpace(mesh)
+        face_space = Spaces.FaceFiniteDifferenceSpace(cent_space)
+        z = Fields.coordinate_field(cent_space).z
+        O = ones(FT, face_space)
+
+        # Initial condition
+        q_init = pulse.(z, 0.0, z₀, zₕ, z₁)
+        y = Fields.FieldVector(q = q_init)
+
+        # Unitary, constant advective velocity
+        w = Geometry.WVector.(speed .* O)
+
+        # Solve the ODE
+        parameters = (; w, Δt, limiter_method)
+        prob = ODEProblem(
+            ClimaODEFunction(; T_exp! = tendency!),
+            y,
+            (t₀, t₁),
+            parameters,
+        )
+        sol = solve(
+            prob,
+            ExplicitAlgorithm(SSP33ShuOsher()),
+            dt = Δt,
+            saveat = Δt,
+        )
+
+        q_final = sol.u[end].q
+        q_analytic = pulse.(z, t₁, z₀, zₕ, z₁)
+        err = norm(q_final .- q_analytic)
+        rel_mass_err = norm((sum(q_final) - sum(q_init)) / sum(q_init))
+
+
+        plot(q_final)
+        Plots.png(
+            Plots.plot!(q_analytic, title = "$(typeof(limiter_method))"),
+            joinpath(
+                path,
+                "exact_and_computed_advected_square_wave_TVDSlopeLimitedFlux_" *
+                "$(typeof(limiter_method))_" *
+                plot_string[i] *
+                ".png",
+            ),
+        )
+        @test err ≤ 0.25
+        @test rel_mass_err ≤ 10eps()
+    end
+end

examples/hybrid/sphere/deformation_flow.jl

@@ -83,6 +83,11 @@ const FCTZalesak = Operators.FCTZalesak(
     bottom = Operators.FirstOrderOneSided(),
     top = Operators.FirstOrderOneSided(),
 )
+const SlopeLimitedFlux = Operators.TVDSlopeLimitedFlux(
+    bottom = Operators.FirstOrderOneSided(),
+    top = Operators.FirstOrderOneSided(),
+    method = Operators.MinModLimiter(),
+)
 const FCTBorisBook = Operators.FCTBorisBook(
     bottom = Operators.FirstOrderOneSided(),
     top = Operators.FirstOrderOneSided(),
@@ -179,6 +184,15 @@ function vertical_tendency!(Yₜ, Y, cache, t)
                     )
                 ),
             )
+        elseif fct_op == SlopeLimitedFlux
+            @. ρqₜ_n -= vdivf2c(
+                ᶠwinterp(ᶜJ, Y.c.ρ) * (
+                    upwind1(face_uᵥ, q_n) + SlopeLimitedFlux(
+                        upwind3(face_uᵥ, q_n) - upwind1(face_uᵥ, q_n),
+                        Y.c.ρ * q_n / dt,
+                    )
+                ),
+            )
         else
             error("unrecognized FCT operator $fct_op")
         end
@@ -306,10 +320,16 @@ tracer_ranges(sol) =
         return maximum(q_n) - minimum(q_n)
     end
 
+@info "Slope Limited Solutions"
+tvd_sol = run_deformation_flow(false, SlopeLimitedFlux, _dt)
+lim_tvd_sol = run_deformation_flow(true, SlopeLimitedFlux, _dt)
+@info "Third-Order Upwind Solutions"
 third_upwind_sol = run_deformation_flow(false, upwind3, _dt)
-fct_sol = run_deformation_flow(false, FCTZalesak, _dt)
 lim_third_upwind_sol = run_deformation_flow(true, upwind3, _dt)
+@info "Zalesak Flux-Corrected Transport Solutions"
+fct_sol = run_deformation_flow(false, FCTZalesak, _dt)
 lim_fct_sol = run_deformation_flow(true, FCTZalesak, _dt)
+@info "First-Order Upwind Solutions"
 lim_first_upwind_sol = run_deformation_flow(true, upwind1, _dt)
 lim_centered_sol = run_deformation_flow(true, nothing, _dt)
 
@@ -386,8 +406,10 @@ for (sol, suffix) in (
     (lim_first_upwind_sol, "_lim_first_upwind"),
     (third_upwind_sol, "_third_upwind"),
     (fct_sol, "_fct"),
+    (tvd_sol, "_tvd"),
     (lim_third_upwind_sol, "_lim_third_upwind"),
     (lim_fct_sol, "_lim_fct"),
+    (lim_tvd_sol, "_lim_tvd"),
 )
     for (sol_index, day) in ((1, 6), (2, 12))
         Plots.png(
@@ -400,3 +422,29 @@ for (sol, suffix) in (
         )
     end
 end
+
+for (sol, suffix) in (
+    (lim_centered_sol, "_lim_centered"),
+    (lim_first_upwind_sol, "_lim_first_upwind"),
+    (third_upwind_sol, "_third_upwind"),
+    (fct_sol, "_fct"),
+    (tvd_sol, "_tvd"),
+    (lim_fct_sol, "_lim_fct"),
+    (lim_tvd_sol, "_lim_tvd"),
+)
+    for (sol_index, day) in ((1, 6), (2, 12))
+        Plots.png(
+            Plots.plot(
+                (
+                    ((sol.u[sol_index].c.ρq.:3) ./ sol.u[sol_index].c.ρ) .- (
+                        lim_third_upwind_sol[sol_index].c.ρq.:3 ./
+                        lim_third_upwind_sol[sol_index].c.ρ
+                    )
+                ),
+                level = 15,
+                clim = (-1, 1),
+            ),
+            joinpath(path, "q3_day_diff_$day$suffix.png"),
+        )
+    end
+end