Add rho*a and u_{3,m} updraft jacobian terms

Author: costachris

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-225
+226
 
 # **README**
 #
@@ -21,6 +21,9 @@
 
 #=
 
+226
+- Add updraft rho*a and u_{3,m} jacobian terms
+
 225
 - Move nonhydrostatic pressure drag calculation to precomputed quantities and 
   remove one reproducibility job

src/prognostic_equations/implicit/implicit_solver.jl

@@ -328,6 +328,14 @@ function ImplicitEquationJacobian(
                     similar(Y.c, TridiagonalRow),
                 (@name(c.ρq_tot), @name(c.sgsʲs.:(1).q_tot)) =>
                     similar(Y.c, TridiagonalRow),
+                (@name(c.ρe_tot), @name(f.sgsʲs.:(1).u₃)) =>
+                    similar(Y.c, BidiagonalRow_ACT3),
+                (@name(c.ρq_tot), @name(f.sgsʲs.:(1).u₃)) =>
+                    similar(Y.c, BidiagonalRow_ACT3),
+                (@name(c.ρe_tot), @name(c.sgsʲs.:(1).ρa)) =>
+                    similar(Y.c, TridiagonalRow),
+                (@name(c.ρq_tot), @name(c.sgsʲs.:(1).ρa)) =>
+                    similar(Y.c, TridiagonalRow),
             )
         else
             ()
@@ -525,6 +533,11 @@ NVTX.@annotate function Wfact!(A, Y, p, dtγ, t)
             p.atmos.turbconv_model isa PrognosticEDMFX ?
             (; p.precomputed.ᶠu₃⁰,) : (;)
         )...,
+        (
+            use_derivative(A.sgs_mass_flux_flag) &&
+            p.atmos.turbconv_model isa PrognosticEDMFX ?
+            (; p.precomputed.ᶜKʲs) : (;)
+        )...,
         p.core.ᶜΦ,
         p.core.ᶠgradᵥ_ᶜΦ,
         p.scratch.ᶜtemp_scalar,
@@ -864,6 +877,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
     end
 
     if p.atmos.turbconv_model isa PrognosticEDMFX
