Navier–Stokes (2D)

F2 Auto-generated 2026-10-30 15:14 UTC 13 plots

Optimize the inflow of a 2D channel flow to minimize drag on a cylinder, differentiating through an incompressible fluid solver.

Fluid flow in 2D. We simulate an incompressible fluid (think of stirred water) and differentiate through the simulation. The headline question is whether each solver returns a trustworthy gradient as the flow becomes chaotic.

Formally, we solve the 2D 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 the kinematic viscosity \(\nu\) as the primary control parameter. The nonlinear advection term \(\nabla\cdot(\mathbf{u}\otimes\mathbf{u})\) transfers energy across scales; at low \(\nu\) the flow develops turbulent cascades and a positive Lyapunov exponent, so long-horizon gradients grow exponentially sensitive to perturbations. This is the regime where backprop-through-simulation is hardest.

How these results were produced

These are example results, produced automatically on GitHub Actions 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, INS.jl, deal.II, FEniCS, Firedrake, TopOpt.jl) on a separate CPU node, whose processor model can differ between runs. Accuracy and gradient metrics are hardware-independent and reproducible. Wall-clock numbers reflect commodity cloud hardware and can vary by 10–15% between runs, so read them for relative scaling between solvers rather than as absolute timings. For numbers that reflect your setup, run the benchmarks yourself on your target hardware.

Boundary conditions

Doubly-periodic square domain \([0, 2\pi]^2\) (the flow wraps around at every edge), with incompressibility \(\nabla\cdot\mathbf{u}=0\) enforced by a pressure projection at each time step. The cylinder experiments use the Schäfer–Turek DFG 2D-1 channel instead: \([0, 2.2] \times [0, 0.41]\) with a cylinder of diameter 0.1 at \((0.2, 0.2)\), a parabolic inflow of mean speed 0.2 on the left, a zero-gradient outflow at \(p=0\) on the right, and no-slip walls; \(\nu = 10^{-3}\), so \(Re = 20\).

Checks at a glance

Periodic box (Taylor–Green, multimode)

Navier–Stokes (2D)

Solver	Forward	Gradient	Forward cost	VJP cost	Reliability
GPU solvers · cost relative to the fastest GPU solver at N=256
jax-cfd	Pass	Pass	1× (fastest)	1× (fastest)	19/19 ok
Warp-NS	Pass	Pass	3.2×	5.0×	19/19 ok
XLB	Pass	Pass	8.7×	18×	19/19 ok
PhiFlow	Pass	Pass	1289×	1186×	19/19 ok
PICT	Borderline	Pass	161×	129×	18/19 ok1 flagged
CPU solvers · cost relative to the fastest CPU solver at N=256
INS.jl	Pass	Pass	1× (fastest)	1× (fastest)	19/19 ok
OpenFOAM	Pass	Not supported	33×	—	9/9 ok

Passes every check: INS.jl, jax-cfd, Warp-NS, XLB, PhiFlow.

Why cells are flagged (1)

PICT · forward/baseline · flagged — error converges at order 0.35 in N, below min_order=0.8 (explained: PICT runs PISO with its default two pressure correctors, and its momentum predictor carries no pressure gradient. The splitting error this leaves scales with ν·dt/dx², which grows as the grid is refined at fixed CFL, so the error levels off on the finest grids and the observed order drops. More correctors remove the floor at a multiple of the wall time and reverse-mode memory)

Not supported (1 solver)

OpenFOAM · 10 experiments — OpenFOAM has no VJP to benchmark; OpenFOAM is non-differentiable (C++, no AD path); DAFoam/OpenFOAM-AD exist but are not deployed in this tesseract

Channel flow past a cylinder

Navier–Stokes (2D)

Optimization task: inflow optimization for drag

Solver	Forward	Optimization	Reliability
GPU solvers
PhiFlow	Pass	Pass	2/2 ok
Warp-NS	Pass	Pass	2/2 ok
XLB	Pass	Pass	2/2 ok
PICT	Pass	Not run	1/2 ok1 not run
jax-cfd	Not run	Not run	0/2 ok2 not run
CPU solvers
OpenFOAM	Pass	Not supported	1/1 ok
INS.jl	Not run	Not run	0/2 ok2 not run

Passes every check: PhiFlow, Warp-NS, XLB.

Why cells are flagged (5)

PICT · optimization/drag_opt · not run — each drag_opt iteration runs a full forward and VJP through PICT's unpreconditioned pressure solves on the body-fitted channel mesh, so the experiment takes hours on the benchmark GPU
INS.jl · forward/cylinder · not run — this tesseract implements only periodic flows; INS.jl itself supports channel BCs (Dirichlet inflow, pressure outflow) and obstacles via Brinkman forcing, which are not wired up here
INS.jl · optimization/drag_opt · not run — this tesseract implements only periodic flows; INS.jl itself supports channel BCs (Dirichlet inflow, pressure outflow) and obstacles via Brinkman forcing, which are not wired up here
jax-cfd · forward/cylinder · not run — this tesseract implements only periodic flows; jax-cfd itself has non-periodic pressure solvers for channel flows, which are not wired up here
jax-cfd · optimization/drag_opt · not run — this tesseract implements only periodic flows; jax-cfd itself has non-periodic pressure solvers for channel flows, which are not wired up here

