feat: chebyshev-gauss quadrature instead of quadgk

Author: haakon-e

docs/src/API.md

@@ -198,8 +198,8 @@ P3Scheme.∫liquid_ice_collisions
 ### Supporting integral methods
 
 ```@docs
-P3Scheme.∫fdD
-P3Scheme.∫fdD_error
+P3Scheme.integrate
+P3Scheme.ChebyshevGauss
 P3Scheme.integral_bounds
 ```
 

src/P3_integral_properties.jl

@@ -1,86 +1,143 @@
 
 """
-    ∫fdD(f, dist; [p = 1e-6], kwargs...)
+    integrate(f, a, b, quad = ChebyshevGauss(100))
+
+ Approximate the definite integral ∫ₐᵇ f(x) dx using Chebyshev-Gauss quadrature of the first kind.
+
+# Mathematical Background
+
+ This method transforms the integral to the standard interval [-1, 1] and applies 
+ Chebyshev-Gauss quadrature. The transformation is:
+ 
+     y = (2x - (a+b)) / (b-a)    →    x = (b-a)y/2 + (a+b)/2
+ 
+ with Jacobian dx/dy = (b-a)/2.
+ 
+ The integral becomes:
+
+     ∫ₐᵇ f(x) dx = (b-a)/2 ∫₋₁¹ f(x(y)) dy
+ 
+ Using the Chebyshev-Gauss quadrature identity:
+
+     ∫₋₁¹ g(y) dy = ∫₋₁¹ [g(y)√(1-y²)] / √(1-y²) dy ≈ (π/n) ∑ᵢ₌₁ⁿ g(yᵢ)√(1-yᵢ²)
+ 
+ where yᵢ = cos((2i-1)π/(2n)) are the Chebyshev nodes of the first kind.
+ 
+ # Final Formula
+
+     ∫ₐᵇ f(x) dx ≈ (b-a)π/(2n) ∑ᵢ₌₁ⁿ f(xᵢ)√(1-yᵢ²)
+ 
+ where:
+ - yᵢ = cos((2i-1)π/(2n)) for i = 1, ..., n
+ - xᵢ = (b-a)yᵢ/2 + (a+b)/2 
+ - All quadrature weights are equal: π/n
 
-Integrate the function `f` over the size distribution `dist`
+# Arguments
+ - `f`: Function to integrate
+ - `a`, `b`: Integration bounds
+
+# Keyword arguments
+ - `quad`: Quadrature scheme, default: `ChebyshevGauss(100)`
+
+# Returns
+ Approximation to the definite integral ∫ₐᵇ f(x) dx
+
+# Notes
+ This method achieves spectral convergence for smooth functions and is particularly 
+ effective for analytic functions. The √(1-y²) weighting factor helps handle 
+ functions with mild singularities at the interval endpoints.
+ Ref: https://en.wikipedia.org/wiki/Chebyshev–Gauss_quadrature
+"""
+function integrate(f, a, b; quad = ChebyshevGauss(100))
+    FT = eltype(float(a))
+    # Pre-compute transformation parameters
+    scale_factor = (b - a) / 2
+    shift = (a + b) / 2
+
+    # Compute integral using Chebyshev-Gauss quadrature
+    result = zero(f(a))
+    a == b && return result  # return early if a == b
+    (; n) = quad
+    for i in 1:n
+        # Node on [-1, 1] interval
+        y = node(quad, FT(i), n)
+
+        # Node on [a, b] interval
+        x = scale_factor * y + shift
+
+        # Total weight √(1 - y²) * wᵢ
+        w = inv_weight_fun(quad, y) * weight(quad, FT(i), n)
+
+        # Accumulate: f(x) * √(1 - y²) * wᵢ
+        result += f(x) * w
+    end
+
+    return scale_factor * result
+end
 
