Add new slope-limited transport operator

NEWS.md

@@ -9,6 +9,12 @@ v0.28.1
 -------
 ### Features
 
+### Add van Leer class operator
+
+Added a new vertical transport option `vanleer_limiter` (for tracer and energy variables)
+which uses methods described in Lin et al. (1994) to apply slope-limited upwinding. Adds 
+operator 
+
 ### Read initial conditions from NetCDF files
 
 Added functionality to allow initial conditions to be overwritten by

config/default_configs/default_config.yml

@@ -206,11 +206,11 @@ viscous_sponge:
   help: "Viscous sponge [`true`, `false` (default)]"
   value: false
 tracer_upwinding:
-  help: "Tracer upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`]"
-  value: none
+  help: "Tracer upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`, `vanleer_limiter`]"
+  value: vanleer_limiter
 energy_upwinding:
-  help: "Energy upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`]"
-  value: none
+  help: "Energy upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`, `vanleer_limiter`]"
+  value: vanleer_limiter
 orographic_gravity_wave:
   help: "Orographic drag on horizontal mean flow [`nothing` (default), `gfdl_restart`, `raw_topo`]"
   value: ~

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-194
+195
 
 # **README**
 #
@@ -21,6 +21,9 @@
 
 #=
 
+195
+- Use `vanleer_limiter` as default.
+
 194
 - Reproducibility infrastructure fixes.
 

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -2,6 +2,7 @@
 ##### Implicit tendencies
 #####
 
+import ClimaCore
 import ClimaCore: Fields, Geometry
 
 NVTX.@annotate function implicit_tendency!(Yₜ, Y, p, t)
@@ -73,6 +74,20 @@ vertical_transport!(
     ::Val{:first_order},
     ᶜdivᵥ,
 ) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)))
+@static if pkgversion(ClimaCore) ≥ v"0.14.22"
+    vertical_transport!(
+        coeff,
+        ᶜρχₜ,
+        ᶜJ,
+        ᶜρ,
+        ᶠu³,
+        ᶜχ,
+        dt,
+        ::Val{:vanleer_limiter},
+        ᶜdivᵥ,
+    ) = @. ᶜρχₜ +=
+        -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠlin_vanleer(ᶠu³, ᶜχ, dt)))
+end
 vertical_transport!(
     coeff,
     ᶜρχₜ,

src/solver/type_getters.jl

@@ -1,6 +1,7 @@
 using Dates: DateTime, @dateformat_str
 import Interpolations
 import NCDatasets
+import ClimaCore
 import ClimaUtilities.OutputPathGenerator
 import ClimaCore: InputOutput, Meshes, Spaces, Quadratures
 import ClimaAtmos.RRTMGPInterface as RRTMGPI
@@ -103,6 +104,19 @@ function get_numerics(parsed_args)
 
     energy_upwinding = Val(Symbol(parsed_args["energy_upwinding"]))
     tracer_upwinding = Val(Symbol(parsed_args["tracer_upwinding"]))
+
+    # Compat
+    if !(pkgversion(ClimaCore) ≥ v"0.14.22") &&
+       energy_upwinding == Val(:vanleer_limiter)
+        energy_upwinding = Val(:none)
+        @warn "energy_upwinding=vanleer_limiter is not supported for ClimaCore $(pkgversion(ClimaCore)), please upgrade. Setting energy_upwinding to :none"
+    end
+    if !(pkgversion(ClimaCore) ≥ v"0.14.22") &&
+       tracer_upwinding == Val(:vanleer_limiter)
+        tracer_upwinding = Val(:none)
+        @warn "tracer_upwinding=vanleer_limiter is not supported for ClimaCore $(pkgversion(ClimaCore)), please upgrade. Setting tracer_upwinding to :none"
+    end
+
     edmfx_upwinding = Val(Symbol(parsed_args["edmfx_upwinding"]))
     edmfx_sgsflux_upwinding =
         Val(Symbol(parsed_args["edmfx_sgsflux_upwinding"]))

src/utils/abbreviations.jl

@@ -1,4 +1,5 @@
 using ClimaCore: Geometry, Operators, MatrixFields
+import ClimaCore
 
 # Alternatively, we could use Vec₁₂₃, Vec³, etc., if that is more readable.
 const C1 = Geometry.Covariant1Vector
@@ -85,6 +86,13 @@ const ᶠfct_zalesak = Operators.FCTZalesak(
     bottom = Operators.FirstOrderOneSided(),
     top = Operators.FirstOrderOneSided(),
 )
+@static if pkgversion(ClimaCore) ≥ v"0.14.22"
+    const ᶠlin_vanleer = Operators.LinVanLeerC2F(
+        bottom = Operators.FirstOrderOneSided(),
+        top = Operators.FirstOrderOneSided(),
+        constraint = Operators.MonotoneLocalExtrema(), # (Mono5)
+    )
+end
 
 const ᶜinterp_matrix = MatrixFields.operator_matrix(ᶜinterp)
 const ᶜleft_bias_matrix = MatrixFields.operator_matrix(ᶜleft_bias)