Not supported (1 solver)

OpenFOAM · 1 experiments — OpenFOAM is a non-differentiable, forward-only solver

Pass / borderline / fail per suite, from the same results as the sections below. See How to read the checks for what each state means and how cost is compared.

Initial Conditions

Visualisation of each initial condition (the starting field a run is launched from) available for this problem. IC plots are generated without running any solver.

Settings

Dfg Inflow

{
  "resolution": 20
}

Flat Inflow

{
  "N": 64,
  "U": 0.5
}

Multimode

{
  "N": 64
}

Tgv

{
  "N": 64
}

Uniform

{
  "N": 64,
  "U": 1.0
}

Forward

Is the prediction right? Forward-pass benchmarks check each solver’s output against a trusted reference (and an analytic solution where one exists): inter-solver agreement, field-level diagnostics, and long-run stability.

Baseline

Relative error vs grid resolution \(N\) against the analytic TGV solution 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 ∈ {16, 32, 64, 128}

{
  "ic": {
    "name": "tgv",
    "seed": 0
  },
  "physics": {
    "N": 16,
    "nu": 0.05,
    "dt": 0.01,
    "steps": 128,
    "cfl_N_base": 128
  },
  "sweep": {
    "key": "N",
    "values": [
      16,
      32,
      64,
      128
    ]
  }
}

Cylinder

Steady flow past a cylinder in a channel (Schäfer–Turek DFG 2D-1, Re = 20) at 10 and 20 cells per diameter. Each solver’s drag coefficient at 20 is scored against the reference Cd = 5.5795. The pressure difference across the cylinder (reference 0.1175) and the lift coefficient are also reported. The vorticity fields are for visual comparison.

Settings

Sweeps resolution ∈ {10, 20}

{
  "ic": {
    "name": "dfg_inflow",
    "seed": 0
  },
  "physics": {
    "resolution": 10
  },
  "sweep": {
    "key": "resolution",
    "values": [
      10,
      20
    ]
  }
}

Tgv Nu Sweep

Relative error vs viscosity \(\nu\) for each solver at a fixed TGV initial condition.

Settings

Sweeps nu ∈ {0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05}

{
  "ic": {
    "name": "tgv",
    "seed": 42
  },
  "physics": {
    "N": 64,
    "dt": 0.05,
    "steps": 20,
    "nu": 0.0001
  },
  "sweep": {
    "key": "nu",
    "values": [
      0.0001,
      0.0005,
      0.001,
      0.005,
      0.01,
      0.05
    ]
  }
}

Results by solver

Solver	Mean rel. error
Warp-NS	6.20e-04
XLB	6.57e-04
PhiFlow	2.39e-03
jax-cfd	2.39e-03
INS.jl	2.39e-03
PICT	3.53e-03
OpenFOAM	3.88e-03

Sorted by mean relative error against the reference solution (lower is closer to the reference; whether a difference matters depends on the task, see the checks at the top of the page).

Cost

What does it cost? Wall-clock scaling of the forward and VJP passes with problem size \(N\) and the number of integration steps. Timings come from dedicated runners with no concurrent workloads; see the reliability note at the top of the page before reading absolute numbers.

Settings

Spatial Cost

Sweeps N ∈ {64, 128, 192, 256}

{
  "physics": {
    "nu": 0.01,
    "dt": 0.01,
    "steps": 100,
    "N": 64
  },
  "cost": {
    "n_trials": 3
  },
  "sweep": {
    "key": "N",
    "values": [
      64,
      128,
      192,
      256
    ]
  }
}

Temporal Cost

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

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

Results by solver

Solver	Device	Forward time	VJP time
jax-cfd	GPU	0.0311 s @ N=256	0.0817 s @ N=256
Warp-NS	GPU	0.0998 s @ N=256	0.407 s @ N=256
XLB	GPU	0.272 s @ N=256	1.51 s @ N=256
PICT	GPU	5 s @ N=256	10.5 s @ N=256
PhiFlow	GPU	40.1 s @ N=256	96.9 s @ N=256
INS.jl	CPU	0.832 s @ N=256	7.3 s @ N=256
OpenFOAM	CPU	27.4 s @ N=256	—

Forward and VJP wall-clock time, each shown at the largest problem size N the solver completed for that pass; GPU solvers first, each group sorted by forward time. The VJP time is one gradient from scratch, a forward plus an adjoint pass. Forward-only solvers have no VJP entry. See the reliability note above before comparing across devices.

Gradient

Is the gradient right? Gradient benchmarks compare each solver’s AD/adjoint gradient against a finite-difference ground truth. We report magnitude error (relative \(L^2\)) and direction agreement (cosine similarity) across parameter, resolution, and horizon sweeps. The horizon sweep in particular exposes how gradients degrade as the rollout lengthens.