-!!! note "Usage"
-    This function is useful for integrating functions over the size distribution `dist`.
-    It is a light wrapper around `QGK.quadgk` that automatically inserts appropriate
-    size distribution thresholds as integration limits.
+"""
+    integrate(f, bnds...; quad = ChebyshevGauss(100))
 
-This method calls [`∫fdD_error`](@ref), which returns both the value of the integral 
-    and the estimated error. Since the error is typically not of interest, this method
-    only returns the value of the integral.
+Integrate the function `f` over each subinterval of the integration bounds, `bnds`.
 
 # Arguments
-- `f`: The function to integrate
-- `dist`: The distribution object, passed to [`integral_bounds`](@ref) to set the integration limits.
-- `p`: The integration bounds are set to the `p`-th and `1-p`-th quantiles of the size distribution.
-    Default: `p = 1e-6` (i.e. 99.9998% of the size distribution is integrated).
-- `kwargs`: Additional optional keyword arguments to pass to [`QGK.quadgk`](https://juliamath.github.io/QuadGK.jl/stable/api/#QuadGK.quadgk)
-    - `rtol`: The relative tolerance for the integration, default: `rtol = sqrt(eps(FT))`
-    - `atol`: The absolute tolerance for the integration, default: `atol = 0`
-    - `maxevals`: The maximum number of function evaluations, default: `maxevals = 10^7`
-    - `order`: The order of the quadrature rule, default: `order = 7`
+ - `f`: Function to integrate
+ - `bnds`: A tuple of bounds, `(a, b, c, d, ...)`
+ - `quad`: Quadrature scheme, default: `ChebyshevGauss(100)`
 
-# Returns
-- `value`: The value of the integral
+ The integral is computed as the sum of the integrals over each subinterval,
+ `(a, b), (b, c), (c, d), ...`.
+"""
+function integrate(f, bnds...; quad = ChebyshevGauss(100))
+    # compute integral over each subinterval (a, b), (b, c), (c, d), ...
+    return sum(integrate(f, a, b; quad) for (a, b) in zip(Base.front(bnds), Base.tail(bnds)))
+end
 
-!!! note "Integral accuracy"
-    To achieve highest accuracy, which can be challenging when integrating over 
-    the size distribution, it is recommended to increase the `order` of the 
-    quadrature rule and set `rtol = 0`. Experimentally, `order = 44` has been found
-    to be sufficient for most cases. 
+"""
+    ChebyshevGauss(n)
 
-    For convenience, passing `accurate = true` will set `rtol = 0` and `order = 44`.
+Quadrature scheme for Chebyshev-Gauss quadrature of the first kind.
 
-# Examples
+# Arguments
+ - `n`: Number of quadrature points
 
-```jldoctest
-julia> import CloudMicrophysics.Parameters as CMP,
-              CloudMicrophysics.P3Scheme   as P3
+# Available methods
 
-julia> params = CMP.ParametersP3(Float64);
+    node(::ChebyshevGauss, i::FT, n) where {FT}
+    weight(::ChebyshevGauss, i::FT, n) where {FT}
+    inv_weight_fun(::ChebyshevGauss, y)
 
-julia> state = P3.get_state(params; F_rim = 0.0, ρ_rim = 400.0, L_ice = 0.002, N_ice = 1000.0);
+- `node(quad, i, n)`: Return the `i`-th node of the `n`-point Chebyshev-Gauss quadrature scheme.
+- `weight(quad, i, n)`: Return the `i`-th weight of the `n`-point Chebyshev-Gauss quadrature scheme.
+- `inv_weight_fun(quad, y)`: Return the inverse of the weight function `w(x)`.
 
-julia> logλ = P3.get_distribution_logλ(state);
 
-julia> f(D) = D^3 * P3.size_distribution(state, logλ)(D);  # Define a function to integrate
 
-julia> P3.∫fdD(f, state, logλ; p = 0.01)  # Integrate the function
-0.0008519464332296608
+# Mathematical Background
 
-julia> P3.∫fdD(state, logλ; p = 0.01) do D  # Integrate with a `do`-block
-           P3.ice_mass(state, D) * P3.size_distribution(state, logλ)(D)
-       end
-0.0017027833723511712
-```
-"""
-function ∫fdD(f, state::P3State, logλ; p = 1e-6, moment_order = 0, kwargs...)
-    return ∫fdD_error(f, state, logλ; p, moment_order, kwargs...)[1]
-end
+A method to approximate the value of integrals of the kind
 