+
         if use_derivative(sgs_advection_flag)
             (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = p
             (; bdmr_l, bdmr_r, bdmr) = p
@@ -1068,7 +1082,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
             # add updraft mass flux contributions to grid-mean
             if use_derivative(sgs_mass_flux_flag)
 
-                (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs, ᶠu³) = p
+                (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs, ᶠu³, ᶜKʲs) = p
                 (; bdmr_l, bdmr_r, bdmr) = p
                 is_third_order = edmfx_upwinding == Val(:third_order)
                 ᶠupwind = is_third_order ? ᶠupwind3 : ᶠupwind1
@@ -1088,8 +1102,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
 
                 # Jacobian contributions of updraft massflux to grid-mean
 
-                ∂ᶜupdraft_mass_flux_∂ᶜh_tot = ᶠbidiagonal_matrix_ct3
-                @. ∂ᶜupdraft_mass_flux_∂ᶜh_tot =
+                ∂ᶜupdraft_mass_flux_∂ᶜscalar = ᶠbidiagonal_matrix_ct3
+                @. ∂ᶜupdraft_mass_flux_∂ᶜscalar =
                     DiagonalMatrixRow(
                         (ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) * (ᶠu³ʲs.:(1) - ᶠu³),
                     ) ⋅ ᶠinterp_matrix() ⋅
@@ -1098,7 +1112,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
                 # Derivative of total energy tendency with respect to updraft MSE
                 ## grid-mean ρe_tot
                 @. ∂ᶜρe_tot_err_∂ᶜρ +=
-                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜh_tot ⋅
+                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
                     DiagonalMatrixRow(
                         (
                             -(1 + ᶜkappa_m) * ᶜspecific.e_tot -
@@ -1107,34 +1121,76 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
                     )
 
                 @. ∂ᶜρe_tot_err_∂ᶜρq_tot +=
-                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜh_tot ⋅
+                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
                     DiagonalMatrixRow(ᶜkappa_m * ∂e_int_∂q_tot / ᶜρ)
 
                 @. ∂ᶜρe_tot_err_∂ᶜρe_tot +=
-                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜh_tot ⋅
+                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
                     DiagonalMatrixRow((1 + ᶜkappa_m) / ᶜρ)
 
                 ∂ᶜρe_tot_err_∂ᶜmseʲ =
                     matrix[@name(c.ρe_tot), @name(c.sgsʲs.:(1).mse)]
                 @. ∂ᶜρe_tot_err_∂ᶜmseʲ =
-                    -(dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜh_tot)
+                    -(dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar)
 
                 ## grid-mean ρq_tot
-                ∂ᶜupdraft_mass_flux_∂ᶜqtot = ∂ᶜupdraft_mass_flux_∂ᶜh_tot
                 ∂ᶜρq_tot_err_∂ᶜρ = matrix[@name(c.ρq_tot), @name(c.ρ)]
                 @. ∂ᶜρq_tot_err_∂ᶜρ +=
-                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜqtot ⋅
+                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
                     DiagonalMatrixRow(-(ᶜspecific.q_tot) / ᶜρ)
 
                 ∂ᶜρq_tot_err_∂ᶜρq_tot = matrix[@name(c.ρq_tot), @name(c.ρq_tot)]
                 @. ∂ᶜρq_tot_err_∂ᶜρq_tot +=
-                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜqtot ⋅
+                    dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
                     DiagonalMatrixRow(1 / ᶜρ)
 
                 ∂ᶜρq_tot_err_∂ᶜq_totʲ =
                     matrix[@name(c.ρq_tot), @name(c.sgsʲs.:(1).q_tot)]
                 @. ∂ᶜρq_tot_err_∂ᶜq_totʲ =
-                    -(dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜqtot)
+                    -(dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar)
+
+                # grid-mean ∂/∂(u₃ʲ)
+                ∂ᶜρe_tot_err_∂ᶠu₃ʲ =
+                    matrix[@name(c.ρe_tot), @name(f.sgsʲs.:(1).u₃)]
+                @. ∂ᶜρe_tot_err_∂ᶠu₃ʲ =
+                    dtγ * -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                        ᶠinterp(
+                            (Y.c.sgsʲs.:(1).mse + ᶜKʲs.:(1) - ᶜh_tot) *
+                            ᶜρʲs.:(1) *
+                            ᶜJ *
+                            draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
+                        ) / ᶠJ * (g³³(ᶠgⁱʲ)),
+                    )
+
+
+                ∂ᶜρq_tot_err_∂ᶠu₃ʲ =
+                    matrix[@name(c.ρq_tot), @name(f.sgsʲs.:(1).u₃)]
+                @. ∂ᶜρq_tot_err_∂ᶠu₃ʲ =
+                    dtγ * -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                        ᶠinterp(
+                            (Y.c.sgsʲs.:(1).q_tot - ᶜspecific.q_tot) *
+                            ᶜρʲs.:(1) *
+                            ᶜJ *
+                            draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
+                        ) / ᶠJ * (g³³(ᶠgⁱʲ)),
+                    )
+
+                # grid-mean ∂/∂(rho*a)
+                ∂ᶜρe_tot_err_∂ᶜρa =
+                    matrix[@name(c.ρe_tot), @name(c.sgsʲs.:(1).ρa)]
+                @. ∂ᶜρe_tot_err_∂ᶜρa =
+                    dtγ * -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                        (ᶠu³ʲs.:(1) - ᶠu³) *
+                        ᶠinterp((Y.c.sgsʲs.:(1).mse + ᶜKʲs.:(1) - ᶜh_tot)) / ᶠJ,
+                    ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(ᶜJ)
+
+                ∂ᶜρq_tot_err_∂ᶜρa =
+                    matrix[@name(c.ρq_tot), @name(c.sgsʲs.:(1).ρa)]
+                @. ∂ᶜρq_tot_err_∂ᶜρa =
+                    dtγ * -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                        (ᶠu³ʲs.:(1) - ᶠu³) *
+                        ᶠinterp((Y.c.sgsʲs.:(1).q_tot - ᶜspecific.q_tot)) / ᶠJ,
+                    ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(ᶜJ)
 
             end