From f1afa084a12f6fc9c6582a18d7eb1e292396a0e9 Mon Sep 17 00:00:00 2001 From: akshaysridhar Date: Fri, 5 Apr 2024 10:01:22 -0700 Subject: [PATCH] Introduce abstractions for tvd limiters Update flux correction stencils Try new TVD test Fix vanLeer --- examples/column/tvd_fct_advection.jl | 132 ++++++++++++++++++++++++ src/Operators/finitedifference.jl | 149 +++++++++++++++++++++++++++ 2 files changed, 281 insertions(+) create mode 100644 examples/column/tvd_fct_advection.jl diff --git a/examples/column/tvd_fct_advection.jl b/examples/column/tvd_fct_advection.jl new file mode 100644 index 0000000000..33d4880427 --- /dev/null +++ b/examples/column/tvd_fct_advection.jl @@ -0,0 +1,132 @@ +using Test +using LinearAlgebra +using OrdinaryDiffEq: ODEProblem, solve +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_limiters" +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) = parameters + FT = Spaces.undertype(axes(y.q)) + divf2c = Operators.DivergenceF2C( + bottom = Operators.SetValue(Geometry.WVector(FT(0))), + 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(), + ) + TVDLimitedFlux = Operators.TVDLimitedFlux( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + ) + @. yₜ.q = + -divf2c( + upwind1(w, y.q) + TVDLimitedFlux( + upwind3(w, y.q) - upwind1(w, y.q), + y.q, # Δt <-- Check factor ? + ), + ) +end + +# Define a pulse wave or square wave +pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀ + +FT = Float64 +t₀ = FT(0.0) +Δt = 0.0001 +t₁ = 10000Δt +z₀ = FT(0.0) +zₕ = FT(1.0) +z₁ = FT(1.0) +speed = FT(1.0) + +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) + 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) + 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)) + + @test err ≤ 0.25 + @test rel_mass_err ≤ 10eps() + + plot(q_final) + Plots.png( + Plots.plot!(q_analytic, title = "TVD Limited Flux"), + joinpath( + path, + "exact_and_computed_advected_square_wave_TVDLimitedFlux_" * + plot_string[i] * + ".png", + ), + ) +end diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 2a2043dd27..061db35826 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -1795,6 +1795,155 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) end +######Limited Flux Methods###### +""" + U = TVDLimitedFlux(;boundaries) + U.(A, Φ) +""" +### Here we define abstract types and specific implementations +### TVD limiters (see, for instance, Durran's notes) +abstract type AbstractTVDLimiter end +struct RZeroLimiter <: AbstractTVDLimiter end +struct RHalfLimiter <: AbstractTVDLimiter end +struct RMaxLimiter <: AbstractTVDLimiter end +struct MinModLimiter <: AbstractTVDLimiter end +struct KorenLimiter <: AbstractTVDLimiter end +struct SuperbeeLimiter <: AbstractTVDLimiter end +struct MonotonizedCentralLimiter <: AbstractTVDLimiter end +struct VanLeerLimiter <: AbstractTVDLimiter end + +@inline function compute_limiter_coeff(r, ::RZeroLimiter) + return zero(eltype(r)) +end + +@inline function compute_limiter_coeff(r, ::RHalfLimiter) + return one(eltype(r)) * 1 / 2 +end + +@inline function compute_limiter_coeff(r, ::RMaxLimiter) + return one(eltype(r)) +end + +@inline function compute_limiter_coeff(r, ::MinModLimiter) + return max(zero(eltype(r)), min(one(eltype(r)), r)) +end + +@inline function compute_limiter_coeff(r, ::KorenLimiter) + return max(zero(eltype(r)), min(2r, min(1 + 2r / 3, 2), 2)) +end + +@inline function compute_limiter_coeff(r, ::SuperbeeLimiter) + return max(zero(eltype(r)), min(one(eltype(r)), r), min(2, r)) +end + +@inline function compute_limiter_coeff(r, ::MonotonizedCentralLimiter) + return max(zero(eltype(r)), min(2r, (1 + r) / 2, 2)) +end + +@inline function compute_limiter_coeff(r, ::VanLeerLimiter) + return (r + abs(r)) / (1 + abs(r) + eps(eltype(r))) +end +# ??? Do we want to allow flux method types to be determined here? +struct TVDLimitedFlux{BCS} <: AdvectionOperator + bcs::BCS +end + +TVDLimitedFlux(; kwargs...) = TVDLimitedFlux(NamedTuple(kwargs)) + +return_eltype(::TVDLimitedFlux, A, Φ) = + Geometry.Contravariant3Vector{eltype(eltype(A))} + +return_space( + ::TVDLimitedFlux, + A_space::AllFaceFiniteDifferenceSpace, + Φ_space::AllCenterFiniteDifferenceSpace, +) = A_space + +function fct_tvd(Aⱼ₋₁₂, Aⱼ₊₁₂, Aⱼ₊₃₂, ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, rⱼ₊₁₂) + # 1/dt is in ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, ϕⱼ₋₁ᵗᵈ, ϕⱼᵗᵈ, ϕⱼ₊₁ᵗᵈ, ϕⱼ₊₂ᵗᵈ + stable_zero = zero(eltype(Aⱼ₊₁₂)) + stable_one = one(eltype(Aⱼ₊₁₂)) + # Test with various limiter methods + # Aⱼ₊₁₂ is the antidiffusive flux (see Durran textbook for notation) + Cⱼ₊₁₂ = compute_limiter_coeff(rⱼ₊₁₂, RZeroLimiter()) + return Cⱼ₊₁₂ * Aⱼ₊₁₂ +end + +stencil_interior_width(::TVDLimitedFlux, A_space, Φ_space) = + ((-1, 1), (-half - 1, half + 1)) + +Base.@propagate_inbounds function stencil_interior( + ::TVDLimitedFlux, + loc, + space, + idx, + hidx, + A_field, + Φ_field, +) + # cell center variables + ϕⱼ₋₁ = getidx(space, Φ_field, loc, idx - half - 1, hidx) + ϕⱼ = getidx(space, Φ_field, loc, idx - half, hidx) + ϕⱼ₊₁ = getidx(space, Φ_field, loc, idx + half, hidx) + ϕⱼ₊₂ = getidx(space, Φ_field, loc, idx + half + 1, hidx) + # cell face variables + Aⱼ₊₁₂ = Geometry.contravariant3( + getidx(space, A_field, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + Aⱼ₋₁₂ = Geometry.contravariant3( + getidx(space, A_field, loc, idx - 1, hidx), + Geometry.LocalGeometry(space, idx - 1, hidx), + ) + Aⱼ₊₃₂ = Geometry.contravariant3( + getidx(space, A_field, loc, idx + 1, hidx), + Geometry.LocalGeometry(space, idx + 1, hidx), + ) + # See filter options below + rⱼ₊₁₂ = compute_slope_ratio(ϕⱼ, ϕⱼ₋₁, ϕⱼ₊₁) + + return Geometry.Contravariant3Vector( + fct_tvd(Aⱼ₋₁₂, Aⱼ₊₁₂, Aⱼ₊₃₂, ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, rⱼ₊₁₂), + ) +end + +@inline function compute_slope_ratio(ϕⱼ, ϕⱼ₋₁, ϕⱼ₊₁) + return (ϕⱼ - ϕⱼ₋₁) / (ϕⱼ₊₁ - ϕⱼ + eps(eltype(ϕⱼ))) +end + + +boundary_width(::TVDLimitedFlux, ::AbstractBoundaryCondition) = 2 + +Base.@propagate_inbounds function stencil_left_boundary( + ::TVDLimitedFlux, + bc::FirstOrderOneSided, + loc, + space, + idx, + hidx, + A_field, + Φ_field, +) + @assert idx <= left_face_boundary_idx(space) + 1 + + return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) +end + +Base.@propagate_inbounds function stencil_right_boundary( + ::TVDLimitedFlux, + bc::FirstOrderOneSided, + loc, + space, + idx, + hidx, + A_field, + Φ_field, +) + @assert idx <= right_face_boundary_idx(space) - 1 + + return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) +end + """