Skip to content
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
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Thermodynamics"
uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
authors = ["Climate Modeling Alliance"]
version = "0.15.4"
version = "0.15.5"

[deps]
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
Expand Down
33 changes: 10 additions & 23 deletions docs/src/HowToGuide.md
Original file line number Diff line number Diff line change
Expand Up @@ -76,15 +76,12 @@ e_int = -7.0e4
q_tot = 0.01

# Solve for phase equilibrium
# We select a root-finding method from RootSolvers.jl (e.g., SecantMethod, NewtonsMethod)
# We specify the formulation (TD.ρe() means inputs are Density & Internal Energy)
# We use the convenience method which handles defaults automatically
# For ρe(), this defaults to the optimized fixed-iteration solver
sol = TD.saturation_adjustment(
RS.SecantMethod, # Method
params, # Parameter set
TD.ρe(), # Formulation
ρ, e_int, q_tot, # Inputs
10, # Max iterations
1e-3 # Tolerance
ρ, e_int, q_tot # Inputs
)

T_equil = sol.T
Expand Down Expand Up @@ -193,30 +190,23 @@ params_f32 = TD.Parameters.ThermodynamicsParameters(FT)
For GPU performance, avoiding branch divergence is important. For the `ρe` formulation, use the fixed-iteration path by passing `forced_fixed_iters=true` as the last positional argument.

```julia
# GPU-optimized broadcasting
# GPU-optimized broadcasting using the convenience method
# This automatically uses the fast, branch-free fixed-iteration solver
sol = TD.saturation_adjustment.(
RS.NewtonsMethod,
Ref(params_f32),
Ref(TD.ρe()),
ρ_gpu, e_int_gpu, q_tot_gpu,
Ref(3), # maxiter (3 iterations → ~0.1 K accuracy)
Ref(1e-4), # tol (ignored when forced_fixed_iters=true)
nothing, # T_guess (ignored when forced_fixed_iters=true)
true, # forced_fixed_iters
ρ_gpu, e_int_gpu, q_tot_gpu
)
```

For CPU single calls, you can omit the last two arguments to use the standard solver:

```julia
# CPU single-call (standard solver)
# CPU single-call (uses defaults)
sol = TD.saturation_adjustment(
RS.NewtonsMethod,
params_f32,
TD.ρe(),
ρ_val, e_int_val, q_tot_val,
10, # maxiter
1e-4, # tolerance
ρ_val, e_int_val, q_tot_val
)
```

