Merge #889

889: Examples: Add BB FCT to deformation flow r=valeriabarra a=valeriabarra Adding a WIP PR for now to gather feedback. This adds the BB FCT to the deformation flow example. - [x] Code follows the [style guidelines](https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/) OR N/A. - [x] Unit tests are included OR N/A. - [x] Code is exercised in an integration test OR N/A. - [x] Documentation has been added/updated OR N/A. Co-authored-by: Daniel-Huang <[email protected]> Co-authored-by: Valeria Barra <[email protected]>
CliMA · Aug 23, 2022 · b44786c · b44786c
2 parents f352ec1 + 4dc0aaa
commit b44786c
Show file tree

Hide file tree

Showing 2 changed files with 129 additions and 50 deletions.
diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl
@@ -43,12 +43,13 @@ const z_c = 5.0e3 # initial altitude of tracers
 const R_t = R / 2 # horizontal half-width of tracers
 const Z_t = 1000.0 # vertical half-width of tracers
 const κ₄ = 1.0e16 # hyperviscosity
-const C = 1.0 # flux-correction coefficient
 
 # time constants
 T = 86400.0 * 12.0
 dt = 60.0 * 60.0
 
+FT = Float64
+
 # set up function space
 function sphere_3D(
     R = 6.37122e6,
@@ -57,7 +58,6 @@ function sphere_3D(
     zelem = 36,
     npoly = 4,
 )
-    FT = Float64
     vertdomain = Domains.IntervalDomain(
         Geometry.ZPoint{FT}(zlim[1]),
         Geometry.ZPoint{FT}(zlim[2]);
@@ -121,7 +121,17 @@ y0 = map(coords) do coord
     ρq3 = ρ_ref(z) * q3
     ρq4 = ρ_ref(z) * q4
 
-    return (ρ = ρ_ref(z), ρq1 = ρq1, ρq2 = ρq2, ρq3 = ρq3, ρq4 = ρq4)
+    return (
+        ρ = ρ_ref(z),
+        ρq1 = ρq1,
+        ρq2 = ρq2,
+        ρq3 = ρq3,
+        ρq4 = ρq4,
+        ρq1_td = ρq1,
+        ρq2_td = ρq2,
+        ρq3_td = ρq3,
+        ρq4_td = ρq4,
+    )
 end
 
 y0 = Fields.FieldVector(
@@ -130,11 +140,15 @@ y0 = Fields.FieldVector(
     ρq2 = y0.ρq2,
     ρq3 = y0.ρq3,
     ρq4 = y0.ρq4,
+    ρq1_td = y0.ρq1,
+    ρq2_td = y0.ρq2,
+    ρq3_td = y0.ρq3,
+    ρq4_td = y0.ρq4,
 )
 
 function rhs!(dydt, y, parameters, t, alpha, beta)
 
-    (; τ, coords, face_coords, ystar) = parameters
+    (; τ, coords, face_coords, ystar, y_td) = parameters
 
     ϕ = coords.lat
     λ = coords.long
@@ -165,18 +179,23 @@ function rhs!(dydt, y, parameters, t, alpha, beta)
     uₕ = Geometry.Covariant12Vector.(Geometry.UVVector.(uu, uv))
     w = Geometry.Covariant3Vector.(Geometry.WVector.(uw))
 
+    # aliases
     ρ = y.ρ
     ρq1 = y.ρq1
     ρq2 = y.ρq2
     ρq3 = y.ρq3
     ρq4 = y.ρq4
 
-    dρ = dydt.ρ
     dρq1 = dydt.ρq1
     dρq2 = dydt.ρq2
     dρq3 = dydt.ρq3
     dρq4 = dydt.ρq4
 
+    ρq1_td = y_td.ρq1
+    ρq2_td = y_td.ρq2
+    ρq3_td = y_td.ρq3
+    ρq4_td = y_td.ρq4
+
     # No change in density: divergence-free flow
     @. dydt.ρ = beta * dydt.ρ + 0 .* y.ρ
 
@@ -197,6 +216,10 @@ function rhs!(dydt, y, parameters, t, alpha, beta)
         bottom = Operators.ThirdOrderOneSided(),
         top = Operators.ThirdOrderOneSided(),
     )
+    FCTBB = Operators.FCTBorisBook(
+        bottom = Operators.FirstOrderOneSided(),
+        top = Operators.FirstOrderOneSided(),
+    )
     hdiv = Operators.Divergence()
     hwdiv = Operators.WeakDivergence()
     hgrad = Operators.Gradient()
@@ -221,36 +244,93 @@ function rhs!(dydt, y, parameters, t, alpha, beta)
     cw = If2c.(w)
     cuvw = Geometry.Covariant123Vector.(uₕ) .+ Geometry.Covariant123Vector.(cw)
 
-    # The following, simple Flux-corrected Transport regression test (falling back to the 3rd-order upwinding, with C=1) is implemented following
-    # Ref: https://link.springer.com/book/10.1007/978-1-4419-6412-0 , Sec. 5.4 (p. 221)
-    corrected_antidiff_flux = similar(dρq1)
-    @. dρq1 = beta * dρq1 - alpha * hdiv(cuvw * ρq1) + alpha * ystar.ρq1
-    @. corrected_antidiff_flux = vdivf2c(
-        Ic2f(ρ) *
-        C *
-        (
-            third_order_upwind_c2f(w, ρq1 / ρ) -
-            first_order_upwind_c2f(w, ρq1 / ρ)
-        ),
-    )
-    @. dρq1 -=
-        alpha * (
-            vdivf2c(Ic2f(ρ) * first_order_upwind_c2f(w, ρq1 / ρ)) +
-            corrected_antidiff_flux
+    # 1) Vertical transport for ρq1:
+    # 1.1) vertical advection by vertical velocity, corrected by BB FCT:
+    @. ρq1_td =
+        beta * ρq1 -
+        alpha * vdivf2c(Ic2f(ρ) * first_order_upwind_c2f(w, ρq1 / ρ))
+    @. dρq1 =
+        beta * (dρq1 - ρq1) + ρq1_td -
+        alpha * vdivf2c(
+            Ic2f(ρ) * FCTBB(
+                (
+                    third_order_upwind_c2f(w, ρq1 / ρ) -
+                    first_order_upwind_c2f(w, ρq1 / ρ)
+                ),
+                ρq1_td / ρ / alpha,
+            ),
         )
+
+    # 1.2) vertical advection by horizontal velocity:
     @. dρq1 -= alpha * vdivf2c(Ic2f(uₕ * ρq1))
+    # 2) Horizontal transport for ρq1 (includes both horizontal and vertical velocity components)
+    @. dρq1 -= alpha * hdiv(cuvw * ρq1) - alpha * ystar.ρq1
+
+    # 1) Vertical transport for ρq2:
+    # 1.1) vertical advection by vertical velocity, corrected by BB FCT:
+    @. ρq2_td =
+        beta * ρq2 -
+        alpha * vdivf2c(Ic2f(ρ) * first_order_upwind_c2f(w, ρq2 / ρ))
+    @. dρq2 =
+        beta * (dρq2 - ρq2) + ρq2_td -
+        alpha * vdivf2c(
+            Ic2f(ρ) * FCTBB(
+                (
+                    third_order_upwind_c2f(w, ρq2 / ρ) -
+                    first_order_upwind_c2f(w, ρq2 / ρ)
+                ),
+                ρq2_td / ρ / alpha,
+            ),
+        )
 
-    @. dρq2 = beta * dρq2 - alpha * hdiv(cuvw * ρq2) + alpha * ystar.ρq2
-    @. dρq2 -= alpha * vdivf2c(Ic2f(ρ) * third_order_upwind_c2f.(w, ρq2 ./ ρ))
+    # 1.2) vertical advection by horizontal velocity:
     @. dρq2 -= alpha * vdivf2c(Ic2f(uₕ * ρq2))
+    # 2) Horizontal transport for ρq2 (includes both horizontal and vertical velocity components)
+    @. dρq2 -= alpha * hdiv(cuvw * ρq2) - alpha * ystar.ρq2
+
+    # 1) Vertical transport for ρq3:
+    # 1.1) vertical advection by vertical velocity, corrected by BB FCT:
+    @. ρq3_td =
+        beta * ρq3 -
+        alpha * vdivf2c(Ic2f(ρ) * first_order_upwind_c2f(w, ρq3 / ρ))
+    @. dρq3 =
+        beta * (dρq3 - ρq3) + ρq3_td -
+        alpha * vdivf2c(
+            Ic2f(ρ) * FCTBB(
+                (
+                    third_order_upwind_c2f(w, ρq3 / ρ) -
+                    first_order_upwind_c2f(w, ρq3 / ρ)
+                ),
+                ρq3_td / ρ / alpha,
+            ),
+        )
 
-    @. dρq3 = beta * dρq3 - alpha * hdiv(cuvw * ρq3) + alpha * ystar.ρq3
-    @. dρq3 -= alpha * vdivf2c(Ic2f(ρ) * third_order_upwind_c2f.(w, ρq3 ./ ρ))
+    # 1.2) vertical advection by horizontal velocity:
     @. dρq3 -= alpha * vdivf2c(Ic2f(uₕ * ρq3))
+    # 2) Horizontal transport for ρq3 (includes both horizontal and vertical velocity components)
+    @. dρq3 -= alpha * hdiv(cuvw * ρq3) - alpha * ystar.ρq3
+
+    # 1) Vertical transport for ρq4:
+    # 1.1) vertical advection by vertical velocity, corrected by BB FCT:
+    @. ρq4_td =
+        beta * ρq4 -
+        alpha * vdivf2c(Ic2f(ρ) * first_order_upwind_c2f(w, ρq4 / ρ))
+    @. dρq4 =
+        beta * (dρq4 - ρq4) + ρq4_td -
+        alpha * vdivf2c(
+            Ic2f(ρ) * FCTBB(
+                (
+                    third_order_upwind_c2f(w, ρq4 / ρ) -
+                    first_order_upwind_c2f(w, ρq4 / ρ)
+                ),
+                ρq4_td / ρ / alpha,
+            ),
+        )
 
-    @. dρq4 = beta * dρq4 - alpha * hdiv(cuvw * ρq4) + alpha * ystar.ρq4
-    @. dρq4 -= alpha * vdivf2c(Ic2f(ρ) * third_order_upwind_c2f.(w, ρq4 ./ ρ))
+    # 1.2) vertical advection by horizontal velocity:
     @. dρq4 -= alpha * vdivf2c(Ic2f(uₕ * ρq4))
+    # 2) Horizontal transport for ρq4 (includes both horizontal and vertical velocity components)
+    @. dρq4 -= alpha * hdiv(cuvw * ρq4) - alpha * ystar.ρq4
 
     Spaces.weighted_dss!(dydt.ρ)
     Spaces.weighted_dss!(dydt.ρq1)
@@ -263,7 +343,9 @@ end
 
 # Set up RHS function
 ystar = copy(y0)
-parameters = (; τ, coords, face_coords, ystar)
+y_td = copy(y0)
+
+parameters = (; τ, coords, face_coords, ystar, y_td)
 
 rhs!(ystar, y0, parameters, 0.0, dt, 1)
 
@@ -287,7 +369,7 @@ q1_error =
 q2_error =
     norm(sol.u[end].ρq2 ./ ρ_ref.(coords.z) .- y0.ρq2 ./ ρ_ref.(coords.z)) /
     norm(y0.ρq2 ./ ρ_ref.(coords.z))
-@test q2_error ≈ 0.0 atol = 0.03
+@test q2_error ≈ 0.0 atol = 0.032
 
 q3_error =
     norm(sol.u[end].ρq3 ./ ρ_ref.(coords.z) .- y0.ρq3 ./ ρ_ref.(coords.z)) /
@@ -297,28 +379,28 @@ q3_error =
 q4_error =
     norm(sol.u[end].ρq4 ./ ρ_ref.(coords.z) .- y0.ρq4 ./ ρ_ref.(coords.z)) /
     norm(y0.ρq4 ./ ρ_ref.(coords.z))
-@test q4_error ≈ 0.0 atol = 0.03
+@test q4_error ≈ 0.0 atol = 0.025
 
 # Tracer mass conservation checks
-q1_initial_mass = sum(y0.ρq1 ./ ρ_ref.(coords.z))
-q1_mass = sum(sol.u[end].ρq1 ./ ρ_ref.(coords.z))
+q1_initial_mass = sum(y0.ρq1)
+q1_mass = sum(sol.u[end].ρq1)
 q1_rel_mass_err = norm((q1_mass - q1_initial_mass) / q1_initial_mass)
-@test q1_rel_mass_err ≈ 0.0 atol = 1.1e-5
+@test q1_rel_mass_err ≈ 0.0 atol = 8e1eps(FT)
 
-q2_initial_mass = sum(y0.ρq2 ./ ρ_ref.(coords.z))
-q2_mass = sum(sol.u[end].ρq2 ./ ρ_ref.(coords.z))
+q2_initial_mass = sum(y0.ρq2)
+q2_mass = sum(sol.u[end].ρq2)
 q2_rel_mass_err = norm((q2_mass - q2_initial_mass) / q2_initial_mass)
-@test q2_rel_mass_err ≈ 0.0 atol = 2.2e-6
+@test q2_rel_mass_err ≈ 0.0 atol = 8e1eps(FT)
 
-q3_initial_mass = sum(y0.ρq3 ./ ρ_ref.(coords.z))
-q3_mass = sum(sol.u[end].ρq3 ./ ρ_ref.(coords.z))
+q3_initial_mass = sum(y0.ρq3)
+q3_mass = sum(sol.u[end].ρq3)
 q3_rel_mass_err = norm((q3_mass - q3_initial_mass) / q3_initial_mass)
-@test q3_rel_mass_err ≈ 0.0 atol = 2.5e-6
+@test q3_rel_mass_err ≈ 0.0 atol = 2e2eps(FT)
 
-q4_initial_mass = sum(y0.ρq4 ./ ρ_ref.(coords.z))
-q4_mass = sum(sol.u[end].ρq4 ./ ρ_ref.(coords.z))
+q4_initial_mass = sum(y0.ρq4)
+q4_mass = sum(sol.u[end].ρq4)
 q4_rel_mass_err = norm((q4_mass - q4_initial_mass) / q4_initial_mass)
-@test q4_rel_mass_err ≈ 0.0 atol = 2.2e-6
+@test q4_rel_mass_err ≈ 0.0 atol = 8e1eps(FT)
 
 # visualization artifacts
 ENV["GKSwstype"] = "nul"

diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl
@@ -1502,7 +1502,7 @@ return_space(
     arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
 ) = velocity_space
 
-@inline function fct_boris_book(v, a⁻⁻, a⁻, a⁺, a⁺⁺, step)
+@inline function fct_boris_book(v, a⁻⁻, a⁻, a⁺, a⁺⁺)
     if v != zero(eltype(v))
         sign(v) ⊠ (RecursiveApply.rmap(
             max,
@@ -1512,8 +1512,8 @@ return_space(
                 RecursiveApply.rmap(abs, v),
                 RecursiveApply.rmap(
                     min,
-                    sign(v) ⊠ (a⁺⁺ - a⁺) ⊠ step,
-                    sign(v) ⊠ (a⁻ - a⁻⁻) ⊠ step,
+                    sign(v) ⊠ (a⁺⁺ - a⁺),
+                    sign(v) ⊠ (a⁻ - a⁻⁻),
                 ),
             ),
         ))
@@ -1524,7 +1524,7 @@ return_space(
             RecursiveApply.rmap(
                 min,
                 v,
-                RecursiveApply.rmap(min, (a⁺⁺ - a⁺) ⊠ step, (a⁻ - a⁻⁻) ⊠ step),
+                RecursiveApply.rmap(min, (a⁺⁺ - a⁺), (a⁻ - a⁻⁻)),
             ),
         )
     end
@@ -1543,10 +1543,7 @@ stencil_interior_width(::FCTBorisBook, velocity, arg) =
         getidx(velocity, loc, idx, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    step = Geometry.LocalGeometry(space, idx, hidx).∂x∂ξ[1]
-    return Geometry.Contravariant3Vector(
-        fct_boris_book(vᶠ, a⁻⁻, a⁻, a⁺, a⁺⁺, step),
-    )
+    return Geometry.Contravariant3Vector(fct_boris_book(vᶠ, a⁻⁻, a⁻, a⁺, a⁺⁺))
 end
 
 boundary_width(::FCTBorisBook, ::FirstOrderOneSided, velocity, arg) = 2