Tutorial: 3D Turbulence Simulation

This tutorial demonstrates setting up and running a 3D turbulent flow simulation using Tarang.jl.

Physical Problem

We simulate decaying turbulence in a triply-periodic box, governed by the incompressible Navier-Stokes equations:

\[\begin{aligned} \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} &= -\nabla p + \nu \nabla^2 \mathbf{u} \\ \nabla \cdot \mathbf{u} &= 0 \end{aligned}\]

where $\mathbf{u} = (u, v, w)$ is the velocity, $p$ the pressure (divided by density), and $\nu$ the kinematic viscosity.

Formulation: vorticity and vector potential

Tarang's pure-Fourier time-stepping path integrates each field with an implicit linear operator that is diagonal in Fourier space (a Laplacian, a hyper/fractional Laplacian, a constant damping, or a derivative of the field itself). The pressure-velocity saddle-point system is not of that form: an implicit $\nabla p$ couples a different field into the momentum equation, and the solver rejects it (under MPI it raises an error; in serial the coupled system is singular and the implicit solve is silently dropped, leaving $p = 0$ and a velocity that is no longer divergence-free).

In a triply-periodic box the pressure can be removed analytically, which is what we do here. Taking the curl of the momentum equation gives the vorticity equation, and the velocity is recovered from a vector potential $\mathbf{A}$ in the Coulomb gauge:

\[\begin{aligned} \frac{\partial \boldsymbol{\omega}}{\partial t} - \nu \nabla^2 \boldsymbol{\omega} &= -(\mathbf{u}\cdot\nabla)\boldsymbol{\omega} + (\boldsymbol{\omega}\cdot\nabla)\mathbf{u} \\ \nabla^2 \mathbf{A} &= -\boldsymbol{\omega}, \qquad \mathbf{u} = \nabla \times \mathbf{A} \end{aligned}\]

Every implicit term is now a Laplacian (diagonal), every nonlinear term is explicit and differentiated along Fourier axes only, and the two constraint equations are exactly the Poisson / algebraic-substitution pattern the solver refreshes each step.

A velocity built as $\mathbf{u} = \nabla\times\mathbf{A}$ is divergence-free whatever $\mathbf{A}$ is, so $\max|\nabla\cdot\mathbf{u}| \approx 8\times10^{-15}$ (measured on the $32^3$ run below) confirms nothing on its own. The check that actually tests the Poisson solve is the reconstructed velocity itself: after one step max|u1| = 0.999940, matching the analytic viscous decay of the Taylor-Green amplitude, $e^{-3\nu\,\Delta t} = 0.999940$, to seven digits.

Domain Setup

Coordinates and Distribution

using Tarang
using MPI
using Printf

MPI.Init()

coords = CartesianCoordinates("x", "y", "z")
dist   = Distributor(coords; dtype=Float64, device=CPU())

Leave mesh unset unless you have a reason: the distributor picks a valid process mesh for you ((1, 2) on 2 ranks, (2, 2) on 4 ranks).

Spectral Bases

All directions are periodic, so we use Fourier bases:

L  = 2π       # domain size
N  = 32       # N^3 grid points
Re = 100.0
nu = 1.0 / Re

x_basis = RealFourier(coords["x"]; size=N, bounds=(0.0, L), dealias=3/2)
y_basis = RealFourier(coords["y"]; size=N, bounds=(0.0, L), dealias=3/2)
z_basis = RealFourier(coords["z"]; size=N, bounds=(0.0, L), dealias=3/2)

domain = Domain(dist, (x_basis, y_basis, z_basis))

Fields and Problem Definition

The unknowns are the three vorticity components (time-stepped), the three vector-potential components and the three velocity components (both solved from constraints each step), plus one gauge constant per Poisson equation:

w1 = ScalarField(domain, "w1"); w2 = ScalarField(domain, "w2"); w3 = ScalarField(domain, "w3")
A1 = ScalarField(domain, "A1"); A2 = ScalarField(domain, "A2"); A3 = ScalarField(domain, "A3")
u1 = ScalarField(domain, "u1"); u2 = ScalarField(domain, "u2"); u3 = ScalarField(domain, "u3")

tau_A1 = ScalarField(dist, "tau_A1", (), Float64)
tau_A2 = ScalarField(dist, "tau_A2", (), Float64)
tau_A3 = ScalarField(dist, "tau_A3", (), Float64)

problem = IVP([w1, w2, w3, A1, A2, A3, u1, u2, u3, tau_A1, tau_A2, tau_A3])
add_parameters!(problem, nu=nu)

# Vorticity transport: advection + vortex stretching on the explicit side,
# viscous diffusion on the implicit side.
for (w, u) in (("w1", "u1"), ("w2", "u2"), ("w3", "u3"))
    add_equation!(problem,
        "∂t($w) - nu*Δ($w) = -(u1*∂x($w) + u2*∂y($w) + u3*∂z($w))" *
        " + (w1*∂x($u) + w2*∂y($u) + w3*∂z($u))")
