Solver Reference

Each solver card shows its numerical scheme, AD strategy, and Tesseract schema, and lists the experiments it is excluded from in each problem config. Each card also says what it took to run the solver here: the effort per setup (the highest level among the interventions it needs, L0 documented use to L3 changes to the solver), and under Interventions and deviations every intervention with its reason, evidence and upstream issue, the setup deviations, and the interventions nobody has a reason for yet. Click Details on any solver to see all of this and its field tables. The ∂ column marks fields that support automatic differentiation (VJP/JVP).

Legend. Each card shows a Numerics: line — the discretization followed by the numerical method — and three badge categories: language (backend / runtime), AD: strategy (autodiff / adjoint / hybrid / forward-only), and GPU on GPU-capable solvers.

Checks. Where benchmark results are available, the card footer shows how many of the forward, gradient and optimization checks the solver passes on each setup it runs: H heat transfer, S structural mechanics, F2/F3 Navier–Stokes 2D/3D. Checks the method cannot support are left out. See Results for what each check means.

Navier–Stokes (Grid)

Navier–Stokes (Grid)

Exponax

Numerics: Spectral, ETDRK2

jax AD: autodiff GPU

Exponential time-differencing Runge–Kutta (ETDRK) spectral integrator for 2D and 3D Navier–Stokes, enforcing incompressibility to machine precision by construction. Gradients via JAX source-transformation AD

F2 periodic · 1/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://fkoehler.site/exponax/

Path: mosaic/tesseracts/navier-stokes-grid/exponax

Image: exponax_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward borderline, gradient pass; F2 Channel flow past a cylinder: forward not supported, optimization not supported; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Excluded experiments

• ns-grid · forward/cylinder: a Fourier pseudo-spectral method represents only periodic fields, so it cannot express a channel with walls, an inflow and an outflow
• ns-grid · optimization/drag_opt: a Fourier pseudo-spectral method represents only periodic fields, so it cannot express a channel with walls, an inflow and an outflow

Effort: Periodic box (2D) L2 · Forced Kolmogorov flow (2D) L3 — 2 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
vorticity-velocity-conversion: The 2-D path converts the velocity to Exponax’s vorticity variable (q = -curl u) with a spectral curl, runs Exponax’s vorticity stepper, and converts back with a spectral Biot-Savart inversion. The mean velocity, which the vorticity form drops, is carried separately, under the mean force and drag when there are any.	L2	required	Periodic box (2D), Forced Kolmogorov flow (2D) (cannot-express): The benchmark’s state is velocity, and Exponax 0.2.0’s velocity steppers (NavierStokesVelocity, KolmogorovFlowVelocity) accept only three spatial dimensions; in 2-D it offers the vorticity form only. Evidence: https://github.com/Ceyron/exponax/blob/v0.2.0/exponax/stepper/_navier_stokes.py (num_spatial_dims checks).	—
canonical-force-injection: With a body force, the 2-D path sets the injection term inside Exponax’s KolmogorovFlowVorticity stepper (the private _integrator._nonlinear_fun.injection) to the spectral curl of the force, as complex coefficients.	L3	required	Forced Kolmogorov flow (2D) (cannot-express): Exponax 0.2.0 builds the injection from injection_mode and injection_scale as a real coefficient on one mode, which samples sin(κy) at y = j·h. The canonical force is sampled at (j + ½)h, so the public options put it half a cell off, a phase error of κh/2 or about 20% of the laminar Kolmogorov solution at N = 64. Evidence: https://github.com/Ceyron/exponax/blob/v0.2.0/exponax/nonlin_fun/_vorticity_convection.py (the injection term); tests/test_ns_grid_forcing.py::test_laminar_kolmogorov_transient; #272.	—

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]	✓	Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None	✓	Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.
drag	Array[(1,), Float32]	✓	Linear drag coefficient
order	int		ETDRK time integration order (0-4)
kolmogorov_forcing	bool		Enable Kolmogorov sinusoidal body forcing for sustained turbulence
injection_mode	int		Wavenumber at which energy is injected (Kolmogorov forcing)
injection_scale	Array[(1,), Float32]	✓	Amplitude of the Kolmogorov body forcing

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None		Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

INS.jl

Numerics: Finite Volume, Projection, RK4

julia AD: autodiff GPU

Julia staggered finite-volume solver with pressure projection and RK4 time integration, on the GPU through its CUDA backend (KernelAbstractions.jl). Periodic boxes use a spectral pressure solve; the channel a direct one, with the cylinder as a Brinkman penalization. Differentiates through the time loop via Zygote.jl reverse-mode AD.

F2 periodic · 2/2 pass F2 cylinder · 2/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://agdestein.github.io/IncompressibleNavierStokes.jl/dev/

Path: mosaic/tesseracts/navier-stokes-grid/ins-jl

Image: ins_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward pass, gradient pass; F2 Channel flow past a cylinder: forward pass, optimization pass; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Effort: Periodic box (2D) none · Channel flow past a cylinder L2 · Forced Kolmogorov flow (2D) none · Triply periodic box (3D) none — 7 interventions, 1 with unknown reason — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
hand-written-rk4-rollout: The VJP differentiates a hand-written, non-mutating RK4 loop over create_right_hand_side with Zygote.pullback. On periodic boxes apply runs the library’s in-place solve_unsteady with RKMethods.RK44 instead, except with a body force, a drag or snapshots, where it runs the hand-written loop too, so that with a force both integrate the same right-hand side. No comparison with the library’s force! hook (examples/Kolmogorov2D.jl) is recorded. In the channel apply runs the hand-written loop too, with the channel’s right-hand side, and records the force on the cylinder after every step.	L2	choice	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D): Zygote cannot differentiate array mutation, so the VJP needs a non-mutating rollout; INS.jl’s documented way is a loop over its non-mutating timestep. No reason is recorded for writing RK4 by hand instead. Keeping the in-place stepper for apply was faster on the CPU, before the move to the GPU (#307): 0.56 s against 2.49 s for the non-mutating rollout at N = 256, 100 steps. The two paths agreed to 2.2e-4 (2-D) and 1.0e-4 (3-D) in float32. Evidence: https://github.com/agdestein/IncompressibleNavierStokes.jl/blob/v5.0.0/examples/Differentiability.jl; #103. Channel flow past a cylinder: The VJP needs the non-mutating rollout for the same reason as on the periodic boxes. apply shares it, so the drag gradient check compares the VJP with exactly the rollout it differentiates, and the per-step force comes from the same loop. solve_unsteady was not tried for the channel’s apply. Evidence: tests/test_ins_jl_cylinder.py::test_drag_gradient_matches_fd.	—
brinkman-cylinder: Represents the cylinder as a Brinkman penalization force −χu/η in the momentum right-hand side, with χ = 1 on every face touching a cell of the shared staircase (solid_cells) and η = dt. The force on the cylinder is the penalization’s integral, the sum of χu/η Δx² over the faces, after every step.	L2	required	Channel flow past a cylinder (cannot-express): INS.jl 5.0.0 has no obstacles or immersed boundaries. Its documented hook for extra momentum terms is a body force (examples/Actuator2D.jl), and the penalization is one. The force is integrated explicitly: at η = dt / 2.5 a run at 20 cells per diameter went non-finite. Halving η from dt moved Cd by 0.1% at 10 cells per diameter. Evidence: #325, notes.	—
differentiable-inlet: Gives the channel a no-slip DirichletBC at the inlet and adds the inflow profile to the inlet faces of the state after INS.jl applies the boundary conditions, in every right-hand side.	L2	required	Channel flow past a cylinder (no-gradient): INS.jl’s rrule for apply_bc_u returns no tangent for the setup, which holds a DirichletBC’s value, and building a Setup inside Zygote.pullback fails with “Mutating arrays is not supported – called push!”. The inflow profile, which drag_opt optimizes, would get no gradient. Added after the BCs, the drag gradient matched a central difference to 1e-5 (prototype, 40 steps at 10 cells per diameter) and 3e-5 over drag_opt’s horizon. Evidence: src/boundary_conditions.jl of IncompressibleNavierStokes.jl v5.0.0 (the rrule of apply_bc_u); tests/test_ins_jl_cylinder.py::test_drag_gradient_matches_fd.	—
separable-pressure-solver: Solves the channel’s pressure Poisson equation by diagonalizing INS.jl’s Laplacian, which on the uniform grid is the Kronecker sum of an x and a y operator read off laplacian_mat(setup), instead of with the psolver_direct that default_psolver picks. The transforms run in Float64.	L2	required	Channel flow past a cylinder (crash, too-slow): On CuArrays psolver_direct needs INS.jl’s CUDSS extension, which imports CUDSS.CUDA. CUDSS 0.7 and later, the only versions compatible with CUDA.jl 6, no longer define it, so the extension fails to load, Pkg.precompile() exits 1 and psolver_direct raises “Unsupported array type”. psolver_transform asserts that every boundary is periodic or Dirichlet, which the PressureBC outlet is not. psolver_cg on the GPU was stopped after 600 s at 10 cells per diameter, where this solver’s run takes 2.7 s. Evidence: ext/CUDSSExt.jl and src/pressure.jl of IncompressibleNavierStokes.jl v5.0.0; CUDSS.jl’s registry compat (CUDSS 0.7 and later require CUDACore 6); #325, notes.	Not reported. INS.jl’s ext/CUDSSExt.jl would need to import CUDA (or CUDACore) directly instead of through CUDSS.
float32-setup: Builds the Setup from Float32 grid vectors and keeps literals eltype-preserving, so the solve and its pullback run in float32.	L0	choice	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D), Channel flow past a cylinder: With Float64 grid vectors the state and the Zygote tape were promoted to Float64 after the first projection, so INS.jl paid Float64 cost with float32 inputs and outputs. Float32 matches the other float32 solvers; #103 offered keeping Float64 and disclosing it instead. No cost factor was measured separately, on the CPU or on the GPU. Evidence: #103, item 3; INS.jl’s manual on floating point precision (docs/src/manual/backends.md at v5.0.0).	—
cuda-backend (tuned): With a GPU present, builds the Setup with backend = CUDABackend() and moves the initial state to the device, so the solve and its Zygote pullback run on INS.jl’s CUDA backend instead of the default CPU backend.	L0	choice	Periodic box (2D), Forced Kolmogorov flow (2D) (too-slow): On the CPU, vjp_cost at N = 256 took 9.20 s against 0.96 to 2.59 s on a shared A100. Small cases got slower: spatial_cost at N = 64 went from 0.085 s to 0.16 to 0.64 s. Every status and forward error is unchanged. Evidence: #307, notes. Triply periodic box (3D) (too-slow): On the CPU, vjp_cost at N = 64 took 54.6 s against 1.58 to 3.72 s, and 100 steps 47.2 s against 2.38 to 3.15 s. The 5120-step horizon_sweep_limits VJP needs 28.3 GB on the device and runs out of memory on the T4. Evidence: #307, notes. Channel flow past a cylinder: With one Julia thread on the CPU, a forward run to t = 8 took 3.4 s at 10 cells per diameter and 25.5 s at 20, against 2.7 s and 5.4 s on an A100. drag_opt’s VJP (10 cells per diameter, t = 4) took 8.2 s against 3.9 s. Warm calls; the A100 was not shared. Evidence: #325, notes.	—
julia-cpu-target-generic: The image sets JULIA_CPU_TARGET=generic, so the package images Julia precompiles at build time target a generic x86-64 CPU rather than the build machine’s. Code compiled at run time is unaffected.	L0	unknown	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D), Channel flow past a cylinder: No reason recorded; present since the initial commit. Julia’s docs recommend it where images run on CPUs other than the one that built them. The solver’s kernels run on the GPU (#307) and are compiled at run time, so only host-side code from the package images could be slower; not measured.	—

Reason unknown. Nobody has measured whether these are needed; a measurement, or removing them, would settle it:

• julia-cpu-target-generic (L0): The image sets JULIA_CPU_TARGET=generic, so the package images Julia precompiles at build time target a generic x86-64 CPU rather than the build machine’s. Code compiled at run time is unaffected.

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]	✓	Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None	✓	Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None	✓	Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

