Scientist: denario-6 Date: 2026-05-08
This dataset contains 1001 snapshots from a 3D isothermal hydrodynamic driven turbulence simulation run with AthenaK. The scientific goal is to identify and track vortices across time, compute the trajectory of each vortex's center of vorticity, and determine whether vortex motion follows a Gaussian random walk or a Lévy flight.
All VTK snapshots are located at:
/home/node/work/projects/ns_turbulence_vortex/data/Turb.hydro_w.NNNNN.vtk
where NNNNN ranges from 18903 to 19903 (inclusive), giving 1001 time snapshots.
Download status: files are being downloaded from HuggingFace (pedrota2000/NS_simulation). Before analysis, check how many files are available:
import os, glob
files = sorted(glob.glob('/home/node/work/projects/ns_turbulence_vortex/data/Turb.hydro_w.*.vtk'))
print(f"Available files: {len(files)}, indices {int(files[0].split('.')[-2])} to {int(files[-1].split('.')[-2])}")Use however many files are available at analysis time (at least ~300 are likely ready).
- Format: Legacy VTK binary (DATASET STRUCTURED_POINTS), readable with
pyvista - Grid: 128×128×128 cells (Fortran/C ordering: shape is (128, 128, 128) after reshape)
- Domain: [-0.5, 0.5]³ (unit periodic box)
- Cell spacing: Δx = Δy = Δz = 1/128 ≈ 0.0078125
- Origin: (-0.5, -0.5, -0.5)
Reading example:
import pyvista as pv
import numpy as np
mesh = pv.read('/home/node/work/projects/ns_turbulence_vortex/data/Turb.hydro_w.18903.vtk')
dens = mesh['dens'].reshape(128, 128, 128)
velx = mesh['velx'].reshape(128, 128, 128)
vely = mesh['vely'].reshape(128, 128, 128)
velz = mesh['velz'].reshape(128, 128, 128)| Field | Description | Observed range (t=189.03) |
|---|---|---|
| dens | Mass density ρ | [0.985, 1.007] |
| velx | Velocity component v_x | [-0.865, +0.691] |
| vely | Velocity component v_y | [-0.687, +0.717] |
| velz | Velocity component v_z | [-0.795, +0.853] |
| s_00 | Passive scalar tracer (uniform) | ~0.047 (constant) |
The density is nearly uniform (isothermal EOS), so dynamical structure is encoded in velocity.
| Parameter | Value |
|---|---|
| Code | AthenaK |
| Problem | 3D driven turbulence |
| Grid resolution | 128 × 128 × 128 |
| Domain | [-0.5, 0.5]³ (periodic) |
| EOS | Isothermal, sound speed c_s = 5.0 |
| Time integrator | RK2, CFL = 0.3 |
| Spatial reconstruction | PLM (piecewise linear) |
| Riemann solver | HLLE |
| Driving type | Continuously forced, solenoidal |
| Energy injection rate | dedt = 1×10⁻⁴ |
| Correlation time | tcorr = 5.0 |
| Driving wavenumber | nlow=1, nhigh=3 (peak at n=2) |
| Sol fraction | 1.0 (purely solenoidal) |
- Output cadence: Δt = 0.01 simulation time units
- Time range: t ≈ 189.03 to 199.03
- File index encoding: output_index = NNNNN; sim_time ≈ NNNNN / 100
Vorticity ω = ∇ × v has three components:
dx = 1.0 / 128
omega_x = np.gradient(velz, dx, axis=1) - np.gradient(vely, dx, axis=2)
omega_y = np.gradient(velx, dx, axis=2) - np.gradient(velz, dx, axis=0)
omega_z = np.gradient(vely, dx, axis=0) - np.gradient(velx, dx, axis=1)
omega_mag = np.sqrt(omega_x**2 + omega_y**2 + omega_z**2)Observed vorticity magnitude: [0.048, 28.4], mean ≈ 6.4 (at t=189.03).
-
Vortex detection: At each time snapshot, identify 3D vortices using the Q-criterion, λ₂-criterion, or local vorticity magnitude thresholding. Q > 0 regions (where rotation dominates strain) are standard vortex cores. Alternatively, threshold on |ω| > threshold (e.g., mean + 2σ).
-
Vortex labeling: Connected-component labeling (scipy.ndimage.label or similar) on the vortex mask to identify individual vortices. Account for periodic boundary conditions.
-
Center of vorticity: For each labeled vortex, compute the vorticity-weighted centroid:
x_c = Σ (|ω(r)| · x) / Σ |ω(r)|(sum over voxels in the vortex region). This is the "center of vorticity."
-
Vortex tracking: Link vortices across consecutive time snapshots (e.g., nearest-neighbor in centroid space, Hungarian algorithm, or overlap-based matching). Build trajectories for each tracked vortex.
-
Statistical analysis of trajectories:
- Compute displacement increments Δr_i = r(t_{i+1}) - r(t_i) for each tracked vortex
- Test whether increments follow a Gaussian distribution (Gaussian random walk) or a heavy-tailed / power-law distribution (Lévy flight)
- Compute the mean squared displacement (MSD) as a function of lag time τ: MSD(τ) = <|r(t+τ) - r(t)|²>
- Gaussian RW: MSD ∝ τ (diffusive)
- Lévy flight: MSD ∝ τ^α with α > 1 (superdiffusive) or anomalous scaling
- Fit the step-size distribution P(Δr) to both Gaussian and Lévy-stable distributions; compare via AIC/BIC or KS test
- Compute the kurtosis of displacement increments (Gaussian: κ = 3; heavy-tailed: κ >> 3)
- Plot vortex trajectories in 3D and projected onto 2D planes
- Plot the displacement distribution and compare to Gaussian/Lévy fits
-
Produce many plots: vorticity field snapshots, vortex identification at multiple times, individual vortex trajectory plots, MSD vs lag time (log-log), step-size PDF with Gaussian and Lévy fits, kurtosis evolution, etc.
- The domain is periodic: vortex coordinates may wrap around the boundary. Use periodic distance metrics and handle centroid wrapping carefully (e.g., check if vortex spans the boundary and unwrap coordinates).
- Memory usage: loading all 1001 files at once (~42 GB) is impractical. Process files sequentially or in batches, keeping only the vortex centroid data (not full fields) in memory.
- Parallel processing: use up to 8 workers for VTK reading (I/O bound). For vorticity computation (CPU bound), use joblib or multiprocessing with 8–16 workers.
- The simulation is subsonic (Mach number ≈ velocity / c_s ≈ 0.8/5.0 ≈ 0.16), so compressibility effects are small.
- The turbulence is statistically stationary (driven) — no transient startup to worry about.
- Vortices in 3D turbulence are elongated filament-like structures (tubes), not point-like. The centroid is a meaningful summary statistic.