###### 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,\]the gradient will be

\[\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.

##### Least Squares

TODO

comments powered by Disqus