Vortex Particle Flow Map (VPFM) Simulations

This repository implements the Hybrid Eulerian-Lagrangian Vortex Particle Flow Maps method (arXiv:2505.21946). It is a high-performance fluid solver designed to capture intricate turbulent details with zero numerical dissipation.

🚀 Efficiency & Performance Principles

Traditional fluid solvers face a dilemma:

• Eulerian (Grid-based): Fast $O(N)$, but suffers from Numerical Dissipation. Vortices blur and die out over time due to advection errors.
• Lagrangian (Vortex Methods): Preserves vorticity perfectly, but calculating velocity requires the Biot-Savart Law, which is Extremely Slow $O(N^2)$.

The VPFM Solution

This Hybrid method combines the best of both worlds to achieve Real-Time Performance on consumer GPUs:

1. Particle-in-Cell (PIC) Hybrid: We carry vorticity ($\omega$) on particles to prevent dissipation, but we transfer it to a grid to solve for velocity.
   • Grid Solver: Solving the Poisson equation $\nabla^2 \psi = -\omega$ is $O(N \log N)$ or even $O(N)$ with Multigrid/FFT methods, compared to $O(N^2)$ for direct interaction.
2. GPU Acceleration (Taichi): The entire pipeline fits in GPU memory. Memory access is coalesced using Taichi's localized data structures (ti.root.dense).
3. Sparse Computation: We only simualte where the vorticity is. Ambient space costs almost nothing.

---

🛠️ Implementation Tricks & Methods

1. The Flow Map Jacobian ($F$)

Instead of just advecting vorticity, we track the deformation of space itself. Each particle carries a tensor $F$ (Flow Map Jacobian) and $H$ (Hessian).

• The Trick: $\omega(t) = F \cdot \omega_0$.
• By evolving $F$ rather than $\omega$ directly, we capture the stretching and tilting of vortex tubes analytically. This allows us to resolve sub-grid features that would otherwise be lost.

2. P2G / G2P Operations

Transferring data between particles and grid is the source of most errors.

• Method: We use Quadratic B-Spline Kernels (spanning 3x3x3 cells) for transfer.
• Normalization: A critical fix implemented here is Weight Normalization. $\omega_{grid} = \sum (\omega_p \cdot W) / \sum W$. Without this, particles bunching up would artificially create "super-vortices" and explode the simulation.

3. Reinitialization Cycles

The Flow Map $F$ eventually becomes too distorted.

• The Trick: We periodically "bake" the current vorticity state into the particles and reset $F$ to the Identity Matrix ($I$).
• Short/Long Cycles: We use two frequencies. A short cycle ($N_S \approx 5$) updates the gradient, and a long cycle ($N_L \approx 10$) resets the map entirely. This keeps the Jacobian well-conditioned.

4. Boundary Repulsion (Artificial)

Pure VPFM doesn't handle wall collisions well (particles sticks).

• Constraint: We implement a 5-cell "Repulsion Zone" near boundaries. Any particle entering this zone receives an immediate velocity impulse normal to the wall, simulating pressure rebound without solving a complex pressure projection.

---

📂 Project Structure

File	Description
vortex_chaos_3d.py	3D Simulation of 4 interacting vortex filaments (Instability/Chaos).
vortex_collision_3d.py	3D Simulation of two colliding vortex rings (Annihilation).
vortex_separation_3d.py	3D Simulation of two repelling vortex rings.
vortex_trefoil_3d.py	3D Simulation of a single Trefoil Knot vortex.
vortex_flow_map_3d.py	3D Base Simulation (Coaxial Leapfrogging Rings).
vortex_flow_map_2d.py	2D Simulation (Kelvin-Helmholtz Instability & Vortex Dipoles).
vortex_gravity_2d.py	2D Simulation with Gravity (Dam Break / Rayleigh-Taylor).
docs/	Contains GIF visualizations for the README.

---

📂 Simulation Details

1. vortex_chaos_3d.py (Vortex Tangle)

Scenario: Pure turbulence study.

• Logic: Initializes 4 parallel vertical vortex tubes. A sinusoidal deviation is applied to their positions ($x + \sin(y)$).
• Dynamics: The slight perturbation triggers the Crow Instability. The tubes induce velocity on each other, braiding together before reconnecting in a violent topological change.
• Visual: Uses a high-pass filter (w_mag > 2.0) to render only the "bones" of the turbulence.

