# Perspectives on Incompressibility

10 Jan 2017

###### Notation

We primarily denote derivatives by some variable primarily through subscripts

$\frac{\partial f}{\partial t} = f_t$

but because temporal derivatives are so common we will also denote temporal derivatives by dot notation

$\dot f = f_t.$

$q$ will be used to represent Euclidean positions and $x$ will repsresnt general coordinates. $q_t$ and $u$ will both be the velocity field of particles $q$.

### Incompressibility

For a diffeomorphism of a domain $\Omega$ by some function $\Phi(q;t) = \Phi^t: \Omega: \mathbb{R}^n \times \mathbb{R} \rightarrow \mathbb{R}^n$ we define incompressibility as

$\| \Phi^t(S)\| = \|S\| \forall S \subset \Omega$

for volume defined in the usual way by

$\|S\| = \int_{S} dq.$

The volume of $S$ after being advected is given by

$\| \Phi^t(S)\| = \int_{\Phi^t(S)} dq = \int_{S} d(\Phi^t(q)) = \int_{S} \nabla \cdot \Phi^t(q) dq.$

Obviously if incompressibility is maintained

$0 = \frac{d}{dt}\| \Phi^t(S)\| = \frac{d}{dt}\int_{S} \nabla \cdot \Phi^t(q) dq. = \int_{S} \nabla \cdot \Phi^t(q) dq.$

The differential form of this statement is that

$\nabla \cdot \Phi^t(q;t) = 0\forall q \in \Omega$.

Usually one starts with some PDE

$q^t_t = f(q^t, q^t_t,...)$

and incompressibility is added by the gradient of some scalar field $p$ that guarantees $\nabla q_t = 0$:

$q^t_t = f(q^t, q^t_t,...) + \nabla p$

There are a few perspectives for what pressure $p$ is that are valuable to understand:

#### $L^2$ Projection

This is perhaps the most fundamental geometric perspective for incompressibility. Ignoring the fancy notation, all this really depends on is the fact that differentiation is a linear map and that there’s a nice orthogonal decomposition of arbitrary (but sufficiently smooth) vector fields using different differential operators. What this boils down to is that making a vector field incompressible is the removal of the part of it that has divergence, which is a projection.

Recall the Hodge-deRham decomposition of a $1$-form into three orthogonal components

$\omega = d\alpha + \delta \beta + \gamma$ for $\alpha$ a 0-form, $\beta$ a 2-form, and $\gamma$ harmonic.

Now note that $\delta \delta = 0$, $\delta \gamma = 0$, and $\left(\delta u^{\flat}\right)^{\sharp} = \nabla \cdot u$.

Therefore, so long as for the decomposition

$u^{\flat} = d\alpha + \delta \beta + \gamma,$

$\left(\delta u^{\flat} \right)^{\sharp}= \left(\delta d\alpha + \delta \delta \beta + \delta \gamma\right)^{\sharp} = \left(\delta d\alpha \right)^{\sharp}$

so long as $\delta d\alpha = 0$ the vector field is divergence free. By the definition of harmonic vector fields we see that $\alpha$ is orthogonal to the kernel of $\delta d$ $\alpha$ and so being incompressible is simply setting $\alpha = 0$. The required pressure is therefore just $-(d\alpha)^\sharp$ and pressure projection is

$\left(u^{\flat} - d\alpha\right)^{\sharp}$

where $\alpha$ can be computed by

$\delta u^{\flat} = \delta d\alpha = \Delta \alpha.$

#### Variational energy minimization

The variational perspective for pressure is that its the vector field that its the arugment for the minimal kinetic energy necessary to be divergence-free:

$p = \arg\min_p\int_\Omega \| u - \nabla p \|^2 dx.$

This idea comes from orthogonoality in the sense that it’s finding a $p$ that removes as much $\nabla \cdot$ from $u$ as possible. The key advantage of this formulation is that it does not depend on splitting as projection does. That is, the previous formulation is more amenable for solving for $u$ and $p$ simultaneously (as well as other things if there is a need).

By doing some arithmetic we see

$\frac{1}{2}\frac{\partial}{\partial p} \left(\int_\Omega \|u\|^2 + \|\nabla p\|^2 - 2u\cdot \nabla p dx\right)$ $= \int_\Omega \nabla^T \nabla p - \nabla \cdot u dx$

#### Augmented Lagrangian

Consider the Lagrangian mechanics perspective of dynamics where the general positions and general velocities $x,\dot x$ are defined by minimizing some Lagrangian $\mathcal{L}(x,\dot x)$ and the Euler-Lagrange equations

$\mathcal{L}_x(x,\dot x) - \frac{d}{dt}\mathcal{L}_{\dot x}(x,\dot x) = 0.$

