⚗️ QMOM-Intro: An introduction to the quadrature method of moments

A robust Python implementation for solving Population Balance Equations (PBEs) using the Quadrature Method of Moments (QMOM).

This repository provides an easy-to-use implementation of QMOM to simulate the time evolution of the number density function (NDF) of a general phase-space variable $\xi$, by tracking the statistical moments of the NDF. It benchmarks the QMOM solution against Monte Carlo (MC) simulations to validate accuracy.

---

📝 Authors and License

This code has been developed by Mohamed Amine Bouguezzoul during the course of his graduate degree at Polytechnique Montréal, under the supervision of Bianca Viggiano, Bruno Blais, and Fabian Denner. This code is under the copyright of its developers and made available as open-source software under the terms of the MIT License.

---

📂 Repository Structure

The codebase is organized into the following modules:

File Name	Type	Key Responsibility
QMOM.ipynb	📓 Driver	The main simulation hub. Considering the generic single phase-space dimension $\xi$, it defines the pertinent physics ($\dot{\xi}$), sets initial conditions, runs the time-loop, and plots QMOM vs. MC results.
Wheeler.py	🧠 Core	Implements the Adaptive Wheeler Algorithm. Converts statistical moments into quadrature nodes and weights via the Jacobi matrix.
ssp_rk_solver.py	⏱️ Solver	A Strong Stability Preserving (SSP) Runge-Kutta solver. Ensures moment realizability using adaptive time-stepping and log-convexity checks.
flux_calculator.py	⚙️ Utility	Computes the numerical flux (RHS of the transport equation) for both growth (0th order) and aggregation (1st order) processes.
testing_utils.py	🧪 Test	Generates "ground truth" data using Monte Carlo sampling from normal distributions.

---

🧮 Mathematical Insights

1. The Moment Transport Equation

Instead of tracking the full number density function $n(\xi, t)$, QMOM tracks its moments $M_k = \int \xi^k n(\xi) \ \text{d}\xi$. The general evolution equation is:

$$ \frac{\partial M_k}{\partial t} = \underbrace{\int_0^{\infty} k \xi^{k-1} \dot{\xi} n(\xi) \ \text{d}\xi}_{\text{Growth}} + \underbrace{S_k}_{\text{Source Terms}} + \underbrace{\delta_{k,0} J}_{\text{Nucleation}} $$

Tip

The Closure Problem: The growth term often depends on higher moments (e.g., if $\dot{\xi} \propto \xi^{2}$, we need $M_{k+1}$). We cannot solve this directly.

The QMOM Solution: We approximate the NDF as a sum of Dirac deltas (quadrature approximation): $$n(\xi) \approx \sum_{i=1}^{N} w_i \delta(\xi - \xi_i)$$ This transforms the integral into a solvable summation: $\sum w_i k \xi_i^{k-1} \dot{\xi}(\xi_i)$.

2. The Wheeler Algorithm (Moment Inversion)

The Wheeler.py module solves the inverse problem: Given moments $M_k$, find weights $w_i$ and nodes $\xi_i$.

It constructs a recursive relationship using the Chebyshev algorithm to build a tridiagonal Jacobi Matrix ($J_n$):

$$ J_n = \begin{pmatrix} a_0 & \sqrt{b_1} & 0 & \dots \\ \sqrt{b_1} & a_1 & \sqrt{b_2} & \dots \\ 0 & \sqrt{b_2} & a_2 & \dots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} $$

• Nodes ($\xi_i$): Eigenvalues of $J_n$.
• Weights ($w_i$): Derived from the first component of the eigenvectors.

3. Stability & Realizability

A major challenge in moment methods is ensuring the moments remain realizable (i.e., they correspond to a valid, non-negative probability distribution).

Important

Realizability Check: The ssp_rk_solver.py module enforces the log-convexity condition during every time step: $$M_k^2 \le M_{k-1} M_{k+1}$$ If a time step produces moments that violate this condition, the time-step size $\Delta t$ is halved immediately.

---

🚀 Usage Guide

1. Setting up the Physics

In QMOM.ipynb, you define the particle growth law inside the particle_dynamics function:

def particle_dynamics(t, nodes_flat):
    # Example: Size-dependent growth rate dx/dt = x^2
    dxidt = xi**2 
    return dxidt

QMOM Implementation Guide

2. Initializing the Distribution

You can set the initial distribution using a multivariate normal approximation in the notebook:

mu = np.ones(dim)       # Mean
cov = np.eye(dim)       # Covariance
initial_moments = ...   # Calculated via MC_multinormal_moments

3. Flux Calculation

The flux_calculator.py allows for different physical processes:

• zeroth_order_local_flux: Used for standard growth, diffusion, and nucleation (processes affecting single particles).
• first_order_local_flux: Used for aggregation or collision (processes involving particle pairs).

4. Running the Solver

The simulation uses an SSP-RK3 scheme. This is slower than standard RK4 but is strictly required to maintain physical moment sets.

# From QMOM.ipynb
sol_y = ssp_rk.ssp_rk_solver(
    lambda t, m: compute_flux(t, m, N), 
    [t, t+dt], 
    current_moments
)

📊 Visual Output

The notebook automatically generates comparison plots between the Monte Carlo (MC) ground truth and the QMOM approximation.

• MC (Green Squares): Computed by tracking 1000+ individual particles.
• QMOM (Blue Triangles): Computed by evolving just the first $2N$ moments.

Results are saved to the QMOM_output directory.

🎥 Use Case Example: Visualizing Convergence

The animation below demonstrates the accuracy of the Quadrature Method of Moments (QMOM) compared to the "ground truth" Monte Carlo (MC) method.

It visualizes the time evolution of the first 12 moments for a particle population undergoing non-linear size-dependent growth, governed by the law:

$$\frac{\text{d}\xi}{\text{d}t} = \xi^2$$

What to look for:

• Green Squares (MC): The reference solution derived from tracking discrete particles.
• Blue Triangles (QMOM): The moment-based approximation.
• The Effect of Nodes ($N$): Watch how the QMOM solution (Blue) snaps onto the Monte Carlo solution (Green) as the number of quadrature nodes ($N$) increases.
  • With low $N$, the approximation struggles to capture higher-order moments.
  • As $N$ increases, the "Blue Triangles" align perfectly with the "Green Squares," validating that increasing the number of nodes allows QMOM to accurately capture the complex evolution of the distribution.

---

Related References

• Marchisio & Fox (2013) — Comprehensive textbook mainly focusing on QMOM and its derivatives
• Nguyen et al. (2016) — The appendix contains the SSP-RK solver adapted to the moment integration