Skip to content

Commit f3b1e5f

Browse files
oalcabesoalcabes
committed
fixing if else statement
1 parent 53e5133 commit f3b1e5f

File tree

1 file changed

+146
-146
lines changed

1 file changed

+146
-146
lines changed

‎src/prognostic_equations/implicit/manual_sparse_jacobian.jl

Lines changed: 146 additions & 146 deletions
Original file line numberDiff line numberDiff line change
@@ -555,184 +555,185 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
555555
DiagonalMatrixRow(-Geometry.WVector(ᶜwₚ) / ᶜρ) - (I,)
556556
end
557557

558-
if p.atmos.moisture_model isa NonEquilMoistModel &&
559-
use_derivative(noneq_cloud_formation_flag)
560-
p_vapₛₗ(tps, ts) = TD.saturation_vapor_pressure(tps, ts, TD.Liquid())
561-
p_vapₛᵢ(tps, ts) = TD.saturation_vapor_pressure(tps, ts, TD.Ice())
562-
563-
function ∂p_vapₛₗ_∂T(tps, ts)
564-
T = TD.air_temperature(tps, ts)
565-
Rᵥ = TD.Parameters.R_v(tps)
566-
Lᵥ = TD.latent_heat_vapor(tps, ts)
567-
return p_vapₛₗ(tps, ts) * Lᵥ / (Rᵥ * T^2)
568-
end
569-
function ∂p_vapₛᵢ_∂T(tps, ts)
570-
T = TD.air_temperature(tps, ts)
571-
Rᵥ = TD.Parameters.R_v(tps)
572-
Lₛ = TD.latent_heat_sublim(tps, ts)
573-
return p_vapₛᵢ(tps, ts) * Lₛ / (Rᵥ * T^2)
574-
end
558+
if p.atmos.moisture_model isa NonEquilMoistModel &&
559+
use_derivative(noneq_cloud_formation_flag)
560+
p_vapₛₗ(tps, ts) = TD.saturation_vapor_pressure(tps, ts, TD.Liquid())
561+
p_vapₛᵢ(tps, ts) = TD.saturation_vapor_pressure(tps, ts, TD.Ice())
562+
563+
function ∂p_vapₛₗ_∂T(tps, ts)
564+
T = TD.air_temperature(tps, ts)
565+
Rᵥ = TD.Parameters.R_v(tps)
566+
Lᵥ = TD.latent_heat_vapor(tps, ts)
567+
return p_vapₛₗ(tps, ts) * Lᵥ / (Rᵥ * T^2)
568+
end
569+
function ∂p_vapₛᵢ_∂T(tps, ts)
570+
T = TD.air_temperature(tps, ts)
571+
Rᵥ = TD.Parameters.R_v(tps)
572+
Lₛ = TD.latent_heat_sublim(tps, ts)
573+
return p_vapₛᵢ(tps, ts) * Lₛ / (Rᵥ * T^2)
574+
end
575575

576-
function ∂qₛₗ_∂T(tps, ts)
577-
T = TD.air_temperature(tps, ts)
578-
Rᵥ = TD.Parameters.R_v(tps)
579-
Lᵥ = TD.latent_heat_vapor(tps, ts)
580-
qᵥ_sat_liq = TD.q_vap_saturation_liquid(tps, ts)
581-
return qᵥ_sat_liq * (Lᵥ / (Rᵥ * T^2) - 1 / T)
582-
end
583-
function ∂qₛᵢ_∂T(tps, ts)
584-
T = TD.air_temperature(tps, ts)
585-
Rᵥ = TD.Parameters.R_v(tps)
586-
Lₛ = TD.latent_heat_sublim(tps, ts)
587-
qᵥ_sat_ice = TD.q_vap_saturation_ice(tps, ts)
588-
return qᵥ_sat_ice * (Lₛ / (Rᵥ * T^2) - 1 / T)
589-
end
576+
function ∂qₛₗ_∂T(tps, ts)
577+
T = TD.air_temperature(tps, ts)
578+
Rᵥ = TD.Parameters.R_v(tps)
579+
Lᵥ = TD.latent_heat_vapor(tps, ts)
580+
qᵥ_sat_liq = TD.q_vap_saturation_liquid(tps, ts)
581+
return qᵥ_sat_liq * (Lᵥ / (Rᵥ * T^2) - 1 / T)
582+
end
583+
function ∂qₛᵢ_∂T(tps, ts)
584+
T = TD.air_temperature(tps, ts)
585+
Rᵥ = TD.Parameters.R_v(tps)
586+
Lₛ = TD.latent_heat_sublim(tps, ts)
587+
qᵥ_sat_ice = TD.q_vap_saturation_ice(tps, ts)
588+
return qᵥ_sat_ice * (Lₛ / (Rᵥ * T^2) - 1 / T)
589+
end
590590

