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.

翻译：暂无翻译