@@ -600,15 +600,22 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
600600 # ),
601601 # ) - (I,)
602602
603+ # This scratch variable computation could be skipped if no tracers are present
604+ @. p. scratch. ᶜbidiagonal_adjoint_matrix_c3 =
605+ dtγ * (- ClimaAtmos. ᶜprecipdivᵥ_matrix ()) ⋅
606+ DiagonalMatrixRow (ClimaAtmos. ᶠinterp (ᶜρ * ᶜJ) / ᶠJ)
607+
603608 MatrixFields. unrolled_foreach (tracer_info) do (ρχₚ_name, wₚ_name, _)
604609 MatrixFields. has_field (Y, ρχₚ_name) || return
605610 ∂ᶜρχₚ_err_∂ᶜρχₚ = matrix[ρχₚ_name, ρχₚ_name]
606611 ᶜwₚ = MatrixFields. get_field (p. precomputed, wₚ_name)
612+ # TODO : come up with read-able names for the intermediate computations...
613+ @. p. scratch. ᶠband_matrix_wvec =
614+ ClimaAtmos. ᶠright_bias_matrix () ⋅
615+ DiagonalMatrixRow (ClimaCore. Geometry. WVector (- (ᶜwₚ) / ᶜρ))
607616 @. ∂ᶜρχₚ_err_∂ᶜρχₚ =
608- dtγ * - (ᶜprecipdivᵥ_matrix ()) ⋅
609- DiagonalMatrixRow (ᶠinterp (ᶜρ * ᶜJ) / ᶠJ) ⋅
610- ᶠright_bias_matrix () ⋅
611- DiagonalMatrixRow (- Geometry. WVector (ᶜwₚ) / ᶜρ) - (I,)
617+ p. scratch. ᶜbidiagonal_adjoint_matrix_c3 ⋅
618+ p. scratch. ᶠband_matrix_wvec - (I,)
612619 end
613620
614621 end
@@ -634,24 +641,31 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
634641 ᶜtke⁰ = @. lazy (specific (Y. c. sgs⁰. ρatke, Y. c. ρ))
635642 ᶜmixing_length_field = p. scratch. ᶜtemp_scalar_3
636643 ᶜmixing_length_field .= ᶜmixing_length (Y, p)
637- ᶜK_u = @. lazy (eddy_viscosity (turbconv_params, ᶜtke⁰, ᶜmixing_length_field))
644+ ᶜK_u = p. scratch. ᶜtemp_scalar_4
645+ @. ᶜK_u = eddy_viscosity (turbconv_params, ᶜtke⁰, ᶜmixing_length_field)
638646 ᶜprandtl_nvec = @. lazy (
639647 turbulent_prandtl_number (params, ᶜlinear_buoygrad, ᶜstrain_rate_norm),
640648 )
641- ᶜK_h = @. lazy (eddy_diffusivity (ᶜK_u, ᶜprandtl_nvec))
649+ ᶜK_h = p. scratch. ᶜtemp_scalar_6
650+ @. ᶜK_h = eddy_diffusivity (ᶜK_u, ᶜprandtl_nvec)
642651 end
643652
653+
654+ @. p. scratch. ᶜbidiagonal_adjoint_matrix_c3 =
655+ ᶜadvdivᵥ_matrix () ⋅ DiagonalMatrixRow (ᶠinterp (ᶜρ) * ᶠinterp (ᶜK_h))
644656 @. ᶜdiffusion_h_matrix =
645- ᶜadvdivᵥ_matrix () ⋅ DiagonalMatrixRow ( ᶠinterp (ᶜρ) * ᶠinterp (ᶜK_h)) ⋅
646- ᶠgradᵥ_matrix ()
657+ p . scratch . ᶜbidiagonal_adjoint_matrix_c3 ⋅ ᶠgradᵥ_matrix ()
658+
647659 if (
648660 MatrixFields. has_field (Y, @name (c. sgs⁰. ρatke)) ||
649661 ! isnothing (p. atmos. turbconv_model) ||
650662 ! disable_momentum_vertical_diffusion (p. atmos. vertical_diffusion)
651663 )
652- @. ᶜdiffusion_u_matrix =
664+ @. p . scratch . ᶜbidiagonal_adjoint_matrix_c3 =
653665 ᶜadvdivᵥ_matrix () ⋅
654- DiagonalMatrixRow (ᶠinterp (ᶜρ) * ᶠinterp (ᶜK_u)) ⋅ ᶠgradᵥ_matrix ()
666+ DiagonalMatrixRow (ᶠinterp (ᶜρ) * ᶠinterp (ᶜK_u))
667+ @. ᶜdiffusion_u_matrix =
668+ p. scratch. ᶜbidiagonal_adjoint_matrix_c3 ⋅ ᶠgradᵥ_matrix ()
655669 end
656670
657671 ∂ᶜρe_tot_err_∂ᶜρ = matrix[@name (c. ρe_tot), @name (c. ρ)]
@@ -805,15 +819,14 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
805819 ) - (I,)
806820 ∂ᶜq_totʲ_err_∂ᶠu₃ʲ =
807821 matrix[@name (c. sgsʲs.:(1 ). q_tot), @name (f. sgsʲs.:(1 ). u₃)]
822+ @. p. scratch. ᶜbidiagonal_adjoint_matrix_c3 =
823+ - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
824+ ᶠset_upwind_bcs (
825+ ᶠupwind (CT3 (sign (ᶠu³ʲ_data)), Y. c. sgsʲs.:(1 ). q_tot),
826+ ) * adjoint (C3 (sign (ᶠu³ʲ_data))),
827+ ) + DiagonalMatrixRow (Y. c. sgsʲs.:(1 ). q_tot) ⋅ ᶜadvdivᵥ_matrix ()
808828 @. ∂ᶜq_totʲ_err_∂ᶠu₃ʲ =
809- dtγ * (
810- - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
811- ᶠset_upwind_bcs (
812- ᶠupwind (CT3 (sign (ᶠu³ʲ_data)), Y. c. sgsʲs.:(1 ). q_tot),
813- ) * adjoint (C3 (sign (ᶠu³ʲ_data))),
814- ) +
815- DiagonalMatrixRow (Y. c. sgsʲs.:(1 ). q_tot) ⋅ ᶜadvdivᵥ_matrix ()
816- ) ⋅ DiagonalMatrixRow (g³³ (ᶠgⁱʲ))
829+ dtγ * p. scratch. ᶜbidiagonal_adjoint_matrix_c3 ⋅ DiagonalMatrixRow (g³³ (ᶠgⁱʲ))
817830
818831 if p. atmos. moisture_model isa NonEquilMoistModel && (
819832 p. atmos. microphysics_model isa Microphysics1Moment ||
@@ -938,15 +951,14 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
938951 ) - (I,)
939952 ∂ᶜmseʲ_err_∂ᶠu₃ʲ =
940953 matrix[@name (c. sgsʲs.:(1 ). mse), @name (f. sgsʲs.:(1 ). u₃)]
954+ @. p. scratch. ᶜbidiagonal_adjoint_matrix_c3 =
955+ - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
956+ ᶠset_upwind_bcs (
957+ ᶠupwind (CT3 (sign (ᶠu³ʲ_data)), Y. c. sgsʲs.:(1 ). mse),
958+ ) * adjoint (C3 (sign (ᶠu³ʲ_data))),
959+ ) + DiagonalMatrixRow (Y. c. sgsʲs.:(1 ). mse) ⋅ ᶜadvdivᵥ_matrix ()
941960 @. ∂ᶜmseʲ_err_∂ᶠu₃ʲ =
942- dtγ * (
943- - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
944- ᶠset_upwind_bcs (
945- ᶠupwind (CT3 (sign (ᶠu³ʲ_data)), Y. c. sgsʲs.:(1 ). mse),
946- ) * adjoint (C3 (sign (ᶠu³ʲ_data))),
947- ) +
948- DiagonalMatrixRow (Y. c. sgsʲs.:(1 ). mse) ⋅ ᶜadvdivᵥ_matrix ()
949- ) ⋅ DiagonalMatrixRow (g³³ (ᶠgⁱʲ))
961+ dtγ * p. scratch. ᶜbidiagonal_adjoint_matrix_c3 ⋅ DiagonalMatrixRow (g³³ (ᶠgⁱʲ))
950962
951963 ∂ᶜρaʲ_err_∂ᶜq_totʲ =
952964 matrix[@name (c. sgsʲs.:(1 ). ρa), @name (c. sgsʲs.:(1 ). q_tot)]
@@ -1081,10 +1093,11 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
10811093 # vertical diffusion of updrafts
10821094 if use_derivative (sgs_vertdiff_flag)
10831095 α_vert_diff_tracer = CAP. α_vert_diff_tracer (params)
1084- @. ᶜdiffusion_h_matrix =
1096+ @. p . scratch . ᶜbidiagonal_adjoint_matrix_c3 =
10851097 ᶜadvdivᵥ_matrix () ⋅
1086- DiagonalMatrixRow (ᶠinterp (ᶜρʲs.:(1 )) * ᶠinterp (ᶜK_h)) ⋅
1087- ᶠgradᵥ_matrix ()
1098+ DiagonalMatrixRow (ᶠinterp (ᶜρʲs.:(1 )) * ᶠinterp (ᶜK_h))
1099+ @. ᶜdiffusion_h_matrix =
1100+ p. scratch. ᶜbidiagonal_adjoint_matrix_c3 ⋅ ᶠgradᵥ_matrix ()
10881101
10891102 @. ∂ᶜmseʲ_err_∂ᶜmseʲ +=
10901103 dtγ * DiagonalMatrixRow (1 / ᶜρʲs.:(1 )) ⋅ ᶜdiffusion_h_matrix
@@ -1151,6 +1164,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
11511164
11521165 # non-hydrostatic pressure drag
11531166 # (quadratic drag term treated implicitly, buoyancy term explicitly)
1167+
11541168 if use_derivative (sgs_nh_pressure_flag)
11551169 (; ᶠu₃⁰) = p. precomputed
11561170 α_d = CAP. pressure_normalmode_drag_coeff (turbconv_params)
@@ -1173,6 +1187,8 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
11731187 (ᶠinterp (ᶜρ * ᶜJ) / ᶠJ) * (ᶠu³ʲs.:(1 ) - ᶠu³),
11741188 ) ⋅ ᶠinterp_matrix () ⋅
11751189 DiagonalMatrixRow (Y. c. sgsʲs.:(1 ). ρa / ᶜρʲs.:(1 ))
1190+ @. p. scratch. ᶜtridiagonal_matrix_scalar =
1191+ dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar
11761192
11771193 # Derivative of total energy tendency with respect to updraft MSE
11781194 # # grid-mean ρe_tot
@@ -1189,9 +1205,8 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
11891205 ) / abs2 (TD. cv_m (thermo_params, ᶜts))
11901206
11911207 ᶜq_tot = @. lazy (specific (Y. c. ρq_tot, Y. c. ρ))
1192-
11931208 @. ∂ᶜρe_tot_err_∂ᶜρ +=
1194- dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
1209+ p . scratch . ᶜtridiagonal_matrix_scalar ⋅
11951210 DiagonalMatrixRow (
11961211 (
11971212 - (1 + ᶜkappa_m) * ᶜe_tot -
@@ -1200,7 +1215,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
12001215 )
12011216
12021217 @. ∂ᶜρe_tot_err_∂ᶜρq_tot +=
1203- dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
1218+ p . scratch . ᶜtridiagonal_matrix_scalar ⋅
12041219 DiagonalMatrixRow ((
12051220 ᶜkappa_m * ∂e_int_∂q_tot / ᶜρ +
12061221 ᶜ∂kappa_m∂q_tot * (
@@ -1210,27 +1225,27 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
12101225 ))
12111226
12121227 @. ∂ᶜρe_tot_err_∂ᶜρe_tot +=
1213- dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
1228+ p . scratch . ᶜtridiagonal_matrix_scalar ⋅
12141229 DiagonalMatrixRow ((1 + ᶜkappa_m) / ᶜρ)
12151230
12161231 ∂ᶜρe_tot_err_∂ᶜmseʲ =
12171232 matrix[@name (c. ρe_tot), @name (c. sgsʲs.:(1 ). mse)]
12181233 @. ∂ᶜρe_tot_err_∂ᶜmseʲ =
1219- - (dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar )
1234+ - (p . scratch . ᶜtridiagonal_matrix_scalar )
12201235
12211236 # # grid-mean ρq_tot
12221237 @. ∂ᶜρq_tot_err_∂ᶜρ +=
1223- dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
1238+ p . scratch . ᶜtridiagonal_matrix_scalar ⋅
12241239 DiagonalMatrixRow (- (ᶜq_tot) / ᶜρ)
12251240
12261241 @. ∂ᶜρq_tot_err_∂ᶜρq_tot +=
1227- dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
1242+ p . scratch . ᶜtridiagonal_matrix_scalar ⋅
12281243 DiagonalMatrixRow (1 / ᶜρ)
12291244
12301245 ∂ᶜρq_tot_err_∂ᶜq_totʲ =
12311246 matrix[@name (c. ρq_tot), @name (c. sgsʲs.:(1 ). q_tot)]
12321247 @. ∂ᶜρq_tot_err_∂ᶜq_totʲ =
1233- - (dtγ * ᶜadvdivᵥ_matrix () ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar )
1248+ - (p . scratch . ᶜtridiagonal_matrix_scalar )
12341249
12351250 # grid-mean ∂/∂(u₃ʲ)
12361251 ∂ᶜρe_tot_err_∂ᶠu₃ = matrix[@name (c. ρe_tot), @name (f. u₃)]
@@ -1256,45 +1271,46 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
12561271 ) / ᶠJ * (g³³ (ᶠgⁱʲ)),
12571272 )
12581273
1274+ @. p. scratch. ᶠdiagonal_matrix_ct3xct3 = DiagonalMatrixRow (
1275+ ᶠinterp (
1276+ (Y. c. sgsʲs.:(1 ). q_tot - ᶜq_tot) *
1277+ ᶜρʲs.:(1 ) *
1278+ ᶜJ *
1279+ draft_area (Y. c. sgsʲs.:(1 ). ρa, ᶜρʲs.:(1 )),
1280+ ) / ᶠJ * (g³³ (ᶠgⁱʲ)),
1281+ )
1282+
12591283 ∂ᶜρq_tot_err_∂ᶠu₃ = matrix[@name (c. ρq_tot), @name (f. u₃)]
12601284 @. ∂ᶜρq_tot_err_∂ᶠu₃ +=
1261- dtγ * ᶜadvdivᵥ_matrix () ⋅ DiagonalMatrixRow (
1262- ᶠinterp (
1263- (Y. c. sgsʲs.:(1 ). q_tot - ᶜq_tot) *
1264- ᶜρʲs.:(1 ) *
1265- ᶜJ *
1266- draft_area (Y. c. sgsʲs.:(1 ). ρa, ᶜρʲs.:(1 )),
1267- ) / ᶠJ * (g³³ (ᶠgⁱʲ)),
1268- )
1285+ dtγ * ᶜadvdivᵥ_matrix () ⋅ p. scratch. ᶠdiagonal_matrix_ct3xct3
12691286
12701287 ∂ᶜρq_tot_err_∂ᶠu₃ʲ =
12711288 matrix[@name (c. ρq_tot), @name (f. sgsʲs.:(1 ). u₃)]
12721289 @. ∂ᶜρq_tot_err_∂ᶠu₃ʲ =
1273- dtγ * - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
1274- ᶠinterp (
1275- (Y. c. sgsʲs.:(1 ). q_tot - ᶜq_tot) *
1276- ᶜρʲs.:(1 ) *
1277- ᶜJ *
1278- draft_area (Y. c. sgsʲs.:(1 ). ρa, ᶜρʲs.:(1 )),
1279- ) / ᶠJ * (g³³ (ᶠgⁱʲ)),
1280- )
1290+ dtγ * - (ᶜadvdivᵥ_matrix ()) ⋅ p. scratch. ᶠdiagonal_matrix_ct3xct3
12811291
12821292 # grid-mean ∂/∂(rho*a)
12831293 ∂ᶜρe_tot_err_∂ᶜρa =
12841294 matrix[@name (c. ρe_tot), @name (c. sgsʲs.:(1 ). ρa)]
1295+ @. p. scratch. ᶠtemp_CT3_2 =
1296+ (ᶠu³ʲs.:(1 ) - ᶠu³) *
1297+ ᶠinterp ((Y. c. sgsʲs.:(1 ). mse + ᶜKʲs.:(1 ) - ᶜh_tot)) / ᶠJ
1298+ @. p. scratch. ᶜbidiagonal_matrix_scalar =
1299+ dtγ * - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (p. scratch. ᶠtemp_CT3_2)
12851300 @. ∂ᶜρe_tot_err_∂ᶜρa =
1286- dtγ * - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
1287- (ᶠu³ʲs.:(1 ) - ᶠu³) *
1288- ᶠinterp ((Y. c. sgsʲs.:(1 ). mse + ᶜKʲs.:(1 ) - ᶜh_tot)) / ᶠJ,
1289- ) ⋅ ᶠinterp_matrix () ⋅ DiagonalMatrixRow (ᶜJ)
1301+ p. scratch. ᶜbidiagonal_matrix_scalar ⋅ ᶠinterp_matrix () ⋅
1302+ DiagonalMatrixRow (ᶜJ)
12901303
12911304 ∂ᶜρq_tot_err_∂ᶜρa =
12921305 matrix[@name (c. ρq_tot), @name (c. sgsʲs.:(1 ). ρa)]
1306+ @. p. scratch. ᶠtemp_CT3_2 =
1307+ (ᶠu³ʲs.:(1 ) - ᶠu³) *
1308+ ᶠinterp ((Y. c. sgsʲs.:(1 ). q_tot - ᶜq_tot)) / ᶠJ
1309+ @. p. scratch. ᶜbidiagonal_matrix_scalar =
1310+ dtγ * - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (p. scratch. ᶠtemp_CT3_2)
12931311 @. ∂ᶜρq_tot_err_∂ᶜρa =
1294- dtγ * - (ᶜadvdivᵥ_matrix ()) ⋅ DiagonalMatrixRow (
1295- (ᶠu³ʲs.:(1 ) - ᶠu³) *
1296- ᶠinterp ((Y. c. sgsʲs.:(1 ). q_tot - ᶜq_tot)) / ᶠJ,
1297- ) ⋅ ᶠinterp_matrix () ⋅ DiagonalMatrixRow (ᶜJ)
1312+ p. scratch. ᶜbidiagonal_matrix_scalar ⋅ ᶠinterp_matrix () ⋅
1313+ DiagonalMatrixRow (ᶜJ)
12981314
12991315 # grid-mean tracers
13001316 if p. atmos. moisture_model isa NonEquilMoistModel && (
0 commit comments