591-
function Γₗ(tps, ts)
592-
cₚ_air = TD.cp_m(tps, ts)
593-
Lᵥ = TD.latent_heat_vapor(tps, ts)
594-
return 1 + (Lᵥ / cₚ_air) * ∂qₛₗ_∂T(tps, ts)
595-
end
596-
function Γᵢ(tps, ts)
597-
cₚ_air = TD.cp_m(tps, ts)
598-
Lₛ = TD.latent_heat_sublim(tps, ts)
599-
return 1 + (Lₛ / cₚ_air) * ∂qₛᵢ_∂T(tps, ts)
600-
end
591+
function Γₗ(tps, ts)
592+
cₚ_air = TD.cp_m(tps, ts)
593+
Lᵥ = TD.latent_heat_vapor(tps, ts)
594+
return 1 + (Lᵥ / cₚ_air) * ∂qₛₗ_∂T(tps, ts)
595+
end
596+
function Γᵢ(tps, ts)
597+
cₚ_air = TD.cp_m(tps, ts)
598+
Lₛ = TD.latent_heat_sublim(tps, ts)
599+
return 1 + (Lₛ / cₚ_air) * ∂qₛᵢ_∂T(tps, ts)
600+
end
601601

602-
cmc = CAP.microphysics_cloud_params(params)
603-
τₗ = cmc.liquid.τ_relax
604-
τᵢ = cmc.ice.τ_relax
605-
function limit(q, dt, n::Int)
606-
return q / float(dt) / n
607-
end
608-
609-
function ∂ρqₗ_err_∂ρqᵪ(tps, ts, cmc, dt, deriv, limit_deriv)
610-
FT_inner = eltype(tps)
611-
q = TD.PhasePartition(tps, ts)
612-
ρ = TD.air_density(tps, ts)
613-
614-
S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cmc.liquid, tps, q, ρ, Tₐ(tps, ts))
615-
616-
if S > FT_inner(0)
617-
if S <= limit(TD.vapor_specific_humidity(q), dt, 2)
618-
if TD.vapor_specific_humidity(q) + TD.liquid_specific_humidity(q) > FT_inner(0)
619-
return deriv
602+
cmc = CAP.microphysics_cloud_params(params)
603+
τₗ = cmc.liquid.τ_relax
604+
τᵢ = cmc.ice.τ_relax
605+
function limit(q, dt, n::Int)
606+
return q / float(dt) / n
607+
end
608+
609+
function ∂ρqₗ_err_∂ρqᵪ(tps, ts, cmc, dt, deriv, limit_deriv)
610+
FT_inner = eltype(tps)
611+
q = TD.PhasePartition(tps, ts)
612+
ρ = TD.air_density(tps, ts)
613+
614+
S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cmc.liquid, tps, q, ρ, Tₐ(tps, ts))
615+
616+
if S > FT_inner(0)
617+
if S <= limit(TD.vapor_specific_humidity(q), dt, 2)
618+
if TD.vapor_specific_humidity(q) + TD.liquid_specific_humidity(q) > FT_inner(0)
619+
return deriv
620+
else
621+
return FT_inner(0)
622+
end
620623
else
621-
return FT_inner(0)
624+
return -limit_deriv
622625
end
623626
else
624-
return -limit_deriv
625-
end
626-
else
627-
if abs(S) <= limit(TD.liquid_specific_humidity(q), dt, 2)
628-
if TD.vapor_specific_humidity(q) + TD.liquid_specific_humidity(q) > FT_inner(0)
629-
return -deriv
627+
if abs(S) <= limit(TD.liquid_specific_humidity(q), dt, 2)
628+
if TD.vapor_specific_humidity(q) + TD.liquid_specific_humidity(q) > FT_inner(0)
629+
return -deriv
630+
else
631+
return FT_inner(0)
632+
end
630633
else
631-
return FT_inner(0)
634+
return -limit_deriv
632635
end
633-
else
634-
return -limit_deriv
635636
end
636637
end
637-
end
638638

