Make the quadrature scheme the configuration choice

The P3 size-distribution quadrature was configured through a single
quadrature_order knob, with build_quadrature silently mapping a magic
set of orders {16, 32, 40, 64} to GaussLegendre and everything else to
ChebyshevGauss. The scheme — with its parameters, notably the order —
is now the explicit choice: pass a scheme instance,

    Microphysics2MParams(FT; with_ice = true,
        quadrature = Quadrature.GaussLegendre(40))

and build_quadrature only materializes it in the working float type
(GaussLegendre nodes are rebuilt in FT; ChebyshevGauss is closed-form
and passes through). The default ChebyshevGauss(100) preserves the
previous default behavior.

Author: haakon-e

src/Quadrature.jl

@@ -192,7 +192,7 @@ nodes/weights as `SVector{N, FT}`, so per-`integrate` access is a static lookup
 a GPU kernel.
 
 Arbitrary orders `n ≥ 1` are supported (the order is the type parameter `N`).
-`40` is the ClimaAtmos production `quadrature_order`.
+`40` is the ClimaAtmos production order.
 
 # GPU / type-stability
 
@@ -254,26 +254,24 @@ end
 GaussLegendre(n::Int) = GaussLegendre(Float64, n)
 
 """
-    build_quadrature(FT, quadrature_order)
-
-Select and **construct** the quadrature rule for the P3 size-distribution
-integrals from a single `quadrature_order` knob, in element type `FT`. This is
-the host-side, one-shot builder; the returned object is `isbits` and stored on a
+    build_quadrature(FT, scheme::QuadratureRule)
+
+Materialize the quadrature `scheme` for the P3 size-distribution integrals in
+element type `FT`. The scheme — with its parameters, notably the order — is
+the caller's choice. [`ChebyshevGauss`](@ref) has closed-form nodes and passes
+through unchanged; [`GaussLegendre`](@ref) is rebuilt so its stored
+nodes/weights adopt `FT` (a Float64 rule would leak Float64 into Float32
+integrals). Host-side and one-shot; the result is `isbits` and stored on a
 parameter struct for reuse in the (GPU) hot loop.
 
-Gauss-Legendre is preferred for the orders where it is meaningfully more
-accurate than Chebyshev-Gauss on the smooth P3 integrands (≈20× lower error on
-the dominant ice-rain collision integral at matched `n`; see [`GaussLegendre`](@ref)),
-namely `quadrature_order ∈ {16, 32, 40, 64}` (incl. the ClimaAtmos production
-order 40). Any other order falls back to [`ChebyshevGauss`](@ref), preserving the
-default behaviour for non-preferred orders.
+Scheme guidance: at matched order, Gauss-Legendre is substantially more
+accurate on the smooth P3 integrands (≈20× lower error on the dominant
+ice-rain collision integral; see [`GaussLegendre`](@ref)); Chebyshev-Gauss
+remains the conservative default (the cusp-limited `ice_self_collection`
+diagonal is quadrature-limited under both schemes).
 """
-function build_quadrature(::Type{FT}, quadrature_order::Int) where {FT}
-    return if quadrature_order in (16, 32, 40, 64)
-        GaussLegendre(FT, quadrature_order)
-    else
-        ChebyshevGauss(quadrature_order)
-    end
-end
+build_quadrature(::Type{FT}, scheme::ChebyshevGauss) where {FT} = scheme
+build_quadrature(::Type{FT}, scheme::GaussLegendre) where {FT} =
+    GaussLegendre(FT, scheme.n)
 
 end # module Quadrature

src/parameters/Microphysics2MParams.jl

@@ -71,24 +71,23 @@ which constructs the parameterization with components:
     immersion-cap rates. (A prognostic activation-memory model is deferred
     to a follow-up PR.)"
     inp_depletion_model::INPDM = NIceProxyDepletion()
