Information geometric regularization

The first inviscid regularization of high-speed gas dynamics

The problem: Shock formation in gas dynamics

This post summarizes joint work with Ruijia Cao.

Shock formation

The density and velocity profiles of compressible and (nearly) inviscid gases steepen over time, as faster gas particles catch up to slower ones. Eventually, the profiles become almost discontinuous. These shock waves are crucial phenomena in many domains, such as astrophysics and aerospace applications.

Shock waves occur in supernovae and atmospheric reentry (copyright by Nasa).

Shocks— a multiscale problem

Shock waves appear discontinuous on macroscopic scales. But on sufficiently small scales the viscosity of the gas preserves their smoothness. This usually happens on the scale of the mean free path length — the average distance that gas particles travel between collisions. In atmosspheric conditions, this distance is on the order of microns, vastly smaller than the scales relevant to aerospatial simulations and impossible to resolve in numerical simulations.

Shocks are macroscopically discontinuous...

The multiscale nature of shocks. Shocks appear as macroscopic discontinuities (left image, taken from Takayama 2019). But they are smooth on the microscopic mean free scale (right image, taken from Koreeda et al. 1995).

Numerical challenges due to shocks

Most numerical methods use smooth functions, such as polynomial, trigonometric, neural network, or radial basis functions, to approximate the solution. If the solution is discontinuous, the Gibbs-Runge oscillations can cause the simulation to break down. Artifically adding viscosity avoids oscillations but introduces excessive energy dissipation.

Directly using higher-order methods on problems with shocks causes oscillatory artifacts.

Simply adding viscosity oversmoothes the solution.

Oscillation and dissipation. We plot the velocity $u$ and density $\rho$ of a unidimensional problem. Left: Using higher order methods without additional viscosity results causes Gibbs-Runge oscillations to spoil the simulation. Right: Adding viscosity avoids this, but smears out the solution over time.

Previous approach: Viscous regularization

The field of shock capturing designs discrete or continuous models that are well-behaved on grid-resolvable scales but correctly account for the effects of the microscopic shock, without excessive energy dissipation. Explicitly or implicitly, virtually all attempts add viscosity locally near the shock to regularize it without spurious dissipation in smooth regions of the flow. This delicate balance is a major obstacle in simulating flow involving complex interactions of shock waves with physical oscillations due to acoustic or turbulent effects. The video below shows the behavior of a typical approach, on a shock and a high-frequency pressure wave, for varying regularization strengths. Weak regularization results in near-discontinuous solutions, while strong regularization damps out the sound wave. Furthermore, the shock profile’s derivative always has a near-singularity, limiting smoothness-based approximations.

Viscosity-based regularization struggles to regularize shocks without dissipating physical oscillations.

Regularization with local viscosity faces a tradeoff. Strong regularization is needed to expand shocks on grid-resolvable scales. But this creates excessive dissipation of physical oscillatory structures due to acoustic or turbulent effects. Even for strong regularization, the derivative is near-discontinuous, which creates Gibbs-Runge oscillations in strong shock waves. We plot the pressure $p$ of the nominal and regularized solution.

Our solution: Information geometric regularization (IGR)