jax-cfd

Numerics: Finite Difference, Projection, RK4

jax AD: autodiff GPU

JAX-native incompressible Navier–Stokes solver on a staggered MAC grid with finite-difference advection and spectral FFT pressure projection. Gradients via JAX source-transformation AD through the time loop.

F2 periodic · 2/2 pass F2 cylinder · 0/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://github.com/google/jax-cfd

Path: mosaic/tesseracts/navier-stokes-grid/jax-cfd

Image: jax_cfd_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward pass, gradient pass; F2 Channel flow past a cylinder: forward blocked, optimization blocked; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Excluded experiments

• ns-grid · forward/cylinder: jax-cfd 0.2.1 cannot express the channel through its API. Its advection raises NotImplementedError for any inflow (only no-penetration and periodic boundaries are supported), it pads a zero-gradient boundary only for cell-centered fields, so the outlet’s face-centered u_x cannot have one, its boundary values are scalars, so the inlet cannot carry a profile, and it has no obstacles. Running it would take a hand-written convection, inlet and outlet conditions and obstacle on jax-cfd’s grid primitives, which this tesseract does not have
• ns-grid · optimization/drag_opt: jax-cfd 0.2.1 cannot express the channel through its API. Its advection raises NotImplementedError for any inflow (only no-penetration and periodic boundaries are supported), it pads a zero-gradient boundary only for cell-centered fields, so the outlet’s face-centered u_x cannot have one, its boundary values are scalars, so the inlet cannot carry a profile, and it has no obstacles. Running it would take a hand-written convection, inlet and outlet conditions and obstacle on jax-cfd’s grid primitives, which this tesseract does not have

Effort: Periodic box (2D) none · Forced Kolmogorov flow (2D) L2 · Triply periodic box (3D) none — 2 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
central-advection-rk4 (tuned): Passes convect=advection.convect_linear and time_stepper=time_stepping.classic_rk4 to semi_implicit_navier_stokes instead of its defaults, van Leer advection with forward Euler.	L0	choice	Periodic box (2D), Forced Kolmogorov flow (2D) (inaccurate-forward): With the default scheme the Taylor-Green errors were 1.3e-2 to 1.5e-2, 5 to 6 times the peer median (borderline), the convergence order was 1.25 and the fd_check error 5.5e-4. With this scheme: 2.3e-3 to 2.4e-3 (1.0 times the median), order 2.00 and 9.3e-6. Central advection is unstable under forward Euler and RK2 in the PR’s probes, hence RK4. The forward pass is 2 to 3 times slower. Evidence: #225, before/after and scheme tables (merged through #232). Triply periodic box (3D): Not compared with the default on this setup; jax-cfd joined ns-3d-grid with this scheme, and every cell was ok. Evidence: #283, notes.	https://github.com/google/jax-cfd/issues/259 asked which advection scheme the maintainers would pick; no reply, and the project says it is no longer maintained.
summed-forcing: Adds the body force and the linear drag in one forcing function passed to semi_implicit_navier_stokes, in place of jax-cfd’s sum_forcings.	L2	required	Forced Kolmogorov flow (2D) (crash): sum_forcings adds the terms with equations.sum_fields, which calls jax.tree_map. The pinned JAX 0.10.0 no longer has jax.tree_map, so combining the force and the drag through the library would raise an AttributeError. This comes from reading both sources and was not run. Evidence: https://github.com/google/jax-cfd/blob/0.2.1/jax_cfd/base/equations.py (sum_fields); #272.	No issue filed: the project says it is no longer maintained.

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]	✓	Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None		Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.
density	Array[(1,), Float32]	✓	Density of the fluid
inner_steps	int		Time steps of dt per output step, so a run advances steps * inner_steps * dt

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None		Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

OpenFOAM

Numerics: Finite Volume, PISO/PIMPLE, Euler

cpp AD: forward-only

OpenFOAM icoFoam PISO solver run as a forward-only reference baseline (C++ subprocess). No reverse-mode AD; used to anchor cross-solver comparisons.

F2 periodic · 0/1 pass F2 cylinder · 1/1 pass F3 periodic · 1/1 pass Details

Upstream docs: https://www.openfoam.com/documentation/overview

Path: mosaic/tesseracts/navier-stokes-grid/openfoam

Image: openfoam_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward borderline, gradient not supported; F2 Channel flow past a cylinder: forward pass, optimization not supported; F3 Triply periodic box (Taylor–Green): forward pass, gradient not supported, optimization not supported — see Results.

Excluded experiments

• ns-3d-grid · cost/vjp_cost: standard icoFoam has no VJP to benchmark
• ns-3d-grid · gradient: standard icoFoam is non-differentiable (C++, no AD path); DAFoam/OpenFOAM-AD exist but are not deployed in this tesseract
• ns-3d-grid · optimization: standard icoFoam is a non-differentiable, forward-only solver
• ns-grid · forward/kolmogorov_laminar: this tesseract runs icoFoam, which has no momentum source term, so it cannot apply the body force or the drag. OpenFOAM’s foamRun incompressibleFluid solver supports both through an fvModel, but this tesseract does not use it
• ns-grid · forward/kolmogorov: this tesseract runs icoFoam, which has no momentum source term, so it cannot apply the body force or the drag. OpenFOAM’s foamRun incompressibleFluid solver supports both through an fvModel, but this tesseract does not use it
• ns-grid · cost/vjp_cost: OpenFOAM has no VJP to benchmark
• ns-grid · gradient: OpenFOAM is non-differentiable (C++, no AD path); DAFoam/OpenFOAM-AD exist but are not deployed in this tesseract
• ns-grid · optimization: OpenFOAM is a non-differentiable, forward-only solver

Effort: Periodic box (2D) L0 · Channel flow past a cylinder L0 · Triply periodic box (3D) none — 6 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: application. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
momentum-tolerance (tuned): The periodic icoFoam case solves the momentum predictor to an absolute tolerance of 1e-8 instead of the cavity tutorial’s 1e-5. The 3-D periodic case uses the same file.	L0	required	Periodic box (2D) (inaccurate-forward): At 1e-5 the 2-D Taylor-Green error stalls at 2.7e-4 between N = 64 and N = 128, which made forward/baseline borderline. At 1e-8 it reaches 5.9e-5 at N = 128. Evidence: #222, table of U tolerances from 1e-5 to 1e-9 (merged through #232); 3-D results were unchanged.	—
periodic-pressure-solver (tuned): The periodic boxes solve the pressure with GAMG and the DICGaussSeidel smoother instead of the cavity tutorial’s PCG with DIC, at the tutorial’s tolerances (1e-6, relTol 0.05).	L0	choice	Periodic box (2D) (too-slow): With PCG/DIC, spatial_cost at N = 256 takes 37.4 s against 16.4 s, and temporal_cost at 1000 steps 44.7 s against 27.8 s. Every status is unchanged; forward errors move in the fifth significant digit. Evidence: #306, before/after table. Triply periodic box (3D): Cost moves by less than 8% either way (spatial_cost at N = 64: 28.7 s with PCG/DIC against 27.1 s); the 3-D box shares the 2-D box’s fvSolution. Evidence: #306, before/after table.	—
periodic-mpi-decomposition (tuned): The periodic boxes are decomposed with scotch and run under mpirun, one rank per physical core and at least 10,000 cells per rank. Smaller cases run serially, as does the channel.	L0	choice	Periodic box (2D) (too-slow): Serial, spatial_cost at N = 256 takes 9.90 s against 6.39 s decomposed. Statuses are unchanged. With the CPUs loaded by other solvers, decomposed runs were slower than serial ones (#306); the measurement assumes CI’s cost jobs run alone. Evidence: #314, before/after table; #306, notes. Triply periodic box (3D) (too-slow): Serial, spatial_cost at N = 64 takes 25.6 s against 12.9 s decomposed, and 100 steps at N = 48 take 18.2 s against 7.24 s. At 10 steps decomposing and reconstructing outweigh the solve (3.92 s serial against 4.33 s). On a loaded CPU, N = 64 took 150 to 280 s at 24 ranks against 29 s serial (#306). Evidence: #314, before/after table; #306, notes.	—
channel-momentum-tolerance: The cylinder channel solves the momentum predictor to 1e-8 instead of the tutorial’s 1e-5.	L0	choice	Channel flow past a cylinder: None measured: at 1e-5, forward/cylinder gives the same Cd to five digits with no speed-up. drag_opt was not run at 1e-5. Evidence: #269, notes.	—
channel-non-orthogonal-corrector: One non-orthogonal corrector (nNonOrthogonalCorrectors 1) on the body-fitted O-grid, where the offsetCylinder tutorial uses 0. The periodic case uses 0, as its tutorial does.	L0	required	Channel flow past a cylinder (crash): Without it, the run at 20 cells per diameter fails (maximum Courant number 11.9). Evidence: #269, before/after table.	—
graded-channel-mesh (tuned): The O-grid cells shrink radially to half the size of the box cells at the cylinder wall, so the smallest cell is about 0.35 dx, and the wake block coarsens to 10 dx at the outlet.	L0	choice	Channel flow past a cylinder (too-slow): A uniform mesh takes 111 s for forward/cylinder against 82 s, both measured with the tighter pressure tolerance and MPI decomposition the channel used before #269. Its Cd (5.567318 and 5.578718 at 10 and 20 cells per diameter) is closer to the reference 5.5795 than the graded mesh’s (5.563948 and 5.576608). Evidence: #269, before/after table; smallest cell from #263.	—

Inputs

Field	Type	∂	Description
v0	GridVectorField		Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]		Kinematic viscosity ν.
dt	Array[(1,), Float32]		Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None		Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.

Outputs

Field	Type	∂	Description
result	GridVectorField		Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None		Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

PhiFlow

Numerics: Finite Difference, 2nd-order central FD under incompressible RK4 (periodic); semi-Lagrangian advection, explicit diffusion (channel); CG pressure projection

jax AD: hybrid GPU

