DoPhys Applied Mathematics & Kernel Group

Technical Reports & Mathematical Monographs

Complete variational formulations, spectral stability proofs, error bounds, and GPU PTX micro-architectural analyses underlying the DoPhys Core engine.

TR-2025-04 Published: October 18, 2025 • 19 Pages
DOI: 10.5281/dophys.2025.04

Eliminating Non-Associative Floating-Point Drift in Multi-GPU Krylov Solvers

H. Chen (DoPhys Zurich), M. Lindqvist (ETH Zurich / DoPhys), S. K. Venkatesh (DoPhys Boston)

Abstract: In distributed-memory Krylov subspace solvers ($\text{FGMRES}$, $\text{BiCGSTAB}$), global inner products $\langle \mathbf{u}, \mathbf{v} \rangle = \sum_{i=1}^N u_i v_i$ are evaluated via warp-level shuffle reductions followed by MPI $\text{AllReduce}$ trees. Because IEEE-754 floating-point addition $\mathbf{fl}(a \oplus b)$ is non-associative, thread scheduling jitter and changes in domain decomposition $P$ introduce $\mathcal{O}(N \epsilon_{\text{mach}})$ rounding discrepancies. In chaotic or stiffly coupled multiphysics systems, these perturbations grow as $\exp(\lambda_{\max} t)$, violating strict bit-reproducibility requirements in aerospace hardware-in-the-loop ($\text{HIL}$) certification.

Mathematical Formulation & Super-Accumulator Architecture: Every IEEE-754 binary64 product $p_i = u_i \cdot v_i$ (computed exactly via fused multiply-add splitting $p_i = h_i + l_i$ with $h_i = \mathbf{fl}(u_i v_i)$ and $l_i = \mathbf{fma}(u_i, v_i, -h_i)$) has an exact binary representation within the exponent range $[2^{-1074}, 2^{1024})$. We map each exact product into a thread-block-shared Kulisch fixed-point accumulator $\mathcal{A} \in \mathbb{Z}$ comprising $W = 34$ unsigned 64-bit words ($2176\text{ bits}$ total, including $80\text{ carry-guard bits}$):

Theorem 2.1 (Partition-Invariant Exact Inner Product & CGS2 Stability)

Let $\Omega = \bigsqcup_{k=1}^{P} \Omega_k$ be an arbitrary disjoint partitioning of degrees of freedom across $P$ MPI ranks. Each product $u_i v_i = (-1)^{s_i} m_i 2^{e_i - 52}$ is bit-shifted into word index $w_i = \lfloor (e_i + 1074)/64 \rfloor$ with bit offset $b_i = (e_i + 1074) \bmod 64$:

$$\mathcal{A}(\mathbf{u}, \mathbf{v}) = \sum_{k=1}^{P} \sum_{i \in \Omega_k} (-1)^{s_i} m_i \, 2^{e_i - 52} = \sum_{w=0}^{33} A_w \, 2^{64w - 1074}, \qquad \mathbf{RN}\big(\mathcal{A}(\mathbf{u}, \mathbf{v})\big) = \mathbf{fl}_{\text{exact}}\!\left(\sum_{i=1}^{N} u_i v_i\right)$$

When embedded inside Classical Gram-Schmidt with re-orthogonalization ($\text{CGS2}$), the Arnoldi basis $\mathbf{V}_m = [\mathbf{v}_1, \dots, \mathbf{v}_m]$ satisfies the partition-invariant loss-of-orthogonality bound:

$$\left\|\mathbf{I}_m - \mathbf{V}_m^T \mathbf{V}_m\right\|_2 \le \frac{c_1 \, m^{3/2} \, \epsilon_{\text{mach}}}{1 - c_2 \, \kappa_2(\mathbf{K}_m) \, \epsilon_{\text{mach}}} \approx 1.11 \times 10^{-16}, \qquad \left\|\mathbf{r}_m^{(P_1)} - \mathbf{r}_m^{(P_2)}\right\|_\infty \equiv 0$$
Artifact Evaluated: 0-ULP Bit-Identical (1 to 64 GPUs) Overhead vs CUB: +3.1%
TR-2025-03 Published: August 29, 2025 • 27 Pages
DOI: 10.5281/dophys.2025.03

Preconditioned Implicit Adjoints for Stiff Thermo-Mechanical Contact

A. Subramanian (DoPhys Boston), K. Vance (DoPhys Zurich)

Abstract: Gradient-based topology and shape optimization of high-heat-flux cold plates and micro-actuators requires differentiating through coupled nonlinear hyperelasticity, Fourier heat conduction, and unilateral frictional contact. Standard non-smooth penalty or active-set contact formulations exhibit discontinuous Jacobians across stick-slip and gap-closure transitions, causing adjoint Krylov solvers to stagnate or diverge.

Variational Formulation & Transposed Schur Preconditioning: We replace non-smooth Signorini-Coulomb conditions on the mortar contact boundary $\Gamma_c$ with a $C^2$-continuous log-barrier normal potential $b(d_n; \hat{d})$ and a regularized tangential friction dissipation potential $D_\epsilon(\mathbf{v}_\tau)$:

