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

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.

• The paper demonstrates a matrix-free Galerkin multigrid solver that reduces iteration counts and wall-time for 3D SIMP linear systems on a single GPU.
• It employs a mixed-precision schedule using BF16, FP32, and FP64 to balance computational speed with numerical stability.
• The study highlights a failure-mode diagnostic that identifies robustness limits in moderate-contrast topology optimization scenarios.

Matrix-Free Galerkin Multigrid Hierarchies for 3D SIMP Linear Systems

Introduction

This work addresses the computational bottleneck in large-scale 3D SIMP-based topology optimization, specifically the repeated solution of large, density-dependent linear elasticity systems on a single GPU. The core contribution is a matrix-free Galerkin geometric multigrid (GMG) solver stack that builds a preconditioning hierarchy around a matrix-free operator in a manner compatible with current GPU architectures, systematically evaluates mixed-precision (BF16/FP32/FP64) schedules where feasible, and provides an empirical diagnostic for failure modes. This GMG strategy is positioned as an architectural and algorithmic advance compared to flat Jacobi-preconditioned CG solvers, particularly in regimes with moderate density contrast.

Methodology

The solver integrates a fully matrix-free fine-level operator directly leveraging fused gather--GEMM--scatter kernels. The Galerkin multigrid hierarchy is constructed such that:

• The fine-level operator is never assembled, avoiding sparse storage bottlenecks;
• Level-1 operator is assembled by efficient, local matrix-free Galerkin aggregation, requiring only elementwise modulus data;
• Deeper levels utilize sparse-matrix triple products, enabling consistency with the filtered fine-level operator.

The precise placement of this GMG stack within the SIMP optimization loop is depicted in (Figure 1).

Figure 1: Solver placement inside one SIMP iteration. The density update loop remains unchanged; the main change is the introduction of the GMG preconditioning stack between the current density and sensitivity update.

Mixed-precision is handled through a precision-descent schedule: the expensive fine-level smoother utilizes BF16 tensor cores, Level-1 employs FP32, and deeper levels revert to FP64 for numerical stability, while the outer Krylov solver and residual computations operate in FP64 throughout. This asymmetric allocation is justified by the cost profile and error-propagation properties of the hierarchy (Figure 2).

Figure 2: Precision descent across the GMG hierarchy, with BF16 at the fine-level, FP32 at Level-1, and FP64 on deeper levels, while residuals and the Krylov basis are held in FP64.

A degree-ν Chebyshev-Jacobi smoother is integrated by default, with ablations provided against Jacobi smoothers. Transfer between levels is performed by structured $2{:}1$ coarsening and prolongation with careful Dirichlet boundary masking.

An empirical spectral-proxy ($$\varepsilon_{\mathrm{BF16}}\cdot\kappa_{\mathrm{eff}}$$) is adopted as a diagnostic screen for expected BF16 stability, motivated by standard mixed-precision numerical linear algebra analysis: reliable low-precision operation is only expected when the effective condition number is sufficiently small after multigrid preconditioning.

Numerical Results

Several key experimental findings are reported for 3D elasticity problems with up to one million elements and 3.09M DOFs on an NVIDIA RTX 4090.

Outer Iteration Counts and Convergence

On heterogeneous 3D cantilever benchmarks, the FP64-GMG hierarchy substantially reduces iteration counts versus Jacobi-preconditioned CG, but robustness diminishes with increasing mesh size and density contrast. Pass rates under a 200-iteration cap are 7/9 at 64k, 4/9 at 216k, and 1/9 at 512k across diverse binary-contrast configurations (Figure 3).

Figure 3: GMG outer PCG iteration count versus mesh size over 27 heterogeneous test cases.

Wall-Time Acceleration

On uniform-density problems ($\rho=0.5$, $p=3$), FP32-GMG yields mean per-linear-solve wall-time ratios of $1.62\times$, $1.75\times$, and $3.12\times$ versus Jacobi-PCG at 64k, 216k, and 512k, respectively. Notably, these speedups are achieved versus a capped, non-converged reference—Jacobi-PCG does not reach the desired residual tolerance in any trial at these scales. BF16-GMG is not faster than FP32-GMG, and its wall-time ratio versus Jacobi-PCG only exceeds unity at larger problems (Figures 4, 5).

Figure 4: Per-solve wall time scaling across mesh sizes for Jacobi-PCG, FP32-GMG, and BF16-GMG.

Figure 5: Wall time ratios of FP32-GMG and BF16-GMG versus capped Jacobi-PCG.

GMG preconditioning changes the outer-solver regime, producing rapid residual reductions where Jacobi-PCG stagnates (Figure 6).

Figure 6: Residual histories at 64k, 216k, 512k. GMG preconditioning delivers rapid convergence while Jacobi-PCG stagnates.

Kernel Throughput and Roofline Analysis