Periodic boxes: 2nd-order central differences under phiflow’s incompressible RK4 with CG pressure projection. The cylinder channel: semi-Lagrangian advection with operator splitting, phiflow’s wind-tunnel boundaries and a staircase cylinder. Gradients via JAX autodiff with adjoint CG solves.

F2 periodic · 2/2 pass F2 cylinder · 2/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://tum-pbs.github.io/PhiFlow/

Path: mosaic/tesseracts/navier-stokes-grid/phiflow

Image: phiflow_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward pass, gradient pass; F2 Channel flow past a cylinder: forward pass, optimization pass; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Effort: Periodic box (2D) L3 · Channel flow past a cylinder L2 · Forced Kolmogorov flow (2D) L3 · Triply periodic box (3D) L3 — 9 interventions, 3 with unknown reason — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
int64-sparse-indices: Monkeypatches phiml’s private NumPyBackend.ravel_multi_index to widen flat indices to int64 when the index space exceeds int32.	L3	required	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D) (crash): phiml casts sparse-matrix flat indices back to int32; they wrap at 2-D N >= 216 and 3-D N >= 36, and building the pressure matrix fails. In the benchmark that failed 12 cost cells (ns-grid N = 256, ns-3d-grid N >= 48). With the patch they complete and other results are bitwise identical. Evidence: https://github.com/pasteurlabs/mosaic/issues/155; #232, table; docstring of _ravel_multi_index_int64.	https://github.com/tum-pbs/PhiML/pull/16 (open)
explicit-cg-solve: The channel’s pressure solve passes math.Solve(“CG”, …) instead of the “auto” method of make_incompressible’s default Solve(). The periodic boxes’ solve uses incompressible_rk4’s default method and tolerances.	L0	unknown	Channel flow past a cylinder: No reason recorded.	—
periodic-cg-rank-deficiency: The periodic boxes’ pressure solve passes rank_deficiency=0, so phiml runs CG on the singular pressure system, whose right-hand side make_incompressible balances to zero mean, instead of adding a Lagrange multiplier that makes the system indefinite.	L0	required	Forced Kolmogorov flow (2D) (non-finite, incorrect-gradient, too-slow): CG is not guaranteed to converge on an indefinite system, and under jit phiml returns an unconverged solve without raising. At 256², 852 of 6400 solves to t = 4 stopped unconverged (128 of them as diverged, with residuals up to 10⁴ times the right-hand side), and two of three states went non-finite on CI. The chain gradient was non-finite on CI and off by 12% at the best ε on the A100. With rank_deficiency=0 every solve converged, in about 375 iterations against 570, both checks passed, and forward/kolmogorov took 936 s against 5,585 s. Evidence: forward/kolmogorov and gradient/chain_fd_check of #272; PR description. Periodic box (2D) (too-slow): No check changed. The cost experiments got 1.3x to 7.5x faster, growing with N, and all four PhiFlow cost anomalies went away. Evidence: #272, PR description. Triply periodic box (3D) (too-slow): No check changed. The cost experiments got 1.4x to 16x faster, growing with N (15.7x for the VJP at N = 64). Evidence: #272, PR description.	—
periodic-adjoint-relative-tolerance: The periodic boxes’ adjoint pressure solve (the Solve’s gradient_solve) stops at a relative residual of 1e-5 only, instead of inheriting the forward solve’s absolute tolerance of 1e-5.	L0	choice	Forced Kolmogorov flow (2D) (incorrect-gradient): The adjoint solve’s right-hand side scales with the cotangent. In the chain gradient check, whose cotangents are small, the absolute tolerance stopped the adjoint solves after 0 to 17 iterations, against about 100 at the relative tolerance, and the gradient missed the finite differences by 3.6e-3 at the best ε, against 1.1e-4 with the relative tolerance; the peers miss by 7e-5 to 2e-3. The check passes either way. Evidence: gradient/chain_fd_check of Periodic box (2D), Triply periodic box (3D): Not measured separately; no check changed.	—
channel-cg-tolerance (tuned): The channel’s pressure CG stops at a relative and absolute residual of 1e-2 instead of phiml’s float32 default of 1e-5.	L0	choice	Channel flow past a cylinder (too-slow): At 1e-5, forward/cylinder took 293 s against 43 s and drag_opt 10,782 s against 1,271 s, on a shared GPU. Against 1e-7, Cd, Cl and the pressure drop move by at most 0.07%, and drag_opt reaches a drag reduction of 47.60% against 47.66%; every status is unchanged. At 1e-1 the Cd error moves by 4%. Tuned with the channel-pressure-warm-start correction in place, so the tolerance applies to the correction’s residual. Evidence: #308, notes.	—
rk4-central-differences: Periodic boxes advect and diffuse with 2nd-order central differences under fluid.incompressible_rk4 (pressure_order=2, four pressure projections per step) instead of the default semi-Lagrangian advection with operator splitting. The channel keeps the default, since incompressible_rk4 takes no obstacles.	L0	choice	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D) (inaccurate-forward): PhiFlow’s docs point users who want accurate periodic flow to this pairing and label semi-Lagrangian advection first order. Against the previous central-difference Euler step it changes errors by a few percent and costs 2.4 to 4.4x more, but stays stable where that step went NaN (2-D, N = 32, nu = 0, 800 steps). The semi-Lagrangian default was measured about 10x less accurate on the periodic box in the cleanup notes; that table is not published. pressure_order=2 follows the comparison notebook; the function default (4) was not measured. Evidence: #288; https://tum-pbs.github.io/PhiFlow/Fluid_Simulation.html (“Higher-Order Simulations”).	—
inlined-obstacle-projection: Writes the channel’s obstacle projection out instead of calling fluid.make_incompressible with an obstacle: face masks from the shared staircase solid_cells, solid faces zeroed after advection and after diffusion, a masked_laplace CG solve over the fluid cells and a masked green-gauss pressure gradient.	L2	unknown	Channel flow past a cylinder (cannot-express): No reason recorded. A plausible one, not checked: make_incompressible rasterises its own obstacle geometry, so it cannot take the staircase mask the other grid solvers share, and it does not return the momentum removed at solid faces, which body_force needs.	—
channel-pressure-warm-start: Solves each step’s channel pressure as a correction to the previous step’s pressure, from a zero initial guess, instead of passing the previous pressure as the solve’s x0.	L2	unknown	Channel flow past a cylinder (too-slow): Per the code comment, it keeps the CG iterations low in the near-steady flow; not measured for PhiFlow, and the channel’s CG tolerance was tuned with it in place. Passing the previous pressure as x0 fails under reverse mode with an UnexpectedTracerError inside lax.scan. Evidence: Comment in channel_fwd.step (#239).	https://github.com/tum-pbs/PhiML/issues/17 (x0 tracer error)
body-force-output: Computes the per-step force on the cylinder as the momentum the solid faces remove in a step (advection and diffusion) plus the pressure on the faces between fluid and solid cells.	L2	required	Channel flow past a cylinder (cannot-express): PhiFlow has no force-on-obstacle routine (checked for 3.2.0, not re-checked for the pinned 3.4.0), so without this there is no body_force output for forward/cylinder and drag_opt to score. Evidence: #239 makes body_force a required output.	—

Reason unknown. Nobody has measured whether these are needed; a measurement, or removing them, would settle it:

• explicit-cg-solve (L0): The channel’s pressure solve passes math.Solve(“CG”, …) instead of the “auto” method of make_incompressible’s default Solve(). The periodic boxes’ solve uses incompressible_rk4’s default method and tolerances.
• inlined-obstacle-projection (L2): Writes the channel’s obstacle projection out instead of calling fluid.make_incompressible with an obstacle: face masks from the shared staircase solid_cells, solid faces zeroed after advection and after diffusion, a masked_laplace CG solve over the fluid cells and a masked green-gauss pressure gradient.
• channel-pressure-warm-start (L2): Solves each step’s channel pressure as a correction to the previous step’s pressure, from a zero initial guess, instead of passing the previous pressure as the solve’s x0.

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]	✓	Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None	✓	Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None	✓	Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

PICT

Numerics: Finite Volume, PISO, BDF1

pytorch AD: autodiff GPU

GPU-accelerated differentiable PISO solver built on PyTorch with custom CUDA kernels (PISOtorch). VJP w.r.t. v0, inflow_profile, and viscosity flows through the PISO time loop via PISOtorch_diff (tracks VELOCITY, BOUNDARY_VELOCITY, and VISCOSITY gradients). dt is a Python scalar consumed at construction time; zero gradients are returned for it. The VJP scales the cotangents to a largest entry of 1 before PICT’s backward pass and scales the gradients back, because PICT’s linear solves stop at an absolute residual and would cut the adjoint solves short for small cotangents. PICT has no releases; this runs commit 72a6e86.

F2 periodic · 1/2 pass F2 cylinder · 1/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://github.com/tum-pbs/PICT

Path: mosaic/tesseracts/navier-stokes-grid/pict

Image: pict_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward borderline, gradient pass; F2 Channel flow past a cylinder: forward pass, optimization blocked; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Excluded experiments

• ns-grid · optimization/drag_opt: 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

