Merge pull request #3853 from CliMA/ck/slim_cache2

Remove more dependency on specific quantities

Author: charleskawczynski

perf/jet_test_nfailures.jl

@@ -38,7 +38,7 @@ using Test
     # inference. By increasing this counter, we acknowledge that
     # we have introduced an inference failure. We hope to drive
     # this number down to 0.
-    n_allowed_failures = 249
+    n_allowed_failures = 253
     @show n
     if n < n_allowed_failures
         @info "Please update the n-failures to $n"

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -961,7 +961,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures!
     (; params) = p
     (; dt) = p
     (; ᶜp, ᶜu, ᶜts) = p.precomputed
-    (; q_tot) = p.precomputed.ᶜspecific
     (; ustar, obukhov_length) = p.precomputed.sfc_conditions
     (; ᶜtke⁰) = p.precomputed
     (;
@@ -1079,14 +1078,13 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
     microphys_0m_params = CAP.microphysics_0m_params(p.params)
     (; dt) = p
     (; ᶜts, ᶜSqₜᵖ⁰) = p.precomputed
-    (; q_tot) = p.precomputed.ᶜspecific
 
     # Environment precipitation sources (to be applied to grid mean)
     @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources(
         thermo_params,
         microphys_0m_params,
         dt,
-        q_tot,
+        specific(Y.c.ρq_tot, Y.c.ρ),
         ᶜts,
     )
     return nothing
@@ -1102,7 +1100,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
     #microphys_1m_params = CAP.microphysics_1m_params(p.params)
 
     #(; ᶜts, ᶜSqₜᵖ⁰, ᶜSeₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed
-    #(; q_tot) = p.precomputed.ᶜspecific
     #(; ᶜqᵣ, ᶜqₛ) = p.precomputed
 
     #ᶜSᵖ = p.scratch.ᶜtemp_scalar

src/cache/precipitation_precomputed_quantities.jl

@@ -453,39 +453,29 @@ function set_precipitation_surface_fluxes!(
 )
     (; surface_rain_flux, surface_snow_flux) = p.precomputed
     (; col_integrated_precip_energy_tendency,) = p.conservation_check
-    (; ᶜwᵣ, ᶜwₛ, ᶜwₗ, ᶜwᵢ, ᶜspecific) = p.precomputed
+    (; ᶜwᵣ, ᶜwₛ, ᶜwₗ, ᶜwᵢ) = p.precomputed
     ᶜJ = Fields.local_geometry_field(Y.c).J
     ᶠJ = Fields.local_geometry_field(Y.f).J
     sfc_J = Fields.level(ᶠJ, Fields.half)
     sfc_space = axes(sfc_J)
 
     # Jacobian-weighted extrapolation from interior to surface, consistent with
     # the reconstruction of density on cell faces, ᶠρ = ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ
-    int_J = Fields.Field(Fields.field_values(Fields.level(ᶜJ, 1)), sfc_space)
-    int_ρ = Fields.Field(Fields.field_values(Fields.level(Y.c.ρ, 1)), sfc_space)
+    sfc_lev(x) =
+        Fields.Field(Fields.field_values(Fields.level(x, 1)), sfc_space)
+    int_J = sfc_lev(ᶜJ)
+    int_ρ = sfc_lev(Y.c.ρ)
     sfc_ρ = @. lazy(int_ρ * int_J / sfc_J)
 
     # Constant extrapolation to surface, consistent with simple downwinding
-    sfc_qᵣ = Fields.Field(
-        Fields.field_values(Fields.level(ᶜspecific.q_rai, 1)),
-        sfc_space,
-    )
-    sfc_qₛ = Fields.Field(
-        Fields.field_values(Fields.level(ᶜspecific.q_sno, 1)),
-        sfc_space,
-    )
-    sfc_qₗ = Fields.Field(
-        Fields.field_values(Fields.level(ᶜspecific.q_liq, 1)),
-        sfc_space,
-    )
-    sfc_qᵢ = Fields.Field(
-        Fields.field_values(Fields.level(ᶜspecific.q_ice, 1)),
-        sfc_space,
-    )
-    sfc_wᵣ = Fields.Field(Fields.field_values(Fields.level(ᶜwᵣ, 1)), sfc_space)
-    sfc_wₛ = Fields.Field(Fields.field_values(Fields.level(ᶜwₛ, 1)), sfc_space)
-    sfc_wₗ = Fields.Field(Fields.field_values(Fields.level(ᶜwₗ, 1)), sfc_space)
-    sfc_wᵢ = Fields.Field(Fields.field_values(Fields.level(ᶜwᵢ, 1)), sfc_space)
+    sfc_wᵣ = sfc_lev(ᶜwᵣ)
+    sfc_wₛ = sfc_lev(ᶜwₛ)
+    sfc_wₗ = sfc_lev(ᶜwₗ)
+    sfc_wᵢ = sfc_lev(ᶜwᵢ)
+    sfc_qᵣ = lazy.(specific.(sfc_lev(Y.c.ρq_rai), sfc_ρ))
+    sfc_qₛ = lazy.(specific.(sfc_lev(Y.c.ρq_sno), sfc_ρ))
+    sfc_qₗ = lazy.(specific.(sfc_lev(Y.c.ρq_liq), sfc_ρ))
+    sfc_qᵢ = lazy.(specific.(sfc_lev(Y.c.ρq_ice), sfc_ρ))
 
     @. surface_rain_flux = sfc_ρ * (sfc_qᵣ * (-sfc_wᵣ) + sfc_qₗ * (-sfc_wₗ))
     @. surface_snow_flux = sfc_ρ * (sfc_qₛ * (-sfc_wₛ) + sfc_qᵢ * (-sfc_wᵢ))

src/parameterized_tendencies/radiation/radiation.jl

@@ -425,12 +425,16 @@ end
 
 function radiation_model_cache(Y, radiation_mode::RadiationDYCOMS)
     FT = Spaces.undertype(axes(Y.c))
-    NT = NamedTuple{(:z, :ρ, :q_tot), NTuple{3, FT}}
+    # The NT type is needed for the `column_reduce!` call below because
+    # `column_reduce!` computes the output, `isoline_z_ρ_ρq`, in-place. `Y.c.ρ`
+    # and `Y.c.ρq_tot` are input arguments, and `nt1` / `nt2` are output entry
+    # values at each point along the vertical column.
+    NT = NamedTuple{(:z, :ρ, :ρq_tot), NTuple{3, FT}}
     return (;
         ᶜκρq = similar(Y.c, FT),
         ∫_0_∞_κρq = similar(Spaces.level(Y.c, 1), FT),
         ᶠ∫_0_z_κρq = similar(Y.f, FT),
-        isoline_z_ρ_q = similar(Spaces.level(Y.c, 1), NT),
+        isoline_z_ρ_ρq = similar(Spaces.level(Y.c, 1), NT),
         ᶠradiation_flux = similar(Y.f, Geometry.WVector{FT}),
         net_energy_flux_toa = [Geometry.WVector(FT(0))],
         net_energy_flux_sfc = [Geometry.WVector(FT(0))],
@@ -440,13 +444,13 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS)
     @assert !(p.atmos.moisture_model isa DryModel)
 
     (; params) = p
-    (; ᶜspecific, ᶜts) = p.precomputed
-    (; ᶜκρq, ∫_0_∞_κρq, ᶠ∫_0_z_κρq, isoline_z_ρ_q, ᶠradiation_flux) =
+    (; ᶜκρq, ∫_0_∞_κρq, ᶠ∫_0_z_κρq, isoline_z_ρ_ρq, ᶠradiation_flux) =
         p.radiation
+    (; ᶜts) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
     cp_d = CAP.cp_d(params)
     FT = Spaces.undertype(axes(Y.c))
-    NT = NamedTuple{(:z, :ρ, :q_tot), NTuple{3, FT}}
+    NT = NamedTuple{(:z, :ρ, :ρq_tot), NTuple{3, FT}}
     ᶜz = Fields.coordinate_field(Y.c).z
     ᶠz = Fields.coordinate_field(Y.f).z
 
@@ -467,14 +471,14 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS)
     q_tot_isoline = FT(0.008)
     Operators.column_reduce!(
         (nt1, nt2) ->
-            abs(nt1.q_tot - q_tot_isoline) < abs(nt2.q_tot - q_tot_isoline) ?
-            nt1 : nt2,
-        isoline_z_ρ_q,
-        Base.broadcasted(NT ∘ tuple, ᶜz, Y.c.ρ, ᶜspecific.q_tot),
+            abs(specific(nt1.ρq_tot, nt1.ρ) - q_tot_isoline) <
+            abs(specific(nt2.ρq_tot, nt2.ρ) - q_tot_isoline) ? nt1 : nt2,
+        isoline_z_ρ_ρq,
+        Base.broadcasted(NT ∘ tuple, ᶜz, Y.c.ρ, Y.c.ρq_tot),
     )
 