Coupled Thermo-Contact Adjoint System & Schur Complement Factorization
$$f_\epsilon(\|\mathbf{v}_\tau\|) = \mu_f \, p_n(d_n) \left\{ \begin{array}{ll} \frac{\|\mathbf{v}_\tau\|}{\epsilon_v} \left(2 - \frac{\|\mathbf{v}_\tau\|}{\epsilon_v}\right), & \|\mathbf{v}_\tau\| < \epsilon_v \\[6pt] 1, & \|\mathbf{v}_\tau\| \ge \epsilon_v \end{array} \right., \qquad p_n(d_n) = -\frac{\partial b(d_n; \hat{d})}{\partial d_n}$$

The discrete adjoint vector $\boldsymbol{\lambda} = [\boldsymbol{\lambda}_u^T, \boldsymbol{\lambda}_T^T]^T$ satisfies the transposed block-coupled system $\mathcal{J}^T \boldsymbol{\lambda} = -\nabla_{\mathbf{u}, T} \mathcal{L}^T$, which we precondition via an upper block-triangular approximate Schur complement $\mathcal{P}_{\text{adj}}^{-T}$:

$$\begin{bmatrix} \mathbf{K}_{uu} + \mathbf{H}_{\text{contact}}^{(n,\tau)} & \mathbf{K}_{Tu}^T \\[4pt] \mathbf{K}_{uT}^T & \mathbf{M}_T / \Delta t + \mathbf{K}_{TT} \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda}_u \\[4pt] \boldsymbol{\lambda}_T \end{bmatrix} = -\begin{bmatrix} \partial \mathcal{J} / \partial \mathbf{u} \\[4pt] \partial \mathcal{J} / \partial \mathbf{T} \end{bmatrix}, \qquad \widetilde{\mathbf{S}}_T^T = \left(\frac{\mathbf{M}_T}{\Delta t} + \mathbf{K}_{TT}\right) - \mathbf{K}_{uT}^T \, \text{diag}(\mathbf{K}_{uu})^{-1} \mathbf{K}_{Tu}^T$$
Complex-Step FD Verification: $|\nabla_{\text{adj}} - \nabla_{\text{CSFD}}| / |\nabla_{\text{CSFD}}| < 2.4 \times 10^{-11}$
TR-2025-01 Published: May 12, 2025 • 15 Pages
DOI: 10.5281/dophys.2025.01

Zero-Copy Surrogate Training: Bridging C++20 Field Tensors and JAX XLA

H. Chen (DoPhys Zurich), R. Takahashi (DoPhys Singapore)

Abstract: Purely data-driven neural PDE surrogates (FNO, GraphCast, MeshGraphNets) frequently violate conservation of mass $\nabla \cdot \mathbf{u} = 0$ and accumulate unphysical energy at high wavenumbers during autoregressive rollouts. Embedding a classical Poisson projection step inside the training loop stabilizes rollouts, but copying unstructured mesh tensors between Python ML frameworks and C++ solvers over PCIe degrades throughput by an order of magnitude.

Self-Adjoint Discrete Leray-Helmholtz Projection via XLA TypedFFI: We expose DoPhys's matrix-free AMG-preconditioned Poisson solver directly inside XLA HLO fusion graphs using the XLA CustomCall v4 (TypedFFI) interface and zero-copy DLPack device pointers. Given an intermediate neural velocity prediction $\mathbf{u}^* \in \mathbb{R}^{3N_c}$, the discrete Leray projection operator $\mathbb{P}_h$ enforces exact discrete solenoidality:

Self-Adjoint Discrete Hodge Decomposition & Custom VJP Rule
$$\mathbf{u}_{\text{div-free}} = \mathbb{P}_h(\mathbf{u}^*) = \mathbf{u}^* - \mathbf{M}_u^{-1} \mathbf{G}_h \left(\mathbf{D}_h \mathbf{M}_u^{-1} \mathbf{G}_h\right)^{\!\dagger} \mathbf{D}_h \mathbf{u}^*, \qquad \mathbf{D}_h = -\mathbf{G}_h^T$$

Because the discrete divergence $\mathbf{D}_h$ and gradient $\mathbf{G}_h$ are negative adjoints with respect to the inner product weighted by cell volumes $\mathbf{M}_u$, the projection operator $\mathbb{P}_h$ is an orthogonal projector ($\mathbb{P}_h^2 = \mathbb{P}_h, \ \mathbf{M}_u \mathbb{P}_h = \mathbb{P}_h^T \mathbf{M}_u$). Consequently, the backward vector-Jacobian product ($\text{VJP}$) for cotangent vector $\bar{\mathbf{u}}$ is simply:

$$\text{vjp}\big(\mathbb{P}_h, \mathbf{u}^*\big)(\bar{\mathbf{u}}) = \mathbb{P}_h^*\big(\bar{\mathbf{u}}\big) = \mathbf{M}_u^{-1} \mathbb{P}_h\big(\mathbf{M}_u \bar{\mathbf{u}}\big), \qquad \|\mathbf{D}_h \mathbf{u}_{\text{div-free}}\|_\infty \le 10^{-14}$$
JAX 0.4.30+ / XLA Custom Call API v4 • $4.8\times$ End-to-End Speedup on $8\times\text{H100}$