Triply periodic box

F3 Navier–Stokes (3D) Generated 2026-11-06 17:21 UTC 10 plots

Recover the initial velocity field of a 3D flow from a later snapshot, differentiating through the Navier–Stokes rollout.

Fluid flow in 3D. The 3D counterpart of the 2D fluid problem, and the harder stress test for differentiable simulation: real turbulence lives in 3D, and the gradients are correspondingly more delicate.

We solve the 3D incompressible Navier–Stokes equations \(\partial_t \mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \nu\,\nabla^2\mathbf{u}\), \(\nabla\cdot\mathbf{u}=0\), with viscosity \(\nu\) as the primary control parameter. Unlike 2D, the 3D equations admit vortex stretching, which amplifies vorticity and brings on chaos far sooner. The chaos horizon is \(T^\ast \approx 8\text{–}16\) versus \(T^\ast > 64\) in 2D (at \(\nu=10^{-3}\), \(N=16\)). As a result gradient norms grow along the rollout here rather than decaying as they do in 2D.

How these results were produced

These are example results, produced automatically on GitHub-hosted runners and refreshed on every release. Each solver runs on its intended device: GPU-capable solvers on a Tesla T4 GPU node, CPU-only solvers (OpenFOAM, deal.II, FEniCS, Firedrake, TopOpt.jl) on a separate CPU node. Each benchmark job has its virtual machine to itself, but the hardware and the load on the host can change from run to run, so wall-clock times vary. Read them for relative scaling between solvers rather than as absolute timings. Accuracy and gradient metrics change little across hardware, although float32 arithmetic is not identical across devices and chaotic rollouts amplify small discrepancies. For numbers that reflect your setup, run the benchmark yourself on your target hardware.

Boundary conditions

Triply periodic cubic domain \([0, 2\pi]^3\) (the flow wraps around on all three axes). Incompressibility \(\nabla\cdot\mathbf{u}=0\) is enforced by a pressure projection at each time step. No walls or inflow/outflow boundaries.

Checks

The results overview explains the columns and states.

Solver	Forward	Gradient	Optimization	Forward cost	VJP cost	Peak memory	Effort
GPU solversForward and VJP at N=64, multiples of the fastest GPU solver, measured on Tesla T4
PICT	Pass	Pass	Pass	1.28 sfastest	3.31 sfastest	14 GiB	L0
Exponax	Pass	Pass	Pass	—	—	—	L0
INS.jl	Pass	Pass	Pass	—	—	—	L0
jax-cfd	Pass	Pass	Pass	—	—	—	L0
PhiFlow	Pass	Pass	Pass	—	—	—	L3
Warp-NS	Pass1 of 5 run	Pass	Pass	—	—	—	L2
XLB modified setup	Pass	Pass	Pass	—	—	—	L2
CPU solversForward and VJP, multiples of the fastest CPU solver
OpenFOAM	Pass1 of 5 run	Not supported	Not supported	—	—	—	L0

Solvers differ in how general they are. Exponax and jax-cfd do not run Channel flow past a cylinder. Exponax, INS.jl, OpenFOAM, PhiFlow, Warp-NS, XLB and jax-cfd do not run Forced Kolmogorov flow.

Modified setups (1)

These solvers run a substitute for part of the setup, so their checks hold for the modified setup.

XLB: Sub-steps each step to bound the Mach number. Effect: LBM ties the time step to the grid through the lattice Mach number, so XLB runs k times the prescribed number of steps (k = 9 at N = 256 on the 2-D box), and its cost includes them. Without sub-stepping the Taylor-Green case ran at Ma of about 0.88, with errors 84 to 114x the peer median.

Flags and exclusions (41)

