Skip to content

Commit 1d2a30d

Browse files
szy21szy21
authored
Merge pull request #3747 from CliMA/zs/jac_u3_qt
modify the derivative of p and \rho with respect to qt
2 parents 7954db9 + 4ff8085 commit 1d2a30d

File tree

3 files changed

+69
-14
lines changed

3 files changed

+69
-14
lines changed

‎.buildkite/pipeline.yml

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -816,7 +816,6 @@ steps:
816816
artifact_paths: "prognostic_edmfx_trmm_column/output_active/*"
817817
agents:
818818
slurm_mem: 20GB
819-
soft_fail: true
820819

821820
- label: ":genie: Prognostic EDMFX TRMM with 0-moment in a column"
822821
command: >

‎reproducibility_tests/ref_counter.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
232
1+
233
22

33
# **README**
44
#
@@ -21,6 +21,9 @@
2121

2222
#=
2323
24+
233
25+
- modify the derivative of p and \rho with respect to qt
26+
2427
232
2528
- Use lazy broadcasting instead of temp scalar in implicit solver for kappa_m vars,
2629
which fixes a bug that the temp scalar is updated before it is reused.

‎src/prognostic_equations/implicit/implicit_solver.jl

Lines changed: 65 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -606,7 +606,9 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
606606
Δcv_v = FT(CAP.cv_v(params)) - cv_d
607607
T_0 = FT(CAP.T_0(params))
608608
R_d = FT(CAP.R_d(params))
609+
ΔR_v = FT(CAP.R_v(params)) - R_d
609610
cp_d = FT(CAP.cp_d(params))
611+
Δcp_v = FT(CAP.cp_v(params)) - cp_d
610612
# This term appears a few times in the Jacobian, and is technically
611613
# minus ∂e_int_∂q_tot
612614
∂e_int_∂q_tot = T_0 * (Δcv_v - R_d) - FT(CAP.e_int_v0(params))
@@ -625,6 +627,14 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
625627
TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts),
626628
)
627629

630+
# Using abs2 because ^2 results in allocation
631+
ᶜ∂kappa_m∂q_tot = @. lazy(
632+
(
633+
ΔR_v * TD.cv_m(thermo_params, ᶜts) -
634+
Δcv_v * TD.gas_constant_air(thermo_params, ᶜts)
635+
) / abs2(TD.cv_m(thermo_params, ᶜts)),
636+
)
637+
628638
if use_derivative(topography_flag)
629639
@. ∂ᶜK_∂ᶜuₕ = DiagonalMatrixRow(
630640
adjoint(CTh(ᶜuₕ)) + adjoint(ᶜinterp(ᶠu₃)) * g³ʰ(ᶜgⁱʲ),
@@ -681,7 +691,13 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
681691
if MatrixFields.has_field(Y, @name(c.ρq_tot))
682692
∂ᶠu₃_err_∂ᶜρq_tot = matrix[@name(f.u₃), @name(c.ρq_tot)]
683693
@. ∂ᶠu₃_err_∂ᶜρq_tot =
684-
dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m * ∂e_int_∂q_tot)
694+
dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow((
695+
ᶜkappa_m * ∂e_int_∂q_tot +
696+
ᶜ∂kappa_m∂q_tot * (
697+
cp_d * T_0 + ᶜspecific.e_tot - ᶜK - ᶜΦ +
698+
∂e_int_∂q_tot * ᶜspecific.q_tot
699+
)
700+
))
685701
end
686702

687703
∂ᶠu₃_err_∂ᶜuₕ = matrix[@name(f.u₃), @name(c.uₕ)]
@@ -806,8 +822,13 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
806822
∂ᶜρe_tot_err_∂ᶜρq_tot = matrix[@name(c.ρe_tot), @name(c.ρq_tot)]
807823
∂ᶜρq_tot_err_∂ᶜρ = matrix[@name(c.ρq_tot), @name(c.ρ)]
808824
@. ∂ᶜρe_tot_err_∂ᶜρq_tot +=
809-
dtγ * ᶜdiffusion_h_matrix ⋅
810-
DiagonalMatrixRow(ᶜkappa_m * ∂e_int_∂q_tot / ᶜρ)
825+
dtγ * ᶜdiffusion_h_matrix ⋅ DiagonalMatrixRow((
826+
ᶜkappa_m * ∂e_int_∂q_tot / ᶜρ +
827+
ᶜ∂kappa_m∂q_tot * (
828+
cp_d * T_0 + ᶜspecific.e_tot - ᶜK - ᶜΦ +
829+
∂e_int_∂q_tot * ᶜspecific.q_tot
830+
)
831+
))
811832
@. ∂ᶜρq_tot_err_∂ᶜρ =
812833
dtγ * ᶜdiffusion_h_matrix ⋅
813834
DiagonalMatrixRow(-(ᶜspecific.q_tot) / ᶜρ)
@@ -904,6 +925,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
904925
TD.cv_m(thermo_params, ᶜtsʲs.:(1)),
905926
)
906927

