Skip to content

Commit

Permalink
PhaseEquil outer constructor -> PhaseEquil_ρeq
Browse files Browse the repository at this point in the history
  • Loading branch information
charleskawczynski committed Sep 13, 2021
1 parent 505890d commit b78b479
Show file tree
Hide file tree
Showing 8 changed files with 117 additions and 78 deletions.
4 changes: 2 additions & 2 deletions docs/ThreeDimensionalInput.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,11 @@ TS_no_err = Array{ThermodynamicState}(undef, prod(dims));
q_tot_all[p] = q_tot[k]

Thermodynamics.error_on_non_convergence() = false
TS_no_err[p] = PhaseEquil(param_set, e_int[j], ρ[i], q_tot[k])
TS_no_err[p] = PhaseEquil_ρeq(param_set, ρ[i], e_int[j], q_tot[k])
Thermodynamics.error_on_non_convergence() = true
# @show p/prod(linear_indices.indices)*100
try
TS[p] = PhaseEquil(param_set, e_int[j], ρ[i], q_tot[k])
TS[p] = PhaseEquil_ρeq(param_set, ρ[i], e_int[j], q_tot[k])
catch
TS[p] = nothing
end
Expand Down
4 changes: 2 additions & 2 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ Users are encouraged to first establish a thermodynamic state with one of our
a moist thermodynamic state using

```julia
ts = PhaseEquil(param_set, e_int, ρ, q_tot);
ts = PhaseEquil_ρeq(param_set, ρ, e_int, q_tot);
```

