Mathematics Behind Particle Flow Maps

Navier-Stokes Equation (Governing Equations for Velocity)

We begin with the incompressible Navier-Stokes equations, which govern the fluid motion.

\[ \left\{ \begin{aligned} &\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \Delta \mathbf{u} \\ &\nabla \cdot \mathbf{u} = 0 \end{aligned} \right. \]

The left-hand side of the momentum equation can also be written as the material derivative of the velocity field, denoted as \( \frac{D \mathbf{u}}{D t} \), which captures the total rate of change experienced by a fluid particle moving with the flow.

Governing Equations for Impulse

To understand Particle Flow Maps, we first introduce some gauge variables, i.e., the impulse \(\mathbf{m}\), and the vorticity \(\boldsymbol{\omega}\). The evolution of \(\mathbf{m}\) follows:

\[ \left\{ \begin{aligned} &\frac{D \mathbf{m}}{D t} = -(\nabla \mathbf{u})^T \mathbf{m} + \nu \Delta \mathbf{m} \\ &\mathbf{u} = \mathbf{m} - \nabla \left( \Delta^{-1} \left( \nabla \cdot \mathbf{m} \right) \right) \end{aligned} \right. \]

Governing Equations for Vorticity

We now introduce the vorticity variable \( \boldsymbol{\omega} = \nabla \times \mathbf{u} \), which characterizes local rotational motion of the fluid. Its evolution is governed by:

\[ \left\{ \begin{aligned} &\frac{D \boldsymbol{\omega}}{D t} = (\boldsymbol{\omega} \cdot \nabla)\mathbf{u} + \nu \Delta \boldsymbol{\omega} \\ &\mathbf{u} = -\nabla \times \Delta^{-1} \boldsymbol{\omega} \end{aligned} \right. \]

Flow Map Preliminaries

We consider a velocity field \( \mathbf{u}(\mathbf{x}, t) \) defined over the fluid domain \( \Omega \), which specifies the fluid velocity at each point \( \mathbf{x} \) and time \( t \). Let \( \mathbf{X} \in \Omega \) denote the initial position of a fluid particle at time \( t = 0 \).

The forward flow map \( \boldsymbol{\phi}(\cdot, t): \Omega \rightarrow \Omega \) describes the trajectory of the particle as it moves with the flow:

\[ \left\{ \begin{aligned} &\frac{\partial \boldsymbol{\phi}(\mathbf{X}, t)}{\partial t} = \mathbf{u}(\boldsymbol{\phi}(\mathbf{X}, t), t) \\ &\boldsymbol{\phi}(\mathbf{X}, 0) = \mathbf{X} \\ &\boldsymbol{\phi}(\mathbf{X}, t) = \mathbf{x} \end{aligned} \right. \]

This map tracks the position \( \mathbf{x} \) of a fluid particle that starts from \( \mathbf{X} \) at \( t = 0 \) and arrives at time \( t \). The corresponding backward flow map \( \boldsymbol{\psi}(\cdot, t): \Omega \rightarrow \Omega \) retrieves the original position from the current one:

\[ \left\{ \begin{aligned} &\boldsymbol{\psi}(\mathbf{x}, 0) = \mathbf{x} \\ &\boldsymbol{\psi}(\mathbf{x}, t) = \mathbf{X} \end{aligned} \right. \]

To capture how the flow map and its inverse deform infinitesimal regions in the fluid, we define their Jacobian matrices:

\[ \left\{ \begin{aligned} &\mathcal{F}(\mathbf{X}, t) = \frac{\partial \boldsymbol{\phi}(\mathbf{X}, t)}{\partial \mathbf{X}} \\ &\mathcal{T}(\mathbf{x}, t) = \frac{\partial \boldsymbol{\psi}(\mathbf{x}, t)}{\partial \mathbf{x}} \end{aligned} \right. \]

Under the Lagrangian framework, the evolution of these Jacobian matrices follows:

\[ \left\{ \begin{aligned} &\frac{D \mathcal{F}}{D t} = \nabla \mathbf{u} \cdot \mathcal{F} \\ &\frac{D \mathcal{T}}{D t} = -\mathcal{T} \cdot \nabla \mathbf{u} \end{aligned} \right. \]

For further details, we refer readers to the works of Cortez (1995).

Physical Model of Particle Flow Maps

The Particle Flow Map (PFM) method leverages the observation that the trajectory of a particle naturally defines a perfect flow map. This insight enables the construction of accurate, bidirectional flow maps using particles as carriers.

Given a particle trajectory spanning from time \( a \) to \( c \), any intermediate time \( t \in [a, c] \) can be treated as the origin of a flow map. The destination is fixed at time \( c \), regardless of which intermediate \( t \) is chosen.

During advection, particles carry various physical quantities—including \( \boldsymbol{\omega} \), \( \mathbf{m} \), \( \nabla \boldsymbol{\omega} \), and \( \nabla \mathbf{m} \)—as well as flow map information like \( \mathcal{F} \), \( \mathcal{T} \), \( \nabla \mathcal{F} \), and \( \nabla \mathcal{T} \). These quantities evolve along the particle trajectory, enabling the mapping of fluid variables from time \( t \) to \( c \).

Unlike earlier approaches that estimate stretching effects using finite difference approximations, the PFM method directly computes stretching terms using the Jacobians of the flow maps, offering more accurate and stable modeling of fluid transport from the initial state.

Impulse Transport via Particle Flow Maps

As shown in the works of Cortez (1995), the evolution of the impulse variable \( \mathbf{m} \) can be expressed using the backward flow map Jacobian \( \mathcal{T}_t \):

\[ \mathbf{m}(\mathbf{x}, t) = \mathcal{T}_t^\top(\mathbf{x}) \, \mathbf{m}(\boldsymbol{\psi}(\mathbf{x}), 0) \]

This equation maps the initial impulse at the material point \( \boldsymbol{\psi}(\mathbf{x}) \) to the current location \( \mathbf{x} \), using the transpose of the backward flow map Jacobian. The Jacobian \( \mathcal{T}_t(\mathbf{x}) \) is defined as \( \frac{\partial \boldsymbol{\psi}(\mathbf{x}, t)}{\partial \mathbf{x}} \).

For particle-to-grid (P2G) transfers later, it is also necessary to evaluate the gradient of the impulse variable, \( \nabla \mathbf{m} \), evolves as:

\[ \nabla \mathbf{m}(\mathbf{x}, t) = \mathcal{T}_t^\top \, \nabla_{\boldsymbol{\psi}} \mathbf{m}(\boldsymbol{\psi}(\mathbf{x}), 0) \, \mathcal{T}_t + \nabla \mathcal{T}_t^\top \, \mathbf{m}(\boldsymbol{\psi}(\mathbf{x}), 0) \]

Here, \( \nabla_{\boldsymbol{\psi}} \) denotes the gradient with respect to the initial coordinates, and \( \nabla \mathcal{T}_t \) refers to the backward flow map Hessian.

Vorticity Transport via Particle Flow Maps

Similarly, the evolution of vorticity on particle flow maps can be described analytically based on the formulation by Cortez (1995). The flow map Jacobian \( \mathcal{F}_t \) propagates the initial vorticity field \( \boldsymbol{\omega}_0 \) to any later time \( t \) as follows:

\[ \boldsymbol{\omega}(\mathbf{x}, t) = \mathcal{F}_t(\mathbf{x}) \, \boldsymbol{\omega}(\boldsymbol{\psi}(\mathbf{x}), 0) \]

This expression indicates that the vorticity at a current spatial position \( \mathbf{x} \) and time \( t \) is determined by transforming the initial vorticity at the corresponding material point \( \boldsymbol{\psi}(\mathbf{x}) \) through the forward flow map Jacobian \( \mathcal{F}_t \).

For particle-to-grid (P2G) transfers later, it is also necessary to evaluate the gradient of vorticity, \( \nabla \boldsymbol{\omega} \). Taking spatial gradients of the above expression, we obtain the following evolution equation for \( \nabla \boldsymbol{\omega} \) under the flow map framework:

\[ \nabla \boldsymbol{\omega}(\mathbf{x}, t) = \mathcal{F}_t \, \nabla_{\boldsymbol{\psi}} \boldsymbol{\omega}(\boldsymbol{\psi}(\mathbf{x}), 0) \, \mathcal{T}_t + \nabla \mathcal{F}_t \, \boldsymbol{\omega}(\boldsymbol{\psi}(\mathbf{x}), 0) \]

This relation captures how both the Jacobian and its gradient contribute to the transport of spatial vorticity variation during particle advection.

P2G Transfer

After advecting the impulse \( \mathbf{m} \) and vorticity \( \boldsymbol{\omega} \) fields from their initial states using particle flow maps, we interpolate these quantities from particles back to the Eulerian grid in order to compute the grid-based velocity field.

For impulse, once the convected impulse \( \mathbf{m}_c \) and its gradient \( \nabla \mathbf{m}_c \) are updated on particles, we perform a particle-to-grid (P2G) transfer to obtain the grid impulse field \( \mathbf{m}_i \). The transfer is performed via APIC:

\[ \mathbf{m}_i \gets \frac{ \displaystyle \sum_p w_{ip} \left( \mathbf{m}_c^p + \nabla \mathbf{m}_c^p \cdot (\mathbf{x}_i - \mathbf{x}_p) \right) }{ \displaystyle \sum_p w_{ip} } \]

This grid-based impulse field \( \mathbf{m}_i \) is then used to solve a Poisson equation to remove its divergence component, ultimately yielding the incompressible velocity field \( \mathbf{u}_i \) on the grid.

Similarly, the vorticity field \( \boldsymbol{\omega} \) and its gradient \( \nabla \boldsymbol{\omega} \) are transferred from particles to grid nodes to compute the grid vorticity \( \boldsymbol{\omega}_{i} \) as:

\[ \boldsymbol{\omega}_{i} = \frac{ \displaystyle \sum_p w_{ip} \left( \boldsymbol{\omega}^p_c + \nabla \boldsymbol{\omega}^p_c \cdot (\mathbf{x}_{i} - \mathbf{x}_p) \right) }{ \displaystyle \sum_p w_{ip} } \]

This grid-based vorticity field \( \boldsymbol{\omega}_i \) is then used to reconstruct the velocity field on the grid.