🔢 mathematics

A Nodal Discontinuous Galerkin Method with Rank-Adaptive Velocity Space Representation for the Multiscale BGK Model

This paper presents a novel hybrid algorithm for the Boltzmann-BGK equation that combines a nodal discontinuous Galerkin method in physical space with a rank-adaptive decomposition in velocity space, utilizing multiscale implicit integrators to achieve robust, high-order accuracy and reduced computational complexity across various gas kinetics regimes.

Original authors: Andres Galindo-Olarte, Joseph Nakao, Mirjeta Pasha, Jing-Mei Qiu, William Taitano

Published 2026-06-09

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Andres Galindo-Olarte, Joseph Nakao, Mirjeta Pasha, Jing-Mei Qiu, William Taitano

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict how a giant cloud of gas behaves. This isn't just a simple puff of smoke; it's a complex system where billions of tiny particles are zipping around, bumping into each other, and reacting to their surroundings. In physics, we use a massive, complicated equation called the Boltzmann-BGK equation to describe this dance.

The problem is that this equation is incredibly hard to solve on a computer. It's like trying to track the position and speed of every single grain of sand on a beach simultaneously. The more detail you want, the more computer power you need, until eventually, the task becomes impossible. This is known as the "curse of dimensionality."

This paper introduces a clever new way to solve this puzzle. Here is the breakdown of their method using simple analogies:

1. The Two-Part Strategy: The Map and the Crowd

The authors split the problem into two different worlds:

The Physical World (Where the gas is): Think of this as a map of a city. The city might have weird shapes, holes, or complex buildings (complex geometry). To handle this, they use a method called Nodal Discontinuous Galerkin (DG). Imagine this as a high-resolution, flexible grid that can wrap perfectly around any shape, like a custom-made suit of armor. It's very accurate but can be heavy to carry if you try to make it too detailed everywhere.
The Velocity World (How fast the gas is moving): Think of this as a crowd of people running in different directions. Usually, even in a chaotic crowd, there are patterns. Most people might be walking in a similar direction, or the crowd might be mostly calm. The authors realized that in the "velocity" part of the equation, the data is often low-rank.

The "Low-Rank" Analogy:
Imagine you have a photo of a crowd. A "full-rank" approach tries to save every single pixel individually. A "low-rank" approach is like realizing the crowd is mostly wearing the same color shirt. Instead of saving the color for every person, you just save the "average color" and a few notes on who is wearing something different. This shrinks the file size massively without losing the main picture.

The paper's big innovation is a hybrid approach: They keep the "map" (physical space) detailed and full-rank to handle complex shapes, but they compress the "crowd" (velocity space) using this low-rank trick. This saves a huge amount of computer memory and time.

2. The Time Machine: IMEX Runge-Kutta

The gas particles move in two ways: they drift along (transport) and they crash into each other (collisions).

Drifting is easy to predict.
Colliding is tricky. If the gas is very dense, particles crash into each other billions of times per second. Simulating every single crash would freeze the computer.

The authors use a Time-Stepping Method (IMEX) that acts like a smart time machine. It takes big, fast steps for the easy drifting parts, but it switches to a slow, careful, "implicit" mode when it detects a collision-heavy moment. This allows the simulation to run fast without missing the important physics.

3. Keeping the Balance: Conservation and Limits

When you compress data (like the low-rank trick), you risk losing important rules, like the law that says "mass cannot be created or destroyed."

The Fix: They invented a special "correction" step (Quadrature Corrected Moments). Think of it like a referee in a game who constantly checks the score. If the math starts to drift and lose track of the total mass or energy, the referee instantly corrects it to ensure the laws of physics are obeyed.
The Shock Absorber: Sometimes, gas creates sharp shockwaves (like a sonic boom). Computers often get confused by these sharp edges and start vibrating wildly (spurious oscillations). The authors added a slope limiter, which acts like a shock absorber on a car. When the solution gets too bumpy, the limiter smooths it out just enough to keep the simulation stable without ruining the accuracy.

