From 643a85476a156b644043e7a3e646dc0f032ba7af Mon Sep 17 00:00:00 2001 From: akshaysridhar Date: Fri, 5 Apr 2024 10:01:22 -0700 Subject: [PATCH] Introduce abstractions for TVD slope limiter functions (Durran 1999) and van Leer limiters as in Lin(1994) Update numerical flux stencils to use tvd limiters Update column examples and references Update deformation flow example to use limiters Co-authored-by: Charles Kawczynski --- .buildkite/pipeline.yml | 14 + docs/refs.bib | 14 + docs/src/operators.md | 1 + examples/column/tvd_advection.jl | 192 +++++++++ examples/column/vanleer_advection.jl | 150 +++++++ examples/hybrid/sphere/deformation_flow.jl | 65 ++- src/Operators/finitedifference.jl | 396 +++++++++++++++++- .../finitedifference/convergence_column.jl | 167 ++++++++ 8 files changed, 993 insertions(+), 6 deletions(-) create mode 100644 examples/column/tvd_advection.jl create mode 100644 examples/column/vanleer_advection.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index ebbe121849..f5976713dd 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -1450,6 +1450,20 @@ 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 Lin vanLeer Limiter Advection Eq" + key: "cpu_lvl_column_advect" + command: + - "julia --color=yes --project=.buildkite examples/column/vanleer_advection.jl" + artifact_paths: + - "examples/column/output/vanleer_advection/*" - label: ":computer: Column BB FCT Advection Eq" key: "cpu_bb_fct_column_advect" diff --git a/docs/refs.bib b/docs/refs.bib index 63d542a994..fb5bab7d2c 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -124,6 +124,20 @@ @article{GubaOpt2014 doi = {10.1016/j.jcp.2014.02.029} } +@article{Lin1994, + author = {Shian-Jiann Lin and Winston C. Chao and Y. C. Sud and G. K. Walker}, + title = {A Class of the van Leer-type Transport Schemes and Its Application to the Moisture Transport in a General Circulation Model}, + journal = {Monthly Weather Review}, + year = {1994}, + publisher = {American Meteorological Society}, + address = {Boston MA, USA}, + volume = {122}, + number = {7}, + doi = {10.1175/1520-0493(1994)122<1575:ACOTVL>2.0.CO;2}, + pages= {1575 - 1593}, + url = {https://journals.ametsoc.org/view/journals/mwre/122/7/1520-0493_1994_122_1575_acotvl_2_0_co_2.xml} +} + @article{Nair2005, title = {A Discontinuous Galerkin Transport Scheme on the Cubed Sphere}, author = {Nair, Ramachandran D and Thomas, Stephen J and Loft, Richard D}, diff --git a/docs/src/operators.md b/docs/src/operators.md index 559c28aec7..c112837dc9 100644 --- a/docs/src/operators.md +++ b/docs/src/operators.md @@ -64,6 +64,7 @@ LeftBiasedC2F RightBiasedC2F LeftBiasedF2C RightBiasedF2C +AbstractTVDSlopeLimiter ``` ### Derivative operators diff --git a/examples/column/tvd_advection.jl b/examples/column/tvd_advection.jl new file mode 100644 index 0000000000..bc8a9c041c --- /dev/null +++ b/examples/column/tvd_advection.jl @@ -0,0 +1,192 @@ +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 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(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(), + ) + FCTZalesak = Operators.FCTZalesak( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + ) + TVDSlopeLimited = Operators.TVDSlopeLimitedFlux( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = limiter_method, + ) + + If = Operators.InterpolateC2F() + + if limiter_method == "Zalesak" + @. yₜ.q = + -divf2c( + upwind1(w, y.q) + FCTZalesak( + upwind3(w, y.q) - upwind1(w, y.q), + y.q / Δt, + y.q / Δt - divf2c(upwind1(w, y.q)), + ), + ) + else + Δfluxₕ = @. w * If(y.q) + Δfluxₗ = @. upwind1(w, y.q) + @. yₜ.q = + -divf2c( + upwind1(w, y.q) + + TVDSlopeLimited(upwind3(w, y.q) - upwind1(w, y.q), y.q / Δt, w), + ) + end +end + +# Define a pulse wave or square wave + +FT = Float64 +t₀ = FT(0.0) +t₁ = FT(6) +z₀ = FT(0.0) +zₕ = FT(2π) +z₁ = FT(1.0) +speed = FT(-1.0) +pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀ + +n = 2 .^ 8 +elemlist = 2 .^ [3, 4, 5, 6, 7, 8, 9, 10] +Δt = FT(0.3) * (20π / n) +@info "Timestep Δt[s]: $(Δt)" + +domain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(-10π), + Geometry.ZPoint{FT}(10π); + boundary_names = (:bottom, :top), +) + +stretch_fns = [Meshes.Uniform()] +plot_string = ["uniform"] + +for (i, stretch_fn) in enumerate(stretch_fns) + limiter_methods = ( + Operators.RZeroLimiter(), + Operators.RMaxLimiter(), + Operators.KorenLimiter(), + Operators.SuperbeeLimiter(), + Operators.MonotonizedCentralLimiter(), + "Zalesak", + ) + for (j, limiter_method) in enumerate(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)) + + if j == 1 + fig = Plots.plot(q_analytic; label = "Exact", color = :red) + end + linstyl = [:dash, :dot, :dashdot, :dashdotdot, :dash, :solid] + clrs = [:orange, :gray, :green, :maroon, :pink, :blue] + if limiter_method == "Zalesak" + fig = plot!( + q_final; + label = "Zalesak", + linestyle = linstyl[j], + color = clrs[j], + dpi = 400, + xlim = (-0.5, 1.1), + ylim = (-15, 10), + ) + else + fig = plot!( + q_final; + label = "$(typeof(limiter_method))"[21:end], + linestyle = linstyl[j], + color = clrs[j], + dpi = 400, + xlim = (-0.5, 1.1), + ylim = (-15, 10), + ) + end + fig = plot!(legend = :outerbottom, legendcolumns = 2) + if j == length(limiter_methods) + Plots.png( + fig, + joinpath( + path, + "SlopeLimitedFluxSolution_" * + "$(typeof(limiter_method))"[21:end] * + ".png", + ), + ) + end + @test err ≤ 0.25 + @test rel_mass_err ≤ 10eps() + end +end diff --git a/examples/column/vanleer_advection.jl b/examples/column/vanleer_advection.jl new file mode 100644 index 0000000000..060945ccd6 --- /dev/null +++ b/examples/column/vanleer_advection.jl @@ -0,0 +1,150 @@ +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 = "vanleer_advection" +path = joinpath(@__DIR__, "output", dir) +mkpath(path) + + +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(0))), + top = Operators.SetValue(Geometry.WVector(FT(0))), + ) + VanLeerMethod = Operators.LinVanLeerC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = limiter_method, + ) + + If = Operators.InterpolateC2F() + + @. yₜ.q = -divf2c(VanLeerMethod(w, y.q, Δt)) +end + +# Define a pulse wave or square wave + +FT = Float64 +t₀ = FT(0.0) +t₁ = FT(6) +z₀ = FT(0.0) +zₕ = FT(2π) +z₁ = FT(1.0) +speed = FT(-1.0) +pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀ + +n = 2 .^ 8 +elemlist = 2 .^ [3, 4, 5, 6, 7, 8, 9, 10] +Δt = FT(0.3) * (20π / n) +@info "Timestep Δt[s]: $(Δt)" + +domain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(-10π), + Geometry.ZPoint{FT}(10π); + 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.AlgebraicMean(), + Operators.PosDef(), + Operators.Mono4(), + Operators.Mono5(), + ) + for (j, limiter_method) in enumerate(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)) + + if j == 1 + fig = Plots.plot(q_analytic; label = "Exact", color = :red) + end + linstyl = [:solid, :dot, :dashdot, :dash] + clrs = [:orange, :blue, :green, :black] + fig = plot!( + q_final; + label = "$(typeof(limiter_method))"[21:end], + linestyle = linstyl[j], + color = clrs[j], + dpi = 400, + xlim = (-0.5, 1.1), + ylim = (-17, 12), + ) + fig = plot!(legend = :outerbottom, legendcolumns = 2) + if j == length(limiter_methods) + Plots.png( + fig, + joinpath( + path, + "LinVanLeerFluxLimiter_" * + "$(typeof(limiter_method))"[21:end] * + plot_string[i] * + ".png", + ), + ) + end + @test err ≤ 0.25 + @test rel_mass_err ≤ 10eps() + end +end diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl index 9b4ce6e5c8..d35a9d60e5 100644 --- a/examples/hybrid/sphere/deformation_flow.jl +++ b/examples/hybrid/sphere/deformation_flow.jl @@ -83,6 +83,16 @@ const FCTZalesak = Operators.FCTZalesak( bottom = Operators.FirstOrderOneSided(), top = Operators.FirstOrderOneSided(), ) +const SlopeLimitedFlux = Operators.TVDSlopeLimitedFlux( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = Operators.MinModLimiter(), +) +const LinVanLeerFlux = Operators.LinVanLeerC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = Operators.Mono5(), +) const FCTBorisBook = Operators.FCTBorisBook( bottom = Operators.FirstOrderOneSided(), top = Operators.FirstOrderOneSided(), @@ -179,6 +189,19 @@ 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), + q_n / dt, + face_uᵥ, + ) + ), + ) + elseif fct_op == LinVanLeerFlux + @. ρqₜ_n -= + vdivf2c(ᶠwinterp(ᶜJ, Y.c.ρ) * LinVanLeerFlux(face_uᵥ, q_n, dt)) else error("unrecognized FCT operator $fct_op") end @@ -306,10 +329,19 @@ 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 "vanLeer Flux Solutions" +lvl_sol = run_deformation_flow(false, LinVanLeerFlux, _dt) +lim_lvl_sol = run_deformation_flow(true, LinVanLeerFlux, _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 +418,12 @@ for (sol, suffix) in ( (lim_first_upwind_sol, "_lim_first_upwind"), (third_upwind_sol, "_third_upwind"), (fct_sol, "_fct"), + (tvd_sol, "_tvd"), + (lvl_sol, "_lvl"), (lim_third_upwind_sol, "_lim_third_upwind"), (lim_fct_sol, "_lim_fct"), + (lim_tvd_sol, "_lim_tvd"), + (lim_lvl_sol, "_lim_lvl"), ) for (sol_index, day) in ((1, 6), (2, 12)) Plots.png( @@ -400,3 +436,30 @@ 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"), + (lvl_sol, "_lvl"), + (lim_fct_sol, "_lim_fct"), + (lim_lvl_sol, "_lim_lvl"), +) + 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 diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 12740e2b28..dea0035129 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -1320,6 +1320,202 @@ Base.@propagate_inbounds function stencil_right_boundary( stencil_interior(op, loc, space, idx - 1, hidx, velocity, arg) end +struct LinVanLeerC2F{BCS} <: AdvectionOperator + bcs::BCS +end + +""" + LimiterConstraint +Following the van Leer class of limiters as noted in +[Lin1994](@cite), four limiter constraint options are provided +- AlgebraicMean: Algebraic mean, this guarantees neither positivity nor monotonicity (eq 2) +- PosDef: Positive-definite with implicit diffusion based on local stencil extrema (eq 3b, 3c, 5a, 5b) +- Mono4: Monotonicity preserving harmonic mean, this implies a strong monotonicity constraint (eq 4) +- Mono5: Monotonicity preserving, with extrema bounded by the edge cells in the stencil (eq 5) + +The diffusion implied by these methods is proportional to the local upwind CFL number. +The `mismatch` Δ𝜙 = 0 returns the first-order upwind method. Special cases (discussed in Lin et al (1994)) include setting the 𝜙_min = 0 or 𝜙_max = saturation mixing ratio for water vapor +are not considered here in favour of the generalized local extrema in equation (5a, 5b). +""" +abstract type LimiterConstraint end +struct AlgebraicMean <: LimiterConstraint end +struct PosDef <: LimiterConstraint end +struct Mono4 <: LimiterConstraint end +struct Mono5 <: LimiterConstraint end + +LinVanLeerC2F(; kwargs...) = LinVanLeerC2F(NamedTuple(kwargs)) + +return_eltype(::LinVanLeerC2F, V, A, dt) = + Geometry.Contravariant3Vector{eltype(eltype(V))} + +return_space( + ::LinVanLeerC2F, + velocity_space::AllFaceFiniteDifferenceSpace, + arg_space::AllCenterFiniteDifferenceSpace, + dt, +) = velocity_space + +function compute_Δ𝛼_linvanleer(a⁻, a⁰, a⁺, v, dt, ::Mono5) + Δ𝜙_avg = ((a⁰ - a⁻) + (a⁺ - a⁰)) / 2 + min𝜙 = min(a⁻, a⁰, a⁺) + max𝜙 = max(a⁻, a⁰, a⁺) + 𝛼 = min(abs(Δ𝜙_avg), 2 * (a⁰ - min𝜙), 2 * (max𝜙 - a⁰)) + Δ𝛼 = sign(Δ𝜙_avg) * 𝛼 * (1 - sign(v) * v * dt) +end + +function compute_Δ𝛼_linvanleer(a⁻, a⁰, a⁺, v, dt, ::Mono4) + Δ𝜙_avg = ((a⁰ - a⁻) + (a⁺ - a⁰)) / 2 + if sign(a⁰ - a⁻) == sign(a⁺ - a⁰) + return ((a⁰ - a⁻) * (a⁺ - a⁰)) / (Δ𝜙_avg + eps(eltype(v))) + else + return eltype(v)(0) + end +end + +function posdiff(x, y) + ifelse(x - y >= eltype(x)(0), x - y, eltype(x)(0)) +end +function compute_Δ𝛼_linvanleer(a⁻, a⁰, a⁺, v, dt, ::PosDef) + Δ𝜙_avg = ((a⁰ - a⁻) + (a⁺ - a⁰)) / 2 + min𝜙 = min(a⁻, a⁰, a⁺) + max𝜙 = max(a⁻, a⁰, a⁺) + return sign(Δ𝜙_avg) * + min(abs(Δ𝜙_avg), 2 * posdiff(a⁺, min𝜙), 2 * posdiff(max𝜙, a⁺)) +end + +function compute_Δ𝛼_linvanleer(a⁻, a⁰, a⁺, v, dt, ::AlgebraicMean) + return ((a⁰ - a⁻) + (a⁺ - a⁰)) / 2 +end + +function slope_limited_product(v, a⁻, a⁻⁻, a⁺, a⁺⁺, dt, method) + # Following Lin et al. (1994) + # https://doi.org/10.1175/1520-0493(1994)122<1575:ACOTVL>2.0.CO;2 + if v >= 0 + # Eqn (2,5a,5b,5c) + Δ𝛼 = compute_Δ𝛼_linvanleer(a⁻⁻, a⁻, a⁺, v, dt, method) + return v ⊠ (a⁻ ⊞ RecursiveApply.rdiv(Δ𝛼, 2)) + else + # Eqn (2,5a,5b,5c) + Δ𝛼 = compute_Δ𝛼_linvanleer(a⁻, a⁺, a⁺⁺, v, dt, method) + return v ⊠ (a⁺ ⊟ RecursiveApply.rdiv(Δ𝛼, 2)) + end +end + +stencil_interior_width(::LinVanLeerC2F, velocity, arg, dt) = + ((0, 0), (-half - 1, half + 1), (0, 0)) + +Base.@propagate_inbounds function stencil_interior( + ℱ::LinVanLeerC2F, + loc, + space, + idx, + hidx, + velocity, + arg, + dt, +) + a⁻ = getidx(space, arg, loc, idx - half, hidx) + a⁻⁻ = getidx(space, arg, loc, idx - half - 1, hidx) + a⁺ = getidx(space, arg, loc, idx + half, hidx) + a⁺⁺ = getidx(space, arg, loc, idx + half + 1, hidx) + vᶠ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + return Geometry.Contravariant3Vector( + slope_limited_product(vᶠ, a⁻, a⁻⁻, a⁺, a⁺⁺, dt, ℱ.bcs.method), + ) +end + +boundary_width(::LinVanLeerC2F, ::AbstractBoundaryCondition) = 2 + +Base.@propagate_inbounds function stencil_left_boundary( + ::LinVanLeerC2F, + bc::FirstOrderOneSided, + loc, + space, + idx, + hidx, + velocity, + arg, + dt, +) + @assert idx <= left_face_boundary_idx(space) + 1 + v = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg) + a⁺ = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺)) +end + +Base.@propagate_inbounds function stencil_right_boundary( + ::LinVanLeerC2F, + bc::FirstOrderOneSided, + loc, + space, + idx, + hidx, + velocity, + arg, + dt, +) + @assert idx >= right_face_boundary_idx(space) - 1 + v = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁻ = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg) + return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺)) + +end + +Base.@propagate_inbounds function stencil_left_boundary( + ℱ::LinVanLeerC2F, + bc::ThirdOrderOneSided, + loc, + space, + idx, + hidx, + velocity, + arg, + dt, +) + @assert idx <= left_face_boundary_idx(space) + 1 + + vᶠ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + + return Geometry.Contravariant3Vector(vᶠ * a) +end + +Base.@propagate_inbounds function stencil_right_boundary( + ℱ::LinVanLeerC2F, + bc::ThirdOrderOneSided, + loc, + space, + idx, + hidx, + velocity, + arg, + dt, +) + @assert idx <= right_face_boundary_idx(space) - 1 + + vᶠ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + + return Geometry.Contravariant3Vector(vᶠ * a) +end + """ U = Upwind3rdOrderBiasedProductC2F(;boundaries) U.(v, x) @@ -1606,9 +1802,6 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(vᶠ))) end - - -######################### """ U = FCTZalesak(;boundaries) U.(A, Φ, Φᵗᵈ) @@ -1667,6 +1860,10 @@ function fct_zalesak( stable_zero = zero(eltype(Aⱼ₊₁₂)) stable_one = one(eltype(Aⱼ₊₁₂)) + # 𝒮5.4.2 (1) Durran (5.32) Zalesak's cosmetic correction + # which is usually omitted but used in Durran's textbook + # implementation of the flux corrected transport method. + # (Textbook suggests mixed results in 3 reported scenarios) if ( Aⱼ₊₁₂ * (ϕⱼ₊₁ᵗᵈ - ϕⱼᵗᵈ) < stable_zero && ( Aⱼ₊₁₂ * (ϕⱼ₊₂ᵗᵈ - ϕⱼ₊₁ᵗᵈ) < stable_zero || @@ -1675,15 +1872,19 @@ function fct_zalesak( ) Aⱼ₊₁₂ = stable_zero end + + # 𝒮5.4.2 (2) + # If flow is nondivergent, ϕᵗᵈ are not needed in the formulae below ϕⱼᵐᵃˣ = max(ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₋₁ᵗᵈ, ϕⱼᵗᵈ, ϕⱼ₊₁ᵗᵈ) ϕⱼᵐⁱⁿ = min(ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₋₁ᵗᵈ, ϕⱼᵗᵈ, ϕⱼ₊₁ᵗᵈ) Pⱼ⁺ = max(stable_zero, Aⱼ₋₁₂) - min(stable_zero, Aⱼ₊₁₂) + # Zalesak also requires, in equation (5.33) Δx/Δt, which for the + # reference element we may assume Δζ = 1 between interfaces Qⱼ⁺ = (ϕⱼᵐᵃˣ - ϕⱼᵗᵈ) Rⱼ⁺ = (Pⱼ⁺ > stable_zero ? min(stable_one, Qⱼ⁺ / Pⱼ⁺) : stable_zero) Pⱼ⁻ = max(stable_zero, Aⱼ₊₁₂) - min(stable_zero, Aⱼ₋₁₂) Qⱼ⁻ = (ϕⱼᵗᵈ - ϕⱼᵐⁱⁿ) Rⱼ⁻ = (Pⱼ⁻ > stable_zero ? min(stable_one, Qⱼ⁻ / Pⱼ⁻) : stable_zero) - ϕⱼ₊₁ᵐᵃˣ = max(ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, ϕⱼᵗᵈ, ϕⱼ₊₁ᵗᵈ, ϕⱼ₊₂ᵗᵈ) ϕⱼ₊₁ᵐⁱⁿ = min(ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, ϕⱼᵗᵈ, ϕⱼ₊₁ᵗᵈ, ϕⱼ₊₂ᵗᵈ) Pⱼ₊₁⁺ = max(stable_zero, Aⱼ₊₁₂) - min(stable_zero, Aⱼ₊₃₂) @@ -1696,7 +1897,6 @@ function fct_zalesak( Cⱼ₊₁₂ = (Aⱼ₊₁₂ ≥ stable_zero ? min(Rⱼ₊₁⁺, Rⱼ⁻) : min(Rⱼ⁺, Rⱼ₊₁⁻)) return Cⱼ₊₁₂ * Aⱼ₊₁₂ - end stencil_interior_width(::FCTZalesak, A_space, Φ_space, Φᵗᵈ_space) = @@ -1787,7 +1987,182 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) end +""" + U = TVDSlopeLimitedFlux(;boundaries) + U.(𝒜, Φ, 𝓊) +𝒜, following the notation of Durran (Numerical Methods for Fluid +Dynamics, 2ⁿᵈ ed.) is the antidiffusive flux given by +𝒜 = ℱʰ - ℱˡ +where h and l superscripts represent the high and lower order (monotone) +fluxes respectively. The effect of the TVD limiters is then to +adjust the flux +F_{j+1/2} = F^{l}_{j+1/2} + C_{j+1/2}(F^{h}_{j+1/2} - F^{l}_{j+1/2}) +where C_{j+1/2} is the multiplicative limiter which is a function of +the ratio of the slope of the solution across a cell interface. +C=1 recovers the high order flux. +C=0 recovers the low order flux. + +Supported limiter types are +- RZeroLimiter (returns low order flux) +- RHalfLimiter (flux multiplier == 1/2) +- RMaxLimiter (returns high order flux) +- MinModLimiter +- KorenLimiter +- SuperbeeLimiter +- MonotonizedCentralLimiter +""" +abstract type AbstractTVDSlopeLimiter end +struct RZeroLimiter <: AbstractTVDSlopeLimiter end +struct RHalfLimiter <: AbstractTVDSlopeLimiter end +struct RMaxLimiter <: AbstractTVDSlopeLimiter end +struct MinModLimiter <: AbstractTVDSlopeLimiter end +struct KorenLimiter <: AbstractTVDSlopeLimiter end +struct SuperbeeLimiter <: AbstractTVDSlopeLimiter end +struct MonotonizedCentralLimiter <: AbstractTVDSlopeLimiter 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 / 3 + 2r / 3, 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 + +struct TVDSlopeLimitedFlux{BCS} <: AdvectionOperator + bcs::BCS +end + +TVDSlopeLimitedFlux(; method, kwargs...) = + TVDSlopeLimitedFlux((; method, kwargs...)) + +return_eltype(::TVDSlopeLimitedFlux, A, Φ, 𝓊) = + Geometry.Contravariant3Vector{eltype(eltype(A))} + +return_space( + ::TVDSlopeLimitedFlux, + A_space::AllFaceFiniteDifferenceSpace, + Φ_space::AllCenterFiniteDifferenceSpace, + 𝓊_space::AllFaceFiniteDifferenceSpace, +) = A_space + +function tvd_limited_flux(Aⱼ₋₁₂, Aⱼ₊₁₂, ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, rⱼ₊₁₂, method) + stable_zero = zero(eltype(Aⱼ₊₁₂)) + stable_one = one(eltype(Aⱼ₊₁₂)) + Cⱼ₊₁₂ = compute_limiter_coeff(rⱼ₊₁₂, method) + @assert Cⱼ₊₁₂ <= eltype(Aⱼ₊₁₂)(2) + @assert Cⱼ₊₁₂ >= eltype(Aⱼ₊₁₂)(0) + return Cⱼ₊₁₂ * Aⱼ₊₁₂ +end + +stencil_interior_width(::TVDSlopeLimitedFlux, A_space, Φ_space, 𝓊_space) = + ((-1, 1), (-half - 1, half + 1), (-1, +1)) + +Base.@propagate_inbounds function stencil_interior( + ℱ::TVDSlopeLimitedFlux, + loc, + space, + idx, + hidx, + A_field, + Φ_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) + 𝓊 = Geometry.contravariant3( + getidx(space, 𝓊_field, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, 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), + ) + # See filter options below + rⱼ₊₁₂ = compute_slope_ratio(ϕⱼ, ϕⱼ₋₁, ϕⱼ₊₁, ϕⱼ₊₂, 𝓊) + + return Geometry.Contravariant3Vector( + tvd_limited_flux( + Aⱼ₋₁₂, + Aⱼ₊₁₂, + ϕⱼ₋₁, + ϕⱼ, + ϕⱼ₊₁, + ϕⱼ₊₂, + rⱼ₊₁₂, + ℱ.bcs.method, + ), + ) +end + +@inline function compute_slope_ratio(ϕⱼ, ϕⱼ₋₁, ϕⱼ₊₁, ϕⱼ₊₂, 𝓊) + if 𝓊 >= 0 + return (ϕⱼ - ϕⱼ₋₁) / (ϕⱼ₊₁ - ϕⱼ + eps(eltype(ϕⱼ))) + else + return (ϕⱼ₊₂ - ϕⱼ₊₁) / (ϕⱼ₊₁ - ϕⱼ + eps(eltype(ϕⱼ))) + end +end + +boundary_width(::TVDSlopeLimitedFlux, ::AbstractBoundaryCondition) = 2 + +Base.@propagate_inbounds function stencil_left_boundary( + ::TVDSlopeLimitedFlux, + bc::FirstOrderOneSided, + loc, + space, + idx, + hidx, + A_field, + Φ_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( + ::TVDSlopeLimitedFlux, + bc::FirstOrderOneSided, + loc, + space, + idx, + hidx, + A_field, + Φ_field, + 𝓊_field, +) + @assert idx <= right_face_boundary_idx(space) - 1 + return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) +end """ A = AdvectionF2F(;boundaries) @@ -3441,3 +3816,14 @@ Base.@propagate_inbounds function apply_stencil!( end return field_out end +# Compute slope ratio 𝜃 and limiter coefficient 𝜙 +#𝜃 = compute_slope_ratio(a⁻, a⁻⁻, a⁺, a⁺⁺, v) +#𝜙 = compute_limiter_coeff(𝜃, method) + + +#@assert 0 <= 𝜙 <= 2 +#if v >= 0 +# return v ⊠ (a⁻ ⊞ RecursiveApply.rdiv((a⁺ - a⁻) ⊠ 𝜙 ,2)) +#else +# return v ⊠ (a⁺ ⊟ RecursiveApply.rdiv((a⁺ - a⁻) ⊠ 𝜙 ,2)) # Current working solution +#end diff --git a/test/Operators/finitedifference/convergence_column.jl b/test/Operators/finitedifference/convergence_column.jl index 58ac0c078b..54370762da 100644 --- a/test/Operators/finitedifference/convergence_column.jl +++ b/test/Operators/finitedifference/convergence_column.jl @@ -688,6 +688,173 @@ end @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2] end +@testset "Lin et al. (1994) van Leer class limiter (Mono5)" begin + FT = Float64 + n_elems_seq = 2 .^ (5, 6, 7, 8, 9, 10) + + err_adv_wc = zeros(FT, length(n_elems_seq)) + + Δh = zeros(FT, length(n_elems_seq)) + device = ClimaComms.device() + + 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) + + cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) + fs = Spaces.FaceFiniteDifferenceSpace(cs) + + centers = getproperty(Fields.coordinate_field(cs), :z) + + # Unitary, constant advective velocity + w = Geometry.WVector.(ones(fs)) + # c = sin(z), scalar field defined at the centers + c = sin.(centers) + + SLMethod = Operators.LinVanLeerC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = Operators.Mono5(), + ) + + divf2c = Operators.DivergenceF2C() + flux = SLMethod.(w, c, FT(0)) + adv_wc = divf2c.(flux) + + Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)] + + # Error + err_adv_wc[k] = norm(adv_wc .- cos.(centers)) + end + + # Check convergence rate + conv_adv_wc = convergence_rate(err_adv_wc, Δh) + + # LinVanLeer limited flux conv, with f(z) = sin(z) + @test conv_adv_wc[1] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[2] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[3] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[4] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[5] ≈ 1.5 atol = 0.01 + +end + +@testset "Lin et al. (1994) van Leer class limiter (Mono4)" begin + FT = Float64 + n_elems_seq = 2 .^ (5, 6, 7, 8, 9, 10) + + err_adv_wc = zeros(FT, length(n_elems_seq)) + + Δh = zeros(FT, length(n_elems_seq)) + device = ClimaComms.device() + + 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) + + 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) + + # Unitary, constant advective velocity + w = Geometry.WVector.(ones(fs)) + # c = sin(z), scalar field defined at the centers + c = sin.(centers) + + SLMethod = Operators.LinVanLeerC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = Operators.Mono4(), + ) + + divf2c = Operators.DivergenceF2C() + flux = SLMethod.(w, c, FT(0)) + adv_wc = divf2c.(flux) + + Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)] + + # Error + err_adv_wc[k] = norm(adv_wc .- cos.(centers)) + end + + # Check convergence rate + conv_adv_wc = convergence_rate(err_adv_wc, Δh) + + # LinVanLeer limited flux conv, with f(z) = sin(z) + @test conv_adv_wc[1] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[2] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[3] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[4] ≈ 1.5 atol = 0.01 + @test conv_adv_wc[5] ≈ 1.5 atol = 0.01 + +end + +@testset "Lin et al. (1994) van Leer class limiter (PosDef)" begin + FT = Float64 + n_elems_seq = 2 .^ (5, 6, 7, 8, 9, 10) + + err_adv_wc = zeros(FT, length(n_elems_seq)) + + Δh = zeros(FT, length(n_elems_seq)) + device = ClimaComms.device() + + 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) + + 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) + + # Unitary, constant advective velocity + w = Geometry.WVector.(ones(fs)) + # c = sin(z), scalar field defined at the centers + c = sin.(centers) + + SLMethod = Operators.LinVanLeerC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = Operators.PosDef(), + ) + + divf2c = Operators.DivergenceF2C() + flux = SLMethod.(w, c, FT(0)) + adv_wc = divf2c.(flux) + + Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)] + + # Error + err_adv_wc[k] = norm(adv_wc .- cos.(centers)) + end + + # Check convergence rate + conv_adv_wc = convergence_rate(err_adv_wc, Δh) + + # LinVanLeer limited flux conv, with f(z) = sin(z) + @test conv_adv_wc[1] ≈ 1.0 atol = 0.01 + @test conv_adv_wc[2] ≈ 1.0 atol = 0.01 + @test conv_adv_wc[3] ≈ 1.0 atol = 0.01 + @test conv_adv_wc[4] ≈ 1.0 atol = 0.01 + @test conv_adv_wc[5] ≈ 1.0 atol = 0.01 + +end + @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)