PICT on cost/vjp_cost/by_steps: failed. container_died: ConnectionError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))
Exponax on cost/spatial_cost: not run
Exponax on cost/temporal_cost: not run
Exponax on cost/vjp_cost/by_N: not run
Exponax on cost/vjp_cost/by_steps: not run
Exponax on gradient/horizon_sweep_limits: not run
INS.jl on cost/spatial_cost: not run
INS.jl on cost/temporal_cost: not run
INS.jl on cost/vjp_cost/by_N: not run
INS.jl on cost/vjp_cost/by_steps: not run
INS.jl on gradient/horizon_sweep_limits: not run
jax-cfd on cost/spatial_cost: not run
jax-cfd on cost/temporal_cost: not run
jax-cfd on cost/vjp_cost/by_N: not run
jax-cfd on cost/vjp_cost/by_steps: not run
jax-cfd on gradient/horizon_sweep_limits: not run
PhiFlow on cost/spatial_cost: not run
PhiFlow on cost/temporal_cost: not run
PhiFlow on cost/vjp_cost/by_N: not run
PhiFlow on cost/vjp_cost/by_steps: not run
PhiFlow on gradient/horizon_sweep_limits: not run
Warp-NS on forward/agreement: not run
Warp-NS on forward/physical_laws/vs_N: not run
Warp-NS on forward/physical_laws/vs_nu: not run
Warp-NS on forward/physical_laws/vs_steps: not run
Warp-NS on cost/spatial_cost: not run
Warp-NS on cost/temporal_cost: not run
Warp-NS on cost/vjp_cost/by_N: not run
Warp-NS on cost/vjp_cost/by_steps: not run
Warp-NS on gradient/horizon_sweep_limits: not run
XLB on cost/spatial_cost: not run
XLB on cost/temporal_cost: not run
XLB on cost/vjp_cost/by_N: not run
XLB on cost/vjp_cost/by_steps: not run
XLB on gradient/horizon_sweep_limits: not run
OpenFOAM on forward/agreement: not run
OpenFOAM on forward/physical_laws/vs_N: not run
OpenFOAM on forward/physical_laws/vs_nu: not run
OpenFOAM on forward/physical_laws/vs_steps: not run
OpenFOAM on cost/spatial_cost: not run
OpenFOAM on cost/temporal_cost: not run

Resource limits (2)

Sweeps a solver completed only up to some size or rollout length. It ran out of memory or time or, in sweeps that look for the limit, failed past the last value it completed. Unless that happened at the first value, the experiment counts as ok, and the limit shows how far the solver goes on this hardware.

PICT on cost/vjp_cost/by_N, from N=96
PICT on gradient/horizon_sweep_limits, from steps=10240

Not supported (1 solver)

OpenFOAM, 9 experiments. standard icoFoam has no VJP to benchmark; standard icoFoam is non-differentiable (C++, no AD path); DAFoam/OpenFOAM-AD exist but are not deployed in this tesseract; standard icoFoam is a non-differentiable, forward-only solver

Initial conditions

The initial conditions and fields that this setup’s experiments start from. They are plotted without running any solver.

Settings

abc

{
  "N": 32
}

rand_div_free

{
  "N": 32
}

tgv3d

{
  "N": 32
}

Forward

Is the prediction right? Forward experiments compare each solver’s output with a trusted reference, an analytic solution where one exists, and the other solvers.

Agreement between solvers

3D velocity magnitude fields and kinetic energy spectra per solver, swept over viscosity ν, compared against a converged spectral reference. (The 3D TGV initial condition is not advection-free, so the linearized analytic decay is only a short-horizon consistency check, not an accuracy reference.) The reference is an Exponax run at N=64, spectrally truncated to N=16. Exponax’s low error here, and its top place in the forward ranking, therefore largely measure agreement with its own higher-resolution solution.

Settings

Sweeps nu ∈ {0.001, 0.01, 0.05}

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 16,
    "nu": 0.001,
    "dt": 0.01,
    "steps": 50
  },
  "sweep": {
    "key": "nu",
    "values": [
      0.001,
      0.01,
      0.05
    ]
  }
}

Agreement with the reference

Relative error vs grid resolution N against a converged spectral reference at t=1.28, refining space and time together at a fixed CFL number. The legend gives each solver’s observed order of accuracy p.

Settings

Sweeps N ∈ {8, 16, 32}

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 8,
    "nu": 0.05,
    "dt": 0.04,
    "steps": 32,
    "cfl_N_base": 32
  },
  "sweep": {
    "key": "N",
    "values": [
      8,
      16,
      32
    ]
  }
}

Results by solver

Solver	Mean rel. error
Exponax	7.26e-04
PICT	7.96e-03
XLB	8.45e-03
PhiFlow	3.98e-02
INS.jl	4.14e-02
jax-cfd	4.14e-02

From Agreement between solvers.

Cost

What does it cost? Wall-clock time and peak memory of the forward and VJP passes as the problem size \(N\) and the number of steps grow. See the note on how these results were produced before reading absolute numbers.

Settings

spatial_cost

Sweeps N ∈ {16, 32, 48, 64, 96, 128}

{
  "physics": {
    "nu": 0.01,
    "dt": 0.01,
    "steps": 50,
    "N": 16
  },
  "cost": {
    "n_trials": 3
  },
  "sweep": {
    "key": "N",
    "values": [
      16,
      32,
      48,
      64,
      96,
      128
    ]
  }
}

temporal_cost