4. What They Tested

They didn't just talk about theory; they ran the code on several test cases:

Smooth Flow: A gentle, flowing gas. Their method was highly accurate and fast.
Standing Shock: A gas hitting a wall and stopping abruptly. Their method handled the sharp edge without crashing, especially when they used the "shock absorber" (limiters).
Sod Shock Tube: A classic explosion test where high-pressure gas rushes into low-pressure gas. They showed that their method could handle everything from very "thin" gas (where particles rarely hit) to very "thick" gas (where they hit constantly).
Variable Conditions: They even tested a scenario where the gas density changed from place to place, proving the method works in mixed environments.

The Bottom Line

The paper claims to have built a new, faster, and more flexible engine for simulating gas. By treating the "where" (space) and the "how fast" (velocity) differently—using a detailed map for the shape and a compressed summary for the speed—they can solve problems that were previously too expensive for computers.

They emphasize that this is a proof-of-concept for a 1D space / 2D velocity setup, but the framework is designed to be the foundation for tackling even more complex, real-world engineering problems in the future. They did not claim to have solved specific clinical or medical applications in this paper, but rather focused on the mathematical and computational breakthrough for gas dynamics.

Technical Summary: A Nodal Discontinuous Galerkin Method with Rank-Adaptive Velocity Space Representation for the Multiscale BGK Model

Problem Statement

The Boltzmann transport equation governs the evolution of particle distribution functions (PDF) in phase space, bridging microscopic and macroscopic descriptions of gas dynamics. However, numerical modeling of this equation faces the "curse of dimensionality," where computational complexity scales as $O(N^6)$ for a 3D physical and 3D velocity space. Additionally, the equation presents significant challenges regarding:

Multiscale Nature: The system spans from kinetic regimes ( $\varepsilon \sim O(1)$ ) to fluid regimes ( $\varepsilon \ll 1$ ), requiring schemes that remain stable and accurate across scales without inverting dense integral operators.
Geometric Complexity: Real-world engineering applications often involve unstructured domains with non-conformal boundaries, holes, or complex geometries, which are incompatible with the structured hypercubes typically required by existing low-rank tensor decomposition methods.
Conservation and Stability: Numerical schemes must strictly conserve mass, momentum, and energy, maintain positivity, and handle the stiffness of the collision operator without generating spurious oscillations near discontinuities.

Methodology

The authors propose a hybrid numerical framework that decouples the treatment of physical space and velocity space to address these challenges simultaneously. The method combines a Nodal Discontinuous Galerkin (NDG) approach in physical space with a rank-adaptive decomposition in velocity space, integrated via an Implicit-Explicit (IMEX) Runge-Kutta scheme.

1. Spatial Discretization (Physical Space)

The physical domain is discretized using a Nodal Discontinuous Galerkin (NDG) method. This choice allows for:

Arbitrary high-order accuracy in space.
Flexibility in handling unstructured meshes and complex geometries (e.g., polyhedral elements).
Use of Lagrange polynomial bases defined on Gauss-Legendre quadrature nodes, simplifying the evaluation of solution coefficients to nodal values.

2. Velocity Space Discretization (Rank-Adaptive)

To combat the curse of dimensionality, the method applies a rank-adaptive decomposition exclusively to the velocity subspace.

Low-Rank Representation: The solution coefficients at each spatial node are approximated as low-rank matrices using Singular Value Decomposition (SVD). The rank is dynamically adjusted based on a tolerance $\vartheta$ .
Algorithmic Efficiency: The summation of low-rank matrices (required for flux and collision terms) is performed using a QR-SVD approach (Algorithm 1) to maintain minimal rank.
Integration: A "Low-Rank Double Integral" (LRDI) function is employed to compute velocity moments, reducing computational complexity from $O(N_v^2)$ to $O(N_v r)$ , where $r$ is the rank.

3. Time Integration and Multiscale Handling