here, `ρ` is the density of the moist air, and the internal energy `e_int =
Expand Down Expand Up @@ -102,7 +102,7 @@ do timestep # timestepping loop
e_int = e_tot - 0.5 * (u^2 + v^2 + w^2) - geopotential

# compute temperature, pressure and condensate specific humidities,
ts = PhaseEquil(param_set, e_int, ρ, q_tot);
ts = PhaseEquil_ρeq(param_set, ρ, e_int, q_tot);
T = air_temperature(ts);
q = PhasePartition(ts);
p = air_pressure(ts);
Expand Down
47 changes: 24 additions & 23 deletions src/relations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1645,10 +1645,10 @@ function saturation_adjustment_given_ρθq(
_T_min::FT = T_min(param_set)
T_1 = max(
_T_min,
air_temperature_given_θρq(
air_temperature_given_ρθq(
param_set,
θ_liq_ice,
ρ,
θ_liq_ice,
PhasePartition(q_tot),
),
) # Assume all vapor
Expand All @@ -1657,10 +1657,10 @@ function saturation_adjustment_given_ρθq(
if unsaturated && T_1 > _T_min
return T_1
else
T_2 = air_temperature_given_θρq(
T_2 = air_temperature_given_ρθq(
param_set,
θ_liq_ice,
ρ,
θ_liq_ice,
PhasePartition(q_tot, FT(0), q_tot),
) # Assume all ice
T_2 = bound_upper_temperature(T_1, T_2)
Expand Down Expand Up @@ -1741,10 +1741,10 @@ function saturation_adjustment_given_pθq(
tol::AbstractTolerance,
) where {FT <: Real}
_T_min::FT = T_min(param_set)
T_1 = air_temperature_given_θpq(
T_1 = air_temperature_given_pθq(
param_set,
θ_liq_ice,
p,
θ_liq_ice,
PhasePartition(q_tot),
) # Assume all vapor
ρ = air_density(param_set, T_1, p, PhasePartition(q_tot))
Expand All @@ -1753,10 +1753,10 @@ function saturation_adjustment_given_pθq(
if unsaturated && T_1 > _T_min
return T_1
else
T_2 = air_temperature_given_θpq(
T_2 = air_temperature_given_pθq(
param_set,
θ_liq_ice,
p,
θ_liq_ice,
PhasePartition(q_tot, FT(0), q_tot),
) # Assume all ice
T_2 = bound_upper_temperature(T_1, T_2)
Expand Down Expand Up @@ -2000,7 +2000,7 @@ function temperature_and_humidity_given_TᵥρRH(
end

"""
air_temperature_given_θρq(param_set, θ_liq_ice, ρ, q::PhasePartition)
air_temperature_given_ρθq(param_set, ρ, θ_liq_ice, q::PhasePartition)
The temperature given
- `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
Expand All @@ -2009,10 +2009,10 @@ The temperature given
and, optionally,
- `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air.
"""
function air_temperature_given_θρq(
function air_temperature_given_ρθq(
param_set::APS,
θ_liq_ice::FT,
ρ::FT,
θ_liq_ice::FT,
q::PhasePartition{FT} = q_pt_0(FT),
) where {FT <: Real}

Expand All @@ -2028,7 +2028,7 @@ function air_temperature_given_θρq(
end

"""
air_temperature_given_θρq_nonlinear(param_set, θ_liq_ice, ρ, q::PhasePartition)
air_temperature_given_ρθq_nonlinear(param_set, ρ, θ_liq_ice, q::PhasePartition)
Computes temperature `T` given
Expand All @@ -2043,15 +2043,15 @@ and, optionally,
- `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air,
by finding the root of
`T - air_temperature_given_θpq(param_set,
θ_liq_ice,
`T - air_temperature_given_pθq(param_set,
air_pressure(param_set, T, ρ, q),
θ_liq_ice,
q) = 0`
"""
function air_temperature_given_θρq_nonlinear(
function air_temperature_given_ρθq_nonlinear(
param_set::APS,
θ_liq_ice::FT,
ρ::FT,
θ_liq_ice::FT,
maxiter::Int,
tol::AbstractTolerance,
q::PhasePartition{FT} = q_pt_0(FT),
Expand All @@ -2060,10 +2060,10 @@ function air_temperature_given_θρq_nonlinear(
_T_max::FT = T_max(param_set)
sol = find_zero(
T ->
T - air_temperature_given_θpq(
T - air_temperature_given_pθq(
param_set,
θ_liq_ice,
air_pressure(param_set, heavisided(T), ρ, q),
θ_liq_ice,
q,
),
SecantMethod(_T_min, _T_max),
Expand All @@ -2074,7 +2074,7 @@ function air_temperature_given_θρq_nonlinear(
if !sol.converged
if print_warning()
@print("-----------------------------------------\n")
@print("maxiter reached in air_temperature_given_θρq_nonlinear:\n")
@print("maxiter reached in air_temperature_given_ρθq_nonlinear:\n")
@print(" Method=SecantMethod")
@print(", θ_liq_ice=", θ_liq_ice)
@print(", ρ=", ρ)
Expand All @@ -2093,10 +2093,11 @@ function air_temperature_given_θρq_nonlinear(
end

"""
air_temperature_given_θpq(
air_temperature_given_pθq(
param_set,
p,
θ_liq_ice,
p[, q::PhasePartition]
[q::PhasePartition]
)
The air temperature where
Expand All @@ -2107,10 +2108,10 @@ The air temperature where
and, optionally,
- `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air.
"""
function air_temperature_given_θpq(
function air_temperature_given_pθq(
param_set::APS,
θ_liq_ice::FT,
p::FT,
θ_liq_ice::FT,
q::PhasePartition{FT} = q_pt_0(FT),
) where {FT <: Real}
return θ_liq_ice * exner_given_pressure(param_set, p, q) +
Expand Down
39 changes: 20 additions & 19 deletions src/states.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ export ThermodynamicState,
PhaseDry_pθ,
PhaseDry_ρp,
PhaseEquil,
PhaseEquil_ρeq,
PhaseEquil_ρTq,
PhaseEquil_pTq,
PhaseEquil_peq,
Expand Down Expand Up @@ -139,7 +140,7 @@ Constructs a [`PhaseDry`](@ref) thermodynamic state from:
- `θ_dry` dry potential temperature
"""
function PhaseDry_ρθ(param_set::APS, ρ::FT, θ_dry::FT) where {FT <: Real}
T = air_temperature_given_θρq(param_set, θ_dry, ρ)
T = air_temperature_given_ρθq(param_set, ρ, θ_dry)
e_int = internal_energy(param_set, T)
return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)
end
Expand Down Expand Up @@ -201,7 +202,7 @@ may be needed).
# Constructors
PhaseEquil(param_set, e_int, ρ, q_tot)
PhaseEquil(param_set, ρ, e_int, q_tot)
# Fields
Expand All @@ -210,10 +211,10 @@ $(DocStringExtensions.FIELDS)
struct PhaseEquil{FT, PS} <: ThermodynamicState{FT}
"parameter set, used to dispatch planet parameter function calls"
param_set::PS
"internal energy"
e_int::FT
"density of air (potentially moist)"
ρ::FT
"internal energy"
e_int::FT
"total specific humidity"
q_tot::FT
"temperature: computed via [`saturation_adjustment`](@ref)"
Expand All @@ -234,10 +235,10 @@ and, optionally
- `sat_adjust_method` the numerical method to use.
See the [`Thermodynamics`](@ref) for options.
"""
function PhaseEquil(
function PhaseEquil_ρeq(
param_set::APS,
e_int::FT,
ρ::FT,
e_int::FT,
q_tot::FT,
maxiter::IT = nothing,
temperature_tol::FTT = nothing,
Expand All @@ -257,7 +258,7 @@ function PhaseEquil(
maxiter,
temperature_tol,
)
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot_safe, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q_tot_safe, T)
end

# Convenience method for comparing Numerical
Expand All @@ -266,17 +267,17 @@ end
# should be in sync with the PhaseEquil(...) constructor
function PhaseEquil_dev_only(
param_set::APS,
e_int::FT,
ρ::FT,
e_int::FT,
q_tot::FT;
maxiter::Int = 8,
temperature_tol::FT = FT(1e-1),
sat_adjust_method::Type{NM} = NewtonsMethod,
) where {FT <: Real, NM}
return PhaseEquil(
return PhaseEquil_ρeq(
param_set,
e_int,
ρ,
e_int,
q_tot,
maxiter,
temperature_tol,
Expand Down Expand Up @@ -319,7 +320,7 @@ function PhaseEquil_ρθq(
)
q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
e_int = internal_energy(param_set, T, q_pt)
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q_tot, T)
end

"""
Expand Down Expand Up @@ -356,7 +357,7 @@ function PhaseEquil_pθq(
ρ = air_density(param_set, T, p, PhasePartition(q_tot))
q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
e_int = internal_energy(param_set, T, q)
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q.tot, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q.tot, T)
end