Effort: Periodic box (2D) L0 · Channel flow past a cylinder L2 · Forced Kolmogorov flow (2D) L0 · Triply periodic box (3D) none — 9 interventions, 4 with unknown reason — 1 setup deviation

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
periodic-solve-tolerance (tuned): The 2-D periodic runs pass 1e-6 as the advection and pressure solve tolerance instead of PICT’s float32 default of 1e-5.	L0	required	Periodic box (2D), Forced Kolmogorov flow (2D) (incorrect-gradient): At 1e-5 the residual error of the forward and adjoint solves caps fd_check at about 1e-3 for every step size: the reported error is 8.8e-4 against 5.7e-5 at 1e-6, which becomes an anomaly against the peer median after #282. Measured without cotangent-scaling, which alone lowers the error at 1e-5 to 2.7e-4 (#294); the two together are not measured. Evidence: #287, before/after table; #294, A/B table.	—
periodic-3d-solve-tolerance: The 3-D periodic box passes the same 1e-6 advection and pressure solve tolerance, through the same code path as the 2-D box.	L0	unknown	Triply periodic box (3D): Not measured on ns-3d-grid: #287 compared tolerances on ns-grid only.	—
cotangent-scaling: The VJP scales the cotangents to a largest entry of 1 before PICT’s backward pass and scales the gradients back, so that the absolute tolerance of PICT’s adjoint solves acts relative to the cotangent.	L2	choice	Periodic box (2D), Forced Kolmogorov flow (2D) (incorrect-gradient): PICT’s adjoint solves stop at an absolute residual, so its VJP is not linear in the cotangent; single calls with O(1) cotangents, which the gradient suite checks, are barely affected. ns-grid fd_check at eps = 1 goes from 1.1e-3 to 2.7e-4 with it; every status stays the same. Through a chain of 16 calls with cotangents of about 1e-4, the gradient was off by up to 73% without it. Measured at the earlier solve tolerance of 1e-5, before periodic-solve-tolerance. Evidence: #294, A/B table; tests/test_pict_vjp.py. Triply periodic box (3D): Statuses unchanged, and fd_check moves within its spread. The recovery’s final IC error is 0.0022 to 0.0040 with it against 0.0013 to 0.0021 without, both far below the failure threshold; #294 has no verified explanation. Evidence: #294, A/B table.	None filed. A relative convergence criterion would need a change to PISOtorch, whose criteria are all absolute.
warm-start-pressure (tuned): Subclasses PISOtorch_simulation.Simulation and overrides linear_solve so that each channel pressure CG solve starts from the same corrector’s pressure one step earlier and stops at 1e-8. In differentiable mode the solve is a torch.autograd.Function on the private PISOtorch_diff._linear_solve_wrapper, with a warm-started adjoint.	L3	choice	Channel flow past a cylinder (too-slow): Without it the channel forward is 2.5x (10 cells per diameter) to 4.8x (20 cells per diameter) slower; with it Cd stays within 0.05%. The code comment gives 2e-5 against a solve from zero at 1e-8, and at 1e-6 Cd at 20 cells per diameter still drifts by 1e-3 over the final 40% of the run. No running experiment exercises the differentiable half: PICT is excluded from drag_opt. Evidence: #239, summary and notes; comment above pressure_tol in _run_channel.	None filed.
pressure-tol-flux-gate: Sets the channel’s pressure_tol to max(1e-6, 2 x the boundary flux imbalance). With the warm start active it only feeds the check at the start of Simulation.run(), which refuses to run when the imbalance exceeds it.	L1	unknown	Channel flow past a cylinder (crash): Per the code comment, the float32 imbalance after balancing is about 4e-8 of the flux (1.4e-6 at 40 cells per diameter), which a fixed 1e-6 would fail. PICT’s float32 default of 1e-5 may already be above it at the benchmark’s 10 and 20 cells per diameter; not checked. Evidence: Comment above pressure_tol in _run_channel.	—
orthogonal-scheme-on-o-grid: Runs the channel with non_orthogonal=False, which drops PICT’s skew correction although the O-grid cells near the box corners are skewed.	L0	required	Channel flow past a cylinder (crash): With the default non_orthogonal=True, CG and float32 BiCGStab break down in the first pressure solve. Whether dropping the correction contributes to PICT’s Cd sitting about 3% high has not been tested. Evidence: #239, notes; comment at the _simulate call in _run_channel.	https://github.com/tum-pbs/PICT/issues/9
channel-in-cell-units: Builds the channel blocks in units of the canonical cell size dx (vertices, viscosity and dt divided by dx), so that PICT’s absolute solver residual is measured on unit cells.	L0	unknown	Channel flow past a cylinder (inaccurate-forward): Not measured. In physical units the same absolute tolerance would mean a different relative accuracy. The periodic box gets unit cells from PICT’s CreateBlockWithSize without asking. Evidence: Header comment of the channel section in tesseract_api.py.	—
inflow-in-autograd-graph: Re-applies the inflow profile to the inlet faces in a PRE callback every step instead of once before the run.	L0	unknown	Channel flow past a cylinder (incorrect-gradient): Per the code comment, a profile set once is not in the autograd graph that PISOtorch_diff records, so the inflow_profile gradient would come back as zeros. Not reproduced; only drag_opt uses this gradient, and PICT is excluded from it. Evidence: Comment in _pre_step in _run_channel.	—
body-force-integration: Computes the force on the cylinder from the cells next to the wall: wall pressure and shear from second-order one-sided fits to the first two cells, integrated over the arc faces in a POST callback every step.	L2	required	Channel flow past a cylinder (cannot-express): PICT has no force or drag output, so body_force, which forward/cylinder scores, could not be returned. #239 reports PICT’s control-volume Cd 7 to 11% off this surface Cd (about 1% for the other solvers); not investigated. Evidence: #239; PISOtorch_simulation.py at PICT 72a6e86 has no force routine.	—

Setup deviations

• advective-outflow (Channel flow past a cylinder): PICT’s advective outflow replaces the prescribed p = 0 outlet: outlet faces start from the inlet profile, boundary fluxes are balanced, the outlet velocity is advected with the mean inflow speed every step, pressure is Neumann there, and the returned pressure is shifted to zero on the outlet face. Effect: Partly measured. In OpenFOAM, an advective velocity outlet with p = 0 leaves Cd, Cl and pressure drop unchanged to six digits, as expected for a steady flow; the missing pressure outlet could not be reproduced there. Its share of PICT’s 3% Cd excess (#239) is unmeasured. Evidence: PICT has no pressure outlet (comment in _run_channel); OpenFOAM A/B of 2026-10-29 in the intervention cleanup notes (policy 3), not yet published.

Reason unknown. Nobody has measured whether these are needed; a measurement, or removing them, would settle it:

• periodic-3d-solve-tolerance (L0): The 3-D periodic box passes the same 1e-6 advection and pressure solve tolerance, through the same code path as the 2-D box.
• pressure-tol-flux-gate (L1): Sets the channel’s pressure_tol to max(1e-6, 2 x the boundary flux imbalance). With the warm start active it only feeds the check at the start of Simulation.run(), which refuses to run when the imbalance exceeds it.
• channel-in-cell-units (L0): Builds the channel blocks in units of the canonical cell size dx (vertices, viscosity and dt divided by dx), so that PICT’s absolute solver residual is measured on unit cells.
• inflow-in-autograd-graph (L0): Re-applies the inflow profile to the inlet faces in a PRE callback every step instead of once before the run.

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]		Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None	✓	Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None	✓	Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

Warp-NS

Numerics: Finite Difference, Projection, Euler

warp AD: autodiff GPU

Chorin projection scheme with explicit Euler time stepping, implemented in NVIDIA Warp CUDA kernels with FFT pressure Poisson on periodic domains. The cylinder channel uses a staggered grid, a staircase cylinder and a conjugate-gradient pressure Poisson solve. Gradients via wp.Tape kernel-level VJP, with adjoint CG solves.

F2 periodic · 2/2 pass F2 cylinder · 2/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://github.com/NVIDIA/warp

Path: mosaic/tesseracts/navier-stokes-grid/warp-ns

Image: warp_ns_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward pass, gradient pass; F2 Channel flow past a cylinder: forward pass, optimization pass; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Excluded experiments

• ns-grid · gradient/chain_fd_check: the solver advances with explicit Euler and central differences for advection, which is linearly unstable at the timestep of the chain (see forward/kolmogorov). Over the chain’s calls the unstable modes grow, so they dominate both the gradient and the finite differences, which then disagree at every ε a float32 run resolves. Where they grow fast enough, the chain from the ground-truth state diverges