The system is evolved using a second-order IMEX Runge-Kutta (RK) method:

Explicit Treatment: The transport term is handled explicitly.
Implicit Treatment: The stiff collision operator is treated implicitly using a Diagonally Implicit RK (DIRK) scheme.
Auxiliary Moment Equations: To facilitate large time steps in the fluid limit ( $\varepsilon \to 0$ ), the method evolves discrete moments of the distribution function. These moments are updated explicitly and used to construct the local Maxwellian for the implicit collision step.

4. Conservation and Stability Mechanisms

Quadrature-Corrected Moments (QCM): To ensure exact discrete conservation of mass, momentum, and energy, the local Maxwellian is parameterized by solving a nonlinear system for its moments at each quadrature node. This forces the discrete moments of the collision operator to vanish exactly.
Slope Limiting: To control spurious oscillations near sharp gradients (shocks), a post-processing slope limiter is applied to the macroscopic moments (not the full kinetic solution) at each RK stage. The paper compares minmod and WENO limiters, noting that WENO provides higher-order accuracy for smooth regions while minmod is sufficient for troubled cell identification.

Key Contributions

Hybrid Framework: The paper introduces a novel formulation that retains full-rank structure in physical space (to handle geometric complexity via NDG) while applying rank-adaptive factorization only in velocity space. This avoids the incompatibility of existing low-rank methods with unstructured physical domains.
Asymptotic Preservation: The scheme is proven to converge to the local Maxwellian distribution in the asymptotic limit $\varepsilon \to 0$ , satisfying the Asymptotic Preserving (AP) property.
Conservative Low-Rank Scheme: The integration of QCM with low-rank decomposition ensures that the reduced-order model strictly conserves macroscopic quantities up to machine precision (or truncation tolerance).
Efficient Multiscale Solver: The combination of IMEX time stepping, auxiliary moment equations, and low-rank velocity integration allows for efficient simulation of both kinetic and fluid regimes.

Numerical Results

The method was validated on a suite of 1d2v (one spatial, two velocity dimensions) BGK problems:

Smooth Flow: Convergence studies demonstrated that the method achieves the expected high-order spatial accuracy (2nd and 3rd order for NDG2 and NDG3, respectively) across various Knudsen numbers ( $\varepsilon \in \{1, 10^{-2}, 10^{-6}\}$ ).
Computational Complexity: CPU runtime tests confirmed that the computational cost scales sub-linearly with respect to the velocity grid size $N_v$ , validating the efficiency of the low-rank approach.
Shock Capturing: In standing shock and Sod shock tube problems, the method successfully captured discontinuities. The application of slope limiters (specifically WENO) effectively suppressed non-physical oscillations in density, velocity, and temperature.
Variable Knudsen Number: The method successfully simulated mixed regimes with spatially dependent $\varepsilon(x)$ , capturing steep gradients and transitions between kinetic and fluid behaviors.
Rank Adaptivity: In shock tube simulations, the solution rank remained low (rank-1) in smooth regions and increased slightly near discontinuities, demonstrating the adaptivity of the representation. In the limit $\varepsilon \to 10^{-13}$ , the solution remained uniformly rank-1, consistent with the Maxwellian equilibrium.

Significance and Claims

The authors claim that this work establishes a foundation for extending modern rank-adaptive techniques to realistic engineering settings where structured representations are infeasible. By decoupling the spatial and velocity treatments, the method offers a pathway to solving the Boltzmann equation in complex geometries while mitigating the curse of dimensionality.

The paper modestly positions itself as a proof-of-principle demonstration on the 1d2v model. The authors explicitly state that the current work serves as a basis for future extensions to multidimensional settings with complex geometries. They identify future directions including rigorous analysis of the AP property with velocity discretization, positivity preservation for low-rank solutions, and the application of these techniques to unstructured meshes in higher dimensions. The method is presented as a unified framework addressing conservation, geometric complexity, and dimensionality, but the authors acknowledge that further investigation is required for full multidimensional generalization.