modified: examples/column/tvd_fct_advection.jl

modified: examples/hybrid/sphere/deformation_flow.jl modified: src/Operators/finitedifference.jl
CliMA · Nov 13, 2024 · e993c97 · e993c97
1 parent f268404
commit e993c97
Show file tree

Hide file tree

Showing 3 changed files with 43 additions and 24 deletions.
diff --git 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
@@ -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
@@ -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