Expand Down Expand Up @@ -257,7 +247,7 @@ If your model handles phase changes (microphysics), you might step `q_liq` and `
ρ = 1.0
e_int = -7.0e4
q_tot = 0.01
sol = TD.saturation_adjustment(RS.NewtonsMethod, params, TD.ρe(), ρ, e_int, q_tot, 15, 1e-4)
sol = TD.saturation_adjustment(params, TD.ρe(), ρ, e_int, q_tot)

# Update thermodynamic variables
T_new = sol.T
Expand Down Expand Up @@ -295,12 +285,9 @@ e_int = TD.internal_energy(params, T_sat, q_tot, q_liq_eq, q_ice_eq)
# Define a function from e_int → q_liq through saturation adjustment
function q_liq_from_e_int(e)
sol = TD.saturation_adjustment(
RS.SecantMethod,
params,
TD.ρe(),
ρ, e, q_tot,
20, # maxiter
1e-5 # tolerance
ρ, e, q_tot
)
return sol.q_liq
end
Expand Down
264 changes: 262 additions & 2 deletions src/saturation_adjustment.jl
Original file line number Diff line number Diff line change
Expand Up @@ -114,14 +114,14 @@ to avoid branch divergence on GPUs.
- `ρ`: Density [kg/m³]
- `e_int`: Specific internal energy [J/kg]
- `q_tot`: Total specific humidity [kg/kg]
- `maxiter`: Number of Newton iterations (3 recommended for ~0.1 K accuracy)
- `maxiter`: Number of Newton iterations (2 recommended for < 0.1 K accuracy)

# Returns
- `NamedTuple` `(; T, q_liq, q_ice, converged)`

# Notes
- The `converged` field is always `true` (no convergence check is performed).
- With `maxiter = 3`, temperature accuracy is better than 0.1 K for typical atmospheric
- With `maxiter = 2`, temperature accuracy is better than 0.1 K for typical atmospheric
conditions (T < 320 K).
- This is an internal helper function. For the public API, use [`saturation_adjustment`](@ref)
with `forced_fixed_iters=true` as a positional argument.
Expand Down Expand Up @@ -672,6 +672,266 @@ function saturation_adjustment(
return (; T, q_liq, q_ice, converged)
end

# ---------------------------------------------
# Convenience methods with reasonable defaults
# ---------------------------------------------

"""
saturation_adjustment(
param_set,
::ρe,
ρ,
e_int,
q_tot;
maxiter::Int = 2,
)

Convenience method for `ρe` formulation with reasonable GPU-optimized defaults.

Uses `RS.NewtonsMethod` with `forced_fixed_iters=true` and `maxiter=2` for fast,
branch-free execution on GPUs. For typical atmospheric conditions (T < 320 K),
this achieves better than 0.1 K accuracy.

For more control over solver parameters, use the full signature with explicit method type.
"""
function saturation_adjustment(
param_set::APS,
::ρe,
ρ,
e_int,
q_tot;
maxiter::Int = 2,
)
return saturation_adjustment_ρe_fixed_iters(
param_set,
ρ,
e_int,
q_tot,
maxiter,
)
end

"""
saturation_adjustment(
param_set,
::pe,
p,
e_int,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)

Convenience method for `pe` formulation with reasonable defaults.

Uses `RS.SecantMethod` with `maxiter=15` and relative tolerance `1e-4`.

For more control over solver parameters, use the full signature with explicit method type.
"""
function saturation_adjustment(
param_set::APS,
::pe,
p,
e_int,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)
FT = eltype(param_set)
_tol = tol === nothing ? RS.RelativeSolutionTolerance(FT(1e-4)) : tol

return saturation_adjustment(
RS.SecantMethod,
param_set,
pe(),
p,
e_int,
q_tot,
maxiter,
_tol,
T_guess,
)
end

"""
saturation_adjustment(
param_set,
::ph,
p,
h,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)

Convenience method for `ph` formulation with reasonable defaults.

Uses `RS.SecantMethod` with `maxiter=15` and relative tolerance `1e-4`.

For more control over solver parameters, use the full signature with explicit method type.
"""
function saturation_adjustment(
param_set::APS,
::ph,
p,
h,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)
FT = eltype(param_set)
_tol = tol === nothing ? RS.RelativeSolutionTolerance(FT(1e-4)) : tol

return saturation_adjustment(
RS.SecantMethod,
param_set,
ph(),
p,
h,
q_tot,
maxiter,
_tol,
T_guess,
)
end

"""
saturation_adjustment(
param_set,
::pθ_li,
p,
θ_li,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)

Convenience method for `pθ_li` formulation with reasonable defaults.

Uses `RS.SecantMethod` with `maxiter=15` and relative tolerance `1e-4`.

For more control over solver parameters, use the full signature with explicit method type.
"""
function saturation_adjustment(
param_set::APS,
::pθ_li,
p,
θ_li,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)
FT = eltype(param_set)
_tol = tol === nothing ? RS.RelativeSolutionTolerance(FT(1e-4)) : tol

return saturation_adjustment(
RS.SecantMethod,
param_set,
pθ_li(),
p,
θ_li,
q_tot,
maxiter,
_tol,
T_guess,
)
end

"""
saturation_adjustment(
param_set,
::ρθ_li,
ρ,
θ_li,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)

Convenience method for `ρθ_li` formulation with reasonable defaults.

Uses `RS.SecantMethod` with `maxiter=15` and relative tolerance `1e-4`.

For more control over solver parameters, use the full signature with explicit method type.
"""
function saturation_adjustment(
param_set::APS,
::ρθ_li,
ρ,
θ_li,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)
FT = eltype(param_set)
_tol = tol === nothing ? RS.RelativeSolutionTolerance(FT(1e-4)) : tol

return saturation_adjustment(
RS.SecantMethod,
param_set,
ρθ_li(),
ρ,
θ_li,
q_tot,
maxiter,
_tol,
T_guess,
)
end

"""
saturation_adjustment(
param_set,
::pρ,
p,
ρ,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)

Convenience method for `pρ` formulation with reasonable defaults.

Uses `RS.SecantMethod` with `maxiter=15` and relative tolerance `1e-4`.

For more control over solver parameters, use the full signature with explicit method type.
"""
function saturation_adjustment(
param_set::APS,
::pρ,
p,
ρ,
q_tot;
maxiter::Int = 15,
tol = nothing,
T_guess = nothing,
)
FT = eltype(param_set)
_tol = tol === nothing ? RS.RelativeSolutionTolerance(FT(1e-4)) : tol

return saturation_adjustment(
RS.SecantMethod,
param_set,
pρ(),
p,
ρ,
q_tot,
maxiter,
_tol,
T_guess,
)
end


# ---------------------------------------------
# Helper functions
# ---------------------------------------------
Expand Down
Loading
Loading