end

# Vector potential: Δ(A) = -ω, one Poisson equation per component.
add_equation!(problem, "Δ(A1) + tau_A1 + w1 = 0")
add_equation!(problem, "Δ(A2) + tau_A2 + w2 = 0")
add_equation!(problem, "Δ(A3) + tau_A3 + w3 = 0")

# Velocity: u = ∇ × A.
add_equation!(problem, "u1 - (∂y(A3) - ∂z(A2)) = 0")
add_equation!(problem, "u2 - (∂z(A1) - ∂x(A3)) = 0")
add_equation!(problem, "u3 - (∂x(A2) - ∂y(A1)) = 0")

# Gauge: fix the mean of each vector-potential component.
add_bc!(problem, "integ(A1) = 0")
add_bc!(problem, "integ(A2) = 0")
add_bc!(problem, "integ(A3) = 0")

Solver Setup

solver = InitialValueSolver(problem, RK443(); dt=2e-3)

# The RHS compiles to a type-specialized plan (no interpreted fallback):
@assert solver.rhs_plan.is_compiled

RK443 is an IMEX Runge-Kutta scheme; RK222, SBDF2 and SBDF3 are also exported and work with this problem. Under MPI the pure-Fourier diagonal-IMEX path takes over automatically — the viscous term stays implicit, the advection and vortex-stretching terms explicit.

┌ Warning: Implicit global matrix: operator Laplacian{...} has no spectral matrix builder and is
│ approximated as IDENTITY-on-operand; the implicit solve will be wrong for this term.

Initial Conditions

Taylor-Green Vortex

The classic test case. The velocity

\[\mathbf{u} = (\sin x \cos y \cos z,\; -\cos x \sin y \cos z,\; 0)\]

has the vorticity $\boldsymbol{\omega} = \nabla\times\mathbf{u}$ we initialize below (the solver reconstructs $\mathbf{A}$ and $\mathbf{u}$ from it on the first step):

x, y, z = local_grids(dist, x_basis, y_basis, z_basis)
Z = reshape(z, 1, 1, :)          # broadcast the third axis

for w in (w1, w2, w3)
    ensure_layout!(w, :g)