-    zi = isoline_z_ρ_q.z
-    ρi = isoline_z_ρ_q.ρ
+    zi = isoline_z_ρ_ρq.z
+    ρi = isoline_z_ρ_ρq.ρ
 
     # TODO: According to the paper, we should remove the ifelse condition that
     # clips the third term to 0 below zi, and we should also replace cp_d with

src/prognostic_equations/advection.jl

@@ -228,7 +228,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     end
 
     if !(p.atmos.moisture_model isa DryModel) && tracer_upwinding != Val(:none)
-        ᶜq_tot = ᶜspecific.q_tot
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
         vtt = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, float(dt), tracer_upwinding)
         vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, float(dt), Val(:none))
         @. Yₜ.c.ρq_tot += vtt - vtt_central

src/prognostic_equations/edmfx_sgs_flux.jl

@@ -40,7 +40,7 @@ function edmfx_sgs_mass_flux_tendency!(
 
     n = n_mass_flux_subdomains(turbconv_model)
     (; edmfx_sgsflux_upwinding) = p.atmos.numerics
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶠu³, ᶜh_tot) = p.precomputed
     (; ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = p.precomputed
     (; ᶜρa⁰, ᶜρ⁰, ᶠu³⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰) = p.precomputed
     if (
@@ -236,7 +236,7 @@ function edmfx_sgs_mass_flux_tendency!(
     a_max = CAP.max_area(turbconv_params)
     n = n_mass_flux_subdomains(turbconv_model)
     (; edmfx_sgsflux_upwinding) = p.atmos.numerics
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶠu³, ᶜh_tot) = p.precomputed
     (; ᶜρaʲs, ᶜρʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs) = p.precomputed
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
@@ -450,7 +450,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; ᶜu, ᶜh_tot, ᶜspecific, ᶜtke⁰, ᶜmixing_length) = p.precomputed
+    (; ᶜu, ᶜh_tot, ᶜtke⁰, ᶜmixing_length) = p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
 

src/prognostic_equations/forcing/external_forcing.jl

@@ -337,7 +337,7 @@ function external_forcing_tendency!(
     # horizontal advection, vertical fluctuation, nudging, subsidence (need to add),
     (; params) = p
     thermo_params = CAP.thermodynamics_params(params)
-    (; ᶜspecific, ᶜts, ᶜh_tot) = p.precomputed
+    (; ᶜts, ᶜh_tot) = p.precomputed
     (;
         ᶜdTdt_fluc,
         ᶜdqtdt_fluc,
@@ -398,11 +398,12 @@ function external_forcing_tendency!(
         ᶜh_tot,
         Val{:first_order}(),
     )
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
     subsidence!(
         Yₜ.c.ρq_tot,
         Y.c.ρ,
         ᶠls_subsidence³,
-        ᶜspecific.q_tot,
+        ᶜq_tot,
         Val{:first_order}(),
     )
 

src/prognostic_equations/forcing/subsidence.jl

@@ -55,7 +55,7 @@ scalar quantity `χ`.
 Arguments:
 - `Yₜ`: The tendency state vector, modified in place.
 - `Y`: The current state vector (used for `Y.c.ρ`).
-- `p`: Cache containing parameters, precomputed fields (`ᶜh_tot`, `ᶜspecific`),
+- `p`: Cache containing parameters, precomputed fields (`ᶜh_tot`),
        atmospheric model configurations (`p.atmos.moisture_model`, `p.atmos.subsidence`),
        and scratch space.
 - `t`: Current simulation time (unused by this specific tendency calculation).
@@ -68,7 +68,7 @@ subsidence_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing    # No subsidence
 function subsidence_tendency!(Yₜ, Y, p, t, ::Subsidence)
     (; moisture_model) = p.atmos
     subsidence_profile = p.atmos.subsidence.prof
-    (; ᶜh_tot, ᶜspecific) = p.precomputed
+    (; ᶜh_tot) = p.precomputed
 
     ᶠz = Fields.coordinate_field(axes(Y.f)).z
     ᶠlg = Fields.local_geometry_field(Y.f)
@@ -78,26 +78,23 @@ function subsidence_tendency!(Yₜ, Y, p, t, ::Subsidence)
 
     # Large-scale subsidence
     subsidence!(Yₜ.c.ρe_tot, Y.c.ρ, ᶠsubsidence³, ᶜh_tot, Val{:first_order}())
-    subsidence!(
-        Yₜ.c.ρq_tot,
-        Y.c.ρ,
-        ᶠsubsidence³,
-        ᶜspecific.q_tot,
-        Val{:first_order}(),
-    )
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+    subsidence!(Yₜ.c.ρq_tot, Y.c.ρ, ᶠsubsidence³, ᶜq_tot, Val{:first_order}())
     if moisture_model isa NonEquilMoistModel
+        ᶜq_liq = @. lazy(specific(Y.c.ρq_liq, Y.c.ρ))
         subsidence!(
             Yₜ.c.ρq_liq,
             Y.c.ρ,
             ᶠsubsidence³,
-            ᶜspecific.q_liq,
+            ᶜq_liq,
             Val{:first_order}(),
         )
+        ᶜq_ice = @. lazy(specific(Y.c.ρq_ice, Y.c.ρ))
         subsidence!(
             Yₜ.c.ρq_ice,
             Y.c.ρ,
             ᶠsubsidence³,
-            ᶜspecific.q_ice,
+            ᶜq_ice,
             Val{:first_order}(),
         )
     end

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -139,15 +139,16 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     ᶜJ = Fields.local_geometry_field(Y.c).J
     ᶠJ = Fields.local_geometry_field(Y.f).J
     (; ᶠgradᵥ_ᶜΦ) = p.core
-    (; ᶜh_tot, ᶜspecific, ᶠu³, ᶜp) = p.precomputed
+    (; ᶜh_tot, ᶠu³, ᶜp) = p.precomputed
 
     @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠu³)
 
     # Central vertical advection of active tracers (e_tot and q_tot)
     vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜh_tot, dt, Val(:none))
     @. Yₜ.c.ρe_tot += vtt
     if !(moisture_model isa DryModel)
-        vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜspecific.q_tot, dt, Val(:none))
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+        vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜq_tot, dt, Val(:none))
         @. Yₜ.c.ρq_tot += vtt
     end
 

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -399,7 +399,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
     (; matrix) = cache
     (; params) = p
     (; ᶜΦ, ᶠgradᵥ_ᶜΦ) = p.core