-    "Number of quadrature nodes used for size-distribution integrals
-    (deposition / sublimation, melting, riming, ice-rain collection,
-    sedimentation). Lower → faster, slightly less accurate. Default 100
-    matches the original P3 paper sensitivity studies; n_elem=128 KiD runs
-    show ~5× speed-up at qorder=40 with negligible bulk error."
-    quadrature_order::Int = 100
-    "Pre-constructed quadrature rule for the size-distribution integrals,
-    built once (host-side) from `quadrature_order` via
-    [`Quadrature.build_quadrature`](@ref) and reused in the (GPU) hot loop.
-    It is `isbits` (`GaussLegendre`/`ChebyshevGauss`), so it ships to device
-    kernels with no per-call construction. See [`Quadrature.GaussLegendre`](@ref)."
-    quad::Q = QUAD.build_quadrature(Float64, quadrature_order)
+    "Quadrature scheme used for the size-distribution integrals (deposition /
+    sublimation, melting, riming, ice-rain collection, sedimentation). The
+    scheme — with its parameters, notably the order — is the choice; pass
+    e.g. `Quadrature.ChebyshevGauss(n)` or `Quadrature.GaussLegendre(FT, n)`
+    (see [`Quadrature.build_quadrature`](@ref) for the element-type
+    materialization the toml constructor applies). Lower order → faster,
+    slightly less accurate; `ChebyshevGauss(100)` matches the original P3
+    paper sensitivity studies, and n_elem=128 KiD runs show ~5× speed-up at
+    `GaussLegendre(40)` with negligible bulk error. The object is `isbits`
+    and reused in the (GPU) hot loop with no per-call construction."
+    quad::Q = QUAD.ChebyshevGauss(100)
 end
 Base.show(io::IO, mime::MIME"text/plain", x::P3IceParams) =
     ShowMethods.verbose_show_type_and_fields(io, mime, x)
 
 P3IceParams(toml_dict::CP.ParamDict;
-    is_limited = true, quadrature_order::Int = 100,
+    is_limited = true, quadrature = QUAD.ChebyshevGauss(100),
     inp_depletion_model = NIceProxyDepletion(τ_act = 300),
 ) = P3IceParams(;
     scheme = ParametersP3(toml_dict),
@@ -98,11 +97,10 @@ P3IceParams(toml_dict::CP.ParamDict;
     ice_nucleation = Frostenberg2023(toml_dict),
     rain_freezing = RainFreezing(toml_dict),
     inp_depletion_model,
-    quadrature_order,
-    # Build the quadrature in the working float type, so its nodes/weights
+    # Materialize the scheme in the working float type, so its nodes/weights
     # adopt the integrand's eltype (a Float64 rule would leak Float64 into the
     # Float32 collision integrals). Construction is host-side and one-shot.
-    quad = QUAD.build_quadrature(CP.float_type(toml_dict), quadrature_order),
+    quad = QUAD.build_quadrature(CP.float_type(toml_dict), quadrature),
 )
 
 """
@@ -147,13 +145,13 @@ Create a `Microphysics2MParams` object from a ClimaParams TOML dictionary.
 """
 Microphysics2MParams(toml_dict::CP.ParamDict;
     with_ice = false, is_limited = true,
-    quadrature_order::Int = 100,
+    quadrature = QUAD.ChebyshevGauss(100),
     inp_depletion_model = NIceProxyDepletion(τ_act = 300),
 ) = Microphysics2MParams(;
     # Warm rain parameters (always present)
     warm_rain = WarmRainParams2M(toml_dict; is_limited),
     # Optional ice phase parameters
     ice = with_ice ?
-          P3IceParams(toml_dict; is_limited, quadrature_order, inp_depletion_model) :
+          P3IceParams(toml_dict; is_limited, quadrature, inp_depletion_model) :
           nothing,
 )

test/bulk_tendencies_quadrature_tests.jl

@@ -3,13 +3,14 @@ using Test
 import ClimaParams as CP
 import CloudMicrophysics as CM
 import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.Quadrature as QUAD
 import CloudMicrophysics.BulkMicrophysicsTendencies as BMT
 import CloudMicrophysics.ThermodynamicsInterface as TDI
 
 """
-Sweep test for the P3-ice quadrature order. `quadrature_order` now lives
-on `P3IceParams`, so the sweep is over `Microphysics2MParams(FT;
-with_ice = true, quadrature_order = n)`.
+Sweep test for the P3-ice quadrature resolution. The quadrature scheme
+(with its order) lives on `P3IceParams`, so the sweep is over
+`Microphysics2MParams(FT; with_ice = true, quadrature = ChebyshevGauss(n))`.
 
 Issue 011 measured the worst-case relative error of one single integral
 (`bulk_liquid_ice_collision_sources`) across 5 states and 7 orders, and
@@ -240,9 +241,9 @@ end
 
 function test_quadrature_order_sweep(FT)
     tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
-    # `quadrature_order` now lives on `P3IceParams`, so building one
+    # the quadrature scheme lives on `P3IceParams`, so building one
     # `Microphysics2MParams` per order keeps the sweep clean.
-    make_mp(n) = CMP.Microphysics2MParams(FT; with_ice = true, quadrature_order = n)
+    make_mp(n) = CMP.Microphysics2MParams(FT; with_ice = true, quadrature = QUAD.ChebyshevGauss(n))
 
     orders_and_tol = [
         (100, FT(2e-3)),
@@ -295,7 +296,23 @@ function test_quadrature_order_sweep(FT)
     end
 end
 
-@testset "BMT quadrature-order kwarg (Float64)" begin
+@testset "BMT quadrature-scheme kwarg (Float64)" begin
     test_quadrature_order_sweep(Float64)
 end
+
+@testset "quadrature scheme choice and materialization" begin
+    for FT in (Float32, Float64)
+        # GaussLegendre: rebuilt in the working float type at the given order
+        mp_gl = CMP.Microphysics2MParams(FT; with_ice = true, quadrature = QUAD.GaussLegendre(40))
+        @test mp_gl.ice.quad isa QUAD.GaussLegendre{FT, 40}
+        # ChebyshevGauss: closed-form nodes, passes through unchanged
+        mp_cg = CMP.Microphysics2MParams(FT; with_ice = true, quadrature = QUAD.ChebyshevGauss(64))
+        @test mp_cg.ice.quad isa QUAD.ChebyshevGauss
+        @test mp_cg.ice.quad.n == 64
+    end
+    # default preserves the previous behavior
+    mp = CMP.Microphysics2MParams(Float64; with_ice = true)
+    @test mp.ice.quad isa QUAD.ChebyshevGauss
+    @test mp.ice.quad.n == 100
+end
 nothing