Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rename PhaseEquil outer constructor to PhaseEquil_ρeq, change arg order #47

Merged
merged 1 commit into from
Sep 13, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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
1 change: 1 addition & 0 deletions docs/src/API.md
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ PhaseDry_ρθ
PhaseDry_ρT
PhaseDry_ρp
PhaseEquil
PhaseEquil_ρeq
PhaseEquil_ρTq
PhaseEquil_pTq
PhaseEquil_pθq
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
43 changes: 22 additions & 21 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,34 +211,34 @@ $(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)"
T::FT
end

"""
PhaseEquil
PhaseEquil_ρeq

Moist thermodynamic phase, given
- `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
- `e_int` internal energy
- `ρ` density
- `e_int` internal energy
- `q_tot` total specific humidity
and, optionally
- `maxiter` maximum iterations for saturation adjustment
- `temperature_tol` temperature tolerance for saturation adjustment
- `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