-"""
-    ∫fdD_error(f, dist; p, [moment_order = 0], [accurate = false], kwargs...)
+    ∫_{-1}^{1} f(x) w(x) dx ≈ ∑_1^n f(x_i) wᵢ(x_i)
 
-Integrate the function `f` over the size distribution `dist`
+where `w(x) = 1 / √(1 - x^2)` is the weight function, `x_i` are the nodes, and `wᵢ` are the weights.
 
-# Returns
-- `value`: The value of the integral
-- `error`: The estimated error of the integral
+If we are interested in only the integral of `f(x)`, we can instead integrate `g(x) = f(x) / w(x)`,
 
-# Notes
-See [`∫fdD`](@ref), which only returns the value of the integral and not the error, for details.
+    ∫_{-1}^{1} g(x) w(x) dx ≈ ∑_1^n g(x_i) wᵢ(x_i) = ∑_1^n f(x_i) wᵢ(x_i) / w(x_i)
+
+# References
+
+- https://en.wikipedia.org/wiki/Chebyshev–Gauss_quadrature
+- https://en.wikipedia.org/wiki/Gaussian_quadrature
 """
-function ∫fdD_error(f, state::P3State, logλ; p, moment_order = 0, accurate = false, kwargs...)
-    # Get integration bounds
-    bnds = integral_bounds(state, logλ; p, moment_order)
-    # Use a more accurate quadrature rule if requested
-    accurate && (kwargs = (; rtol = 0, order = 44, kwargs...))
-    return QGK.quadgk(f, bnds...; kwargs...)
+struct ChebyshevGauss
+    n::Int
 end
+Base.broadcastable(quad::ChebyshevGauss) = (quad,)
+
+@inline node(::ChebyshevGauss, i::FT, n) where {FT} = cospi((2float(FT(i)) - 1) / (2n))
+@inline weight(::ChebyshevGauss, i::FT, n) where {FT} = float(FT(π)) / n
+@inline inv_weight_fun(::ChebyshevGauss, y) = √(1 - y^2)
 
 """
     integral_bounds(state::P3State, logλ; p, moment_order = 0)

src/P3_processes.jl

@@ -60,13 +60,13 @@ end
  - `logλ`: the log of the slope parameter [log(1/m)]
 
 # Keyword arguments
- - `∫kwargs`: Named tuple of keyword arguments passed to [`∫fdD`](@ref)
+ - `∫kwargs...`: Additional keyword arguments passed to the quadrature rule
 
 Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015)).
 """
 function ice_melt(
     velocity_params::CMP.Chen2022VelType, aps::CMP.AirProperties, tps::TDI.PS, Tₐ, ρₐ, dt, state::P3State, logλ;
-    ∫kwargs = (;),
+    ∫kwargs...,
 )
     # Note: process not dependent on `F_liq`
     # (we want ice core shape params)
