From e993c973a13978c38baa5ba74b15e876dde9e6ad Mon Sep 17 00:00:00 2001 From: akshaysridhar Date: Wed, 13 Nov 2024 13:10:20 -0800 Subject: [PATCH] modified: examples/column/tvd_fct_advection.jl modified: examples/hybrid/sphere/deformation_flow.jl modified: src/Operators/finitedifference.jl --- examples/column/tvd_fct_advection.jl | 41 ++++++++++++---------- examples/hybrid/sphere/deformation_flow.jl | 11 ++++-- src/Operators/finitedifference.jl | 15 ++++++-- 3 files changed, 43 insertions(+), 24 deletions(-) diff --git a/examples/column/tvd_fct_advection.jl b/examples/column/tvd_fct_advection.jl index 9e8a10ae35..d380f32ebd 100644 --- a/examples/column/tvd_fct_advection.jl +++ b/examples/column/tvd_fct_advection.jl @@ -62,13 +62,11 @@ function tendency!(yₜ, y, parameters, t) If = Operators.InterpolateC2F() @. yₜ.q = -divf2c( - upwind1(w, y.q) + LimitedFlux( - upwind3(w, y.q) - upwind1(w, y.q), - y.q / Δt, - ), + upwind1(w, y.q) + + LimitedFlux(upwind3(w, y.q) - upwind1(w, y.q), y.q / Δt), ) - # @. yₜ.q = - # -divf2c(w * If(y.q)) + # @. yₜ.q = + # -divf2c(w * If(y.q)) end # Define a pulse wave or square wave @@ -95,15 +93,16 @@ stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(FT(7.0))) plot_string = ["uniform", "stretched"] for (i, stretch_fn) in enumerate(stretch_fns) - limiter_methods = (Operators.RZeroLimiter(), - Operators.RHalfLimiter(), - Operators.RMaxLimiter(), - Operators.MinModLimiter(), - Operators.KorenLimiter(), - Operators.SuperbeeLimiter(), - Operators.MonotonizedCentralLimiter(), - Operators.VanLeerLimiter(), - ) + limiter_methods = ( + Operators.RZeroLimiter(), + Operators.RHalfLimiter(), + Operators.RMaxLimiter(), + Operators.MinModLimiter(), + Operators.KorenLimiter(), + Operators.SuperbeeLimiter(), + Operators.MonotonizedCentralLimiter(), + Operators.VanLeerLimiter(), + ) for limiter_method in limiter_methods @info (limiter_method, stretch_fn) mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n) @@ -127,7 +126,12 @@ for (i, stretch_fn) in enumerate(stretch_fns) (t₀, t₁), parameters, ) - sol = solve(prob, ExplicitAlgorithm(SSP33ShuOsher()), dt = Δt, saveat = Δt) + sol = solve( + prob, + ExplicitAlgorithm(SSP33ShuOsher()), + dt = Δt, + saveat = Δt, + ) q_final = sol.u[end].q q_analytic = pulse.(z, t₁, z₀, zₕ, z₁) @@ -137,10 +141,11 @@ for (i, stretch_fn) in enumerate(stretch_fns) plot(q_final) Plots.png( - Plots.plot!(q_analytic, title = "$(typeof(limiter_method))"), + Plots.plot!(q_analytic, title = "$(typeof(limiter_method))"), joinpath( path, - "exact_and_computed_advected_square_wave_TVDSlopeLimitedFlux_"* "$(typeof(limiter_method))_" * + "exact_and_computed_advected_square_wave_TVDSlopeLimitedFlux_" * + "$(typeof(limiter_method))_" * plot_string[i] * ".png", ), diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl index 830766e78f..b9e3bef55d 100644 --- a/examples/hybrid/sphere/deformation_flow.jl +++ b/examples/hybrid/sphere/deformation_flow.jl @@ -188,8 +188,8 @@ function vertical_tendency!(Yₜ, Y, cache, t) @. ρqₜ_n -= vdivf2c( ᶠwinterp(ᶜJ, Y.c.ρ) * ( upwind1(face_uᵥ, q_n) + SlopeLimitedFlux( - upwind3(face_uᵥ, q_n) - upwind1(face_uᵥ, q_n), - Y.c.ρ * q_n, + upwind3(face_uᵥ, q_n) - upwind1(face_uᵥ, q_n), + Y.c.ρ * q_n, ) ), ) @@ -431,7 +431,12 @@ for (sol, suffix) in ( 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.ρ)), + ( + ((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), ), diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 1e6da794fd..a47f0fa1ba 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -1852,7 +1852,7 @@ end end @inline function compute_limiter_coeff(r, ::KorenLimiter) - return max(zero(eltype(r)), min(2r, min(1/3 + 2r / 3, 2))) + return max(zero(eltype(r)), min(2r, min(1 / 3 + 2r / 3, 2))) end @inline function compute_limiter_coeff(r, ::SuperbeeLimiter) @@ -1888,7 +1888,7 @@ function tvd_limited_flux(Aⱼ₋₁₂, Aⱼ₊₁₂, ϕⱼ₋₁, ϕⱼ, ϕ stable_one = one(eltype(Aⱼ₊₁₂)) # Test with various limiter methods # Aⱼ₊₁₂ is the antidiffusive flux (see Durran textbook for notation) - Cⱼ₊₁₂ = compute_limiter_coeff(rⱼ₊₁₂, method) + Cⱼ₊₁₂ = compute_limiter_coeff(rⱼ₊₁₂, method) # Assert condition for TVD limiter @assert Cⱼ₊₁₂ <= eltype(Aⱼ₊₁₂)(2) @assert Cⱼ₊₁₂ >= eltype(Aⱼ₊₁₂)(0) @@ -1925,7 +1925,16 @@ Base.@propagate_inbounds function stencil_interior( rⱼ₊₁₂ = compute_slope_ratio(ϕⱼ, ϕⱼ₋₁, ϕⱼ₊₁) return Geometry.Contravariant3Vector( - tvd_limited_flux(Aⱼ₋₁₂, Aⱼ₊₁₂, ϕⱼ₋₁, ϕⱼ, ϕⱼ₊₁, ϕⱼ₊₂, rⱼ₊₁₂, ℱ.bcs.method), + tvd_limited_flux( + Aⱼ₋₁₂, + Aⱼ₊₁₂, + ϕⱼ₋₁, + ϕⱼ, + ϕⱼ₊₁, + ϕⱼ₊₂, + rⱼ₊₁₂, + ℱ.bcs.method, + ), ) end