Adds convenience saturation adjustment methods with reasonable defaults (#302)

* Add convenience SA methods with reasonable defaults

* Reduce maxiter to 2 and bump patch version

Author: tapios

Project.toml

@@ -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"

docs/src/HowToGuide.md

@@ -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
@@ -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
 )
 ```
 
@@ -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
@@ -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

src/saturation_adjustment.jl

@@ -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.
@@ -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 
 # ---------------------------------------------