@@ -83,7 +83,8 @@ function ice_melt(
 
     # Integrate
     fac = 4 * K_therm / L_f * (Tₐ - T_freeze)
-    dLdt = fac * ∫fdD(state, logλ; ∫kwargs...) do D
+    bnds = integral_bounds(state, logλ; p = 1e-6)
+    dLdt = fac * integrate(bnds...; ∫kwargs...) do D
         ∂ice_mass_∂D(state, D) * F_v(D) * N′(D) / D
     end
 
@@ -257,7 +258,7 @@ Returns a function `liquid_integrals(Dᵢ)` that computes the liquid particle in
 - `liq_bounds`: integration bounds for liquid particles
 
 # Keyword arguments
-- `∫kwargs...`: Additional keyword arguments passed to `QuadGK.quadgk`
+- `∫kwargs...`: Additional keyword arguments passed to the quadrature rule
 
 # Notes
 The function `liquid_integrals(Dᵢ)` returns a tuple `(∂ₜN_col, ∂ₜM_col, ∂ₜB_col)` 
@@ -268,17 +269,16 @@ The function `liquid_integrals(Dᵢ)` returns a tuple `(∂ₜN_col, ∂ₜM_col
 """
 function get_liquid_integrals(n, ∂ₜV, m_liq, ρ′_rim, liq_bounds; ∫kwargs...)
     function liquid_integrals(Dᵢ)
-        ((∂ₜN_col, ∂ₜM_col, ∂ₜB_col), _) =
-            QGK.quadgk(liq_bounds...; ∫kwargs...) do D
-                return SA.SVector(
-                    # ∂ₜN_col = ∫ ∂ₜV ⋅ n ⋅ dD
-                    ∂ₜV(Dᵢ, D) * n(D),
-                    # ∂ₜM_col = ∫ ∂ₜV ⋅ n ⋅ m_liq ⋅ dD
-                    ∂ₜV(Dᵢ, D) * n(D) * m_liq(D),
-                    # ∂ₜB_col = ∫ ∂ₜV ⋅ n ⋅ m_liq / ρ′_rim ⋅ dD 
-                    ∂ₜV(Dᵢ, D) * n(D) * m_liq(D) / ρ′_rim(Dᵢ, D),
-                )
-            end
+        (∂ₜN_col, ∂ₜM_col, ∂ₜB_col) = integrate(liq_bounds...; ∫kwargs...) do D
+            return SA.SVector(
+                # ∂ₜN_col = ∫ ∂ₜV ⋅ n ⋅ dD
+                ∂ₜV(Dᵢ, D) * n(D),
+                # ∂ₜM_col = ∫ ∂ₜV ⋅ n ⋅ m_liq ⋅ dD
+                ∂ₜV(Dᵢ, D) * n(D) * m_liq(D),
+                # ∂ₜB_col = ∫ ∂ₜV ⋅ n ⋅ m_liq / ρ′_rim ⋅ dD 
+                ∂ₜV(Dᵢ, D) * n(D) * m_liq(D) / ρ′_rim(Dᵢ, D),
+            )
+        end
         return ∂ₜN_col, ∂ₜM_col, ∂ₜB_col
     end
     return liquid_integrals
@@ -299,7 +299,7 @@ Computes the bulk collision rate integrands between ice and liquid particles.
 - `ice_bounds`: integration bounds for ice particles, from [`integral_bounds`](@ref)
 
 # Keyword arguments
-- `∫kwargs...`: Additional keyword arguments passed to `QuadGK.quadgk`
+- `∫kwargs...`: Additional keyword arguments passed to the quadrature rule
 
 # Returns
 A tuple of 8 integrands, see [`∫liquid_ice_collisions`](@ref) for details.
@@ -333,9 +333,7 @@ function ∫liquid_ice_collisions(n_i, ∂ₜM_max, cloud_integrals, rain_integr
             n * 𝟙_wet * ∂ₜM_col,          # ∫𝟙_wet_M_col, wet growth indicator
         )
     end
-
-    (rates, _) = QGK.quadgk(liquid_ice_collisions_integrands, ice_bounds...; ∫kwargs...)
-    return rates
+    return integrate(liquid_ice_collisions_integrands, ice_bounds...; ∫kwargs...)
 end
 
 """
@@ -378,7 +376,7 @@ A tuple `(QCFRZ, QCSHD, NCCOL, QRFRZ, QRSHD, NRCOL, ∫M_col, BCCOL, BRCOL, ∫
 function ∫liquid_ice_collisions(
     state, logλ,
     psd_c, psd_r, L_c, N_c, L_r, N_r,
-    aps, tps, vel, ρₐ, T, m_liq,
+    aps, tps, vel, ρₐ, T, m_liq; ∫kwargs...,
 )
     FT = eltype(state)
 
@@ -388,13 +386,10 @@ function ∫liquid_ice_collisions(
     n_i = DT.size_distribution(state, logλ)               # n_i(Dᵢ)
 
     # Initialize integration buffers by evaluating a representative integral
-    ice_bounds = integral_bounds(state, logλ; p = 0.00001)
-    mm = FT(1e-3)
-    bounds_c = FT[0; 0.01mm; 0.1mm; 1mm; 10mm; 100mm; 1]  # TODO: Replace by quantiles method
-    bounds_r = FT[0; 0.01mm; 0.1mm; 1mm; 10mm; 100mm; 1]  # TODO: Replace by quantiles method
-    segbuf_c = QGK.quadgk_segbuf(n_c, bounds_c...)[3]
-    segbuf_r = QGK.quadgk_segbuf(n_r, bounds_r...)[3]
-    segbuf_ice = QGK.quadgk_segbuf(n_i, ice_bounds...)[3]
+    p = FT(0.00001)
+    ice_bounds = integral_bounds(state, logλ; p)
+    bounds_c = CM2.get_size_distribution_bounds(psd_c, L_c, ρₐ, N_c, p)
+    bounds_r = CM2.get_size_distribution_bounds(psd_r, L_r, ρₐ, N_r, p)
 
     # Integrand components
     # NOTE: We assume collision efficiency, shape (spherical), and terminal velocity is the 
@@ -403,12 +398,10 @@ function ∫liquid_ice_collisions(
     ρ′_rim = compute_local_rime_density(vel, ρₐ, T, state)  # ρ′_rim(Dᵢ, Dₗ)
     ∂ₜM_max = compute_max_freeze_rate(aps, tps, vel, ρₐ, T, state)  # ∂ₜM_max(Dᵢ)
 
-    cloud_integrals = get_liquid_integrals(n_c, ∂ₜV, m_liq, ρ′_rim, bounds_c; eval_segbuf = segbuf_c)  # (∂ₜN_c_col, ∂ₜM_c_col, ∂ₜB_c_col)
-    rain_integrals = get_liquid_integrals(n_r, ∂ₜV, m_liq, ρ′_rim, bounds_r; eval_segbuf = segbuf_r)  # (∂ₜN_r_col, ∂ₜM_r_col, ∂ₜB_r_col)
+    cloud_integrals = get_liquid_integrals(n_c, ∂ₜV, m_liq, ρ′_rim, bounds_c; ∫kwargs...)  # (∂ₜN_c_col, ∂ₜM_c_col, ∂ₜB_c_col)
+    rain_integrals = get_liquid_integrals(n_r, ∂ₜV, m_liq, ρ′_rim, bounds_r; ∫kwargs...)  # (∂ₜN_r_col, ∂ₜM_r_col, ∂ₜB_r_col)
 
-    return ∫liquid_ice_collisions(
-        n_i, ∂ₜM_max, cloud_integrals, rain_integrals, ice_bounds; eval_segbuf = segbuf_ice,
-    )
+    return ∫liquid_ice_collisions(n_i, ∂ₜM_max, cloud_integrals, rain_integrals, ice_bounds; ∫kwargs...)
 end
 
 """
@@ -451,7 +444,7 @@ A `NamedTuple` of `(; ∂ₜq_c, ∂ₜq_r, ∂ₜN_c, ∂ₜN_r, ∂ₜL_rim, 
 function bulk_liquid_ice_collision_sources(
     params, logλ, L_ice, F_rim, ρ_rim,
     psd_c, psd_r, L_c, N_c, L_r, N_r,
-    aps, tps, vel, ρₐ, T,
+    aps, tps, vel, ρₐ, T; ∫kwargs...,
 )
     FT = eltype(params)
     (; τ_wet) = params
@@ -466,7 +459,7 @@ function bulk_liquid_ice_collision_sources(
     (QCFRZ, QCSHD, NCCOL, QRFRZ, QRSHD, NRCOL, ∫∂ₜM_col, BCCOL, BRCOL, ∫𝟙_wet_M_col) = ∫liquid_ice_collisions(
         state, logλ,
         psd_c, psd_r, L_c, N_c, L_r, N_r,
-        aps, tps, vel, ρₐ, T, m_liq,
+        aps, tps, vel, ρₐ, T, m_liq; ∫kwargs...,
     )
 
     # Bulk wet growth fraction