2. vortex_collision_3d.py (Ring Collision)

Scenario: Head-on collision of two vortex rings.

• Logic:
  • Ring 1: $\Gamma = +100$, moving UP.
  • Ring 2: $\Gamma = -100$, moving DOWN.
• Outcome: As they approach, the induced velocity causes them to stretch radially outward. They collide and annihilate/reconnect into a series of smaller secondary rings.
• Color: Rainbow coloring helps identify which parts of the wreckage came from which original ring.

3. vortex_separation_3d.py (Repulsion)

Scenario: Interaction of opposing spins.

• Logic: Two rings placed coaxially.
  • Ring 1: Clockwise spin.
  • Ring 2: Clockwise spin (relative to viewing plane).
• Outcome: Unlike collision, the induced velocity field here points away from the center for both rings. They act like magnets with matching poles, slowing down and reversing direction.

4. vortex_trefoil_3d.py (Trefoil Knot)

Scenario: Topological complexity.

• Logic: Particles are seeded along a parametric curve:
  • $x = \sin(t) + 2\sin(2t)$
  • $y = \cos(t) - 2\cos(2t)$
  • $z = -\sin(3t)$
• Dynamics: The curvature induces self-rotation. The knot rotates as a rigid body initially but eventually breaks down due to self-intersection.

5. vortex_flow_map_3d.py (Leapfrog)

Scenario: The classic textbook example.

• Logic: Two identical rings moving in the same direction.
• Outcome: The rear ring is sucked into the wake of the front ring. It shrinks (speeds up) and passes through the front ring. The front ring expands (slows down) and wraps around the back one. This cycle repeats.

6. vortex_flow_map_2d.py (Two Vortices)

Scenario: Interaction of a Vortex Dipole.

• Logic: Use the "Two Vortices" button in vortex_flow_map_2d.py.
• Setup: Two Gaussian vortices with opposite circulation ($\Gamma = \pm 400$) placed side-by-side.
• Outcome: The induced velocity field effectively couples them, forming a stable Vortex Dipole that propagates linearly through the fluid.
• Visual: Gradient Dye (Coordinate-based coloring).

7. vortex_flow_map_2d.py (Kelvin-Helmholtz Instability)

Scenario: Shear layer instability.

• Logic: The simulation initializes a classic shear flow setup:
  1. Velocity Field: Two horizontal bands of fluid moving in opposite directions ($U_{top} = -1, U_{bottom} = +1$).
  2. Perturbation: A small sinusoidal vertical displacement is applied to the interface particles.
  3. Dynamics: The VPFM solver captures the linear growth of these perturbations. As the shear intensifies, the accumulation of vorticity causes the interface to roll up into stable, rotating structures known as "Kelvin-Helmholtz Cats' Eyes".
• Visual: We use Passive Dye Tracers. A separate set of particles carries a color field (White for top, Blue for bottom) but no vorticity. These tracers are advected by the velocity field to visualize the mixing interface without affecting the physics.

7. vortex_gravity_2d.py (Dam Break / Rayleigh-Taylor)

Scenario: Gravity-driven collapse of a heavy fluid column.

• Logic: This script extends the standard VPFM with Variable Density Physics.
  1. Baroclinic Torque: In a constant density fluid, gravity does not generate vorticity. Here, we implement the full Baroclinic term: $\frac{D\omega}{Dt} = \frac{1}{\rho^2} (\nabla \rho \times \nabla p)$.
  2. Implementation: We approximate this by generating vorticity wherever the density gradient ($\nabla \rho$) is perpendicular to gravity ($g$).
  3. Setup: A rectangular block of High-Density particles ($\rho=2$) is initialized in a Low-Density medium ($\rho=1$).
• Outcome: The misalignment of density and pressure gradients at the vertical wall of the water column generates huge amounts of vorticity. This vorticity drives the "splashing" motion, causing the water to slump, hit the wall, and roll up into turbulent waves.

---

References

• Paper: Hybrid Eulerian-Lagrangian Vortex Particle Flow Maps