Sweeps steps ∈ {10, 50, 100, 500, 1000, 2000}

{
  "physics": {
    "nu": 0.01,
    "dt": 0.01,
    "N": 48,
    "steps": 10
  },
  "cost": {
    "n_trials": 3
  },
  "sweep": {
    "key": "steps",
    "values": [
      10,
      50,
      100,
      500,
      1000,
      2000
    ]
  }
}

Gradient

Is the gradient right? Gradient experiments compare each solver’s AD or adjoint gradient with finite differences, in magnitude (relative \(L^2\) error) and direction (cosine similarity).

Settings

fd_check

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 16,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 10
  },
  "fd": {
    "eps_values": [
      5.0,
      1.0,
      0.1,
      0.01,
      0.001,
      0.0001
    ],
    "n_dirs": 10
  }
}

horizon_sweep_limits

Sweeps steps ∈ {40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480, 40960, 81920}

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 20,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 40
  },
  "sweep": {
    "key": "steps",
    "values": [
      40,
      80,
      160,
      320,
      640,
      1280,
      2560,
      5120,
      10240,
      20480,
      40960,
      81920
    ]
  }
}

jacobian_svd

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 8,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 10
  },
  "jacobian": {
    "n_alphas": 41,
    "alpha_range": 0.3
  }
}

jacobian_svd_nu01

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 8,
    "nu": 0.01,
    "dt": 0.05,
    "steps": 10
  },
  "jacobian": {
    "n_alphas": 41,
    "alpha_range": 0.3
  }
}

jacobian_svd_steps20

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 8,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 20
  },
  "jacobian": {
    "n_alphas": 41,
    "alpha_range": 0.3
  }
}

jacobian_svd_steps40

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 8,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 40
  },
  "jacobian": {
    "n_alphas": 41,
    "alpha_range": 0.3
  }
}

Gradient against finite differences

Finite-difference gradient error U-curves and direction cosine vs perturbation ε for each solver on the 3D Taylor-Green vortex IC.

Settings

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 16,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 10
  },
  "fd": {
    "eps_values": [
      5.0,
      1.0,
      0.1,
      0.01,
      0.001,
      0.0001
    ],
    "n_dirs": 10
  }
}

How far gradients reach in rollout length

Per-solver rollout-limit table reporting step count at first failure, failure type, and wall time per successful step.

Settings

Sweeps steps ∈ {40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480, 40960, 81920}

{
  "ic": {
    "name": "tgv3d",
    "seed": 0
  },
  "physics": {
    "N": 20,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 40
  },
  "sweep": {
    "key": "steps",
    "values": [
      40,
      80,
      160,
      320,
      640,
      1280,
      2560,
      5120,
      10240,
      20480,
      40960,
      81920
    ]
  }
}

Results by solver

Solver	Best-ε FD error	1 − cosine
XLB	4.12e-06	1.33e-11
Warp-NS	1.30e-05	1.53e-10
PhiFlow	1.85e-05	1.48e-10
jax-cfd	1.96e-05	1.15e-10
PICT	2.78e-05	4.05e-10
INS.jl	2.91e-05	1.02e-09
Exponax	1.15e-04	4.87e-09

From Gradient against finite differences.

Optimization

Can you optimize through it? Optimization experiments run a gradient-based optimizer on each solver’s own gradients. A gradient can pass the finite-difference check and still fail to drive an optimization loop.

Initial-condition recovery (L-BFGS)

Final IC recovery error per solver from zero-initialised L-BFGS optimisation with divergence-free projection.

Settings

Sweeps steps ∈ {100}

{
  "ic": {
    "name": "rand_div_free",
    "seed": 0
  },
  "physics": {
    "N": 16,
    "nu": 0.01,
    "dt": 0.02,
    "steps": 100
  },
  "optim": {
    "ic_init_type": "zeros",
    "max_iters": 100,
    "patience": 20,
    "failure_threshold": 2.0,
    "snap_interval": 5,
    "ic_seeds": [
      0,
      1,
      2
    ],
    "record_diagnostics": true
  },
  "sweep": {
    "key": "steps",
    "values": [
      100
    ]
  },
  "optimizer": "bfgs_proj",
  "_exp_key": "recovery_constant_ic_bfgs_proj"
}

Results by solver

Solver	Final IC recovery error	Iterations to plateau
Exponax	9.63e-05	—
INS.jl	2.89e-03	—
PICT	2.97e-03	—
jax-cfd	3.05e-03	—
PhiFlow	3.27e-03	—
Warp-NS	4.48e-03	—
XLB	8.39e-03	—

From Initial-condition recovery (L-BFGS).