To overcome this challenging tradeoff, we recently proposed information geometric regularization (IGR), the first inviscid regularization. It regularizes gas dynamics without adding any viscosity. Instead, it adds to the pressure $p$ the entropic pressure ${\color{#B8420F} \Sigma}$ defined by

\[\rho^{-1} {\color{#B8420F}\Sigma} - {\color{seagreen} \alpha} \operatorname{div}(\rho^{-1} \nabla {\color{#B8420F}\Sigma}) = {\color{#2E6B69}\alpha} \left(\operatorname{trace}^2\left([\mathrm{D} \mathbf{u}]\right) + \operatorname{trace}\left([\mathrm{D} \mathbf{u}]^2\right) \right).\]

Here, $\rho$ is the density, $\mathbf{u}$ the velocity field, and $\mathrm{D} \mathbf{u}$ its Jacobian. The parameter ${\color{seagreen} \alpha}$ controls the strength of the regularization. The regularized shock profiles have thickness $\propto \sqrt{\alpha}$. The video below shows that IGR replaces shocks with smooth profiles without damping acoustic waves.

Information geometric regularization is an inviscid regularization.

Inviscid regularization. Information geometric regularization replaces shocks with smooth profiles, without dissipating physical oscillations. We plot the pressure $p$ of the nominal and regularized solution.

In general relativity, gravitational forces are manifestations of the curvature of spacetime. Similarly the entropic pressure $\Sigma$ modifies the notion of straight line in the space of particle trajectories. In the modified geometry, particle trajectories merge but never cross, preventing the formation of singular shocks while maintaining nominal long-time behavior.

Information geometric regularization modifies the geometry of particle trajectories.

Regularizing particle trajectories. Information geometric regularization modifies the notion of straight line underlying particle trajectories. In the modified geometry, particle paths do not collide but asymptotically approach each other. This avoids singularity formation, without changing the long-time behavior.

To explain this modification in more detail, we start with some background on shock formation.

Background: Four perspectives on shock formation

The geometry of shock formation

In the case of one spatial dimension and vanishing pressure, the geometry of shock formation is especially simple. Gas particles starting at $x$ follow straight lines

\[\Phi_t\left(x\right) = x + t u_0(x),\]

determined by their initial velocity. Equivalently, the deformation map $\Phi_t$ follows a straight line $t \mapsto \Phi_t$ in a suitable function space.

Collison of particles / with the boundary

Shocks form when gas particles collide or equivalently, $\Phi_t$ reaches the boundary of the set of invertible functions (diffeomorphism manifold)

Shocks form when the flow map reaches the boundary of the set of invertible maps.

Shocks from collision. In the pressureless case, Shock formation arises from the collision of particle trajectories (left, center). equivalently, from the flow map $\Phi_t$ reaching the boundary of the set of invertible maps (right).

Collapse and constraint activation

Two lesser-known perspectives arise from mass transport and optimization. The mass density $\rho_t$ at time $t$ is the so-called push-forward $\left(\Phi_t\right)_{\#} \rho_0$ of the initial density $\rho_0$. Shocks form when the mass acumulates on individual points. The map $\Phi_t$ also solves the minimization problem

\[\min \limits_{\partial_x \Phi \geq 0 } - t \int \limits_{\mathbb{R}} u_0(x) \left(\Phi(x) - x\right) \mathrm{d} x + \frac{1}{2}\int \limits_{\mathbb{R}} \left(\Phi(x) - x \right)^2\mathrm{d} x.\]

From this perspective, shocks form when the inequality constraint becomes active.

Shocks form when the mass density collapses on individual points.

Shocks form when the constraint becomes active.

Shocks from mass collapse and constraint satisfaction. The formation of shocks occurs when mass density is compressed to individual points (left) or when the constraint in the minimization form activates (right).

After shock

After shock formation, particle trajectories merge, $\Phi_t$ evolves on the boundary of the set of invertible functions, propagating point masses, and satisfies the constrained problem.

..., the flow map evolves on the boundary, ...

..., and the flow map solves a constrained optimization problem.

After shock formation, particle trajectories merge, flow maps stick to the boundary of the set of invertible functions, atoms of the mass density are propagated, and the flow map solves a constrained optimization problem.

Deriving information geometric regularization

Interior point methods for fluid dynamics

The minimization definition of $\Phi_t$ resembles so-called homotopy interior point methods. To minimize an objective, (like $\int_{\mathbb{R}} u_0(x) \left(\Phi(x) - x\right) \mathrm{d} x$), subject to a constraint (like $\partial_x \Phi \geq 0$), they add the Bregman divergence of a strictly convex function $\psi$ to regularize the problem.

\[\min \limits_{\partial_x \Phi \geq 0 } - t \int \limits_{\mathbb{R}} u_0(x) \left(\Phi(x) - x\right) \mathrm{d} x + \underbrace{ {\color{#B8420F} \psi(\Phi)} {\color{#B0ABA8} -}{\color{#B0ABA8} \psi(\mathrm{Id}) - \left[\mathrm{D}\psi(\mathrm{Id})\right] \left(\Phi - \mathrm{Id}\right)}}_{\text{Bregman divergence } \mathbb{D}\left(\Phi \; \middle\| \; \mathrm{Id}\right) \text{ of } \psi}\]

As the relative weight $t$ of the objective increases, the solution of the regularized problem converges to the minimizer of the original problem. Choosing a so-called barrier $\Psi$ that blows up at the boundary ensures that the constraint never becomes active, greatly simplifying the solution of the regularized problems. The deformation map $\Phi_t$ arises from choosing

\[\psi_0(\Phi) = \frac{1}{2} \int \limits_{\mathbb{R}} \left(\Phi(x)\right)^2 \mathrm{d} x.\]

Shock formation thus arises from the failure of $\psi_{0}$ to keep the regularized solution in the interior of the feasible set — the barrier is too weak! We thus regularize the barrier with a logarithmic term that blows up as we approach the boundary of the constraint set

\[\psi_{\alpha}(\Phi) = \psi_0(\Phi) + \alpha \int \limits_{\mathbb{R}} - \log(\partial_x \Phi(x)) \mathrm{d} x.\]

Upon regularization, solution paths asymptotically approach the boundary.

Barrier functions. In the unregularized problem, the minimizers move on a straight line until they hit the boundary of the constraint set. With the logarithmic barrier, they instead asymptotically approach the domain boundary.

We see that with the new barrier, the solution path asymptotically approaches but never touches the boundary of the feasible set. The choice of $\alpha$ modulates the strength of the regularization. Larger $\alpha$ yield more boundary-avoiding behavior, but for all choices of $\alpha$ the asymptotic behavior is the same as the nominal solution.
Expressing the minimizers as particle trajectories, we see that in the modified trajectories particles approach each other without touching. This avoids formation of singular shocks while maintaining the nominal long-time behavior!

From optimization paths to particle trajectories.

From minimizer paths to particle trajectories. Reexpressing the regularized minimizer paths as particle trajectories, we see that they avoid singular shocks while maintaining the nominal long-time behavior. The size of $\alpha$ determines the strength of the regularization and thus how close the regularized paths are to the nominal ones.

From barriers to geodesics

The minimization form of $\Phi_t$ is highly specific to the unidimensional pressureless case. To obtain a general approach we take a differential geometric perspective. The key idea of differential geometry is to distinguish between points $p \in \mathcal{M}$ on a manifold and infinitesimal directions of displacement living in their tangent space $T_{p}\mathcal{M}$. These two objects are connected by the exponential map $(p, v) \mapsto \operatorname{Exp}_{p}(v)$ that returns the destination reached by traveling from $p$ in the direction of $v \in T_{p} \mathcal{M}$. Linear algebra uses the Euclidean exponential map $\operatorname{Exp}_{p}(v) := p + v$, but in general there are numerous ways to define “traveling on a straigth line in an initial direction”. A strictly convex barrier function defines a so-called dual exponential map $\operatorname{Exp}^{\psi}_{(\cdot)}\left(\cdot\right)$ as the Euclidean exponential map in the space of gradients of $\psi$. The resulting straight lines of the form $t \mapsto \operatorname{Exp}^{\psi}_{p}\left(t v\right)$ are called dual geodesics.

Dual geodesics bend away from the domain boundary.

The geometry of interior point methods. The minimizer paths of homotopy interior point methods are dual geodesics --- a notion of straight line that bends away from the domain boundary.

Crucially for our purposes, the minimizing paths $\Phi_t$ of the last section are dual geodesics generated by $\psi_{\alpha}$. This generalizing our regularization beyond the unidimensional pressureless case by replacing particle trajectories with dual geodesics.

Information geometric regularization and information geometric mechanics

After deriving modified differential equations for $\Phi$, we can recover the resulting velocity $\mathbf{u}$ and density $\rho$, yielding a modifed conservation law with $\Sigma$ added to the pressure. The regularization term $(\psi_{\alpha} - \psi_0) / \alpha$ is equal to the negative Shannon entropy $\int \log(\rho) \rho \mathrm{d} x$ of the mass density. In fact, its dual geodesics are the only notions of straight line that continue indefinitely (geodeiscally completen), are Euclidean in the right coordinate system (flat), and are invariant under sufficient statistics of the density (Markov invariant). Information geometric regularization augments the Euclidean geometry of individual particles with the information geometry of particle distributions. It thus accounts for the dual nature of $\rho$ as being a statistical estimate of physical particles.

IGR in practice

Shock capturing without viscosity and limiters

In practice, information geometric regularization allows dispensing with sophisticated shock capturing approaches of the path and instead deploy high-order linear numerics to the regularized equation. For instance, we can take a finite volume method with fifth-order accurate WENO reconstruction and remove the WENO limiter, using just the fifth order reconstruction applicable to smooth functions. The inviscid nature of IGR allows capturing subtle fine-scale features even when using relatively coarse meshes, whereas the WENO limiter diffuses them.

IGR preserves fine-scale features better than limiters.

Improvements over limiters. IGR solutions are smooth and computable without limiters. We superimpose a 2d Riemann problem with fine-scale entropy fluctuation and show the density of a reference solution, a solution computed using fifth order WENO limiters, and a solution obtained from fifth order reconstruction applied to the IGR equation. Limiters diffuse fine-scale structures and Kelvin-Helmholtz instabilities. IGR preserves them.

IGR at scale

To conclude, I want to point out that IGR is readily applicable to large, three dimensional problems. In recent joint work with Benjamin Willfong, Anand Radhakrishnan, Henry Le Berre, Nikolaos Tselepidis, Benedikt Dorschner, Reuben Budiardja, Brian Cornille, Stephen Abbott, and Spencer Bryngelson we integrated IGR into the exascale-ready MFC Code to conduct the first-ever compressible finite volume simulation with more than 100 trillion grid points! IGR played a key role in this endeavor by simplifying shock treatment, and thus the overall numerics.

Simulating arrays of rocket motors with IGR. These simulations of large arrays of rocket motors allow estimating the amount of hot gas reflected to the base of the rocket due to plume interaction. They were conducted using IGR on supercomputers like OLCF Frontier and CSCS Alps.