-    (; ᶜspecific, ᶠu³, ᶜK, ᶜts, ᶜp, ᶜh_tot) = p.precomputed
+    (; ᶠu³, ᶜK, ᶜts, ᶜp, ᶜh_tot) = p.precomputed
     (;
         ∂ᶜK_∂ᶜuₕ,
         ∂ᶜK_∂ᶠu₃,
@@ -477,20 +477,22 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
     ∂ᶜρ_err_∂ᶠu₃ = matrix[@name(c.ρ), @name(f.u₃)]
     @. ∂ᶜρ_err_∂ᶠu₃ = dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(g³³(ᶠgⁱʲ))
 
-    tracer_info = (
-        (@name(c.ρe_tot), @name(ᶜh_tot)),
-        (@name(c.ρq_tot), @name(ᶜspecific.q_tot)),
-    )
-    MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, χ_name)
+    tracer_info = (@name(c.ρe_tot), @name(c.ρq_tot))
+    MatrixFields.unrolled_foreach(tracer_info) do ρχ_name
         MatrixFields.has_field(Y, ρχ_name) || return
-        ᶜχ = MatrixFields.get_field(p.precomputed, χ_name)
+        ᶜχ = if ρχ_name === @name(c.ρe_tot)
+            p.precomputed.ᶜh_tot
+        else
+            @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+        end
         if use_derivative(topography_flag)
             ∂ᶜρχ_err_∂ᶜuₕ = matrix[ρχ_name, @name(c.uₕ)]