928+
# Note this is the derivative of R_m / cp_m with respect to q_tot
929+
# but we call it ∂kappa_m∂q_totʲ
930+
ᶜ∂kappa_m∂q_totʲ = @. lazy(
931+
(
932+
ΔR_v * TD.cp_m(thermo_params, ᶜtsʲs.:(1)) -
933+
Δcp_v * TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1))
934+
) / abs2(TD.cp_m(thermo_params, ᶜtsʲs.:(1))),
935+
)
936+
907937
∂ᶜq_totʲ_err_∂ᶜq_totʲ =
908938
matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)]
909939
@. ∂ᶜq_totʲ_err_∂ᶜq_totʲ =
@@ -932,10 +962,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
932962
@. ∂ᶜmseʲ_err_∂ᶜq_totʲ =
933963
dtγ * (
934964
-DiagonalMatrixRow(
935-
adjoint(ᶜinterp(ᶠu³ʲs.:(1))) *
936-
ᶜgradᵥ_ᶠΦ *
937-
Y.c.ρ *
938-
ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) * ∂e_int_∂q_tot,
965+
adjoint(ᶜinterp(ᶠu³ʲs.:(1))) * ᶜgradᵥ_ᶠΦ * Y.c.ρ / ᶜp *
966+
(
967+
(ᶜkappa_mʲ / (ᶜkappa_mʲ + 1) * ∂e_int_∂q_tot) +
968+
ᶜ∂kappa_m∂q_totʲ * (
969+
Y.c.sgsʲs.:(1).mse - ᶜΦ +
970+
cp_d * T_0 +
971+
∂e_int_∂q_tot * Y.c.sgsʲs.:(1).q_tot
972+
)
973+
),
939974
)
940975
)
941976

@@ -977,8 +1012,14 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
9771012
),
9781013
) / ᶠJ,
9791014
) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
980-
ᶜJ * ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
981-
∂e_int_∂q_tot,
1015+
ᶜJ * (ᶜρʲs.:(1))^2 / ᶜp * (
1016+
ᶜkappa_mʲ / (ᶜkappa_mʲ + 1) * ∂e_int_∂q_tot +
1017+
ᶜ∂kappa_m∂q_totʲ * (
1018+
Y.c.sgsʲs.:(1).mse - ᶜΦ +
1019+
cp_d * T_0 +
1020+
∂e_int_∂q_tot * Y.c.sgsʲs.:(1).q_tot
1021+
)
1022+
),
9821023
)
9831024
@. ᶠbidiagonal_matrix_ct3_2 =
9841025
DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ) ⋅
@@ -1045,8 +1086,14 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
10451086
dtγ * DiagonalMatrixRow(
10461087
ᶠgradᵥ_ᶜΦ * ᶠinterp(Y.c.ρ) / (ᶠinterp(ᶜρʲs.:(1)))^2,
10471088
) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
1048-
ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
1049-
∂e_int_∂q_tot,
1089+
(ᶜρʲs.:(1))^2 / ᶜp * (
1090+
ᶜkappa_mʲ / (ᶜkappa_mʲ + 1) * ∂e_int_∂q_tot +
1091+
ᶜ∂kappa_m∂q_totʲ * (
1092+
Y.c.sgsʲs.:(1).mse - ᶜΦ +
1093+
cp_d * T_0 +
1094+
∂e_int_∂q_tot * Y.c.sgsʲs.:(1).q_tot
1095+
)
1096+
),
10501097
)
10511098
∂ᶠu₃ʲ_err_∂ᶜmseʲ =
10521099
matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)]
@@ -1154,7 +1201,13 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, t)
11541201

11551202
@. ∂ᶜρe_tot_err_∂ᶜρq_tot +=
11561203
dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅
1157-
DiagonalMatrixRow(ᶜkappa_m * ∂e_int_∂q_tot / ᶜρ)
1204+
DiagonalMatrixRow((
1205+
ᶜkappa_m * ∂e_int_∂q_tot / ᶜρ +
1206+
ᶜ∂kappa_m∂q_tot * (
1207+
cp_d * T_0 + ᶜspecific.e_tot - ᶜK - ᶜΦ +
1208+
∂e_int_∂q_tot * ᶜspecific.q_tot
1209+
)
1210+
))
11581211

11591212
@. ∂ᶜρe_tot_err_∂ᶜρe_tot +=
11601213
dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar ⋅

0 commit comments

Comments
 (0)