639-
function ∂ρqᵢ_err_∂ρqᵪ(tps, ts, cmc, dt, deriv, limit_deriv)
640-
FT_inner = eltype(tps)
641-
q = TD.PhasePartition(tps, ts)
642-
ρ = TD.air_density(tps, ts)
639+
function ∂ρqᵢ_err_∂ρqᵪ(tps, ts, cmc, dt, deriv, limit_deriv)
640+
FT_inner = eltype(tps)
641+
q = TD.PhasePartition(tps, ts)
642+
ρ = TD.air_density(tps, ts)
643643

644-
S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cmc.ice, tps, q, ρ, Tₐ(tps, ts))
644+
S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cmc.ice, tps, q, ρ, Tₐ(tps, ts))
645645

646-
if S > FT_inner(0)
647-
if S <= limit(TD.vapor_specific_humidity(q), dt, 2)
648-
if TD.vapor_specific_humidity(q) + TD.ice_specific_humidity(q) > FT_inner(0)
649-
return deriv
646+
if S > FT_inner(0)
647+
if S <= limit(TD.vapor_specific_humidity(q), dt, 2)
648+
if TD.vapor_specific_humidity(q) + TD.ice_specific_humidity(q) > FT_inner(0)
649+
return deriv
650+
else
651+
return FT_inner(0)
652+
end
650653
else
651-
return FT_inner(0)
654+
return -limit_deriv
652655
end
653656
else
654-
return -limit_deriv
655-
end
656-
else
657-
if abs(S) <= limit(TD.ice_specific_humidity(q), dt, 2)
658-
if TD.vapor_specific_humidity(q) + TD.ice_specific_humidity(q) > FT_inner(0)
659-
return -deriv
657+
if abs(S) <= limit(TD.ice_specific_humidity(q), dt, 2)
658+
if TD.vapor_specific_humidity(q) + TD.ice_specific_humidity(q) > FT_inner(0)
659+
return -deriv
660+
else
661+
return FT_inner(0)
662+
end
660663
else
661-
return FT_inner(0)
664+
return -limit_deriv
662665
end
663-
else
664-
return -limit_deriv
665666
end
666667
end
667-
end
668668

669-
∂ᶜρqₗ_err_∂ᶜρqₗ = matrix[@name(c.ρq_liq), @name(c.ρq_liq)]
670-
∂ᶜρqᵢ_err_∂ᶜρqᵢ = matrix[@name(c.ρq_ice), @name(c.ρq_ice)]
669+
∂ᶜρqₗ_err_∂ᶜρqₗ = matrix[@name(c.ρq_liq), @name(c.ρq_liq)]
670+
∂ᶜρqᵢ_err_∂ᶜρqᵢ = matrix[@name(c.ρq_ice), @name(c.ρq_ice)]
671671

672-
∂ᶜρqₗ_err_∂ᶜρqₜ = matrix[@name(c.ρq_liq), @name(c.ρq_tot)]
673-
∂ᶜρqᵢ_err_∂ᶜρqₜ = matrix[@name(c.ρq_ice), @name(c.ρq_tot)]
672+
∂ᶜρqₗ_err_∂ᶜρqₜ = matrix[@name(c.ρq_liq), @name(c.ρq_tot)]
673+
∂ᶜρqᵢ_err_∂ᶜρqₜ = matrix[@name(c.ρq_ice), @name(c.ρq_tot)]
674674

