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..ac097cbf4b 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.1 + @test conv_adv_wc[2] ≈ 1.5 atol = 0.1 + @test conv_adv_wc[3] ≈ 1.5 atol = 0.1 + @test conv_adv_wc[4] ≈ 1.5 atol = 0.1 + @test conv_adv_wc[5] ≈ 1.5 atol = 0.1 + +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.1 + @test conv_adv_wc[2] ≈ 1.5 atol = 0.1 + @test conv_adv_wc[3] ≈ 1.5 atol = 0.1 + @test conv_adv_wc[4] ≈ 1.5 atol = 0.1 + @test conv_adv_wc[5] ≈ 1.5 atol = 0.1 + +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.1 + @test conv_adv_wc[2] ≈ 1.0 atol = 0.1 + @test conv_adv_wc[3] ≈ 1.0 atol = 0.1 + @test conv_adv_wc[4] ≈ 1.0 atol = 0.1 + @test conv_adv_wc[5] ≈ 1.0 atol = 0.1 + +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)