Settings

Analytic Check

Sweeps nu ∈ {0.001, 0.005, 0.01, 0.02, 0.05}

{
  "ic": {
    "name": "tgv",
    "seed": 0
  },
  "physics": {
    "N": 32,
    "dt": 0.05,
    "steps": 20,
    "nu": 0.001
  },
  "reference": "<callable _tgv_analytic_grad>",
  "sweep": {
    "key": "nu",
    "values": [
      0.001,
      0.005,
      0.01,
      0.02,
      0.05
    ]
  },
  "project": "<callable project_divergence_free>"
}

Fd Check

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

Horizon Sweep

Sweeps steps ∈ {5, 10, 20, 40, 80, 160, 320}

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 16,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 5
  },
  "fd": {
    "eps_values": [
      1.0,
      0.1,
      0.01,
      0.001
    ],
    "n_dirs": 8
  },
  "sweep": {
    "key": "steps",
    "values": [
      5,
      10,
      20,
      40,
      80,
      160,
      320
    ]
  }
}

Jacobian Svd

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 8,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 10
  }
}

Jacobian Svd Nu01

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 8,
    "nu": 0.01,
    "dt": 0.05,
    "steps": 10
  }
}

Jacobian Svd Steps20

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 8,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 20
  }
}

Jacobian Svd Steps40

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 8,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 40
  }
}

Param Sweep

Sweeps nu ∈ {0.05, 0.01, 0.005, 0.001}

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 16,
    "dt": 0.05,
    "steps": 200,
    "nu": 0.05
  },
  "fd": {
    "eps_values": [
      5.0,
      1.0,
      0.1,
      0.01,
      0.001
    ],
    "n_dirs": 15
  },
  "sweep": {
    "key": "nu",
    "values": [
      0.05,
      0.01,
      0.005,
      0.001
    ]
  }
}

Finite-Difference Check

U-curves of finite-difference gradient error vs perturbation size \(\varepsilon\) together with subspace cosine; validates VJP correctness.

Settings

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

Parameter Sweep

Gradient norm, best-\(\varepsilon\) FD error, direction cosine, and U-curves vs the sweep parameter.

Settings

Sweeps nu ∈ {0.05, 0.01, 0.005, 0.001}

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 16,
    "dt": 0.05,
    "steps": 200,
    "nu": 0.05
  },
  "fd": {
    "eps_values": [
      5.0,
      1.0,
      0.1,
      0.01,
      0.001
    ],
    "n_dirs": 15
  },
  "sweep": {
    "key": "nu",
    "values": [
      0.05,
      0.01,
      0.005,
      0.001
    ]
  }
}

Horizon Sweep

Gradient norm, FD error, and direction cosine vs rollout horizon \(T = \mathrm{steps}\cdot dt\).

Settings

Sweeps steps ∈ {5, 10, 20, 40, 80, 160, 320}

{
  "ic": {
    "name": "multimode",
    "seed": 42
  },
  "physics": {
    "N": 16,
    "nu": 0.001,
    "dt": 0.05,
    "steps": 5
  },
  "fd": {
    "eps_values": [
      1.0,
      0.1,
      0.01,
      0.001
    ],
    "n_dirs": 8
  },
  "sweep": {
    "key": "steps",
    "values": [
      5,
      10,
      20,
      40,
      80,
      160,
      320
    ]
  }
}

Results by solver

Solver	Best-ε FD error	1 − cosine
jax-cfd	7.12e-06	3.51e-11
PhiFlow	8.68e-06	8.68e-11
XLB	1.15e-05	1.67e-10
INS.jl	1.53e-05	1.15e-10
Warp-NS	4.07e-05	1.64e-09
PICT	4.62e-05	1.12e-09

Sorted by the best-ε finite-difference error of the gradient (lower is more trustworthy); direction cosine near 1 confirms the gradient points the right way.

Optimization

Can you optimize through it? End-to-end optimization benchmarks run a gradient-based optimizer using each solver’s own gradients: recovery of initial conditions or physical parameters, topology optimization, and drag minimization. This is the ultimate test, since a gradient can pass the finite-difference check yet still fail to drive a full optimization loop.

Drag Optimisation

Drag minimization over the inlet profile of the DFG 2D-1 channel at a fixed flow rate: drag coefficient per iteration, optimized vs initial inflow profiles, and the flow fields before and after.

Settings

{
  "ic": {
    "name": "dfg_inflow",
    "seed": 0
  },
  "physics": {
    "resolution": 10,
    "t_end": 4.0
  },
  "optim": {
    "lr": 0.01,
    "max_iters": 100,
    "patience": 10,
    "min_rel_improvement": 0.001,
    "snap_interval": 10
  }
}

Results by solver

Solver	Final drag	Converged
PhiFlow	2.90e+00	yes
Warp-NS	2.92e+00	yes
XLB	3.15e+00	yes

Sorted by the final objective reached within the iteration budget (lower is better).