+            @. ∂ᶜρχ_err_∂ᶜuₕ =
+                dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ)) ⋅
+                ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ))
         end
+
         ∂ᶜρχ_err_∂ᶠu₃ = matrix[ρχ_name, @name(f.u₃)]
-        use_derivative(topography_flag) && @. ∂ᶜρχ_err_∂ᶜuₕ =
-            dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ)) ⋅
-            ᶠwinterp_matrix(ᶜJ * ᶜρ) ⋅ DiagonalMatrixRow(g³ʰ(ᶜgⁱʲ))
         @. ∂ᶜρχ_err_∂ᶠu₃ =
             dtγ * ᶜadvection_matrix ⋅ DiagonalMatrixRow(ᶠinterp(ᶜχ) * g³³(ᶠgⁱʲ))
     end
@@ -536,12 +538,12 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
     end
 
     tracer_info = (
-        (@name(c.ρq_liq), @name(q_liq), @name(ᶜwₗ)),
-        (@name(c.ρq_ice), @name(q_ice), @name(ᶜwᵢ)),
-        (@name(c.ρq_rai), @name(q_rai), @name(ᶜwᵣ)),
-        (@name(c.ρq_sno), @name(q_sno), @name(ᶜwₛ)),
-        (@name(c.ρn_liq), @name(n_liq), @name(ᶜwnₗ)),
-        (@name(c.ρn_rai), @name(n_rai), @name(ᶜwnᵣ)),
+        (@name(c.ρq_liq), @name(ᶜwₗ)),
+        (@name(c.ρq_ice), @name(ᶜwᵢ)),
+        (@name(c.ρq_rai), @name(ᶜwᵣ)),
+        (@name(c.ρq_sno), @name(ᶜwₛ)),
+        (@name(c.ρn_liq), @name(ᶜwnₗ)),
+        (@name(c.ρn_rai), @name(ᶜwnᵣ)),
     )
     if !(p.atmos.moisture_model isa DryModel) || use_derivative(diffusion_flag)
         ∂ᶜρe_tot_err_∂ᶜρe_tot = matrix[@name(c.ρe_tot), @name(c.ρe_tot)]
