Eliminating Non-Associative Floating-Point Drift in Multi-GPU Krylov Solvers
How we implemented 2048-bit Kulisch fixed-point super-accumulators in shared memory to guarantee $\text{ULP}=0$ bit-identical convergence across MPI ranks.
DoPhys Labs develops high-order finite-volume and symplectic integration engines compiled directly to multi-GPU and bare-metal targets. Designed for robotics state estimation, thermal-fluid digital twins, and inverse material design.
import dophys as dp
from dophys.solvers import ImplicitFVM
# 1. Load unstructured hexahedral mesh (2.4M cells)
mesh = dp.Mesh.from_cgns("stator_cooling_v4.cgns", precision="fp64")
# 2. Configure coupled Navier-Stokes + Fourier heat conduction
solver = ImplicitFVM(
mesh=mesh,
preconditioner="amg_smoothed_aggregation",
max_krylov_subspace=48,
rtol=1e-9
)
# 3. Differentiable adjoint step for boundary optimization
state, grad = solver.step_with_adjoint(
dt=2.5e-4,
inlet_velocity=(14.2, 0.0, 0.0),
heat_flux_w_m2=8.5e4
)
#include <dophys/hil/solver_context.hpp>
int main() {
// Lock 256 MiB static arena; zero heap allocs in loop
dophys::hil::ArenaConfig cfg{
.static_pool_bytes = 256ULL << 20,
.lock_pages = true
};
dophys::hil::SolverContext ctx(cfg, "actuator_rom_v2.dph");
while (ctx.bus_healthy()) {
const auto frame = ctx.read_ethercat_frame();
const auto est = ctx.step_symplectic_ekf(frame, 1e-3);
ctx.publish_torque_command(est.optimal_u);
}
}
import jax
import jax.numpy as jnp
from dophys.xla_ffi import helmholtz_projection_op
@jax.jit
def physics_informed_rollout(u_field, params):
# Neural operator predictor + exact divergence-free FVM projection
u_pred = neural_flux_update(u_field, params)
u_div_free = helmholtz_projection_op(u_pred, rtol=1e-10)
return jnp.mean(jnp.square(u_div_free - u_target))
grads = jax.grad(physics_informed_rollout)(u0, model_params)
Most fast simulators trade physical invariants for frame rate. DoPhys preserves discrete mass, momentum, and symplectic structure down to machine epsilon ($\epsilon_{\text{mach}} \approx 2.22 \times 10^{-16}$) while saturating modern memory bandwidth.
Sum-factorized spectral element and finite-volume kernels that evaluate operators in $\mathcal{O}(p^{d+1})$ complexity without global sparse matrix assembly—achieving 78% of peak roofline memory bandwidth.
Instead of unrolling Krylov iterations on an AD tape, DoPhys solves the discrete adjoint system $[\partial\mathcal{R}/\partial u]^T \lambda = -(\partial\mathcal{J}/\partial u)^T$ via the implicit function theorem—cutting VRAM usage by $40\times$.
custom_vjp bindings for JAX & PyTorchZero dynamic heap allocation after initialization. Designed for deployment inside $1\text{ kHz}$ control loops, aerospace hardware-in-the-loop benches, and $\mathrm{SE}(3)$ state-estimation filters.
Inspect real-world convergence trajectories $\|r_k\|_2 / \|r_0\|_2$ and memory roofline telemetry captured from DoPhys v2.4.2 across three representative multiphysics workloads.
Measured on single-socket AMD EPYC 9354 + NVIDIA RTX 6000 Ada (48GB ECC). Full input decks available in the documentation.
| Problem / Mesh | Degrees of Freedom ($N_{\text{DoF}}$) | Legacy CPU (32-Core MPI) | DoPhys v2.4 (Single GPU) | Adjoint Backward Pass | Speedup |
|---|---|---|---|---|---|
|
3D Lid-Driven Cavity ($Re = 10,000$)
Hexahedral, $p=4$ sum-factorized spectral
|
$4.19 \times 10^6$ | 612.0 ms / step | 11.4 ms / step | 6.8 ms | $53.7\times$ |
|
Nonlinear Hyperelasticity (Neo-Hookean)
Quadratic tetrahedral arm, Newton-Krylov
|
$1.85 \times 10^6$ | 284.5 ms / step | 6.2 ms / step | 3.1 ms | $45.9\times$ |
|
Conjugate Heat Transfer (Cold Plate)
Coupled fluid-solid interface, AMG-PCG
|
$8.92 \times 10^6$ | 1,490.0 ms / step | 29.8 ms / step | 18.5 ms | $50.0\times$ |
|
Articulated $\mathrm{SE}(3)$ Multibody + Frictional Contact
Batched 64-environment RL rollout
|
$64 \times 148\text{ DoF}$ | 19.4 ms / batch | 0.31 ms / batch | 0.22 ms | $62.5\times$ |
Deep dives into numerical stability, GPU kernel fusion, and mathematical proofs in our research archive.
How we implemented 2048-bit Kulisch fixed-point super-accumulators in shared memory to guarantee $\text{ULP}=0$ bit-identical convergence across MPI ranks.
$C^2$-regularized contact barriers and block-triangular Schur complement preconditioners for inverse topology optimization.
Designing custom XLA FFI calls that expose self-adjoint Leray-Helmholtz projection directly inside compiled neural operator pipelines.
The core Python bindings and single-GPU academic license are free for non-commercial research. For multi-node enterprise deployments, source-code escrow, or custom boundary-condition engineering, contact our systems team.