Fused BF16 WMMA and FP32 kernels achieve comparable throughput on the fine-level matrix-free matvecs (3046 vs. 3242 GFLOP/s). Operator throughput is memory-bandwidth bound, and further lowering arithmetic precision to BF16 does not yield significant solver-level speedups with current Q1-hex kernels (Figure 7, Figure 8).

Figure 7: BF16 vs. FP32 fine-level operator throughput shows limited difference due to bandwidth-bound regime.

Figure 8: Roofline analysis—matrix-free smoother is bandwidth-limited, not compute-limited.

Mixed-Precision Failure Screening

Evaluation of the empirical BF16-proxy ($$\varepsilon_{\mathrm{BF16}}\cdot\kappa_{\mathrm{eff}}$$) across 18 heterogeneous cases shows that it is not a reliable convergence classifier: 7/18 cases converge in BF16, matching FP64, but high values or proximity to threshold do not guarantee failure, and low values do not guarantee convergence (Figure 9).

Figure 9: $$\varepsilon_{\mathrm{BF16}}\cdot\kappa_{\mathrm{eff}}$$ across test suite. Most cases fall below the classical threshold, but the screen is not a full BF16 convergence classifier.

Ablations and Smoother Sensitivity

Precision, depth, and smoother policy ablations reveal:

• FP32 smoothing matches FP64 in iteration counts and is modestly faster;
• Additional FP32 coarse levels have modest influence on solve time;
• V-cycles are generally preferable to W-cycles in wall time, even if the latter achieves fewer Krylov steps (Figure 10);
• Degree-2 Jacobi smoothing is optimal in wall time on moderate-condition-number cases, but higher-degree Chebyshev smoothers are beneficial for robustness in high-contrast regimes.

  Figure 10: Ablation study—effect of fine-level precision, GMG depth, cycle type, and smoother degree on wall time and iteration counts.

Large-Scale and External Baseline Scaling

Uniform-modulus 1M-element problems are solved in $1.50\pm0.58$ s and 8.66 GiB VRAM allocation delta, demonstrating practical feasibility on commodity hardware for moderate-contrast fields (Figure 11). CPU PyAMG (assembled-matrix reference) is less efficient than the matrix-free GPU solver in the post-assembly build+solve phase, with a ratio exceeding $$\varepsilon_{\mathrm{BF16}}\cdot\kappa_{\mathrm{eff}}$$0 at 64k, but this gap is much smaller in the hierarchy construction phase (Figure 12).

Figure 11: Setup and solve scaling for increasing problem sizes up to 1M elements.

Figure 12: Post-assembly CPU PyAMG vs. GPU FP32-GMG timing at 64k elements—GPU path is substantially faster overall in solve phase.

Robustness and Failure Modes

Extensive screening across pathological and binarized density configurations exposes the hierarchy’s current limits. The GMG approach remains robust for moderate-contrast or centrally banded states, but fails or stagnates in the presence of checkerboard designs, severe near-void fields, and extreme binary mixtures (Figures 14, 15).

Figure 13: Robustness screen—GMG success (green) and failure (red) for ten representative fixed-seed density configurations.

Figure 14: Robustness basin map across volume fraction and penalization; screening metric correlates weakly with actual convergence.

Qualitative Solutions

A gallery of computed structures across various benchmarks and mesh sizes confirms the solver’s ability to process realistic 3D designs produced by contemporary SIMP-based topology optimization (Figure 15).

Figure 15: Isosurface renders of auxiliary density fields computed with the hierarchy across diverse benchmarks up to 1M elements.

Implications and Future Directions

The results demonstrate that a matrix-free Galerkin GMG hierarchy enables practical, high-throughput linear elasticity solves required in large-scale 3D topology optimization on a single consumer GPU for moderate-contrast problems. BM16 precision can be guarded by monitoring effective spectrum, but is not a robust path to acceleration beyond FP32 for current dense-matrix sizes and bandwidth-limited kernels.

Failure modes are definitively associated with coarse-space inadequacy and near-singular perturbations under extreme material heterogeneity. The introduced spectral-proxy is diagnostic but insufficient—coarse-grid pathology detection, adaptive coarsening, and improved smoothing will be necessary for increased robustness in high-contrast scenarios.

The matrix-free GPU approach is not universally superior; for irregular meshes or multi-GPU scalability, AMG remains competitive. Nonetheless, for structured, single-GPU, matrix-free workflows, the presented GMG approach sets a new engineering baseline and provides a platform for future advances in robust, high-contrast, large-scale topology optimization.

Conclusion

A matrix-free Galerkin multigrid hierarchy, architected for modern GPU hardware and integrated with spectral-proxy diagnostics, enables efficient, scalable solution of large 3D structured SIMP linear systems in moderate-contrast regimes. The approach achieves substantial wall-time reductions and practical solution of million-element problems on consumer hardware but is currently limited by robustness under extreme material heterogeneity and near-singular design fields. The evidence points to the critical interplay between multigrid hierarchy quality, precision scheduling, and density-field characteristics, setting a precise agenda for deeper solver-level advances and broader applicability in future structural optimization frameworks.

---

Reference:

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