@@ -583,13 +585,13 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
         #    DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅
         #    DiagonalMatrixRow(
         #        -1 / ᶜρ * ifelse(
-        #            ᶜspecific.q_tot == 0,
+        #            specific(Y.c.ρq_tot, Y.c.ρ) == 0,
         #            (Geometry.WVector(FT(0)),),
-        #            p.precomputed.ᶜwₜqₜ / ᶜspecific.q_tot,
+        #            p.precomputed.ᶜwₜqₜ / specific(Y.c.ρq_tot, Y.c.ρ),
         #        ),
         #    ) - (I,)
 
-        MatrixFields.unrolled_foreach(tracer_info) do (ρχₚ_name, _, wₚ_name)
+        MatrixFields.unrolled_foreach(tracer_info) do (ρχₚ_name, wₚ_name)
             MatrixFields.has_field(Y, ρχₚ_name) || return
             ∂ᶜρχₚ_err_∂ᶜρχₚ = matrix[ρχₚ_name, ρχₚ_name]
             ᶜwₚ = MatrixFields.get_field(p.precomputed, wₚ_name)
@@ -651,13 +653,14 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
                 dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow(1 / ᶜρ)
         end
 
-        MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, χ_name, _)
+        MatrixFields.unrolled_foreach(tracer_info) do (ρχ_name, _)
             MatrixFields.has_field(Y, ρχ_name) || return
-            ᶜχ = MatrixFields.get_field(ᶜspecific, χ_name)
+            ᶜρχ = MatrixFields.get_field(Y, ρχ_name)
+            ᶜχ = @. lazy(specific(ᶜρχ, Y.c.ρ))
             ∂ᶜρχ_err_∂ᶜρ = matrix[ρχ_name, @name(c.ρ)]
             ∂ᶜρχ_err_∂ᶜρχ = matrix[ρχ_name, ρχ_name]
             ᶜtridiagonal_matrix_scalar = ifelse(
-                χ_name in (@name(q_rai), @name(q_sno), @name(n_rai)),
+                ρχ_name in (@name(c.ρq_rai), @name(c.ρq_sno), @name(c.ρn_rai)),
                 ᶜdiffusion_h_matrix_scaled,
                 ᶜdiffusion_h_matrix,
             )