diff --git a/Project.toml b/Project.toml index 3747cd90..f3657cf2 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/HowToGuide.md b/docs/src/HowToGuide.md index df3d16de..9e1ff14c 100644 --- a/docs/src/HowToGuide.md +++ b/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 diff --git a/src/saturation_adjustment.jl b/src/saturation_adjustment.jl index 1303a047..d46ba319 100644 --- a/src/saturation_adjustment.jl +++ b/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 # --------------------------------------------- diff --git a/test/convergence_saturation_adjustment.jl b/test/convergence_saturation_adjustment.jl index bf4613a4..93f87333 100644 --- a/test/convergence_saturation_adjustment.jl +++ b/test/convergence_saturation_adjustment.jl @@ -222,53 +222,6 @@ const TDTP_SA = TD.TemperatureProfiles end end - # Accuracy check for forced_fixed_iters (GPU optimization) - # We expect better than 0.1 K accuracy with maxiter=3. - @testset "forced_fixed_iters accuracy ($FT)" begin - fixed_maxiter = 3 - fixed_tol = FT(1e-4) # Ignored but passed - max_err = FT(0) - chunk_errs = FT[] - for i in idxs - inp = targets(i) - - # Only applicable to ρe formulation - res = TD.saturation_adjustment( - RS.NewtonsMethod, - param_set, - TD.ρe(), - inp.ρ0, - inp.e_int_ρ, - inp.q0, - fixed_maxiter, - fixed_tol, - nothing, # T_guess (ignored) - true, # forced_fixed_iters - ) - - @test res.converged - err = abs(res.T - inp.T0) - max_err = max(max_err, err) - push!(chunk_errs, err) - - if err > FT(0.1) # Outliers appear to occur, if at all, at unrealistically high temperatures/humidities - @info "Outlier detected" FT i err T_ref = inp.T0 p = inp.p0 q_tot = - inp.q0 ρ = inp.ρ0 - end - - check_partition(res.T, inp.ρ0, inp.q0, res.q_liq, res.q_ice) - end - @info "Max error forced_fixed_iters (maxiter=$fixed_maxiter)" FT max_err - - # Allow a few outliers to exceed strict tolerance - @test max_err < FT(0.5) - - # Check percentage of points with error < 0.1 K - n_strict = count(e -> e < FT(0.1), chunk_errs) - pct_strict = n_strict / length(chunk_errs) - @info "Percentage of points with error < 0.1 K" FT pct_strict - @test pct_strict > 0.98 - end end @testset "ρeq solver variants are timed ($FT)" begin diff --git a/test/default_saturation_adjustment.jl b/test/default_saturation_adjustment.jl new file mode 100644 index 00000000..cfb34b1e --- /dev/null +++ b/test/default_saturation_adjustment.jl @@ -0,0 +1,249 @@ +""" +Tests for default convenience methods of saturation_adjustment. + +This suite verifies that the convenience methods with default settings (without +explicit solver method type) work correctly across all formulations: + +- ρe: Uses RS.NewtonsMethod with forced_fixed_iters=true +- Other formulations: Use RS.SecantMethod +""" + +@testset "Thermodynamics - saturation_adjustment default methods" begin + for FT in (Float32, Float64) + param_set = FT == Float64 ? param_set_Float64 : param_set_Float32 + + # Default error tolerance for all methods (0.1 K) + default_atol_temperature = 0.1 + + approx_tight(a, b) = + isapprox(a, b; rtol = FT(100) * eps(FT), atol = FT(100) * eps(FT)) + + function check_partition(Tsol, ρ_eff, q0, ql, qi) + (ql_exp, qi_exp) = TD.condensate_partition(param_set, Tsol, ρ_eff, q0) + @test approx_tight(ql, ql_exp) + @test approx_tight(qi, qi_exp) + return nothing + end + + @testset "ρe default method accuracy ($FT)" begin + profiles = TestedProfiles.PhaseEquilProfiles(param_set, Array{FT}) + (; T, p, ρ, q_tot) = profiles + + # Sample ~120 points across the full profile grid + idxs = unique(round.(Int, range(1, length(T), length = 250))) + + function targets(i::Int) + T0 = T[i] + p0 = p[i] + ρ0 = ρ[i] + q0 = q_tot[i] + + # ρ-based targets + (q_liq_ρ, q_ice_ρ) = TD.condensate_partition(param_set, T0, ρ0, q0) + e_int_ρ = TD.internal_energy_sat(param_set, T0, ρ0, q0) + + return (; T0, p0, ρ0, q0, e_int_ρ) + end + + # Test default ρe method (forced_fixed_iters=true) + max_err = FT(0) + chunk_errs = FT[] + for i in idxs + inp = targets(i) + + # Use default convenience method (no explicit solver type) + res = TD.saturation_adjustment( + param_set, + TD.ρe(), + inp.ρ0, + inp.e_int_ρ, + inp.q0, + ) + + @test res.converged + err = abs(res.T - inp.T0) + max_err = max(max_err, err) + push!(chunk_errs, err) + + if err > FT(default_atol_temperature) + @info "Outlier detected" FT i err T_ref = inp.T0 p = inp.p0 q_tot = + inp.q0 ρ = inp.ρ0 + end + + check_partition(res.T, inp.ρ0, inp.q0, res.q_liq, res.q_ice) + end + @info "Max error ρe default method" FT max_err + + # Allow a few outliers to exceed strict tolerance + @test max_err < 5 * default_atol_temperature + + # Check percentage of points with error < default_atol_temperature + n_strict = count(e -> e < FT(default_atol_temperature), chunk_errs) + pct_strict = n_strict / length(chunk_errs) * 100 + @info "Percentage of points with error < $(default_atol_temperature) K" FT pct_strict + @test pct_strict > 98 + end + + @testset "Other formulations default methods ($FT)" begin + profiles = TestedProfiles.PhaseEquilProfiles(param_set, Array{FT}) + (; T, p, ρ, q_tot) = profiles + + # Sample fewer points for other formulations (faster test) + idxs = unique(round.(Int, range(1, length(T), length = 100))) + + function targets(i::Int) + T0 = T[i] + p0 = p[i] + ρ0 = ρ[i] + q0 = q_tot[i] + + # ρ-based targets + (q_liq_ρ, q_ice_ρ) = TD.condensate_partition(param_set, T0, ρ0, q0) + e_int_ρ = TD.internal_energy_sat(param_set, T0, ρ0, q0) + p_ρ = TD.air_pressure(param_set, T0, ρ0, q0, q_liq_ρ, q_ice_ρ) + θ_ρ = TD.liquid_ice_pottemp(param_set, T0, ρ0, q0, q_liq_ρ, q_ice_ρ) + + # p-based targets + ρ_p = TD.air_density(param_set, T0, p0, q0) + (q_liq_p, q_ice_p) = TD.condensate_partition(param_set, T0, ρ_p, q0) + e_int_p = TD.internal_energy_sat(param_set, T0, ρ_p, q0) + h_p = TD.enthalpy_sat(param_set, T0, ρ_p, q0) + θ_p = + TD.liquid_ice_pottemp_given_pressure( + param_set, + T0, + p0, + q0, + q_liq_p, + q_ice_p, + ) + + return (; + T0, + p0, + ρ0, + q0, + e_int_ρ, + p_ρ, + θ_ρ, + e_int_p, + h_p, + θ_p, + ) + end + + # Test pe default method + for i in idxs + inp = targets(i) + let (; T, q_liq, q_ice, converged) = TD.saturation_adjustment( + param_set, + TD.pe(), + inp.p0, + inp.e_int_p, + inp.q0, + ) + @test isfinite(T) + @test converged + @test isapprox( + T, + inp.T0; + atol = FT(default_atol_temperature), + rtol = FT(0), + ) + ρ_eff = TD.air_density(param_set, T, inp.p0, inp.q0) + check_partition(T, ρ_eff, inp.q0, q_liq, q_ice) + end + end + + # Test ph default method + for i in idxs + inp = targets(i) + let (; T, q_liq, q_ice, converged) = TD.saturation_adjustment( + param_set, + TD.ph(), + inp.p0, + inp.h_p, + inp.q0, + ) + @test isfinite(T) + @test converged + @test isapprox( + T, + inp.T0; + atol = FT(default_atol_temperature), + rtol = FT(0), + ) + ρ_eff = TD.air_density(param_set, T, inp.p0, inp.q0) + check_partition(T, ρ_eff, inp.q0, q_liq, q_ice) + end + end + + # Test pρ default method + for i in idxs + inp = targets(i) + let (; T, q_liq, q_ice, converged) = TD.saturation_adjustment( + param_set, + TD.pρ(), + inp.p_ρ, + inp.ρ0, + inp.q0, + ) + @test isfinite(T) + @test converged + @test isapprox( + T, + inp.T0; + atol = FT(default_atol_temperature), + rtol = FT(0), + ) + check_partition(T, inp.ρ0, inp.q0, q_liq, q_ice) + end + end + + # Test pθ_li default method + for i in idxs + inp = targets(i) + let (; T, q_liq, q_ice, converged) = TD.saturation_adjustment( + param_set, + TD.pθ_li(), + inp.p0, + inp.θ_p, + inp.q0, + ) + @test isfinite(T) + @test converged + @test isapprox( + T, + inp.T0; + atol = FT(default_atol_temperature), + rtol = FT(0), + ) + ρ_eff = TD.air_density(param_set, T, inp.p0, inp.q0) + check_partition(T, ρ_eff, inp.q0, q_liq, q_ice) + end + end + + # Test ρθ_li default method + for i in idxs + inp = targets(i) + let (; T, q_liq, q_ice, converged) = TD.saturation_adjustment( + param_set, + TD.ρθ_li(), + inp.ρ0, + inp.θ_ρ, + inp.q0, + ) + @test isfinite(T) + @test converged + @test isapprox( + T, + inp.T0; + atol = FT(default_atol_temperature), + rtol = FT(0), + ) + check_partition(T, inp.ρ0, inp.q0, q_liq, q_ice) + end + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 28be22c3..5a6bcead 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,6 +25,7 @@ include("correctness.jl") include("exceptions.jl") include("saturation_adjustment.jl") include("convergence_saturation_adjustment.jl") +include("default_saturation_adjustment.jl") include("type_stability.jl") include("optimization_tests.jl") include("ad_tests.jl")