A Matrix-Free Galerkin Multigrid Solver and Failure-Mode Screen for Single-GPU 3D SIMP Linear Systems

TL;DR

This paper presents a matrix-free geometric multigrid solver for large 3D SIMP linear systems on a single GPU, including a diagnostic failure-mode screening method and mixed-precision evaluation.

Contribution

It introduces a novel matrix-free Galerkin multigrid hierarchy, a failure-mode screening tool, and assesses mixed-precision variants for efficient GPU-based 3D SIMP solves.

Findings

Pass rates vary with problem size and iteration cap, with higher success at smaller sizes.

FP32-GMG outperforms BF16-GMG in wall-time ratios and convergence.

Largest reported solve handles 1 million elements in under 2 seconds.

Abstract

Large 3D SIMP studies require repeated elasticity solves for density-dependent operators whose finest matrices are expensive to assemble and whose conditioning degrades under high contrast. We study this linear-solver layer rather than claiming end-to-end optimization acceleration. The solver builds a matrix-free Galerkin geometric multigrid (GMG) hierarchy around a fused fine operator: the finest level remains matrix-free, the first coarse level is assembled by local Galerkin aggregation, and deeper levels use sparse Galerkin products. The practical default is FP32-GMG; BF16 is evaluated as a guarded mixed-precision variant and diagnostic stress test, not as the main speed mechanism. In a 27-case heterogeneous cantilever sweep, pass rates under a 200-iteration budget are 7/9, 4/9, and 1/9 at 64k, 216k, and 512k elements; converged-only mean iteration counts are about 112, 134, and 146. On uniform rho=0.5, p=3 solves, FP32-GMG gives 1.62x, 1.75x, and 3.12x wall-time ratios relative to the capped flat Jacobi-PCG baseline at the same sizes; that non-converged baseline reaches the 200-iteration cap in all timed trials. BF16-GMG is not faster than FP32-GMG. In 18 fixed-seed heterogeneous BF16 validation cases, 7/18 converge, matching the FP64 count, and 11 cases that pass the spectral screen still fail the 500-iteration cap; the screen is therefore diagnostic rather than a convergence certificate. The largest reported solve is a 1M-element uniform-modulus system solved in 1.50+/-0.58 s with an 8.66 GiB hierarchy-allocation delta during setup, not a peak-memory trace; this point is reported as uniform scaling, not heterogeneous robustness evidence. The contribution is therefore a bounded single-GPU solver result built on an inherited Level 0 matrix-free operator: a Galerkin GMG hierarchy, direct BF16 guard evidence, and an explicit failure-mode screen for structured 3D SIMP linear systems.