Merge pull request #3702 from CliMA/zs/scale_diffusion

add a scaling factor for precipitation diffusion

Author: szy21

.buildkite/Manifest-v1.11.toml

@@ -394,9 +394,9 @@ version = "0.2.12"
 
 [[deps.ClimaParams]]
 deps = ["TOML"]
-git-tree-sha1 = "5f624b24e522ce61bbe0b63ce6fb1c7381b9fa91"
+git-tree-sha1 = "bf8fdee01cf5ea7e5e9fa00ac32367f65448484e"
 uuid = "5c42b081-d73a-476f-9059-fd94b934656c"
-version = "0.10.21"
+version = "0.10.23"
 
 [[deps.ClimaReproducibilityTests]]
 deps = ["OrderedCollections", "PrettyTables"]

Project.toml

@@ -45,7 +45,7 @@ AtmosphericProfilesLibrary = "0.1.7"
 ClimaComms = "0.6.6"
 ClimaCore = "0.14.24"
 ClimaDiagnostics = "0.2.12"
-ClimaParams = "0.10.19"
+ClimaParams = "0.10.23"
 ClimaTimeSteppers = "0.8.2"
 ClimaUtilities = "0.1.22"
 CloudMicrophysics = "0.22.9"

src/cache/temporary_quantities.jl

@@ -65,6 +65,8 @@ function temporary_quantities(Y, atmos)
             BidiagonalMatrixRow{Adjoint{FT, C3{FT}}},
         ),
         ᶜdiffusion_h_matrix = similar(Y.c, TridiagonalMatrixRow{FT}),
+        ᶜdiffusion_h_matrix_scaled = similar(Y.c, TridiagonalMatrixRow{FT}),
+        ᶜtridiagonal_matrix_scalar = similar(Y.c, TridiagonalMatrixRow{FT}),
         ᶜdiffusion_u_matrix = similar(Y.c, TridiagonalMatrixRow{FT}),
         ᶠtridiagonal_matrix_c3 = similar(Y.f, TridiagonalMatrixRow{C3{FT}}),
     )

src/parameters/Parameters.jl

@@ -106,6 +106,10 @@ Base.@kwdef struct ClimaAtmosParameters{
     water_refractive_index::FT
     optics_lookup_temperature_min::FT
     optics_lookup_temperature_max::FT
+    # Hyperdiffusion
+    α_hyperdiff_tracer::FT
+    # Vertical diffusion
+    α_vert_diff_tracer::FT
 end
 
 Base.eltype(::ClimaAtmosParameters{FT}) where {FT} = FT

src/parameters/create_parameters.jl

@@ -105,6 +105,8 @@ atmos_name_map = (;
     :water_refractive_index => :water_refractive_index,
     :optics_lookup_temperature_min => :optics_lookup_temperature_min,
     :optics_lookup_temperature_max => :optics_lookup_temperature_max,
+    :tracer_hyperdiffusion_factor => :α_hyperdiff_tracer,
+    :tracer_vertical_diffusion_factor => :α_vert_diff_tracer,
 )
 
 cloud_parameters(::Type{FT}) where {FT <: AbstractFloat} =

src/prognostic_equations/hyperdiffusion.jl

@@ -230,6 +230,7 @@ NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     (; hyperdiff, turbconv_model) = p.atmos
     isnothing(hyperdiff) && return nothing
 
+    α_hyperdiff_tracer = CAP.α_hyperdiff_tracer(p.params)
     (; ν₄_scalar_coeff) = hyperdiff
     h_space = Spaces.horizontal_space(axes(Y.c))
     h_length_scale = Spaces.node_horizontal_length_scale(h_space) # mean nodal distance
@@ -244,11 +245,13 @@ NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     # TODO: Figure out why caching the duplicated tendencies in ᶜtemp_scalar
     # triggers allocations.
     for (ᶜρχₜ, ᶜ∇²χ, χ_name) in matching_subfields(Yₜ.c, ᶜ∇²specific_tracers)
-        ν₄_scalar = ifelse(χ_name in (:q_rai, :q_sno), 0 * ν₄_scalar, ν₄_scalar)
+        ν₄_scalar = ifelse(
+            χ_name in (:q_rai, :q_sno),
+            α_hyperdiff_tracer * ν₄_scalar,
+            ν₄_scalar,
+        )
         @. ᶜρχₜ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
-        if !(χ_name in (:q_rai, :q_sno))
-            @. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
-        end
+        @. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
     end
     if turbconv_model isa PrognosticEDMFX
         (; ᶜ∇²q_totʲs) = p.hyperdiff

src/prognostic_equations/implicit/implicit_solver.jl

@@ -503,6 +503,8 @@ NVTX.@annotate function Wfact!(A, Y, p, dtγ, t)
         p.scratch.ᶠp_grad_matrix,
         p.scratch.ᶜadvection_matrix,
         p.scratch.ᶜdiffusion_h_matrix,
+        p.scratch.ᶜdiffusion_h_matrix_scaled,
+        p.scratch.ᶜtridiagonal_matrix_scalar,
         p.scratch.ᶜdiffusion_u_matrix,
         p.scratch.ᶠbidiagonal_matrix_ct3,
         p.scratch.ᶠbidiagonal_matrix_ct3_2,
@@ -540,7 +542,13 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
         ᶠbidiagonal_matrix_ct3_2,
         ᶠtridiagonal_matrix_c3,
     ) = p