"""
Expand All @@ -378,7 +379,7 @@ function PhaseEquil_ρTq(
phase_type = PhaseEquil
q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
e_int = internal_energy(param_set, T, q)
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q_tot, T)
end

"""
Expand All @@ -401,7 +402,7 @@ function PhaseEquil_pTq(
ρ = air_density(param_set, T, p, PhasePartition(q_tot))
q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
e_int = internal_energy(param_set, T, q)
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q_tot, T)
end

"""
Expand Down Expand Up @@ -438,7 +439,7 @@ function PhaseEquil_peq(
temperature_tol,
)
ρ = air_density(param_set, T, p, PhasePartition(q_tot_safe))
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot_safe, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q_tot_safe, T)
end

"""
Expand Down Expand Up @@ -483,7 +484,7 @@ function PhaseEquil_ρpq(
T = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt)
e_int = internal_energy(param_set, T, q_pt)
end
return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
return PhaseEquil{FT, typeof(param_set)}(param_set, ρ, e_int, q_tot, T)
end

#####
Expand Down Expand Up @@ -567,10 +568,10 @@ function PhaseNonEquil_ρθq(
) where {FT <: Real}
phase_type = PhaseNonEquil
tol = ResidualTolerance(potential_temperature_tol)
T = air_temperature_given_θρq_nonlinear(
T = air_temperature_given_ρθq_nonlinear(
param_set,
θ_liq_ice,
ρ,
θ_liq_ice,
maxiter,
tol,
q_pt,
Expand All @@ -595,7 +596,7 @@ function PhaseNonEquil_pθq(
θ_liq_ice::FT,
q_pt::PhasePartition{FT},
) where {FT <: Real}
T = air_temperature_given_θpq(param_set, θ_liq_ice, p, q_pt)
T = air_temperature_given_pθq(param_set, p, θ_liq_ice, q_pt)
ρ = air_density(param_set, T, p, q_pt)
e_int = internal_energy(param_set, T, q_pt)
return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt)
Expand Down
4 changes: 2 additions & 2 deletions test/data_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,6 @@ dycoms_dataset_path = get_data_folder(dycoms_dataset)
ρ = Array{FT}(ds_PhaseEquil["ρ"][:])
q_tot = Array{FT}(ds_PhaseEquil["q_tot"][:])

ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 4)
# ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 3) # Fails
ts = PhaseEquil_ρeq.(Ref(param_set), ρ, e_int, q_tot, 4)
# ts = PhaseEquil_ρeq.(Ref(param_set), ρ, e_int, q_tot, 3) # Fails
end
Loading

0 comments on commit b78b479

Please sign in to comment.