From 084f693c51ee504284021fdf466fdd396e1ff45b 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 modified: src/Operators/finitedifference.jl Standard symbols : \scru_space -> u_space Move docstring modified: examples/column/vanleer_advection.jl modified: src/Operators/finitedifference.jl modified: examples/column/tvd_advection.jl modified: examples/hybrid/sphere/deformation_flow.jl modified: src/Operators/finitedifference.jl verbose op+method names modified: examples/column/vanleer_advection.jl Reduce cfl and show eps convergence for mono4, mono5 Remove some eltype conversions Fix name Docs formatting Apply julia formatter Update domain extent Fix names, move constraint outside of BCs Added and updated docs More fixes + info statements method -> constraint ; new kwarg Try Float64 Print info in failing tracer test Update factor in deformation flow tracer condition --- .buildkite/pipeline.yml | 14 + docs/refs.bib | 14 + docs/src/operators.md | 1 + examples/column/tvd_advection.jl | 192 ++++++++ examples/column/vanleer_advection.jl | 168 +++++++ examples/hybrid/sphere/deformation_flow.jl | 122 ++++- src/Operators/finitedifference.jl | 461 +++++++++++++++++- .../finitedifference/convergence_column.jl | 167 +++++++ 8 files changed, 1110 insertions(+), 29 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 ddf9251ce0..ba77571e45 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -1454,6 +1454,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..da86eac1c9 --- /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.TVDLimitedFluxC2F( + 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..3d2105c87a --- /dev/null +++ b/examples/column/vanleer_advection.jl @@ -0,0 +1,168 @@ +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(), + constraint = 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.1) * (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) + @info stretch_fn + limiter_methods = ( + Operators.AlgebraicMean(), + Operators.PositiveDefinite(), + Operators.MonotoneHarmonic(), + Operators.MonotoneLocalExtrema(), + ) + for (j, limiter_method) in enumerate(limiter_methods) + 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 + + @info "Extrema with $(limiter_method), i=$i, j=$j: $(extrema(q_final))" + @show maximum(q_final .- 1) + @show minimum(q_final .- 0) + @show abs(maximum(q_final .- 1)) + monotonicity_preserving = + [Operators.MonotoneHarmonic, Operators.MonotoneLocalExtrema] + if any(x -> limiter_method isa x, monotonicity_preserving) && + stretch_fn == Meshes.Uniform() + @assert abs(maximum(q_final .- 1)) <= eps(FT) + @assert abs(minimum(q_final .- 0)) <= eps(FT) + @assert maximum(q_final) <= FT(1) + elseif limiter_method != Operators.AlgebraicMean() + @assert abs(maximum(q_final .- 1)) <= FT(0.05) + @assert abs(minimum(q_final .- 0)) <= FT(0.05) + @assert maximum(q_final) <= FT(1) + end + + 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 = (-20, 20), + ) + 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..b5adbaeab1 100644 --- a/examples/hybrid/sphere/deformation_flow.jl +++ b/examples/hybrid/sphere/deformation_flow.jl @@ -1,3 +1,7 @@ +#= +julia --project=.buildkite +using Revise; include("examples/hybrid/sphere/deformation_flow.jl") +=# import ClimaComms ClimaComms.@import_required_backends using SciMLBase: ODEProblem, init, solve @@ -26,7 +30,7 @@ const context = ClimaComms.SingletonCommsContext() # http://www-personal.umich.edu/~cjablono/DCMIP-2012_TestCaseDocument_v1.7.pdf, # Section 1.1 -const FT = Float32 # floating point type +const FT = Float64 # floating point type const R = FT(6.37122e6) # radius const grav = FT(9.8) # gravitational constant const R_d = FT(287.058) # R dry (gas constant / mol mass dry air) @@ -83,6 +87,16 @@ const FCTZalesak = Operators.FCTZalesak( bottom = Operators.FirstOrderOneSided(), top = Operators.FirstOrderOneSided(), ) +const SlopeLimitedFlux = Operators.TVDLimitedFluxC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + method = Operators.MinModLimiter(), +) +const LinVanLeerFlux = Operators.LinVanLeerC2F( + bottom = Operators.FirstOrderOneSided(), + top = Operators.FirstOrderOneSided(), + constraint = Operators.MonotoneLocalExtrema(), +) const FCTBorisBook = Operators.FCTBorisBook( bottom = Operators.FirstOrderOneSided(), top = Operators.FirstOrderOneSided(), @@ -179,6 +193,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 +333,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) @@ -353,28 +389,35 @@ max_q5_roundoff_err = 2 * eps(FT) @test lim_centered_ρq_errs[5] ≈ third_upwind_ρ_err atol = max_q5_roundoff_err # Check that the different upwinding modes with the limiter improve the "smoothness" of the tracers. -@test all(tracer_roughnesses(fct_sol) .< tracer_roughnesses(third_upwind_sol)) -@test all( - tracer_roughnesses(lim_third_upwind_sol) .< - 0.9 .* tracer_roughnesses(third_upwind_sol), -) -@test all( - tracer_roughnesses(lim_fct_sol) .< - 0.8 .* tracer_roughnesses(third_upwind_sol), -) -@test all(tracer_ranges(fct_sol) .< tracer_ranges(third_upwind_sol)) -@test all( - tracer_ranges(lim_third_upwind_sol) .< - 0.6 .* tracer_ranges(third_upwind_sol), -) -@test all(tracer_ranges(lim_fct_sol) .< 0.55 .* tracer_ranges(third_upwind_sol)) -@test all( - tracer_ranges(lim_first_upwind_sol) .< - 0.5 .* tracer_ranges(third_upwind_sol), -) -@test all( - tracer_ranges(lim_centered_sol) .< 0.9 .* tracer_ranges(third_upwind_sol), -) +@testset "Test tracer properties" begin + @test all( + tracer_roughnesses(fct_sol) .< tracer_roughnesses(third_upwind_sol), + ) + @test all( + tracer_roughnesses(lim_third_upwind_sol) .< + 0.99 .* tracer_roughnesses(third_upwind_sol), + ) + @test all( + tracer_roughnesses(lim_fct_sol) .< + 0.8 .* tracer_roughnesses(third_upwind_sol), + ) + @test all(tracer_ranges(fct_sol) .< tracer_ranges(third_upwind_sol)) + @test all( + tracer_ranges(lim_third_upwind_sol) .< + 0.6 .* tracer_ranges(third_upwind_sol), + ) + @test all( + tracer_ranges(lim_fct_sol) .< 0.55 .* tracer_ranges(third_upwind_sol), + ) + @test all( + tracer_ranges(lim_first_upwind_sol) .< + 0.5 .* tracer_ranges(third_upwind_sol), + ) + @test all( + tracer_ranges(lim_centered_sol) .< + 0.9 .* tracer_ranges(third_upwind_sol), + ) +end ENV["GKSwstype"] = "nul" using ClimaCorePlots, Plots @@ -386,8 +429,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 +447,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..f9c8e38018 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -1320,6 +1320,212 @@ Base.@propagate_inbounds function stencil_right_boundary( stencil_interior(op, loc, space, idx - 1, hidx, velocity, arg) end +""" + LinVanLeerC2F + +Following the van Leer class of limiters as noted in[Lin1994](@cite), four +limiter constraint options are provided for use with advection operators: + +- `AlgebraicMean`: Algebraic mean, this guarantees neither positivity nor + monotonicity (eq 2, `avg`) +- `PositiveDefinite`: Positive-definite with implicit diffusion based on local + stencil extrema (eq 3b, 3c, 5a, 5b, `posd`) +- `MonotoneHarmonic`: Monotonicity preserving harmonic mean, this implies a strong + monotonicity constraint (eq 4, `mono4`) +- `MonotoneLocalExtrema`: Monotonicity preserving, with extrema bounded by the + edge cells in the stencil (eq 5, `mono5`) + +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). +""" +struct LinVanLeerC2F{BCS, C} <: AdvectionOperator + bcs::BCS + constraint::C +end +abstract type LimiterConstraint end +struct AlgebraicMean <: LimiterConstraint end +struct PositiveDefinite <: LimiterConstraint end +struct MonotoneHarmonic <: LimiterConstraint end +struct MonotoneLocalExtrema <: LimiterConstraint end + +LinVanLeerC2F(; constraint, kwargs...) = + LinVanLeerC2F(NamedTuple(kwargs), constraint) + +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, ::MonotoneLocalExtrema) + Δ𝜙_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, ::MonotoneHarmonic) + Δ𝜙_avg = ((a⁰ - a⁻) + (a⁺ - a⁰)) / 2 + c = sign(v) * v * dt + if sign(a⁰ - a⁻) == sign(a⁺ - a⁰) && Δ𝜙_avg != 0 + return ((a⁰ - a⁻) * (a⁺ - a⁰)) / (Δ𝜙_avg) * (1 - c) + else + return eltype(v)(0) + end +end + +posdiff(x, y) = ifelse(x - y ≥ 0, x - y, eltype(x)(0)) + +function compute_Δ𝛼_linvanleer(a⁻, a⁰, a⁺, v, dt, ::PositiveDefinite) + Δ𝜙_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⁺)) * + (1 - sign(v) * v * dt) +end + +function compute_Δ𝛼_linvanleer(a⁻, a⁰, a⁺, v, dt, ::AlgebraicMean) + return ((a⁰ - a⁻) + (a⁺ - a⁰)) / 2 * (1 - sign(v) * v * dt) +end + +function slope_limited_product(v, a⁻, a⁻⁻, a⁺, a⁺⁺, dt, constraint) + # 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, constraint) + return v ⊠ (a⁻ ⊞ RecursiveApply.rdiv(Δ𝛼, 2)) + else + # Eqn (2,5a,5b,5c) + Δ𝛼 = compute_Δ𝛼_linvanleer(a⁻, a⁺, a⁺⁺, v, dt, constraint) + 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( + op::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, op.constraint), + ) +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( + op::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( + op::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 +1812,6 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(vᶠ))) end - - -######################### """ U = FCTZalesak(;boundaries) U.(A, Φ, Φᵗᵈ) @@ -1667,6 +1870,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 +1882,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 +1907,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 +1997,237 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) end +""" + U = TVDLimitedFluxC2F(;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 + + +""" + U = RZeroLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) + +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +struct RZeroLimiter <: AbstractTVDSlopeLimiter end + +""" + U = RHalfLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +struct RHalfLimiter <: AbstractTVDSlopeLimiter end + +""" + U = RMaxLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) + +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +struct RMaxLimiter <: AbstractTVDSlopeLimiter end + +""" + U = MinModLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) + +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +struct MinModLimiter <: AbstractTVDSlopeLimiter end + +""" + U = KorenLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) + +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +struct KorenLimiter <: AbstractTVDSlopeLimiter end + +""" + U = SuperbeeLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) + +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +struct SuperbeeLimiter <: AbstractTVDSlopeLimiter end + +""" + U = MonotonizedCentralLimiter(;boundaries) + U.(𝒜, Φ, 𝓊) + +A subtype of [`AbstractTVDSlopeLimiter`](@ref) limiter. See +[`AbstractTVDSlopeLimiter`](@ref) for the general formulation. +""" +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 TVDLimitedFluxC2F{BCS, M} <: AdvectionOperator + bcs::BCS + method::M +end + +TVDLimitedFluxC2F(; method, kwargs...) = + TVDLimitedFluxC2F((; kwargs...), method) + +return_eltype(::TVDLimitedFluxC2F, A, Φ, 𝓊) = + Geometry.Contravariant3Vector{eltype(eltype(A))} + +return_space( + ::TVDLimitedFluxC2F, + A_space::AllFaceFiniteDifferenceSpace, + Φ_space::AllCenterFiniteDifferenceSpace, + u_space::AllFaceFiniteDifferenceSpace, +) = A_space + +function tvd_limited_flux(Aⱼ₋₁₂, Aⱼ₊₁₂, ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, rⱼ₊₁₂, constraint) + stable_zero = zero(eltype(Aⱼ₊₁₂)) + stable_one = one(eltype(Aⱼ₊₁₂)) + Cⱼ₊₁₂ = compute_limiter_coeff(rⱼ₊₁₂, constraint) + @assert Cⱼ₊₁₂ <= 2 + @assert Cⱼ₊₁₂ >= 0 + return Cⱼ₊₁₂ * Aⱼ₊₁₂ +end + +stencil_interior_width(::TVDLimitedFluxC2F, A_space, Φ_space, u_space) = + ((-1, 1), (-half - 1, half + 1), (-1, +1)) + +Base.@propagate_inbounds function stencil_interior( + op::TVDLimitedFluxC2F, + 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ⱼ₊₁₂, op.method), + ) +end + +@inline function compute_slope_ratio(ϕⱼ, ϕⱼ₋₁, ϕⱼ₊₁, ϕⱼ₊₂, 𝓊) + if 𝓊 >= 0 + return (ϕⱼ - ϕⱼ₋₁) / (ϕⱼ₊₁ - ϕⱼ + eps(eltype(ϕⱼ))) + else + return (ϕⱼ₊₂ - ϕⱼ₊₁) / (ϕⱼ₊₁ - ϕⱼ + eps(eltype(ϕⱼ))) + end +end + +boundary_width(::TVDLimitedFluxC2F, ::AbstractBoundaryCondition) = 2 + +Base.@propagate_inbounds function stencil_left_boundary( + ::TVDLimitedFluxC2F, + 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( + ::TVDLimitedFluxC2F, + 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 +3881,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..8a647bb67d 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(), + constraint = Operators.MonotoneLocalExtrema(), + ) + + 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(), + constraint = Operators.MonotoneHarmonic(), + ) + + 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(), + constraint = Operators.PositiveDefinite(), + ) + + 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)