Piecewise Constant Velocity Advection

I feel it’s more natural to express position of vortices in barycentric coordinates along with face-id. Then the advection is to be carried on barycentric coordinates rather that the world coordinate of the vortex. Let $\vec{\lambda} = (\lambda_{1}, \lambda_{2}, \lambda_{3})$ where each $\lambda_{i}(\vec{x}, f): \mathbb{R}^3 \times F \to \mathbb{R}$ and $f \in F$ is a face-id.

Then from stream function formulation we know that

$$\vec{u} = J\nabla \psi = \dot{x} \tag 1$$

If we want to perform the advection using barycentric coordinates, we need $\dot{\lambda}$. Further, using multivariate chain-rule we can write (we assume that $f$, the face-id is constant for a time-step)

$$\dot{\lambda_{i}}= \nabla \lambda_{i} \cdot \dot{x} = \nabla \lambda_{i}\cdot \vec{u} \tag 2$$

Now from $(1)$ and $(2)$

$$\dot{\lambda_{i}} = \nabla \lambda_{i} \cdot \vec{u} = \nabla \lambda_{i} \cdot (J \nabla \psi) = \nabla \lambda_{i} \cdot (\hat{n}_{f} \times \nabla \psi) \tag 3$$

where $\hat{n}_{f}$ is a unit face normal for the incident face. We can rearrange $(3)$ using vector triple product identity as:

$$\dot{\lambda_{i}} = \hat{n}_{f} \cdot (\nabla \psi \times \nabla \lambda_{i}) \tag 4$$

Further, the stream function can be written down in terms of barycentric coordinates $\psi = \psi_{1}\lambda_{1} + \psi_{2}\lambda_{2} + \psi_{3}\lambda_{3}$, where $\psi_{i}$ is value of the stream function at $v_{i}$ in face-id $f$. Then,

$$\nabla \psi = \psi_{1}\nabla \lambda_{1} + \psi_{2}\nabla \lambda_{2} + \psi_{3}\nabla \lambda_{3}$$

using this in $(4)$ we get

$$\dot{\lambda_{i}} = \hat{n}_{f} \cdot ((\psi_{1}\nabla \lambda_{1} \times \nabla \lambda_{i}) + (\psi_{2}\nabla \lambda_{2}\times \nabla \lambda_{i}) + (\psi_{3}\nabla \lambda_{3} \times \nabla \lambda_{i})) \tag 5$$

Also, $\nabla \lambda_{i}$ is face-wise constant which is given as:

$$\nabla \lambda_{i} = \frac{\hat{n}_{f} \times \vec{e}_{jk}}{2A_{f}}$$

where $A_{f}$ is area of triangle with face-id $f$, and $\vec{e}_{jk}$ is edge opposite to the vertex $x_{i}$. Using this in Eq. $(5)$ results in lot of cancellation and simplification and finally we get:

$$\begin{align} \dot{\lambda_{1}} = \frac{-\psi_{2} + \psi_{3}}{2A_{f}} \\ \dot{\lambda_{2}} = \frac{\psi_{1} - \psi_{3}}{2A_{f}} \\ \dot{\lambda_{3}} = \frac{-\psi_{1} + \psi_{2}}{2A_{f}} \end{align} \tag 6$$

Note that $\sum_{i=1}^3 \lambda_{i} = 1$, i.e., $\sum_{i=1}^3 \dot{\lambda}_{i} = 0$, and we can verify that Eq. $(6)$ satisfy this constraint.

Edge Crossing

In the reformulation above, we assumed that the face ID $f$ is constant. So what happens when a point vortex crosses an edge? This can be handled by checking for a negative barycentric weight. In that case, we can use a discrete Levi–Civita connection, unfold-translate-rotate-fold adjacent triangles to move the point from one face to another.

Removing Influence of Singular Self Term

After implementing above, the point vortex’s motion was different that expected. After further analysis it turned out to be cause of the self induced velocity due to singular part of a vortex on itself. Looking back a the analytical Green’s Function formulation, this is obvious as well.

$$\vec{u}_{i} = \sum_{j=1 (i \neq j)}^{k} \Gamma_{j} \nabla^{\perp}G(x_{i}, x_{j}) \tag 7$$

Where, $\vec{u}_{i}$ is velocity at the $i^{th}$ point vortex and $G(\cdot, \cdot)$ is Green’s kernel, solution to Poisson equation for delta functions on RHS. Notice the $i \neq j$ part of the sum, this implies that we should exclude a vortex’s influence on itself.

This exclusion seems to be motivated from the fact that evaluation of the Green’s function at the core results in a singular term with infinite velocity. We get rid if this singularity by a barycentric splat of the point vortex to its incident triangle vertices. But after implementing this, it was evident in the simulation that the split of singular part is influencing motion of a point vortex inside triangle. By observation it looked like, after splitting, each triangle vertex acted as a point vortex. This induced a finite, nonzero velocity at the original location of the point vortex. This affected the time evolution and resulted in wrong vortex trajectories.

First fix, I treated each triangle vertex a point vortex with strength $\Gamma \lambda_{i}$. As triangle is a planer domain, I used 2D planer Biot-Sawart Law to calculate the velocity at the original vortex position $\vec{u}_{\text{self}}$. Then the final velocity at the point vortex location is $\vec{u}_{\text{global}} - \vec{u}_{\text{self}}$.

But our mesh domain is a linear approximation and by above approach, we were using wrong model to calculate the $\vec{u}_{self}$. We need to use the same DEC machinery to calculate the self velocity term. One thing to notice is that point vortex’s self velocity (due to the singular part) is a local effect, so we can take a local patch of mesh around the vortex triangle, formulate a Poisson problem with Dirichlet boundary condition. Then the velocity inside the source face would be a true representation for $\vec{u}_{self}$.

But as you could see this is numerically very expensive as we need to setup and solve a mini Poisson problem for each face for every vortex position! But after a lot of pondering I realized, we can actually precompute velocity basis for self velocity per face. Then when needed, just take weighted average of these basis using barycentric weights and vortex strength to get the $\vec{u}_{self}$. This cuts down computation from multiple time per face to once per face. Further, we can do lazy computation to avoid calculating all basis at once.

Below is a detailed sketch of how to get these self velocity basis

Then, given $U_{i}^{f}$ we get self velocity as:

$$\vec{u}_{self} = \Gamma \lambda_{1}U_{1}^{f} + \Gamma \lambda_{2}U_{2}^{f} + \Gamma \lambda_{3}U_{3}^{f} \tag 8$$

Using this singular term removal I was able to get at least a sane behavior, but the path taken by point vortex still doesn’t look natural. This seems be advection artifacts, due to discontinuous velocity field, as each Poisson solve just assumes that the current state is an accurate evolution of previous. But due to artifacts in the mesh, around corners and all the vortice’s relative motion is not accurate enough which is corrupting the next step.

Hamiltonian of the system is fluctuating in range of 0.1 but then starts to grow (mostly due to numerical error) after ~5000 steps.