-    (; ᶜdiffusion_h_matrix, ᶜdiffusion_u_matrix, params) = p
+    (;
+        ᶜdiffusion_h_matrix,
+        ᶜdiffusion_h_matrix_scaled,
+        ᶜtridiagonal_matrix_scalar,
+        ᶜdiffusion_u_matrix,
+        params,
+    ) = p
     (; edmfx_upwinding) = p.atmos.numerics
 
     FT = Spaces.undertype(axes(Y.c))
@@ -716,10 +724,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     end
 
     if use_derivative(diffusion_flag)
+        α_vert_diff_tracer = CAP.α_vert_diff_tracer(params)
         (; ᶜK_h, ᶜK_u) = p
         @. ᶜdiffusion_h_matrix =
             ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * ᶠinterp(ᶜK_h)) ⋅
             ᶠgradᵥ_matrix()
+        @. ᶜdiffusion_h_matrix_scaled =
+            ᶜadvdivᵥ_matrix() ⋅ DiagonalMatrixRow(
+                ᶠinterp(ᶜρ) * ᶠinterp(α_vert_diff_tracer * ᶜK_h),
+            ) ⋅ ᶠgradᵥ_matrix()
         if (
             MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke)) ||
             !isnothing(p.atmos.turbconv_model) ||
@@ -759,10 +772,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
             ᶜq = MatrixFields.get_field(ᶜspecific, q_name)
             ∂ᶜρq_err_∂ᶜρ = matrix[ρq_name, @name(c.ρ)]
             ∂ᶜρq_err_∂ᶜρq = matrix[ρq_name, ρq_name]
+            ᶜtridiagonal_matrix_scalar = ifelse(
+                q_name in (:q_rai, :q_sno),
+                ᶜdiffusion_h_matrix_scaled,
+                ᶜdiffusion_h_matrix,
+            )
             @. ∂ᶜρq_err_∂ᶜρ =
-                dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(-(ᶜq) / ᶜρ)
+                dtγ * ᶜtridiagonal_matrix_scalar ⋅ DiagonalMatrixRow(-(ᶜq) / ᶜρ)
             @. ∂ᶜρq_err_∂ᶜρq +=
-                dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(1 / ᶜρ)
+                dtγ * ᶜtridiagonal_matrix_scalar ⋅ DiagonalMatrixRow(1 / ᶜρ)
         end
 
         if MatrixFields.has_field(Y, @name(c.sgs⁰.ρatke))

src/prognostic_equations/vertical_diffusion_boundary_layer.jl

@@ -18,7 +18,8 @@ function vertical_diffusion_boundary_layer_tendency!(
     ::Union{VerticalDiffusion, DecayWithHeightDiffusion},
 )
     FT = eltype(Y)
-    (; ᶜu, ᶜh_tot, ᶜspecific, ᶜK_u, ᶜK_h, sfc_conditions) = p.precomputed
+    α_vert_diff_tracer = CAP.α_vert_diff_tracer(p.params)
+    (; ᶜu, ᶜh_tot, ᶜspecific, ᶜK_u, ᶜK_h) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
 
     if !disable_momentum_vertical_diffusion(p.atmos.vert_diff)
@@ -38,8 +39,14 @@ function vertical_diffusion_boundary_layer_tendency!(
         ᶜdivᵥ_ρe_tot(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜh_tot)))
 
     ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
+    ᶜK_h_scaled = p.scratch.ᶜtemp_scalar_2
     for (ᶜρχₜ, ᶜχ, χ_name) in matching_subfields(Yₜ.c, ᶜspecific)
         χ_name == :e_tot && continue
+        if χ_name in (:q_rai, :q_sno)
+            @. ᶜK_h_scaled = α_vert_diff_tracer * ᶜK_h
+        else
+            @. ᶜK_h_scaled = ᶜK_h
+        end
         ᶜdivᵥ_ρχ = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
             bottom = Operators.SetValue(C3(FT(0))),