end
get_grid_data(w1) .= -cos.(x) .* sin.(y') .* sin.(Z)
get_grid_data(w2) .= -sin.(x) .* cos.(y') .* sin.(Z)
get_grid_data(w3) .=  2 .* sin.(x) .* sin.(y') .* cos.(Z)
for w in (w1, w2, w3)
    ensure_layout!(w, :c)
end

local_grids returns each rank's local slice of the coordinate axes, so the same code is correct in serial and under MPI. (set!(field, ::Function) builds the global meshgrid and therefore throws a DimensionMismatch when the field is distributed.)

Random Initial Conditions

For developed turbulence, seed the vorticity with noise instead — fill_random! is distributed-safe and reproducible — and band-limit it:

for w in (w1, w2, w3)
    fill_random!(w, "g"; seed=42, distribution="normal", scale=1e-1)
    low_pass_filter!(w; scales=(0.25, 0.25, 0.25))     # keep |k| ≤ N/8
end

The low-pass filter matters. Grid noise puts energy right up to the Nyquist mode, where the spectral derivative is not invertible, and the reconstructed velocity comes out only approximately divergence-free (max|∇·u| ≈ 5e-2 after 20 steps at 16^3). With the noise band-limited, the same run gives max|∇·u| ≈ 3e-13.

Analysis and Output

Diagnostics

Kinetic energy, enstrophy and dissipation are box averages. global_mean reduces the local slab and performs the MPI reduction itself — never wrap it in another MPI.Allreduce:

function kinetic_energy()
    for f in (u1, u2, u3)
        ensure_layout!(f, :g)
    end
    0.5 * global_mean(dist, get_grid_data(u1).^2 .+ get_grid_data(u2).^2 .+ get_grid_data(u3).^2)
end

function enstrophy()
    for f in (w1, w2, w3)
        ensure_layout!(f, :g)
    end
    0.5 * global_mean(dist, get_grid_data(w1).^2 .+ get_grid_data(w2).^2 .+ get_grid_data(w3).^2)
end

dissipation() = 2 * nu * enstrophy()     # ε = 2νΩ

For the Taylor-Green initial condition these give Ω = 0.375 and ε = 7.5e-3 at Re = 100 as soon as the vorticity is set, and KE = 0.125 once the solver has taken a step.

Energy Spectrum

power_spectrum requires every axis to be Fourier (true here) and returns a NamedTuple (k, power, bin_counts, bin_edges) with shell-averaged power at the bin centres k:

function velocity_spectrum()
    s1 = power_spectrum(u1)
    s2 = power_spectrum(u2)
    s3 = power_spectrum(u3)
    return (k = s1.k, E = 0.5 .* (s1.power .+ s2.power .+ s3.power))
end

enstrophy_spectrum and energy_spectrum are also available, but they take a VectorField.

Time Integration

run! drives the loop: it steps the solver, writes any handler created with the solver, fires callbacks at their intervals, and closes the handlers at the end. An Int callback interval fires every N iterations; a Float64 interval fires every Δt of simulation time.

function report(s)
    @root_only @printf("t = %.3f  KE = %.6e  Ω = %.6e  ε = %.6e\n",
                       s.sim_time, kinetic_energy(), enstrophy(), dissipation())
end

@root_only @printf("3D Taylor-Green: Re = %.0f, N = %d^3\n", Re, N)

run!(solver; stop_iteration=100, callbacks=[20 => report], progress=false)

MPI.Finalize()

Output of the run above (serial, 32^3, RK443, dt = 2e-3):

3D Taylor-Green: Re = 100, N = 32^3
t = 0.040  KE = 1.247024e-01  Ω = 3.741384e-01  ε = 7.482767e-03
t = 0.080  KE = 1.244096e-01  Ω = 3.733526e-01  ε = 7.467051e-03
t = 0.120  KE = 1.241216e-01  Ω = 3.726414e-01  ε = 7.452828e-03
t = 0.160  KE = 1.238381e-01  Ω = 3.720037e-01  ε = 7.440074e-03
t = 0.200  KE = 1.235593e-01  Ω = 3.714385e-01  ε = 7.428770e-03

The energy decays at the viscous rate while the vortex tubes are still smooth; run longer (and at higher N) to see the enstrophy rise as the Taylor-Green vortex breaks down.

Saving Snapshots

add_file_handler(path, solver; ...) registers the handler with the solver, so run! writes it automatically — declare it before run!, in place of the bare run! above.

snapshots = add_file_handler("turb3d/turb3d", solver; sim_dt=0.1, max_writes=10)
add_task!(snapshots, w1; name="w1")
add_task!(snapshots, w2; name="w2")
add_task!(snapshots, w3; name="w3")

run!(solver; stop_iteration=100, callbacks=[20 => report], progress=false)

The output is NetCDF-4 with the data in groups (vars, time, grids), one file per rank under MPI:

f = "turb3d/turb3d_s1/turb3d_s1.nc"
Tarang.group_variable_names(f, "vars")       # ["w1", "w2", "w3"]
Tarang.group_ncread(f, "vars", "w3")         # (write, x, y, z)
Tarang.group_ncread(f, "time", "sim_time")   # (write,)

Performance Considerations

Resolution Requirements

MPI Scaling

A 3D domain is decomposed with a 2D process mesh (see the warning above):

Omitting mesh= lets the distributor choose one of these shapes for you.

Memory Usage

Estimate: ~100 bytes per grid point (double precision, velocity + vorticity + vector potential + work arrays)

N_est = 256
bytes_per_point = 100
total_memory = N_est^3 * bytes_per_point / 1e9  # GB
println("Estimated memory: $(total_memory) GB")   # Estimated memory: 1.6777216 GB

Running the Simulation

# Set environment
export OMP_NUM_THREADS=1

# Run with MPI (8 ranks: a 2x4 pencil mesh)
mpiexec -n 8 julia --project turbulence_3d.jl

The script above runs unchanged in serial and under MPI. The time-stepped vorticity — and therefore the enstrophy and dissipation — agrees to roundoff across rank counts: after 20 steps max|w1| is 0.998933308710438 at 1, 2 and 4 ranks and 0.998933308710437 at 8, and the enstrophy agrees to ~15 digits. Do not expect bit-identical results at every rank count — the distributed transforms and reductions reassociate the arithmetic, so the last ulp can move. (Reading the velocity back distributed hits the corruption noted under Diagnostics.)

Visualization

There is no VTK/ParaView writer in Tarang. Read a slice back out of the NetCDF snapshots and hand it to the plotting package of your choice (Plots.jl, Makie.jl — neither is a dependency of Tarang):

f  = "turb3d/turb3d_s1/turb3d_s1.nc"
w3 = Tarang.group_ncread(f, "vars", "w3")    # (write, x, y, z)
xg = Tarang.group_ncread(f, "grids", "x")
yg = Tarang.group_ncread(f, "grids", "y")

slice = w3[end, :, :, 1]                     # last write, z = z[1] plane
# heatmap(xg, yg, slice')

Under MPI each rank writes its own file (turb3d_s1_p0.nc, turb3d_s1_p1.nc, …) containing its slab.

References

1. Taylor, G. I., & Green, A. E. (1937). Mechanism of the production of small eddies from large ones.
2. Pope, S. B. (2000). Turbulent Flows. Cambridge University Press.
3. Canuto, C., et al. (2007). Spectral Methods: Fundamentals in Single Domains.