Effort: Periodic box (2D) L2 · Channel flow past a cylinder L2 · Forced Kolmogorov flow (2D) L2 · Triply periodic box (3D) L2 — 6 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: kernel-toolkit. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
maintainer-advised-rewrite: The periodic solver was restructured following a Warp developer’s advice: native tape adjoints instead of hand-written NumPy adjoints, an on-device tile-FFT Poisson solve, fused pack and unpack kernels, and separate buffer policies for the forward-only and taped paths.	L0	choice	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D) (too-slow): The earlier version was 15 to 30x slower in forward and VJP, measured by direct apply and VJP calls on an A100. Evidence: #132; https://github.com/NVIDIA/warp/discussions/1636.	—
bluestein-non-pow2-fft: Grids whose size is not a power of two run their spectral Poisson solve through a Bluestein chirp-z transform built from power-of-two wp.tile_fft calls. Power-of-two grids keep the native path. The transform is exact to float32 round-off.	L2	required	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D) (crash): wp.tile_fft (cuFFTDx) only compiles power-of-two lengths on the GPU, so the benchmark’s non-power-of-two sizes (N = 192 in 2-D, N = 48 in 3-D) failed at compile time. Evidence: #132; section comment in tesseract_api.py. The scaling of this path was fixed in #232; tests/test_warp_ns_fft.py checks it against NumPy.	Not filed. Warp 1.15.0’s tile_fft documentation lists any size as supported on GPU builds with libmathdx, which is the default, so the failure in #132 contradicts it; it has not been reduced to a standalone reproducer.
channel-cg-adjoint: The backward pass of each channel pressure solve is a hand-written adjoint registered with tape.record_func: a CG solve with the same symmetric operator, started from the previous adjoint solution.	L2	required	Channel flow past a cylinder (no-gradient): The pressure solve (warp.optim.linear.cg, captured as a CUDA graph) is not a Warp kernel, so wp.Tape cannot differentiate it. No failing run is recorded. Evidence: Docstring of channel_solve (#239); tests/test_warp_ns_cylinder.py checks the drag gradient against central finite differences.	—
channel-cg-preconditioner (tuned): Preconditions the channel’s pressure CG with the exact inverse of the same operator without the cylinder, applied through dense 1-D eigenbases along x and y.	L2	choice	Channel flow past a cylinder (too-slow): CG needs about 20 iterations instead of 150 at 10 cells per diameter. Run time halves at 10 cells per diameter (drag_opt), breaks even at 20 (forward/cylinder) and doubles at 40. Measured at the earlier CG tolerance of 1e-5, not at the current 1e-3. Evidence: Docstring of _channel_pressure_solver (#239).	—
channel-cg-graph-capture: On CUDA the channel CG iterations run as a captured conditional graph, replayed per solution buffer, with non-convergence flagged on the device and raised after the run.	L0	choice	Channel flow past a cylinder (too-slow): Avoids a host sync per residual check, per the docstring. Not measured.	—
channel-cg-tolerance (tuned): The channel’s pressure CG stops at a relative residual of 1e-3 (absolute tolerance 0), with a cap of 5000 iterations.	L0	choice	Channel flow past a cylinder (too-slow): At the earlier 1e-5, forward/cylinder took 21.6 s against 15.2 s and drag_opt 407 s against 189 s, on a shared GPU. Against 1e-7, Cd moves by 0.001%, Cl by 0.25% and the pressure drop by 0.015%, and drag_opt’s drag reduction by 0.01 points; every status is unchanged. At 1e-2 Cl moves by 5%. Evidence: #308, notes.	—

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]	✓	Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None	✓	Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None	✓	Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Navier–Stokes (Grid)

XLB

Numerics: Lattice Boltzmann, Streaming

jax AD: autodiff GPU

JAX-accelerated lattice Boltzmann solver (D2Q9 / D3Q27) with BGK collision, running the streaming/collision loop on GPU. Gradients via source-transformation AD; recovers incompressible Navier–Stokes only to O(Ma²) via Chapman–Enskog expansion.

F2 periodic · 2/2 pass F2 cylinder · 2/2 pass F3 periodic · 3/3 pass Details

Upstream docs: https://github.com/Autodesk/XLB

Path: mosaic/tesseracts/navier-stokes-grid/xlb

Image: xlb_navier_stokes_grid

Benchmark checks: F2 Periodic box (Taylor–Green, multimode): forward pass, gradient pass; F2 Channel flow past a cylinder: forward pass, optimization pass; F3 Triply periodic box (Taylor–Green): forward pass, gradient pass, optimization pass — see Results.

Effort: Periodic box (2D) L2 · Channel flow past a cylinder none · Forced Kolmogorov flow (2D) L2 · Triply periodic box (3D) L2 — 9 interventions, 5 with unknown reason — 1 setup deviation

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
float64-everywhere: Enables jax_enable_x64 and runs apply() and the VJP with XLB’s FP64FP64 precision policy, casting outputs to float32. XLB_VJP_FP32=1 switches everything to float32.	L0	unknown	Periodic box (2D), Forced Kolmogorov flow (2D), Channel flow past a cylinder, Triply periodic box (3D) (incorrect-gradient): The code comment says float32 cancellation near omega = 2 corrupts gradients and that float32 gradients fail gradient/fd_check. No float32 fd_check result is recorded, and the comment’s reasoning predates the switch to BGK everywhere (#282). Evidence: Comment above _USE_F64_DIFF; https://github.com/pasteurlabs/mosaic/issues/125 (open) proposes measuring float32 against fd_check.	—
xlb-init-per-velocity-set: Calls xlb.init() again with each operator bundle’s velocity set and precision policy before constructing its operators, because XLB keeps both in global state.	L1	unknown	Periodic box (2D), Forced Kolmogorov flow (2D), Channel flow past a cylinder, Triply periodic box (3D) (crash): The recorded failure, a D3Q27 bundle built after xlb.init() with D2Q9 failing to compile, involved the KBC operator, which #282 removed. Whether the remaining operators need it has not been checked. Evidence: Docstring of _make_ops at 6517f75.	None filed.
sqrt-checkpointing: Runs the time loop as nested jax.lax.scan blocks of about sqrt(steps) steps, with jax.checkpoint on each block and on each step inside it. Reverse mode keeps the block-boundary populations and, while it works through a block, that block’s step inputs.	L0	required	Triply periodic box (3D) (out-of-memory): #135 gives 110 GB at N = 64 for a plain scan, against a 0.3 GB forward. Checkpointing each step but not the blocks still keeps every step’s populations: on an A100 80 GB, the 3-D TGV VJP at N = 48 took 67 GB at 1000 steps and ran out of memory at 2000, against 5.6 GB at both with the blocks checkpointed. At 50 steps, N = 96 took 69 GB against 12 GB, and N = 128 ran out of memory against 38 GB. The block checkpoint makes the VJP 2-14% slower and changes the gradient by at most 5e-9 relative. Checkpointing the blocks without the steps inside them took about twice the memory and was no faster. Evidence: #135 (plain scan); #329 (A100, one container per point). Periodic box (2D), Forced Kolmogorov flow (2D) (too-slow): #135 says it fixes some of XLB’s VJP slowness in 2-D. No number is recorded, and the block checkpoint has not been timed in 2-D. Evidence: #135; https://github.com/pasteurlabs/mosaic/issues/125.	—
channel-checkpointing: The channel runs its time loop through the same block- and step-checkpointed scan as the periodic boxes.	L0	unknown	Channel flow past a cylinder: No reason recorded: the channel reuses the periodic boxes’ helper since #239, and no run without it is recorded. drag_opt differentiates through the whole run.	—
consistent-initial-populations: Periodic runs start from populations consistent with v0: density from a spectral pressure-Poisson solve plus the first-order Chapman-Enskog non-equilibrium part (Mei et al. 2006), instead of the equilibrium at density 1.	L2	required	Periodic box (2D), Forced Kolmogorov flow (2D) (inaccurate-forward): Equilibrium-only populations leave out pressure and viscous stress, which appear as start-up transients. On an offline copy of the solver with sub-stepping, the Taylor-Green error was 0.031 to 0.035 without them and 6e-4 to 8e-4 with them, against a peer median of 2.4e-3. Evidence: https://github.com/pasteurlabs/mosaic/issues/38; #232. Triply periodic box (3D) (inaccurate-forward): The transient made the N = 8 forward/baseline error 0.273, 80x the peers. On the offline copy the 3-D baselines became 0.032, 0.0031 and 0.0025 at N = 8, 16 and 32. Evidence: https://github.com/pasteurlabs/mosaic/issues/38; #232.	—
manual-lbm-loop: The periodic time loop is assembled from XLB’s Stream, Macroscopic, QuadraticEquilibrium and BGK operators inside jax.lax.scan instead of using XLB’s stepper.	L2	unknown	Periodic box (2D), Forced Kolmogorov flow (2D), Triply periodic box (3D) (too-slow): The initial tesseract description gave “O(1) compile time” as the reason. Nothing was measured, and whether XLB’s stepper works under scan and reverse mode was not recorded. Evidence: Initial commit (c47498a) description.	—
exact-difference-forcing: Periodic runs with a body force or drag apply the exact difference method by hand in the manual LBM loop, with the drag evaluated semi-implicitly at the half-force velocity. The initial populations and the output velocity are shifted by half a kick.	L2	required	Forced Kolmogorov flow (2D) (cannot-express): XLB 0.3.1’s ExactDifference operator takes one constant force vector (“currently we are limited to a single force vector not a spatially dependent forcing field”), so it cannot apply the Kolmogorov force sin(κy) or a drag that depends on the velocity. Evidence: https://github.com/Autodesk/XLB/blob/v0.3.1/xlb/operator/force/exact_difference_force.py; tests/test_ns_grid_forcing.py; #272.	—
channel-boundary-conditions: The channel is written by hand on XLB’s operators: boundary link masks for the staircase cylinder, halfway bounce-back on walls and cylinder, velocity bounce-back at the inlet, an anti-bounce-back outlet, the momentum-exchange force, and its own sub-stepped, checkpointed time loop.	L2	choice	Channel flow past a cylinder: No reason recorded. XLB 0.3.1 ships halfway bounce-back, Zou-He, regularized and extrapolation-outflow boundary conditions and a momentum-transfer force operator; whether they work under this tesseract’s scan and reverse mode was not checked. Evidence: https://github.com/Autodesk/XLB/tree/v0.3.1/xlb/operator/boundary_condition	—
channel-characteristic-outlet: The outlet density follows the outgoing characteristic, 1 + (u_x - u_developed)/cs, relative to the fully developed profile with the inlet’s flux, so that p = 0 on the outlet face once the flow is steady.	L2	unknown	Channel flow past a cylinder (inaccurate-forward): Per the docstring, a plain p = 0 face reflects the start-up pressure waves, which then ring between inlet and outlet almost undamped at nu = 1e-3. Not measured. Evidence: Docstring of _channel_run (#239).	—

Setup deviations

• ma-targeted-substepping (Periodic box (2D), Forced Kolmogorov flow (2D), Channel flow past a cylinder, Triply periodic box (3D)): Each prescribed step runs as k = ceil(Ma / Ma_target) LBM steps of dt / k, never coarser than prescribed. Ma_target is 0.1, lowered in proportion to 1/N above N = 64 on periodic boxes. With a body force, Ma uses the larger of max|v0| and a bound on the speed the force can drive, so a flow it spins up from rest stays within the target. 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. Evidence: https://github.com/pasteurlabs/mosaic/issues/38; #232.

Reason unknown. Nobody has measured whether these are needed; a measurement, or removing them, would settle it:

• float64-everywhere (L0): Enables jax_enable_x64 and runs apply() and the VJP with XLB’s FP64FP64 precision policy, casting outputs to float32. XLB_VJP_FP32=1 switches everything to float32.
• xlb-init-per-velocity-set (L1): Calls xlb.init() again with each operator bundle’s velocity set and precision policy before constructing its operators, because XLB keeps both in global state.
• channel-checkpointing (L0): The channel runs its time loop through the same block- and step-checkpointed scan as the periodic boxes.
• manual-lbm-loop (L2): The periodic time loop is assembled from XLB’s Stream, Macroscopic, QuadraticEquilibrium and BGK operators inside jax.lax.scan instead of using XLB’s stepper.
• channel-characteristic-outlet (L2): The outlet density follows the outgoing characteristic, 1 + (u_x - u_developed)/cs, relative to the fully developed profile with the inlet’s flux, so that p = 0 on the outlet face once the flow is steady.

Inputs

Field	Type	∂	Description
v0	GridVectorField	✓	Initial velocity at the cell centers, shape (nx, ny, 1, 2) in 2-D or (nx, ny, nz, 3) in 3-D. Cell (i, j) is centered at ((i + 1/2) dx, (j + 1/2) dx). Default: 64×64 divergence-free random vortex field (seed=42).
viscosity	Array[(1,), Float32]	✓	Kinematic viscosity ν.
dt	Array[(1,), Float32]	✓	Timestep size.
steps	int		Number of timesteps. Total simulated time = steps × dt.
domain_extent	float | list[float]		Physical size of the domain. A scalar L is a square (cube) of side L, and a list gives the length per axis, [Lx, Ly] or [Lx, Ly, Lz]. Periodic boxes take a scalar. The grid spacing must be the same along every axis.
boundary_conditions	GridBC		Per-face boundary conditions. Default: all periodic.
obstacle	GridObstacle | None		Optional solid obstacle embedded in the domain (no-slip walls), in physical coordinates.
inflow_profile	Array[(None,), Float32] | None	✓	Inlet velocity u_x(y) at the x_lo face, sampled at the cell-center heights y_j = (j + 1/2) dx, shape (ny,). u_y is zero on the inlet. Requires x_lo to be ‘dirichlet’ with no value.
forcing	GridVectorField | None		Body force per unit mass (an acceleration) at the cell centers, constant in time, same shape as v0. 2-D periodic boxes only.
linear_drag	Array[(1,), Float32] | None		Coefficient α of a linear drag -αu on the right-hand side of the momentum equation. 2-D periodic boxes only.
snapshot_every	int | None		If set, also return the velocity after every snapshot_every steps as snapshots, so one run yields several times of a trajectory. steps must be a multiple of it. 2-D periodic boxes only.

Outputs

Field	Type	∂	Description
result	GridVectorField	✓	Final velocity at the cell centers, same shape as v0.
pressure	GridScalarField | None		Final kinematic pressure p/ρ at the cell centers, shape (nx, ny, nz). Values inside the obstacle are ignored.
body_force	Array[(None, None), Float32] | None	✓	Force of the fluid on the obstacle per unit depth and density, shape (steps, ndim): row k is (F_x, F_y) at t = (k + 1) dt. Positive along +x and +y, so a flow in +x gives F_x > 0.
snapshots	Array[(None, None, None, None, None), Float32] | None		Velocity at the cell centers after every snapshot_every steps, shape (steps / snapshot_every, nx, ny, 1, 2); the last equals result. Returned only when the inputs set snapshot_every. Not differentiable.

Structural mechanics

Structural mechanics

deal.II

Numerics: Finite Element, Direct (UMFPACK)

cpp AD: forward-only

deal.II Q1 trilinear hexahedral FEM linear-elasticity solver via C++ subprocess (UMFPACK direct). Forward-only — reference implementation for cross-framework validation.

S cantilever · 1/1 pass Details

Upstream docs: https://www.dealii.org/

Path: mosaic/tesseracts/structural-mesh/dealii-structural

Image: dealii_structural_mesh

Benchmark checks: S Linear-elastic cantilever: forward pass, gradient not supported, optimization not supported — see Results.

Excluded experiments

• structural-mesh · cost/vjp_cost: deal.II structural tesseract is forward-only: this tesseract exposes no AD path (no VJP endpoint, no differentiable inputs), so cost/vjp and gradient experiments are categorically out of scope.
• structural-mesh · gradient: deal.II structural tesseract is forward-only: this tesseract exposes no AD path (no VJP endpoint, no differentiable inputs), so cost/vjp and gradient experiments are categorically out of scope.
• structural-mesh · optimization: deal.II structural tesseract has no AD path; topology-optimisation (compliance gradient-based ρ updates) requires a differentiable compliance output.

Effort: Linear-elastic cantilever none — 1 intervention — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
umfpack-direct-solve (tuned): Solves the assembled system with deal.II’s SparseDirectUMFPACK instead of the CG with SSOR preconditioning of the step-8 tutorial the assembly follows.	L0	choice	Linear-elastic cantilever (too-slow): CG with SSOR gives the same compliance but is 1.7x slower at nx = 128. Evidence: #271, notes.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]		Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. Values in [0, 1]: 0 = void, 1 = solid.
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe displacement (values shape (n_groups, 3), None → zero); neumann.mask/values prescribe surface traction (values shape (n_groups, 3)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
E_max	float		Young’s modulus of the fully solid material [MPa].
nu	float		Poisson’s ratio.
xmin	float		Void stiffness ratio (E_min = xmin * E_max).
penal	float		SIMP penalisation exponent (E(ρ) = E_min + (E_max−E_min)·ρ^penal).

Outputs

Field	Type	∂	Description
compliance	Array[(), Float32]		Structural compliance C = F^T U (scalar). Minimising compliance is equivalent to maximising stiffness.

Structural mechanics

FEniCS

Numerics: Finite Element, Direct (MUMPS)

fenics AD: adjoint

FEniCS/DOLFIN CG1 FEM linear-elasticity solver with SIMP penalisation. Gradient ∂C/∂ρ via dolfin-adjoint tape replay (discrete adjoint derived at the variational level via UFL).

S cantilever · 3/3 pass Details

Upstream docs: https://fenicsproject.org/

Path: mosaic/tesseracts/structural-mesh/fenics-structural

Image: fenics_structural_mesh

Benchmark checks: S Linear-elastic cantilever: forward pass, gradient pass, optimization pass — see Results.

Effort: Linear-elastic cantilever none — 2 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
mumps-forward-solver (tuned): The forward solve passes solver_parameters={“linear_solver”: “mumps”} instead of using DOLFIN’s default LU backend (UMFPACK).	L0	choice	Linear-elastic cantilever (too-slow): At nx = 128, apply took 10.1 s before and 8.8 s after the change that introduced it, about 1.15x. Instrumentation put 99.6% of apply in the linear solve, so the gain is attributed to the solver switch. No failure without it is recorded on this setup. Evidence: #182, section 2 (measured together with tape replay).	—
mumps-adjoint-monkeypatch: Replaces fenics_adjoint.types.compat.linalg_solve at import time with a function that always solves with MUMPS, because legacy dolfin-adjoint does not pass the forward solver_parameters to the adjoint solve, which then runs UMFPACK.	L3	choice	Linear-elastic cantilever (too-slow): Without the patch the VJP takes 14.6 s against 11.1 s at nx = 128, and 4.7 s against 3.2 s at nx = 64, with identical gradients. Small meshes show no difference. Evidence: A/B on an A100, 2026-10-26, in the intervention cleanup notes (item 3), not yet published.	Not filed: legacy dolfin-adjoint has had no release since 2019.1.0. Per #182, newer dolfin-adjoint passes linear_solver to the adjoint solve.

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. Values in [0, 1]: 0 = void, 1 = solid.
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe displacement (values shape (n_groups, 3), None → zero); neumann.mask/values prescribe surface traction (values shape (n_groups, 3)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
E_max	float		Young’s modulus of the fully solid material.
nu	float		Poisson’s ratio.
xmin	float		Void stiffness ratio (E_min = xmin * E_max).
penal	float		SIMP penalisation exponent (E(rho) = E_min + (E_max-E_min)*rho^penal).

Outputs

Field	Type	∂	Description
compliance	Array[(), Float32]	✓	Structural compliance C = F^T U (scalar). Minimising compliance is equivalent to maximising stiffness.

Structural mechanics

Firedrake

Numerics: Finite Element, Direct (MUMPS)

firedrake AD: adjoint

Firedrake CG1 FEM linear-elasticity solver with SIMP penalisation. Gradient ∂C/∂ρ via firedrake-adjoint (pyadjoint ReducedFunctional) — an independent tape-based discrete adjoint complement to FEniCS.

S cantilever · 3/3 pass Details

Upstream docs: https://www.firedrakeproject.org/

Path: mosaic/tesseracts/structural-mesh/firedrake-structural

Image: firedrake_structural_mesh

Benchmark checks: S Linear-elastic cantilever: forward pass, gradient pass, optimization pass — see Results.

Effort: nothing beyond ordinary use is recorded

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

No interventions beyond ordinary use.

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. Values in [0, 1]: 0 = void, 1 = solid.
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe displacement (values shape (n_groups, 3), None → zero); neumann.mask/values prescribe surface traction (values shape (n_groups, 3)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
E_max	float		Young’s modulus of the fully solid material.
nu	float		Poisson’s ratio.
xmin	float		Void stiffness ratio (E_min = xmin * E_max).
penal	float		SIMP penalisation exponent (E(rho) = E_min + (E_max-E_min)*rho^penal).

Outputs

Field	Type	∂	Description
compliance	Array[(), Float32]	✓	Structural compliance C = F^T U (scalar). Minimising compliance is equivalent to maximising stiffness.

Structural mechanics

JAX-FEM

Numerics: Finite Element, BiCGSTAB+Jacobi

jax AD: hybrid GPU

JAX-FEM HEX8 linear-elasticity solver with SIMP stiffness E(ρ) = E_min + (E_max − E_min)·ρ^penal. Gradient ∂C/∂ρ via JAX automatic differentiation through the linear FE solve (hybrid implicit-function-theorem adjoint).

S cantilever · 3/3 pass Details

Upstream docs: https://github.com/deepmodeling/jax-fem

Path: mosaic/tesseracts/structural-mesh/jax-fem

Image: jax_fem_structural_mesh

Benchmark checks: S Linear-elastic cantilever: forward pass, gradient pass, optimization pass — see Results.

Effort: Linear-elastic cantilever none — 1 intervention — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
jitted-bicgstab-solver (tuned): Passes JAX-FEM’s own Jacobi-preconditioned BiCGSTAB (the default jax_solver’s algorithm and tolerances), compiled once per matrix size, to ad_wrapper through the custom_solver option for the forward and adjoint solves.	L2	choice	Linear-elastic cantilever (too-slow): The default jax_solver runs bicgstab eagerly, so it compiles 2 XLA programs per forward call and 4 per VJP, a fixed overhead on every call. The umfpack_solver option used before gives identical outputs but is much slower on large meshes. At N = 128, spatial_cost took 11.6 s with umfpack_solver and 1.49 s with this solver, vjp_cost/by_N 23.5 s and 2.91 s. Evidence: #286, before/after table; comment above _bicgstab.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. Values in [0, 1]: 0 = void, 1 = solid.
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe displacement (values shape (n_groups, 3), None → zero); neumann.mask/values prescribe surface traction (values shape (n_groups, 3)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.

Outputs

Field	Type	∂	Description
compliance	Array[(), Float32]	✓	Structural compliance C = F^T U (scalar). Minimising compliance is equivalent to maximising stiffness.

Structural mechanics

TopOpt.jl

Numerics: Finite Element, Direct (CHOLMOD)

julia AD: adjoint

TopOpt.jl HEX8 FEM linear-elasticity solver with SIMP penalisation, Julia + Ferrite.jl backend. Gradient ∂C/∂ρ via analytical SIMP adjoint sensitivity computed alongside the forward solve.

S cantilever · 3/3 pass Details

Upstream docs: https://github.com/JuliaTopOpt/TopOpt.jl

Path: mosaic/tesseracts/structural-mesh/topopt-jl

Image: topopt_jl_structural_grid

Benchmark checks: S Linear-elastic cantilever: forward pass, gradient pass, optimization pass — see Results.

Effort: Linear-elastic cantilever none — 1 intervention, 1 with unknown reason — 1 setup deviation

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
julia-cpu-target-generic: The image sets JULIA_CPU_TARGET=generic, so the package images Julia precompiles at build time target a generic x86-64 CPU rather than the build machine’s. Code compiled at run time is unaffected.	L0	unknown	Linear-elastic cantilever: No reason recorded; added in #6 (“ensure all Tesseracts build on CI”). Julia’s docs recommend it where images run on CPUs other than the one that built them. Whether the precompiled solver code is slower for it has not been measured.	—

Setup deviations

• nodal-load (Linear-elastic cantilever): The harness gives TopOpt.jl the cantilever load as nodal forces (CLOAD) instead of a traction: the shares that integrating the traction over the bilinear face elements gives each node. Effect: None: on the uniform mesh the nodal forces equal the integrated traction. Evidence: TopOpt.jl 0.14.0’s INP reader takes surface loads only as a normal pressure (DLOAD P), and the cantilever load is tangential; _traction_solvers in mosaic/benchmarks/problems/structural_mesh/physics.py.

Reason unknown. Nobody has measured whether these are needed; a measurement, or removing them, would settle it:

• julia-cpu-target-generic (L0): The image sets JULIA_CPU_TARGET=generic, so the package images Julia precompiles at build time target a generic x86-64 CPU rather than the build machine’s. Code compiled at run time is unaffected.

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. Values in [0, 1]: 0 = void, 1 = solid.
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe displacement (values shape (n_groups, 3), None → zero); neumann.mask/values prescribe surface traction (values shape (n_groups, 3)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
E	float		Young’s modulus of the solid material.
nu	float		Poisson’s ratio.
xmin	float		Minimum density (void stiffness) to prevent singular stiffness matrix.

Outputs

Field	Type	∂	Description
compliance	Array[(), Float32]	✓	Structural compliance C = F^T U (scalar). Minimising compliance is equivalent to maximising stiffness.

Structural mechanics

torch-fem

Numerics: Finite Element, CG+Jacobi

pytorch AD: hybrid GPU

torch-fem Solid HEX8 linear-elasticity solver with SIMP stiffness E(ρ) = E_min + (E_max − E_min)·ρ^p — the PyTorch backend in structural-mesh. Gradient ∂C/∂ρ via PyTorch autograd through the FE solve (implicit-function-theorem adjoint); GPU-native.

S cantilever · 3/3 pass Details

Upstream docs: https://github.com/meyer-nils/torch-fem

Path: mosaic/tesseracts/structural-mesh/torch-fem-structural

Image: torch_fem_structural_mesh

Benchmark checks: S Linear-elastic cantilever: forward pass, gradient pass, optimization pass — see Results.

Effort: Linear-elastic cantilever L0 — 5 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
gpu-device: With a GPU present, runs torch-fem on it instead of torch’s default CPU device.	L0	choice	Linear-elastic cantilever: Not compared with the CPU: the tesseract has run on the GPU since its first version.	—
global-default-device: Sets torch’s default device (and dtype, float64) process-wide, so the tensors torch-fem allocates internally land on the GPU with the mesh, material and load tensors.	L0	required	Linear-elastic cantilever (crash): torch-fem allocates auxiliary tensors without device=, so they land on the CPU and collide with the CUDA tensors. No error log is recorded. Evidence: Code comment; torch-fem’s benchmarks/utils.py at v0.12.1 and its current installation docs set the default device the same way.	—
gpu-cg-jacobi: On the GPU, passes method=“cg” and preconditioner=“jacobi” to model.solve for the forward and adjoint solves. On the CPU torch-fem chooses.	L0	choice	Linear-elastic cantilever (too-slow): torch-fem 0.12.1 solves systems below 10,000 DOFs directly with SciPy on the CPU, even when the matrix is on the GPU; larger ones already get CG with Jacobi without AmgX. So the setting keeps the smaller meshes on the device. No isolated measurement. Evidence: #188; resolve_method and resolve_library in torch-fem v0.12.1 src/torchfem/sparse.py.	—
cg-tolerance (tuned): The GPU CG stops at a relative residual of 1e-6 (stol) instead of torch-fem’s 1e-10.	L0	choice	Linear-elastic cantilever (too-slow): At 1e-10, spatial_cost at N = 128 took 0.435 s against 0.338 s and vjp_cost 0.717 s against 0.530 s, on a shared GPU. Every status is unchanged from 1e-12 to 1e-6; forward metrics are bit-identical and gradient checks move by less than 1e-9. Evidence: #301, notes.	—
deterministic-kernels: Enables torch.use_deterministic_algorithms(True, warn_only=True), so torch-fem’s CUDA scatter-adds in assembly sum in a fixed order.	L0	choice	Linear-elastic cantilever: Run-to-run differences were shown on the heat setup only; on this setup the setting is kept for consistency, with no measurable cost. Evidence: #200, notes.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. Values in [0, 1]: 0 = void, 1 = solid.
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe displacement (values shape (n_groups, 3), None → zero); neumann.mask/values prescribe surface traction (values shape (n_groups, 3)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
E_max	float		Young’s modulus of fully solid material (ρ = 1).
nu	float		Poisson’s ratio (density-independent).
xmin	float		Void stiffness ratio E_min / E_max.
p_exp	float		SIMP penalisation exponent p in E(ρ) = E_min + (E_max−E_min)·ρ^p.

Outputs

Field	Type	∂	Description
compliance	Array[(), Float32]	✓	Structural compliance C = F^T U (scalar). Minimising compliance is equivalent to maximising stiffness.

Heat transfer

Heat transfer

deal.II

Numerics: Finite Element, CG+AMG

cpp AD: forward-only

deal.II Q1 trilinear hexahedral FEM heat-conduction solver with SIMP conductivity (C++ subprocess, AMG-preconditioned CG). Forward-only — reference implementation for cross-framework validation.

H slab · 1/1 pass Details

Upstream docs: https://www.dealii.org/

Path: mosaic/tesseracts/thermal-mesh/dealii-heat

Image: dealii_heat_thermal_mesh

Benchmark checks: H Steady conduction in a slab: forward pass, gradient not supported, optimization not supported — see Results.

Excluded experiments

• thermal-mesh · cost/vjp_cost: deal.II thermal tesseract is forward-only: this tesseract exposes no AD path (no VJP endpoint, no differentiable inputs), so cost/vjp and gradient experiments are categorically out of scope.
• thermal-mesh · gradient: deal.II thermal tesseract is forward-only: this tesseract exposes no AD path (no VJP endpoint, no differentiable inputs), so cost/vjp and gradient experiments are categorically out of scope.
• thermal-mesh · optimization: deal.II thermal tesseract exposes no AD path; conductivity recovery (gradient-based ρ updates that minimize the misfit to a target temperature field) requires a differentiable identification-error output.

Effort: Steady conduction in a slab none — 2 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
amg-preconditioner (tuned): Preconditions the CG solve with Trilinos ML algebraic multigrid (PreconditionAMG, elliptic) instead of the SSOR of the step-5 tutorial the assembly follows.	L0	choice	Steady conduction in a slab (too-slow): With SSOR every forward metric is the same, and spatial_cost at nx = 1024 takes 61.1 s against 25.6 s (UMFPACK: 35.6 s). AMG wins from nx = 256 on; at nx = 16 its setup makes it slower, 0.020 s against 0.014 s. Evidence: #305, notes table.	—
cg-tolerance (tuned): The CG solve stops at a residual of 1e-8 times the right-hand-side norm, with a cap of 50000 iterations.	L0	choice	Steady conduction in a slab (too-slow): At the earlier 1e-12 every forward metric is identical and the largest mesh is slower: spatial_cost at nx = 1024 took 60.1 s against 51.1 s. Measured with SSOR; with AMG, #305 checked only that the forward metrics are unchanged at 1e-8. Evidence: #266, before/after table; #305, notes.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]		Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. ρ=1 → fully conducting, ρ=0 → insulating.
source	Array[(None,), Float32]		Per-element volumetric heat source (W/m³), shape (n_cells,). Active slice = source[:hex_mesh.n_faces]. Zero by default (topology-optimisation mode; no body heat source).
target_temperature	Array[(None,), Float32]		Per-node target temperature for the inverse (source-identification) objective, shape (n_nodes,). identification_error =
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe temperature (values shape (n_groups, 1), None → zero); neumann.mask/values prescribe heat flux (values shape (n_groups, 1)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
k_max	float		Maximum thermal conductivity (fully solid material).
p_exp	float		SIMP penalisation exponent p (k(ρ) = k_min + (k_max−k_min)·ρ^p).

Outputs

Field	Type	∂	Description
thermal_compliance	Array[(), Float32]		Total thermal compliance C = ∮_Γ_N q_n · T dΓ (scalar). Equals the work done by the Neumann heat flux on the temperature field.
identification_error	Array[(), Float32]		Source-identification objective:

Heat transfer

FEniCS

Numerics: Finite Element, Direct (MUMPS)

fenics AD: adjoint

FEniCS P1 FEM heat-conduction solver with SIMP penalisation. Gradients (∂C/∂ρ, ∂C/∂source) via dolfin-adjoint tape replay.

H slab · 3/3 pass Details

Upstream docs: https://fenicsproject.org/

Path: mosaic/tesseracts/thermal-mesh/fenics-heat

Image: fenics_heat_thermal_mesh

Benchmark checks: H Steady conduction in a slab: forward pass, gradient pass, optimization pass — see Results.

Effort: Steady conduction in a slab L3 — 3 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
mumps-forward-solver: The forward solve passes solver_parameters={“linear_solver”: “mumps”} instead of using DOLFIN’s default LU backend (UMFPACK).	L0	required	Steady conduction in a slab (out-of-memory): With UMFPACK the forward solve runs out of memory at the largest cost-sweep size (nx = 1024), which showed as an error on the results page. No peak memory was recorded. Evidence: https://github.com/pasteurlabs/mosaic/issues/180 (comment by the FEniCS developer who proposed the change); merged in #182.	—
mumps-adjoint-monkeypatch: Replaces fenics_adjoint.types.compat.linalg_solve at import time with a function that always solves with MUMPS, because legacy dolfin-adjoint does not pass the forward solver_parameters to the adjoint solve, which then runs UMFPACK.	L3	required	Steady conduction in a slab (out-of-memory): With MUMPS in the forward solve only, the gradient computation still fails at nx = 1024, because the adjoint solve runs UMFPACK. Evidence: https://github.com/pasteurlabs/mosaic/issues/180 (comment by the FEniCS developer who proposed the patch); merged in #182.	Not filed: legacy dolfin-adjoint has had no release since 2019.1.0. Per #182, newer dolfin-adjoint passes linear_solver to the adjoint solve; per #180, DOLFINx-adjoint takes adjoint PETSc options.
nodal-misfit-block: A custom pyadjoint Block puts the nodal identification error sum_v mask_v (T_v - T_target,v)^2 on the tape, with hand-written adjoint and tangent-linear components, so pyadjoint differentiates it together with the solve.	L2	required	Steady conduction in a slab (cannot-express, incorrect-gradient): Legacy FEniCS has no vertex quadrature on hexahedra, so the nodal sum cannot be written as a UFL integral. The earlier approximation, an L2 integral rescaled by the node count, gave gradients that did not match the reported value: gradient/source_fd_check flagged FEniCS at 1753x the peer median, and conductivity_recovery_bfgs reached a final/initial misfit of 0.02779 against 0.00020 with the nodal adjoint. Evidence: #167; #184, notes; comment above NodalSquaredErrorBlock.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. ρ=1 → fully conducting, ρ=0 → insulating.
source	Array[(None,), Float32]	✓	Per-element volumetric heat source (W/m³), shape (n_cells,). Active slice = source[:hex_mesh.n_faces]. Zero by default (topology-optimisation mode; no body heat source).
target_temperature	Array[(None,), Float32]		Per-node target temperature for the inverse (source-identification) objective, shape (n_nodes,). identification_error =
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe temperature (values shape (n_groups, 1), None → zero); neumann.mask/values prescribe heat flux (values shape (n_groups, 1)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
k_max	float		Maximum thermal conductivity (fully solid/conducting material).
p_exp	float		SIMP penalisation exponent p (k(ρ) = k_min + (k_max−k_min)·ρ^p).

Outputs

Field	Type	∂	Description
thermal_compliance	Array[(), Float32]	✓	Total thermal compliance C = ∮_Γ_N q_n · T dΓ (scalar). Equals the work done by the Neumann heat flux on the temperature field.
identification_error	Array[(), Float32]	✓	Source-identification objective:

Heat transfer

Firedrake

Numerics: Finite Element, Direct (MUMPS)

firedrake AD: adjoint

Firedrake P1 FEM heat-conduction solver with SIMP penalisation. Gradients (∂C/∂ρ, ∂C/∂source) via firedrake-adjoint, complementing the FEniCS reference.

H slab · 3/3 pass Details

Upstream docs: https://www.firedrakeproject.org/

Path: mosaic/tesseracts/thermal-mesh/firedrake-heat

Image: firedrake_heat_thermal_mesh

Benchmark checks: H Steady conduction in a slab: forward pass, gradient pass, optimization pass — see Results.

Effort: Steady conduction in a slab none — 1 intervention, 1 with unknown reason — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
nodal-misfit-seed: The identification-error gradient seeds the derivative of the temperature functional with a hand-built cofunction dI/dT = 2 (T - T_target) at the nodes that have a target value, because the error is a nodal sum rather than an integral.	L2	unknown	Steady conduction in a slab (incorrect-gradient): The earlier approximation, an L2 integral rescaled by the node count, gave gradients that did not match the reported value: gradient/source_fd_check went from 4.58e-02 to 3.52e-05 with the nodal seed, and the rho gradient was 9% off central differences, with a conductivity_recovery_bfgs final/initial ratio of 0.00313 against 0.00020. Whether a vertex quadrature rule would let the nodal sum be written in UFL, which Firedrake’s custom quadrature may allow, was not checked, so it is not known whether the hand-built seed is needed. Evidence: #162; #167, summary; #184; current form from #229.	—

Reason unknown. Nobody has measured whether these are needed; a measurement, or removing them, would settle it:

• nodal-misfit-seed (L2): The identification-error gradient seeds the derivative of the temperature functional with a hand-built cofunction dI/dT = 2 (T - T_target) at the nodes that have a target value, because the error is a nodal sum rather than an integral.

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. ρ=1 → fully conducting, ρ=0 → insulating.
source	Array[(None,), Float32]	✓	Per-element volumetric heat source (W/m³), shape (n_cells,). Active slice = source[:hex_mesh.n_faces]. Zero by default (topology-optimisation mode; no body heat source).
target_temperature	Array[(None,), Float32]		Per-node target temperature for the inverse (source-identification) objective, shape (n_nodes,). identification_error =
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe temperature (values shape (n_groups, 1), None → zero); neumann.mask/values prescribe heat flux (values shape (n_groups, 1)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
k_max	float		Maximum thermal conductivity (fully solid/conducting material).
p_exp	float		SIMP penalisation exponent p (k(ρ) = k_min + (k_max−k_min)·ρ^p).

Outputs

Field	Type	∂	Description
thermal_compliance	Array[(), Float32]	✓	Total thermal compliance C = ∮_Γ_N q_n · T dΓ (scalar). Equals the work done by the Neumann heat flux on the temperature field.
identification_error	Array[(), Float32]	✓	Source-identification objective:

Heat transfer

JAX-FEM

Numerics: Finite Element, BiCGSTAB+Jacobi

jax AD: hybrid GPU

JAX-FEM HEX8 steady-state heat-conduction solver with SIMP conductivity and optional volumetric source. Gradients ∂C/∂ρ and ∂C/∂source via JAX autodiff through the linear FE solve (hybrid implicit-function-theorem adjoint)

H slab · 3/3 pass Details

Upstream docs: https://github.com/deepmodeling/jax-fem

Path: mosaic/tesseracts/thermal-mesh/jax-fem

Image: jax_fem_thermal_mesh

Benchmark checks: H Steady conduction in a slab: forward pass, gradient pass, optimization pass — see Results.

Effort: Steady conduction in a slab none — 1 intervention — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
jitted-bicgstab-solver (tuned): Passes JAX-FEM’s own Jacobi-preconditioned BiCGSTAB (the default jax_solver’s algorithm and tolerances), compiled once per matrix size, to ad_wrapper through the custom_solver option for the forward and adjoint solves.	L2	choice	Steady conduction in a slab (too-slow): The default jax_solver runs bicgstab eagerly, so it compiles 2 XLA programs per forward call and 4 per VJP, a fixed overhead on every call. The umfpack_solver option used before gives identical outputs but is much slower on large meshes. At N = 1024, spatial_cost took 60.8 s with umfpack_solver and 5.17 s with this solver, vjp_cost/by_N 120 s and 10.8 s. Evidence: #286, before/after table; comment above _bicgstab.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. ρ=1 → fully conducting, ρ=0 → insulating.
source	Array[(None,), Float32]	✓	Per-element volumetric heat source (W/m³), shape (n_cells,). Active slice = source[:hex_mesh.n_faces]. Zero by default (topology-optimisation mode; no body heat source).
target_temperature	Array[(None,), Float32]		Per-node target temperature for the inverse (source-identification) objective, shape (n_nodes,). identification_error =
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe temperature (values shape (n_groups, 1), None → zero); neumann.mask/values prescribe heat flux (values shape (n_groups, 1)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.

Outputs

Field	Type	∂	Description
thermal_compliance	Array[(), Float32]	✓	Total thermal compliance C = ∮_Γ_N q_n · T dΓ (scalar). Equals the work done by the Neumann heat flux on the temperature field.
identification_error	Array[(), Float32]	✓	Source-identification objective:

Heat transfer

torch-fem

Numerics: Finite Element, CG+Jacobi

pytorch AD: hybrid GPU

torch-fem SolidHeat HEX8 heat-conduction solver with IsotropicConductivity3D SIMP material — the first PyTorch backend in thermal-mesh. Gradients via PyTorch autograd through the FE solve; GPU-native.

H slab · 3/3 pass Details

Upstream docs: https://github.com/meyer-nils/torch-fem

Path: mosaic/tesseracts/thermal-mesh/torch-fem-thermal

Image: torch_fem_thermal_mesh

Benchmark checks: H Steady conduction in a slab: forward pass, gradient pass, optimization pass — see Results.

Effort: Steady conduction in a slab L0 — 5 interventions — 0 setup deviations

Interventions and deviations

Abstraction level: framework. Effort counts only required interventions. Each item is an invitation to contribute: an upstream fix or a better configuration moves the solver down the ladder.

Intervention	Level	Necessity	Why, per setup	Upstream
gpu-device: With a GPU present, runs torch-fem on it instead of torch’s default CPU device.	L0	choice	Steady conduction in a slab: Not compared with the CPU: the tesseract has run on the GPU since its first version.	—
global-default-device: Sets torch’s default device (and dtype, float64) process-wide, so the tensors torch-fem allocates internally land on the GPU with the mesh, material and load tensors.	L0	required	Steady conduction in a slab (crash): torch-fem allocates auxiliary tensors without device=, so they land on the CPU and collide with the CUDA tensors. No error log is recorded. Evidence: Code comment; torch-fem’s benchmarks/utils.py at v0.12.1 and its current installation docs set the default device the same way.	—
gpu-cg-jacobi: On the GPU, passes method=“cg” and preconditioner=“jacobi” to model.solve for the forward and adjoint solves. On the CPU torch-fem chooses.	L0	choice	Steady conduction in a slab (too-slow): torch-fem 0.12.1 solves systems below 10,000 DOFs directly with SciPy on the CPU, even when the matrix is on the GPU; larger ones already get CG with Jacobi without AmgX. So the setting keeps the smaller meshes on the device. No isolated measurement. Evidence: #187; resolve_method and resolve_library in torch-fem v0.12.1 src/torchfem/sparse.py.	—
cg-tolerance (tuned): The GPU CG stops at a relative residual of 1e-6 (stol) instead of torch-fem’s 1e-10.	L0	choice	Steady conduction in a slab: No measurable gain here: cost moved by -16% to +11% against 1e-10, within the run-to-run spread on a shared GPU. Every status is unchanged from 1e-12 to 1e-6; forward metrics are bit-identical. Evidence: #301, notes.	—
deterministic-kernels: Enables torch.use_deterministic_algorithms(True, warn_only=True), so torch-fem’s CUDA scatter-adds in assembly sum in a fixed order.	L0	choice	Steady conduction in a slab (non-reproducible): Without it, identical GPU runs differed: the conductivity_recovery_bfgs final loss varied between 241 and 5.7e3, depending on whether the run left a plateau before the iteration cap. With it, reruns are identical (#200, and again at every CG tolerance in #301). It cost about 5 to 15% in #200, on an RTX 2000 Ada. Evidence: #200, notes; #301, notes.	—

Inputs

Field	Type	∂	Description
rho	Array[(None,), Float32]	✓	Per-cell density field, shape (n_max_cells,). Active slice = rho[:hex_mesh.n_faces]. ρ=1 → fully conducting, ρ=0 → insulating.
source	Array[(None,), Float32]	✓	Per-element volumetric heat source (W/m³), shape (n_cells,). Active slice = source[:hex_mesh.n_faces]. Zero by default (topology-optimisation mode; no body heat source).
target_temperature	Array[(None,), Float32]		Per-node target temperature for the inverse (source-identification) objective, shape (n_nodes,). identification_error =
boundary_conditions	MeshBC		Mesh boundary conditions. dirichlet.mask/values prescribe temperature (values shape (n_groups, 1), None → zero); neumann.mask/values prescribe heat flux (values shape (n_groups, 1)).
hex_mesh	HexMesh		Hexahedral mesh (8-node hex elements). n_points / n_faces give the active slice sizes.
k_max	float		Maximum thermal conductivity (fully solid material).
p_exp	float		SIMP penalisation exponent p (k(ρ) = k_min + (k_max−k_min)·ρ^p).

Outputs

Field	Type	∂	Description
thermal_compliance	Array[(), Float32]	✓	Total thermal compliance C = ∮_Γ_N q_n · T dΓ (scalar). Equals the work done by the Neumann heat flux on the temperature field.
identification_error	Array[(), Float32]	✓	Source-identification objective:

Capability matrix

These tables show which benchmark suites each solver runs, based on the exclusions in the problem configs. ✓ means the solver runs every experiment in the suite, — means the method cannot run any of them (not supported, as on the Results pages), and k/n means it runs k of the suite’s n experiments in full; the others are not supported or excluded for a reason that may be lifted. Notes lists the excluded experiment keys by category, and each solver name links to its card, which gives the reasons.

Navier–Stokes (2D)

Solver	Forward	Gradient	Cost	Optimization	Notes
Exponax	9/10	✓	✓	—	not supported: forward/cylinder, optimization/drag_opt
INS.jl	✓	✓	✓	✓
jax-cfd	9/10	✓	✓	0/1	not implemented: forward/cylinder, optimization/drag_opt
OpenFOAM	8/10	—	2/4	—	not implemented: forward/kolmogorov_laminar, forward/kolmogorov; not supported: cost/vjp_cost, gradient, optimization
PhiFlow	✓	✓	✓	✓
PICT	✓	✓	✓	0/1	infeasible: optimization/drag_opt
Warp-NS	✓	9/10	✓	✓	unstable: gradient/chain_fd_check
XLB	✓	✓	✓	✓

Navier–Stokes (3D)

Solver	Forward	Gradient	Cost	Optimization	Notes
Exponax	✓	✓	✓	✓
INS.jl	✓	✓	✓	✓
jax-cfd	✓	✓	✓	✓
OpenFOAM	✓	—	2/4	—	not supported: cost/vjp_cost, gradient, optimization
PhiFlow	✓	✓	✓	✓
PICT	✓	✓	✓	✓
Warp-NS	✓	✓	✓	✓
XLB	✓	✓	✓	✓

Structural mechanics

Solver	Forward	Gradient	Cost	Optimization	Notes
deal.II	✓	—	2/4	—	not supported: cost/vjp_cost, gradient, optimization
FEniCS	✓	✓	✓	✓
Firedrake	✓	✓	✓	✓
JAX-FEM	✓	✓	✓	✓
TopOpt.jl	✓	✓	✓	✓
torch-fem	✓	✓	✓	✓

Heat transfer

Solver	Forward	Gradient	Cost	Optimization	Notes
deal.II	✓	—	2/4	—	not supported: cost/vjp_cost, gradient, optimization
FEniCS	✓	✓	✓	✓
Firedrake	✓	✓	✓	✓
JAX-FEM	✓	✓	✓	✓
torch-fem	✓	✓	✓	✓