Here we assume that $x$ represents positions and densities in a piece of geometry, which implies that $\dot x$ represents a diffeomorphism of the geometry itself.

Incompressibility adding the constraint $\mathbf{D}\dot x = 0$ for some linaer operator $\mathbf{D}$ that satisfies

$\int_{\omega} \mathbf{D}\dot x \approx \int_{\omega} \nabla \cdot \frac{d\Phi_x}{dt} dq$

where $\Phi_x$ is the diffeomorphism induced by the general coordinates $x$.

###### Adjoint of $\mathbf{D}$

Note that the adjoint operator, given the condition that $f=0$ or $\frac{\partial \Phi_x}{\partial t}\cdot N = 0$ for any point along the boundary, can be defined by

$\int_{\omega} f \mathbf{D}\dot x \approx \int_{\omega} f \nabla \cdot \frac{d\Phi_x}{dt} dq = \int_{\omega} \frac{d\Phi_x}{dt} \cdot \nabla f dq \approx \int_{\omega}\dot x \mathbf{D}^* f.$
##### Augmentation

Augmenting this Lagrangian for incompressibility is therefore done by

$\bar {\mathcal{L}}(x,\dot x) = \mathcal{L}(x,\dot x) + \int_\Omega \lambda (\mathbf{D}\dot x - 0) dx.$

Necessary conditions for optimality dictate that

$\bar {\mathcal{L}}_{x}(x,\dot x) = \mathcal{L}_{x}(x,\dot x) = 0,$ $\bar {\mathcal{L}}_{\dot x}(x,\dot x) = \mathcal{L}_{\dot x}(x,\dot x) + \int_\Omega \left[\lambda \mathbf{D}\right]^*dq = 0,$ $\bar {\mathcal{L}}_\lambda(x,\dot x) = \int_\Omega \left(\mathbf{D}\dot x - 0\right) dq = 0,$

where the final condition is simply a restating of the incompresibility condition.

On hte other hand we can the adjoint operator from before

$\int_\Omega \left[\lambda \mathbf{D}\right]^*dq \approx \int_{\Omega} \nabla \lambda dq.$

This $\lambda$ parameter is what is commonly called pressure.

###### An example with viscous fluids

Recall that incompressible Navier-Stokes can be written as

$\frac{\partial u}{\partial t} + u \cdot \nabla u - \mu \Delta u = \nabla p + g$

where $p$ is such that

$\nabla \cdot u = 0.$

If we ignore the term $u \cdot \nabla u$ (which can be dealt with by an ODE integrator) and solve for the future $u$ implicitly we obtain

$\frac{\partial u}{\partial t} - \mu \Delta \left(u + dt \frac{\partial u}{\partial t}\right) = \nabla p + g$ $\left(\mathbf{I} - dt\mu \Delta\right) \frac{\partial u}{\partial t} - \nabla p = g + \mu \Delta u$ $\begin{bmatrix} I - dt\mu \Delta & -\nabla\\ -\nabla & 0 \end{bmatrix} \begin{bmatrix} \frac{\partial u}{\partial t}\\ p \end{bmatrix} \begin{bmatrix} g + \mu \Delta u\\ 0 \end{bmatrix}$

Where the first and second rows represent $\frac{\partial \bar{\mathcal{L}}}{\partial u}$ and $\frac{\partial \bar{\mathcal{L}}}{\partial p}$ respectively for

$\bar {\mathcal{L}}(u,p) = \int_{\Omega}\frac{1}{2} \left[{\frac{\partial u}{\partial t}}\right]^*\left(I - dt\mu \Delta\right)\left[\frac{\partial u}{\partial t}\right] - p \nabla \cdot \frac{\partial u}{\partial t} dx$

which, is the augmented lagrangian for the energy

$\mathcal{L}(u,p) = \int_\Omega\frac{1}{2} \left[{\frac{\partial u}{\partial t}}\right]^*\left(I - dt\mu \Delta\right)\left[\frac{\partial u}{\partial t}\right] dx$

By applying Shur’s complement method to solve for $p$ one obtains from the above matrix

$0 - \nabla \cdot (I - dt\mu \Delta )^{-1} \nabla p = \nabla \cdot -(I - dt\mu \Delta )^{-1} (g + \mu \Delta u)$

which simplifies to

$\nabla \cdot (I - dt\mu \Delta )^{-1} \nabla p = \nabla \cdot (I - dt\mu \Delta )^{-1} (g + \mu \Delta u).$

In the inviscid case, or when we split viscosity and lump it into $g$, this naturally turns into

$\nabla \cdot \nabla p = \nabla \cdot g,$

our standard Poisson problem.

TODO