675-
#if isdefined(Main, :Infiltrator)
676-
# Main.@infiltrate
677-
#end
678-
679-
#@. ∂ᶜρqₗ_err_∂ᶜρqₗ -=
680-
# DiagonalMatrixRow(1 / (τₗ * Γₗ(thermo_params, ᶜts)))
681-
@. ∂ᶜρqₗ_err_∂ᶜρqₗ +=
682-
DiagonalMatrixRow(
683-
∂ρqₗ_err_∂ρqᵪ(
684-
thermo_params, ᶜts, (cmc,), dt, (-1 / (τₗ * Γₗ(thermo_params, ᶜts))), (1/(2*float(dt))),
675+
#if isdefined(Main, :Infiltrator)
676+
# Main.@infiltrate
677+
#end
678+
679+
#@. ∂ᶜρqₗ_err_∂ᶜρqₗ -=
680+
# DiagonalMatrixRow(1 / (τₗ * Γₗ(thermo_params, ᶜts)))
681+
@. ∂ᶜρqₗ_err_∂ᶜρqₗ +=
682+
DiagonalMatrixRow(
683+
∂ρqₗ_err_∂ρqᵪ(
684+
thermo_params, ᶜts, (cmc,), dt, (-1 / (τₗ * Γₗ(thermo_params, ᶜts))), (1/(2*float(dt))),
685+
)
685686
)
686-
)
687-
688-
#@. ∂ᶜρqᵢ_err_∂ᶜρqᵢ -=
689-
# DiagonalMatrixRow(1 / (τᵢ * Γᵢ(thermo_params, ᶜts)))
687+
688+
#@. ∂ᶜρqᵢ_err_∂ᶜρqᵢ -=
689+
# DiagonalMatrixRow(1 / (τᵢ * Γᵢ(thermo_params, ᶜts)))
690690

691-
@. ∂ᶜρqᵢ_err_∂ᶜρqᵢ +=
692-
DiagonalMatrixRow(
693-
∂ρqᵢ_err_∂ρqᵪ(
694-
thermo_params, ᶜts, (cmc,), dt, (-1 / (τᵢ * Γᵢ(thermo_params, ᶜts))), (1/(2*float(dt))),
691+
@. ∂ᶜρqᵢ_err_∂ᶜρqᵢ +=
692+
DiagonalMatrixRow(
693+
∂ρqᵢ_err_∂ρqᵪ(
694+
thermo_params, ᶜts, (cmc,), dt, (-1 / (τᵢ * Γᵢ(thermo_params, ᶜts))), (1/(2*float(dt))),
695+
)
695696
)
696-
)
697697

698-
ᶜp = @. lazy(TD.air_pressure(thermo_params, ᶜts))
699-
ᶜ∂T_∂p = @. lazy(1 / (ᶜρ * TD.gas_constant_air(thermo_params, ᶜts)))
698+
ᶜp = @. lazy(TD.air_pressure(thermo_params, ᶜts))
699+
ᶜ∂T_∂p = @. lazy(1 / (ᶜρ * TD.gas_constant_air(thermo_params, ᶜts)))
700700

701-
# qₛₗ = p_vapₛₗ / p, qₛᵢ = p_vapₛᵢ / p
702-
ᶜ∂qₛₗ_∂p = @. lazy(
703-
-p_vapₛₗ(thermo_params, ᶜts) / ᶜp^2 +
704-
∂p_vapₛₗ_∂T(thermo_params, ᶜts) * ᶜ∂T_∂p / ᶜp,
705-
)
706-
ᶜ∂qₛᵢ_∂p = @. lazy(
707-
-p_vapₛᵢ(thermo_params, ᶜts) / ᶜp^2 +
708-
∂p_vapₛᵢ_∂T(thermo_params, ᶜts) * ᶜ∂T_∂p / ᶜp,
709-
)
701+
# qₛₗ = p_vapₛₗ / p, qₛᵢ = p_vapₛᵢ / p
702+
ᶜ∂qₛₗ_∂p = @. lazy(
703+
-p_vapₛₗ(thermo_params, ᶜts) / ᶜp^2 +
704+
∂p_vapₛₗ_∂T(thermo_params, ᶜts) * ᶜ∂T_∂p / ᶜp,
705+
)
706+
ᶜ∂qₛᵢ_∂p = @. lazy(
707+
-p_vapₛᵢ(thermo_params, ᶜts) / ᶜp^2 +
708+
∂p_vapₛᵢ_∂T(thermo_params, ᶜts) * ᶜ∂T_∂p / ᶜp,
709+
)
710710

711-
ᶜ∂p_∂ρqₜ = @. lazy(
712-
ᶜkappa_m * ∂e_int_∂q_tot +
713-
ᶜ∂kappa_m∂q_tot * (
714-
cp_d * T_0 + ᶜspecific.e_tot - ᶜK - ᶜΦ +
715-
∂e_int_∂q_tot * ᶜspecific.q_tot
716-
),
717-
)
711+
ᶜ∂p_∂ρqₜ = @. lazy(
712+
ᶜkappa_m * ∂e_int_∂q_tot +
713+
ᶜ∂kappa_m∂q_tot * (
714+
cp_d * T_0 + ᶜspecific.e_tot - ᶜK - ᶜΦ +
715+
∂e_int_∂q_tot * ᶜspecific.q_tot
716+
),
717+
)
718718

719-
#@. ∂ᶜρqₗ_err_∂ᶜρqₜ = DiagonalMatrixRow(
720-
# (1 - ᶜρ * ᶜ∂qₛₗ_∂p * ᶜ∂p_∂ρqₜ) / (τₗ * Γₗ(thermo_params, ᶜts)),
721-
#)
722-
@. ∂ᶜρqₗ_err_∂ᶜρqₜ = DiagonalMatrixRow(
723-
∂ρqₗ_err_∂ρqᵪ(
724-
thermo_params, ᶜts, (cmc,), dt, ((1 - ᶜρ * ᶜ∂qₛₗ_∂p * ᶜ∂p_∂ρqₜ) / (τₗ * Γₗ(thermo_params, ᶜts))), FT(0)
719+
#@. ∂ᶜρqₗ_err_∂ᶜρqₜ = DiagonalMatrixRow(
720+
# (1 - ᶜρ * ᶜ∂qₛₗ_∂p * ᶜ∂p_∂ρqₜ) / (τₗ * Γₗ(thermo_params, ᶜts)),
721+
#)
722+
@. ∂ᶜρqₗ_err_∂ᶜρqₜ = DiagonalMatrixRow(
723+
∂ρqₗ_err_∂ρqᵪ(
724+
thermo_params, ᶜts, (cmc,), dt, ((1 - ᶜρ * ᶜ∂qₛₗ_∂p * ᶜ∂p_∂ρqₜ) / (τₗ * Γₗ(thermo_params, ᶜts))), FT(0)
725+
)
725726
)
726-
)
727727

728-
#@. ∂ᶜρqᵢ_err_∂ᶜρqₜ = DiagonalMatrixRow(
729-
# (1 - ᶜρ * ᶜ∂qₛᵢ_∂p * ᶜ∂p_∂ρqₜ) / (τᵢ * Γᵢ(thermo_params, ᶜts)),
730-
#)
731-
@. ∂ᶜρqᵢ_err_∂ᶜρqₜ = DiagonalMatrixRow(
732-
∂ρqᵢ_err_∂ρqᵪ(
733-
thermo_params, ᶜts, (cmc,), dt, ((1 - ᶜρ * ᶜ∂qₛᵢ_∂p * ᶜ∂p_∂ρqₜ) / (τᵢ * Γᵢ(thermo_params, ᶜts))), FT(0)
728+
#@. ∂ᶜρqᵢ_err_∂ᶜρqₜ = DiagonalMatrixRow(
729+
# (1 - ᶜρ * ᶜ∂qₛᵢ_∂p * ᶜ∂p_∂ρqₜ) / (τᵢ * Γᵢ(thermo_params, ᶜts)),
730+
#)
731+
@. ∂ᶜρqᵢ_err_∂ᶜρqₜ = DiagonalMatrixRow(
732+
∂ρqᵢ_err_∂ρqᵪ(
733+
thermo_params, ᶜts, (cmc,), dt, ((1 - ᶜρ * ᶜ∂qₛᵢ_∂p * ᶜ∂p_∂ρqₜ) / (τᵢ * Γᵢ(thermo_params, ᶜts))), FT(0)
734+
)
734735
)
735-
)
736+
end
736737
end
737738

738739
if use_derivative(diffusion_flag)
@@ -1266,4 +1267,3 @@ end
12661267

12671268
invert_jacobian!(::ManualSparseJacobian, cache, ΔY, R) =
12681269
LinearAlgebra.ldiv!(ΔY, cache.matrix, R)
1269-
end

0 commit comments

Comments
 (0)