From 9ccc6bc7055c548e3084b40e3de94f3ebe918f8a Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Tue, 13 Sep 2022 14:28:56 -0700 Subject: [PATCH] Add optional initial temperature guess WIP Revert tolerance change Exercise T_guess in test suite --- src/config_numerical_method.jl | 64 +++++++++++++++++++++++++++++----- src/relations.jl | 40 +++++++++++++++++---- src/states.jl | 36 +++++++++++++++---- test/relations.jl | 39 +++++++++++++++++++++ 4 files changed, 158 insertions(+), 21 deletions(-) diff --git a/src/config_numerical_method.jl b/src/config_numerical_method.jl index b2ab5ba0..3f616e89 100644 --- a/src/config_numerical_method.jl +++ b/src/config_numerical_method.jl @@ -21,6 +21,17 @@ function print_numerical_method( end end +function print_T_guess( + ::Type{sat_adjust_method}, + T_guess::Union{FT, Nothing}, +) where {sat_adjust_method, FT <: Real} + if sat_adjust_method <: RS.NewtonsMethod + KA.@print(", T_guess=", T_guess) + elseif sat_adjust_method <: RS.NewtonsMethodAD + KA.@print(", T_guess=", T_guess) + end +end + ##### ##### Thermodynamic variable inputs: ρ, e_int, q_tot ##### @@ -31,10 +42,14 @@ function sa_numerical_method( e_int::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.NewtonsMethod, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) - T_init = + T_init = if T_guess isa Nothing max(T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor + else + T_guess + end return RS.NewtonsMethod( T_init, T_ -> ∂e_int_∂T(param_set, heavisided(T_), e_int, ρ, q_tot, phase_type), @@ -48,10 +63,14 @@ function sa_numerical_method( e_int::FT, q_tot::FT, ::Type{phase_type}, + T_guess::FT, ) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) - T_init = + T_init = if T_guess isa Nothing max(T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor + else + T_guess + end return RS.NewtonsMethodAD(T_init) end @@ -62,6 +81,7 @@ function sa_numerical_method( e_int::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice @@ -78,6 +98,7 @@ function sa_numerical_method( e_int::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice @@ -98,9 +119,14 @@ function sa_numerical_method_ρpq( p::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil} q_pt = PhasePartition(q_tot) - T_init = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) + T_init = if T_guess isa Nothing + air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) + else + T_guess + end return RS.NewtonsMethodAD(T_init) end @@ -111,6 +137,7 @@ function sa_numerical_method_ρpq( p::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil} q_pt = PhasePartition(q_tot) T_1 = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) - 5 @@ -129,10 +156,14 @@ function sa_numerical_method_peq( e_int::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) - T_init = + T_init = if T_guess isa Nothing max(T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor + else + T_guess + end return RS.NewtonsMethodAD(T_init) end @@ -143,6 +174,7 @@ function sa_numerical_method_peq( e_int::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice @@ -163,12 +195,17 @@ function sa_numerical_method_phq( h::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) - T_init = max( - T_min, - air_temperature_from_enthalpy(param_set, h, PhasePartition(q_tot)), - ) # Assume all vapor + T_init = if T_guess isa Nothing # Assume all vapor + max( + T_min, + air_temperature_from_enthalpy(param_set, h, PhasePartition(q_tot)), + ) + else + T_guess + end return RS.NewtonsMethodAD(T_init) end @@ -179,6 +216,7 @@ function sa_numerical_method_phq( h::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice @@ -198,6 +236,7 @@ function sa_numerical_method_phq( h::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice @@ -221,6 +260,7 @@ function sa_numerical_method_pθq( θ_liq_ice::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil} _T_min::FT = TP.T_min(param_set) _T_max::FT = TP.T_max(param_set) @@ -238,6 +278,7 @@ function sa_numerical_method_pθq( θ_liq_ice::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil} _T_min::FT = TP.T_min(param_set) air_temp(q) = air_temperature_given_pθq(param_set, p, θ_liq_ice, q) @@ -254,9 +295,14 @@ function sa_numerical_method_pθq( θ_liq_ice::FT, q_tot::FT, ::Type{phase_type}, + T_guess::Union{FT, Nothing}, ) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil} T_min::FT = TP.T_min(param_set) air_temp(q) = air_temperature_given_pθq(param_set, p, θ_liq_ice, q) - T_init = max(T_min, air_temp(PhasePartition(q_tot))) # Assume all vapor + T_init = if T_guess isa Nothing + max(T_min, air_temp(PhasePartition(q_tot))) # Assume all vapor + else + T_guess + end return RS.NewtonsMethodAD(T_init) end diff --git a/src/relations.jl b/src/relations.jl index 5846b202..910f65ea 100644 --- a/src/relations.jl +++ b/src/relations.jl @@ -1450,7 +1450,8 @@ end q_tot, phase_type, maxiter, - temperature_tol + temperature_tol, + T_guess, ) Compute the temperature that is consistent with @@ -1464,6 +1465,7 @@ Compute the temperature that is consistent with - `phase_type` a thermodynamic state type - `maxiter` maximum iterations for non-linear equation solve - `temperature_tol` temperature tolerance + - `T_guess` initial temperature guess by finding the root of @@ -1482,6 +1484,7 @@ function saturation_adjustment( ::Type{phase_type}, maxiter::Int, temperature_tol::FT, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil} _T_min::FT = TP.T_min(param_set) cv_d::FT = TP.cv_d(param_set) @@ -1511,6 +1514,7 @@ function saturation_adjustment( e_int, q_tot, phase_type, + T_guess, ), solution_type(), tol, @@ -1522,6 +1526,7 @@ function saturation_adjustment( KA.@print("-----------------------------------------\n") KA.@print("maxiter reached in saturation_adjustment:\n") print_numerical_method(sat_adjust_method) + print_T_guess(sat_adjust_method, T_guess) KA.@print(", e_int=", e_int) KA.@print(", ρ=", ρ) KA.@print(", q_tot=", q_tot) @@ -1545,7 +1550,8 @@ end q_tot, phase_type, maxiter, - temperature_tol + temperature_tol, + T_guess, ) Compute the temperature that is consistent with @@ -1559,6 +1565,7 @@ Compute the temperature that is consistent with - `phase_type` a thermodynamic state type - `maxiter` maximum iterations for non-linear equation solve - `temperature_tol` temperature tolerance + - `T_guess` initial temperature guess by finding the root of @@ -1579,6 +1586,7 @@ function saturation_adjustment_given_peq( ::Type{phase_type}, maxiter::Int, temperature_tol::FT, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil} _T_min::FT = TP.T_min(param_set) cv_d = FT(TP.cv_d(param_set)) @@ -1611,6 +1619,7 @@ function saturation_adjustment_given_peq( e_int, q_tot, phase_type, + T_guess, ), RS.CompactSolution(), tol, @@ -1621,6 +1630,7 @@ function saturation_adjustment_given_peq( KA.@print("-----------------------------------------\n") KA.@print("maxiter reached in saturation_adjustment_peq:\n") print_numerical_method(sat_adjust_method) + print_T_guess(sat_adjust_method, T_guess) KA.@print(", e_int=", e_int) KA.@print(", p=", p) KA.@print(", q_tot=", q_tot) @@ -1646,6 +1656,7 @@ end phase_type, maxiter, temperature_tol + T_guess, ) Compute the temperature that is consistent with @@ -1659,6 +1670,7 @@ Compute the temperature that is consistent with - `phase_type` a thermodynamic state type - `maxiter` maximum iterations for non-linear equation solve - `temperature_tol` temperature tolerance + - `T_guess` initial temperature guess by finding the root of @@ -1679,6 +1691,7 @@ function saturation_adjustment_given_phq( ::Type{phase_type}, maxiter::Int, temperature_tol::FT, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil} _T_min::FT = TP.T_min(param_set) cp_d::FT = TP.cp_d(param_set) @@ -1719,6 +1732,7 @@ function saturation_adjustment_given_phq( h, q_tot, phase_type, + T_guess, ), RS.CompactSolution(), tol, @@ -1729,6 +1743,7 @@ function saturation_adjustment_given_phq( KA.@print("-----------------------------------------\n") KA.@print("maxiter reached in saturation_adjustment_phq:\n") print_numerical_method(sat_adjust_method) + print_T_guess(sat_adjust_method, T_guess) KA.@print(", h=", h) KA.@print(", p=", p) KA.@print(", q_tot=", q_tot) @@ -1752,7 +1767,8 @@ end q_tot, phase_type, maxiter, - temperature_tol + temperature_tol, + T_guess, ) Compute the temperature that is consistent with - `sat_adjust_method` the numerical method to use. @@ -1764,6 +1780,7 @@ Compute the temperature that is consistent with - `phase_type` a thermodynamic state type - `maxiter` maximum iterations for non-linear equation solve - `temperature_tol` temperature tolerance + - `T_guess` initial temperature guess by finding the root of ``` @@ -1786,6 +1803,7 @@ function saturation_adjustment_ρpq( ::Type{phase_type}, maxiter::Int, temperature_tol::FT = sqrt(eps(FT)), + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil} tol = RS.SolutionTolerance(temperature_tol) # Use `oftype` to preserve diagonalized type signatures: @@ -1810,6 +1828,7 @@ function saturation_adjustment_ρpq( p, q_tot, phase_type, + T_guess, ), RS.CompactSolution(), tol, @@ -1820,6 +1839,7 @@ function saturation_adjustment_ρpq( KA.@print("-----------------------------------------\n") KA.@print("maxiter reached in saturation_adjustment_ρpq:\n") print_numerical_method(sat_adjust_method) + print_T_guess(sat_adjust_method, T_guess) KA.@print(", ρ=", ρ) KA.@print(", p=", p) KA.@print(", q_tot=", q_tot) @@ -1867,7 +1887,8 @@ end q_tot, phase_type, maxiter, - tol + tol, + T_guess, ) Compute the temperature `T` that is consistent with @@ -1881,6 +1902,7 @@ Compute the temperature `T` that is consistent with - `tol` absolute tolerance for saturation adjustment iterations. Can be one of: - `SolutionTolerance()` to stop when `|x_2 - x_1| < tol` - `ResidualTolerance()` to stop when `|f(x)| < tol` + - `T_guess` initial temperature guess by finding the root of @@ -1896,6 +1918,7 @@ function saturation_adjustment_given_ρθq( ::Type{phase_type}, maxiter::Int, tol::RS.AbstractTolerance, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, phase_type <: PhaseEquil} _T_min::FT = TP.T_min(param_set) air_temp(q) = air_temperature_given_ρθq(param_set, ρ, θ_liq_ice, q) @@ -1949,7 +1972,8 @@ end q_tot, phase_type, maxiter, - temperature_tol + temperature_tol, + T_guess ) Compute the temperature `T` that is consistent with @@ -1961,7 +1985,8 @@ Compute the temperature `T` that is consistent with - `phase_type` a thermodynamic state type - `temperature_tol` temperature tolerance - `maxiter` maximum iterations for non-linear equation solve -- `sat_adjust_method` the numerical method to use. + - `sat_adjust_method` the numerical method to use. + - `T_guess` initial temperature guess by finding the root of @@ -1978,6 +2003,7 @@ function saturation_adjustment_given_pθq( ::Type{phase_type}, maxiter::Int, temperature_tol::FT, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil} tol = RS.ResidualTolerance(temperature_tol) T_min::FT = TP.T_min(param_set) @@ -2023,6 +2049,7 @@ function saturation_adjustment_given_pθq( θ_liq_ice, q_tot, phase_type, + T_guess, ), RS.CompactSolution(), tol, @@ -2033,6 +2060,7 @@ function saturation_adjustment_given_pθq( KA.@print("-----------------------------------------\n") KA.@print("maxiter reached in saturation_adjustment_given_pθq:\n") print_numerical_method(sat_adjust_method) + print_T_guess(sat_adjust_method, T_guess) KA.@print(", p=", p) KA.@print(", θ_liq_ice=", θ_liq_ice) KA.@print(", q_tot=", q_tot) diff --git a/src/states.jl b/src/states.jl index c1782fbd..badf74ed 100644 --- a/src/states.jl +++ b/src/states.jl @@ -259,7 +259,7 @@ struct PhaseEquil{FT} <: AbstractPhaseEquil{FT} end """ - PhaseEquil_ρeq + PhaseEquil_ρeq(param_set, ρ, e_int, q_tot[, maxiter, temperature_tol, sat_adjust_method, T_guess]) Moist thermodynamic phase, given - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details @@ -271,6 +271,7 @@ and, optionally - `temperature_tol` temperature tolerance for saturation adjustment - `sat_adjust_method` the numerical method to use. See the [`Thermodynamics`](@ref) for options. + - `T_guess` initial guess for temperature in saturation adjustment """ function PhaseEquil_ρeq( param_set::APS, @@ -280,6 +281,7 @@ function PhaseEquil_ρeq( maxiter::IT = nothing, temperature_tol::FTT = nothing, ::Type{sat_adjust_method} = RS.NewtonsMethod, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, IT <: ITERTYPE, FTT <: TOLTYPE(FT)} maxiter === nothing && (maxiter = 8) temperature_tol === nothing && (temperature_tol = FT(1e-1)) @@ -294,6 +296,7 @@ function PhaseEquil_ρeq( phase_type, maxiter, temperature_tol, + T_guess, ) q_pt = PhasePartition_equil(param_set, T, ρ, q_tot_safe, phase_type) p = air_pressure(param_set, T, ρ, q_pt) @@ -325,7 +328,7 @@ function PhaseEquil_dev_only( end """ - PhaseEquil_ρθq(param_set, ρ, θ_liq_ice, q_tot) + PhaseEquil_ρθq(param_set, ρ, θ_liq_ice, q_tot[, maxiter, temperature_tol, sat_adjust_method, T_guess]) Constructs a [`PhaseEquil`](@ref) thermodynamic state from: @@ -335,6 +338,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from: - `q_tot` total specific humidity - `temperature_tol` temperature tolerance for saturation adjustment - `maxiter` maximum iterations for saturation adjustment + - `T_guess` initial guess for temperature in saturation adjustment """ function PhaseEquil_ρθq( param_set::APS, @@ -343,6 +347,7 @@ function PhaseEquil_ρθq( q_tot::FT, maxiter::IT = nothing, temperature_tol::FTT = nothing, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, IT <: ITERTYPE, FTT <: TOLTYPE(FT)} maxiter === nothing && (maxiter = 36) temperature_tol === nothing && (temperature_tol = FT(1e-1)) @@ -356,6 +361,7 @@ function PhaseEquil_ρθq( phase_type, maxiter, tol, + T_guess, ) q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type) p = air_pressure(param_set, T, ρ, q_pt) @@ -411,7 +417,7 @@ function PhaseEquil_pTq( end """ - PhaseEquil_peq(param_set, p, e_int, q_tot) + PhaseEquil_peq(param_set, p, e_int, q_tot[, maxiter, temperature_tol, sat_adjust_method, T_guess]) Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. @@ -419,6 +425,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. - `p` pressure - `e_int` temperature - `q_tot` total specific humidity + - `T_guess` initial guess for temperature in saturation adjustment """ function PhaseEquil_peq( param_set::APS, @@ -428,6 +435,7 @@ function PhaseEquil_peq( maxiter::IT = nothing, temperature_tol::FTT = nothing, ::Type{sat_adjust_method} = RS.SecantMethod, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, IT <: ITERTYPE, FTT <: TOLTYPE(FT)} maxiter === nothing && (maxiter = 40) temperature_tol === nothing && (temperature_tol = FT(1e-2)) @@ -442,6 +450,7 @@ function PhaseEquil_peq( phase_type, maxiter, temperature_tol, + T_guess, ) q_pt = PhasePartition_equil_given_p(param_set, T, p, q_tot_safe, phase_type) ρ = air_density(param_set, T, p, q_pt) @@ -450,7 +459,7 @@ end """ - PhaseEquil_phq(param_set, p, h, q_tot) + PhaseEquil_phq(param_set, p, h, q_tot[, maxiter, temperature_tol, sat_adjust_method, T_guess]) Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. @@ -458,6 +467,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. - `p` pressure - `h` specific enthalpy - `q_tot` total specific humidity + - `T_guess` initial guess for temperature in saturation adjustment """ function PhaseEquil_phq( param_set::APS, @@ -467,6 +477,7 @@ function PhaseEquil_phq( maxiter::IT = nothing, temperature_tol::FTT = nothing, ::Type{sat_adjust_method} = RS.SecantMethod, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method, IT <: ITERTYPE, FTT <: TOLTYPE(FT)} maxiter === nothing && (maxiter = 40) temperature_tol === nothing && (temperature_tol = FT(1e-2)) @@ -481,6 +492,7 @@ function PhaseEquil_phq( phase_type, maxiter, temperature_tol, + T_guess, ) q_pt = PhasePartition_equil_given_p(param_set, T, p, q_tot_safe, phase_type) ρ = air_density(param_set, T, p, q_pt) @@ -489,7 +501,7 @@ function PhaseEquil_phq( end """ - PhaseEquil_ρpq(param_set, ρ, p, q_tot, perform_sat_adjust=true) + PhaseEquil_ρpq(param_set, ρ, p, q_tot[, perform_sat_adjust=true, maxiter, sat_adjust_method, T_guess]) Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. @@ -501,6 +513,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. - `maxiter` maximum number of iterations to perform in saturation adjustment - `sat_adjust_method` the numerical method to use. See the [`Thermodynamics`](@ref) for options. + - `T_guess` initial guess for temperature in saturation adjustment + +TODO: change input argument order: perform_sat_adjust is + unique to this constructor, so it should be last. + (breaking change) """ function PhaseEquil_ρpq( param_set::APS, @@ -509,7 +526,9 @@ function PhaseEquil_ρpq( q_tot::FT, perform_sat_adjust = false, maxiter::Int = 5, + temperature_tol::FT = sqrt(eps(FT)), ::Type{sat_adjust_method} = RS.NewtonsMethodAD, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, sat_adjust_method} phase_type = PhaseEquil{FT} @@ -522,6 +541,8 @@ function PhaseEquil_ρpq( q_tot, phase_type, maxiter, + temperature_tol, + T_guess, ) q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type) e_int = internal_energy(param_set, T, q_pt) @@ -535,7 +556,7 @@ end """ - PhaseEquil_pθq(param_set, θ_liq_ice, q_tot[, maxiter, temperature_tol, sat_adjust_method]) + PhaseEquil_pθq(param_set, θ_liq_ice, q_tot[, maxiter, temperature_tol, sat_adjust_method, T_guess]) Constructs a [`PhaseEquil`](@ref) thermodynamic state from: @@ -546,6 +567,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from: - `temperature_tol` temperature tolerance for saturation adjustment - `maxiter` maximum iterations for saturation adjustment - `sat_adjust_method` the numerical method to use. + - `T_guess` initial guess for temperature in saturation adjustment """ function PhaseEquil_pθq( param_set::APS, @@ -555,6 +577,7 @@ function PhaseEquil_pθq( maxiter::IT = nothing, temperature_tol::FTT = nothing, ::Type{sat_adjust_method} = RS.SecantMethod, + T_guess::Union{FT, Nothing} = nothing, ) where {FT <: Real, IT <: ITERTYPE, FTT <: TOLTYPE(FT), sat_adjust_method} maxiter === nothing && (maxiter = 50) temperature_tol === nothing && (temperature_tol = FT(1e-3)) @@ -569,6 +592,7 @@ function PhaseEquil_pθq( phase_type, maxiter, temperature_tol, + T_guess, ) q_pt = PhasePartition_equil_given_p(param_set, T, p, q_tot_safe, phase_type) ρ = air_density(param_set, T, p, q_pt) diff --git a/test/relations.jl b/test/relations.jl index 49ae5697..a6e31a42 100644 --- a/test/relations.jl +++ b/test/relations.jl @@ -2,6 +2,9 @@ using Test using Thermodynamics using Thermodynamics.TemperatureProfiles using Thermodynamics.TestedProfiles +import Thermodynamics +Thermodynamics.error_on_non_convergence() = true +Thermodynamics.print_warning() = true const TD = Thermodynamics const TP = TD.Parameters @@ -1670,6 +1673,42 @@ end @test all(phase_type .== (nt.phase_type for nt in profiles)) end +@testset "Thermodynamics - test T_guess" begin + ArrayType = Array{Float64} + FT = eltype(ArrayType) + param_set = parameter_set(FT) + profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType) + @unpack p, ρ, e_int, h, θ_liq_ice, q_tot, T, phase_type = profiles + T_guess = T .+ (FT(0.2) .* randn(FT, length(T))) + args = (q_tot, 40, FT(1e-1)) + ts = + PhaseEquil_ρeq.(param_set, ρ, e_int, args..., RS.NewtonsMethod, T_guess) + ts = PhaseEquil_ρθq.(param_set, ρ, θ_liq_ice, args..., T_guess) + ts = PhaseEquil_peq.(param_set, p, e_int, args..., RS.SecantMethod, T_guess) + ts = PhaseEquil_phq.(param_set, p, h, args..., RS.SecantMethod, T_guess) + ts = + PhaseEquil_ρpq.( + param_set, + ρ, + p, + q_tot, + true, + 40, + FT(1e-1), + RS.NewtonsMethodAD, + T_guess, + ) + ts = + PhaseEquil_pθq.( + param_set, + p, + θ_liq_ice, + args..., + RS.SecantMethod, + T_guess, + ) +end + rm(joinpath(@__DIR__, "logfilepath_Float32.toml"); force = true) rm(joinpath(@__DIR__, "logfilepath_Float64.toml"); force = true)