Pith. sign in

REVIEW 2 major objections 3 minor 56 references

A solver-reported 'converged' solve can be 1.35–49.5× over the true residual tolerance; the paper's guarded floor policy rechecks the equilibrium residual before accepting any solution.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

T0 review · deepseek-v4-flash

2026-08-01 17:00 UTC pith:7Y3AGKVB

load-bearing objection Guarded floor-selection is a useful robustness idea, but the FP64 recomputed-residual certificate is unvalidated in exactly the ill-conditioned regime where the paper needs it. the 2 major comments →

arxiv 2607.26382 v1 pith:7Y3AGKVB submitted 2026-07-29 cs.CE cs.NAmath.NA

When a positive SIMP density floor is not enough: solver admissibility and guarded floor selection in matrix-free 3D topology optimization

classification cs.CE cs.NAmath.NA MSC 65F1065N5574P15
keywords topology optimizationSIMPdensity floormatrix-free finite elementsgeometric multigridFGMRESresidual acceptancesolver admissibility
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper shows that in a matrix-free geometric-multigrid FGMRES solver for three-dimensional SIMP topology optimization, the projected residual used by the Krylov method to declare convergence can disagree with a freshly recomputed equilibrium residual. On four of 102 frozen test states the solver reported convergence while ||f−Ku||/||f|| exceeded the 1e-6 tolerance by factors of 1.35 to 49.5; inside an unguarded optimization loop, 22 of 40 state solves reached the iteration cap without any error being raised. The paper proposes a guard: nothing is accepted until the true residual is recomputed, and a 100-iteration probe plus a two-feature rule decides whether to keep the original density floor or escalate to a larger one. The guarded policy preserves the original floor where it is admissible, escalates only under residual evidence, and every selected solve then satisfies the tolerance. This matters because substituting a fixed raised floor changes compliance and sensitivity gradients substantially on hard random states, and a positive density floor alone guarantees neither solver convergence nor trustworthy stopping.

Core claim

The paper's central claim is that solver admissibility of a SIMP density floor is an operational property: a positive floor can yield a nonsingular matrix yet be inadmissible for a particular hierarchy at a requested tolerance. Empirically, the projected-residual stopping estimate diverges from the recomputed true residual on four of 102 held-out states, with the true residual up to 49.5 times tolerance even though the internal flag reads 1e-6. The remedy is a guarded residual-probe floor-selection algorithm: run a 100-iteration probe at the original floor, escalate to 10^-3 or 10^-2 if the residual features indicate difficulty, and never accept a solution until ||f−Ku||/||f|| ≤ 1e-6 is reco

What carries the argument

The load-bearing component is the acceptance guard: a fresh, unpreconditioned evaluation of the relative residual ||f−K(ρ_min)u||/||f|| in FP64 on free degrees of freedom, computed after the solver returns, as the sole certificate that a solution is accepted. Around it sits a decision rule that uses two probe features — the FGMRES residual history at iterations 50 and 100, plus the plateau ratio r100/r50 — to choose whether to try the original floor (10^-12) first or escalate along the ladder {10^-3, 10^-2}. The paper also derives a hierarchy-independent screening rule, ρ_min ≳ c(ρ) ε/τ, from a reduced-operator conditioning indicator, below which no backward-stable FP64 solve can guarantee t

Load-bearing premise

The acceptance guard recomputes the residual 'with the same matrix-free operator path used by the solve,' and the paper does not independently verify that this operator kernel matches an assembled reference matrix; if the operator itself were inconsistent with the intended discretization, the guard would certify equilibrium of the wrong operator.

What would settle it

Pick a held-out or new frozen SIMP state, run the guarded policy, and compare the matrix-free operator K applied to a few vectors against a sparse assembled reference K from the same element stiffness assembly; alternatively, find a state where the recomputed true residual at the selected floor exceeds 1e-6 while the guard reports acceptance — a single such case would refute the 'all selected solves satisfy tolerance' claim.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

If this is right

  • Any Krylov workflow that stops on a projected residual estimate should recompute a true residual before accepting a solution; the guard costs one operator application.
  • A positive SIMP density floor does not guarantee solver admissibility; escalation decisions should be made on residual evidence per state rather than globally.
  • Always using a fixed raised floor alters the operator: on 24 severe random states a fixed 10^-3 floor changes compliance by 31.0% on average and gradient by 0.340, while on optimized designs the changes are 0.48% and 0.008.
  • The guarded policy converts 22-of-40 silent cap-saturating solves in the unguarded cantilever trajectory into 40-of-40 verified solves.
  • At the conventional floor 10^-6, eight of twelve tested states still require escalation, so the phenomenon is not confined to the extreme 10^-12 setting.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same projected-versus-true residual gap may appear in other flexible Krylov solvers with iterate-dependent preconditioners; the paper's fixed-preconditioner control suggests the gap is tied to that variability, not to high contrast alone.
  • The screening rule ρ_min ≳ c(ρ) ε/τ could be tested as a cheap pre-solve check on other matrix-free multigrid stacks and other material interpolation exponents.
  • Reusing the previous floor decision between outer iterations instead of re-probing at every iteration is a plausible cost cut the paper left untested; likewise, a denser floor ladder would trade smaller interventions for more failed attempts.
  • For practitioners who only use the floor as numerical regularization, a fixed raised floor with the same residual guard may be far cheaper than the guarded ladder, at the cost of the fidelity losses the paper quantifies.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

2 major / 3 minor

Summary. The paper studies a matrix-free GMG-FGMRES solver for 3D SIMP topology optimization and documents a failure mode in which the internal projected-residual stopping test reports convergence while a freshly recomputed residual exceeds the target tolerance. On 4 of 102 held-out frozen states, the recomputed residual is 1.35–49.5 times the tolerance (Sec. 4.1.1, Table 4); inside an unguarded 40-iteration optimization trajectory, 22 of 40 state solves hit the 300-iteration cap without raising an error (Sec. 4.1.2, Table 5). The authors propose a guarded floor-selection policy: a 100-iteration probe at the original floor extracts two residual features, a heuristically chosen rule selects the first floor to attempt, and acceptance is always gated on a recomputed true residual ||f−Ku||/||f|| ≤ 10⁻⁶ (Algorithm 1, Sec. 3.5). On the held-out set the probe rule matches 98 of 102 reference classifications, the guard recovers all four missed escalations, and all 102 selected solves pass the recomputed-residual test. The policy preserves the original floor on 24 states and quantifies the fidelity/cost tradeoff: relative to a fixed 10⁻³ floor, it avoids mean compliance changes of 31.0% and gradient changes of 0.340 on severe random states, at 2.5× the mean wall time. A nine-state fixed-preconditioner control supports the attribution of the drift to the iterate-dependent preconditioner.

Significance. If the central claim holds, the paper is a valuable and well-scoped contribution to numerical topology optimization: it gives concrete evidence that a common Krylov stopping criterion can silently certify non-equilibrium states in SIMP loops, and it proposes a simple, parameter-free acceptance safeguard. The evaluation is unusually careful: thresholds were fixed on development states, the held-out set is disjoint, the four missed escalations are reported in detail, and the nine-state control with adaptive components disabled is honestly labeled. The paper also ships open-source code with pinned environments, which substantially strengthens reproducibility. The contribution is pragmatic rather than theoretical — the probe thresholds and ladder are implementation-specific heuristics — but the acceptance guard itself is not circular: it recomputes the residual independently of any fitted quantity, and the probe only controls which floor is tried first. The main unresolved risk is that the acceptance certificate is itself an FP64 residual evaluated in exactly the ill-conditioned regime where the paper shows FP64 conditioning limits; this needs to be addressed before the 'verified' lang

major comments (2)
  1. [Sec. 3.2 / Table 2 and Sec. 4.2] The acceptance certificate η = ||f−Ku||/||f|| is computed in FP64 with the same matrix-free operator path as the solve. The paper nowhere checks the accuracy of this FP64 residual against a higher-precision or independently assembled evaluation. This matters because Sec. 4.2 reports κ up to 1.6e16 at ρ0 = 10⁻¹² and derives the screen ρmin ≳ c(ρ)ε/τ (Eq. 10), below which "no backward-stable FP64 solve can be guaranteed to certify the tolerance," and Sec. 4.4.1 reports two admissible original-floor states with compliance of order 1e11. For such states, rounding in f−Ku can be of order ε·κ relative to ||f||, which can exceed τ = 10⁻⁶. The paper's own screen therefore suggests that the FP64 recomputed residual may not be a reliable certificate at ρ0 = 10⁻¹². Without an extended-precision or compensated-summation cross-check, or a comparison against an assembled sparse K for at least the redu
  2. [Sec. 3.2 and 'Relation to prior implementations'] The recomputed residual uses the same matrix-free operator kernel as the solve, and that kernel is inherited from prior preprints [25, 26] without any independent verification against an assembled reference matrix. If the matrix-free implementation of K in Eq. (5) has a systematic error, then both the projected residual and the acceptance residual are residuals of the same wrong operator, and the guard would certify equilibrium of that wrong operator. This is distinct from the paper's main point about projected-versus-recomputed disagreement, but it is load-bearing for the phrase 'verified solve' and for the claim that the guard is the correctness safeguard. I request a small, concrete validation: compare Kx from the matrix-free operator against a standard assembled sparse stiffness matrix for a representative set of random density fields and boundary conditions (or at least on the reduc
minor comments (3)
  1. [Throughout, esp. Table 2] The word 'true residual' is potentially misleading. The paper's acceptance quantity is a recomputed FP64 residual, not a reference-precision residual. Consider renaming it 'recomputed FP64 residual' in Table 2 and the abstract, reserving 'true' for the exact quantity in Eq. (6).
  2. [Fig. 4b] The caption reports 'fit 0.98, r = 0.99' without specifying which points are included. The text says 27 solves are excluded; please state explicitly that the slope 0.98 and correlation 0.99 are for the remaining 117 solves and clarify whether the regression is on the log-log data only.
  3. [Data Availability] The research-data package is described as 'will be deposited before journal submission,' but no persistent identifier is given. Please add the DOI or repository link, or state clearly that the code repository alone suffices for all claims except the transfer study.

Circularity Check

0 steps flagged

No significant circularity: the acceptance guard is parameter-free and the central claims are empirically evaluated against held-out data; self-citations describe infrastructure rather than the result.

full rationale

The paper's central claims are (1) that the projected FGMRES stopping estimate can fall below 1e-6 while a recomputed f-Ku residual is above tolerance, (2) that a two-feature probe rule matches 98/102 reference classifications on held-out states, and (3) that the guarded policy, which accepts only after recomputing the true residual, recovers the four missed escalations. The acceptance guard is parameter-free and uses no fitted quantity: it directly recomputes ||f-Ku||/||f|| and compares it to tau. The phrase 'all 102 selected solves satisfy the tolerance' is a consequence of Algorithm 1's acceptance gate, not a prediction; the paper transparently states that 'nothing is returned until ||f-Ku||/||f|| ≤ 1e-6 is recomputed' and identifies the guard, not the probe rule, as the load-bearing component. The probe thresholds were fixed on development states before the held-out audit and only influence which floor is tried first, never acceptance. The Sec. 4.2 screening rule fits a per-field constant c on the same reduced direct-solve data and validates conservativeness in-sample; the paper explicitly labels it an empirical diagnostic ('It is not used for acceptance') and a statement about the reduced operator under a backward-stable direct solve, so it does not support the central acceptance claims. The self-citations [25,26] describe the matrix-free operator and multigrid hierarchy used as infrastructure; the present paper claims no new kernel or hierarchy. Whether the inherited operator kernel is faithful to the intended discretization is a correctness/verification limitation, not a circularity, because the projected-vs-recomputed residual comparison is internal to the stated operator path. No derivation in the paper reduces by definition to its own inputs.

Axiom & Free-Parameter Ledger

4 free parameters · 5 axioms · 0 invented entities

No new physical or mathematical entities are introduced. 'Solver admissibility' is a definitional label; the probe features, ladder, and guard are computational procedures, not postulates. The four free parameters listed are all disclosed, and three of them (thresholds, ladder, probe length) affect only the cost/fidelity tradeoff, not the acceptance criterion. The fifth ledger-style quantity, c(rho), is measured rather than fitted but is validated in-sample. The one load-bearing unverified premise is the faithfulness of the matrix-free operator kernel itself.

free parameters (4)
  • Probe rule thresholds (r50 >= 1e-2; r100 >= 1e-4; r100/r50 >= 0.6) = 1e-2 / 1e-4 / 0.6
    Fixed on 16 development states (Sec. 3.4, 4.5.1). They govern which floor is tried first, not acceptance; the paper explicitly declined to retune on held-out outcomes.
  • Floor ladder values L = {1e-3, 1e-2} = 1e-3, 1e-2
    Chosen escalation levels (Sec. 3.3). The paper states the escalated floor is 'the smallest tested floor in L, not the minimal admissible floor.'
  • Probe length (100 iterations) and original floor rho0 = 1e-12 = 100; 1e-12
    Implementation choices setting the cost/fidelity operating point; sensitivity to rho0 is checked on a 12-state subset (Sec. 4.5.1).
  • Per-field conditioning constant c(rho) in Eq. (10) = implicit; 'nearly fixed by the field over six decades'
    Measured from reduced eigensolves of the same 18 fields used to validate the screen rho_min >= c*epsilon/tau; the screening rule is in-sample validated as conservative (Sec. 4.2).
axioms (5)
  • domain assumption The matrix-free operator path used in Sec. 3.2 faithfully evaluates f - K(rho_min)u for the intended discretized elasticity operator, with no independent verification against an assembled K in this paper.
    Both the solver's internal estimate and the acceptance guard call the same matrix-free kernel stack (Table 2, Sec. 3.2); a kernel-level error would invalidate the 'verified' framing.
  • standard math Element stiffness contributions are positive semidefinite and the fully solid constrained stiffness matrix is positive definite on the free degrees of freedom (Proposition 1).
    Needed for the monotone floor-regularization inequality of Eq. (7); standard for ersatz-SIMP with trilinear hexahedra.
  • standard math Classical finite-precision analysis of Krylov stopping estimates (Saad [44]; Carson and Higham [27, 29]) applies to this configuration.
    Used in Sec. 5.3 to interpret the divergence as projected-estimate drift; the FGMRES recurrence is not instrumented in the paper, so the specific path is not identified.
  • domain assumption Empirical convergence of the geometric-multigrid hierarchy is the arbiter of admissibility; no convergence theorem is claimed.
    The paper explicitly states Eq. (7) 'does not imply monotone conditioning improvement' and that hierarchy convergence 'is decided empirically and verified by the acceptance guard' (Sec. 3.3).
  • domain assumption Bernoulli random fields are a representative stress test for solver admissibility.
    Disclosed as 'deliberately harsher than the smooth, filtered iterates an optimizer produces' (Sec. 4); they carry the primary 102-state evaluation.

pith-pipeline@v1.3.0-daily-deepseek · 28213 in / 17882 out tokens · 188361 ms · 2026-08-01T17:00:28.258251+00:00 · methodology

0 comments
read the original abstract

In a matrix-free geometric-multigrid FGMRES solver for three-dimensional SIMP topology optimization, a reported converged solve is not always a converged solve. On four of 102 held-out states, the projected residual used for stopping falls below $10^{-6}$ while a recomputed true residual is 1.35 to 49.5 times the tolerance; in an unguarded optimization trajectory, 22 of 40 state solves reach the iteration cap without raising an error. We formulate floor selection as a verified control problem: probe the frozen state at the original floor, use two residual features to choose the first attempted floor, and accept no solution until $||f-Ku||/||f||\le10^{-6}$ is recomputed. The two-feature rule matches 98 of 102 reference classifications; the residual guard detects the four missed escalations, and all 102 selected solves satisfy the tolerance. Relative to always using a $10^{-3}$ floor, the policy preserves the original operator on 24 admissible states and avoids mean compliance and gradient changes of 31.0% and 0.340 on those severe random states, and 0.48% and 0.008 on seven optimized designs, at 2.5 times the mean wall time. In a 12-state subset of the held-out states, eight still require escalation at the conventional floor $10^{-6}$. In a nine-state control with the preconditioner's adaptive components disabled, every failure is visible and no false acceptance occurs, tying the stopping-estimate drift to the iterate-dependent preconditioner. The recomputed residual is the correctness safeguard; the probe and floor ladder govern an implementation-specific cost-fidelity tradeoff.

Figures

Figures reproduced from arXiv: 2607.26382 by Jun Wang, Shaoliang Yang, Yunsheng Wang.

Figure 1
Figure 1. Figure 1: The two problem families that carry the primary evidence, and what a frozen state [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: The guarded policy and the failure it exists to catch. (a) Control flow: a 100-iteration [PITH_FULL_IMAGE:figures/full_fig_p010_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: The acceptance failure inside a 40-iteration optimization loop. (a, b) Cantilever, 64k [PITH_FULL_IMAGE:figures/full_fig_p011_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Conditioning of the reduced direct problem against the floor for 18 fields at eight floors [PITH_FULL_IMAGE:figures/full_fig_p013_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Held-out evaluation on 102 unseen random states. (a) Decision map in the two probe [PITH_FULL_IMAGE:figures/full_fig_p014_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Cost of the guarded policy. (a) Per-case wall time on the 102 held-out states: bars [PITH_FULL_IMAGE:figures/full_fig_p015_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: What a fixed raised floor changes on states that did not need it. (a) Absolute relative [PITH_FULL_IMAGE:figures/full_fig_p017_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Where escalation is needed, and what the stack contributes. (a–c) Residual histories at [PITH_FULL_IMAGE:figures/full_fig_p018_8.png] view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

56 extracted references · 14 canonical work pages · 2 internal anchors

  1. [1]

    Generating optimal topologies in structural design using a homogenization method.Computer Methods in Applied Mechanics and Engineering, 71(2):197–224, 1988

    Martin Philip Bendsøe and Noboru Kikuchi. Generating optimal topologies in structural design using a homogenization method.Computer Methods in Applied Mechanics and Engineering, 71(2):197–224, 1988. doi: 10.1016/0045-7825(88)90086-2. 27

  2. [2]

    Bendsøe and Ole Sigmund

    Martin P. Bendsøe and Ole Sigmund. Material interpolation schemes in topology optimiza- tion.Archive of Applied Mechanics, 69(9–10):635–654, 1999. doi: 10.1007/s004190050248

  3. [3]

    Bendsøe and Ole Sigmund.Topology Optimization: Theory, Methods, and Applications

    Martin P. Bendsøe and Ole Sigmund.Topology Optimization: Theory, Methods, and Applications. Springer, Berlin, Heidelberg, 2004. doi: 10.1007/978-3-662-05086-6

  4. [4]

    Filters in topology optimization.International Journal for Numerical Methods in Engineering, 50(9):2143–2158, 2001

    Blaise Bourdin. Filters in topology optimization.International Journal for Numerical Methods in Engineering, 50(9):2143–2158, 2001. doi: 10.1002/nme.116

  5. [5]

    Lazarov and Ole Sigmund

    Boyan S. Lazarov and Ole Sigmund. Filters in topology optimization based on Helmholtz- type differential equations.International Journal for Numerical Methods in Engineering, 86 (6):765–781, 2011. doi: 10.1002/nme.3072

  6. [6]

    Guest, Jean H

    James K. Guest, Jean H. Prévost, and Ted Belytschko. Achieving minimum length scale in topology optimization using nodal design variables and projection functions.International Journal for Numerical Methods in Engineering, 61(2):238–254, 2004. doi: 10.1002/nme.1064

  7. [7]

    Lazarov, and Ole Sig- mund

    Erik Andreassen, Anders Clausen, Mattias Schevenels, Boyan S. Lazarov, and Ole Sig- mund. Efficient topology optimization in MATLAB using 88 lines of code.Structural and Multidisciplinary Optimization, 43(1):1–16, 2011. doi: 10.1007/s00158-010-0594-7

  8. [8]

    A new generation 99 line Matlab code for compliance topology optimization and its extension to 3d.Structural and Multidisciplinary Optimization, 62(4):2211–2228, 2020

    Federico Ferrari and Ole Sigmund. A new generation 99 line Matlab code for compliance topology optimization and its extension to 3d.Structural and Multidisciplinary Optimization, 62(4):2211–2228, 2020. doi: 10.1007/s00158-020-02629-w

  9. [9]

    Topology optimization approaches.Structural and Multidis- ciplinary Optimization, 48(6):1031–1055, 2013

    Ole Sigmund and Kurt Maute. Topology optimization approaches.Structural and Multidis- ciplinary Optimization, 48(6):1031–1055, 2013. doi: 10.1007/s00158-013-0978-6

  10. [10]

    Ioannis P. A. Papadopoulos. Numerical analysis of the SIMP model for the topology optimization problem of minimizing compliance in linear elasticity.Numerische Mathematik, 157(1):213–248, 2025. doi: 10.1007/s00211-024-01438-3

  11. [11]

    Shun Wang, Eric de Sturler, and Glaucio H. Paulino. Large-scale topology optimization using preconditioned Krylov subspace methods with recycling.International Journal for Numerical Methods in Engineering, 69(12):2441–2468, 2007. doi: 10.1002/nme.1798

  12. [12]

    Niels Aage and Boyan S. Lazarov. Parallel framework for topology optimization using the method of moving asymptotes.Structural and Multidisciplinary Optimization, 47(4): 493–505, 2013. doi: 10.1007/s00158-012-0869-2

  13. [13]

    Oded Amir, Niels Aage, and Boyan S. Lazarov. On multigrid-CG for efficient topology optimization.Structural and Multidisciplinary Optimization, 49(5):815–829, 2014. doi: 10.1007/s00158-013-1015-5

  14. [14]

    Niels Aage, Erik Andreassen, and Boyan S. Lazarov. Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework.Structural and Multidisciplinary Optimization, 51(3):565–572, 2015. doi: 10.1007/s00158-014-1157-0

  15. [15]

    Lazarov, and Ole Sigmund

    Niels Aage, Erik Andreassen, Boyan S. Lazarov, and Ole Sigmund. Giga-voxel com- putational morphogenesis for structural design.Nature, 550(7674):84–86, 2017. doi: 10.1038/nature23911

  16. [16]

    A system for high-resolution topology optimization.IEEE Transactions on Visualization and Computer Graphics, 22(3):1195–1208,

    Jun Wu, Christian Dick, and Rüdiger Westermann. A system for high-resolution topology optimization.IEEE Transactions on Visualization and Computer Graphics, 22(3):1195–1208,

  17. [17]

    Träff, Anton Rydahl, Sven Karlsson, Ole Sigmund, and Niels Aage

    Erik A. Träff, Anton Rydahl, Sven Karlsson, Ole Sigmund, and Niels Aage. Simple and efficient GPU accelerated topology optimisation: Codes and applications.Computer Methods in Applied Mechanics and Engineering, 410:116043, 2023. doi: 10.1016/j.cma.2023.116043

  18. [18]

    Junpeng Zhao, Tianyuan Qi, and Chunjie Wang. Efficient GPU accelerated topology optimization of composite structures with spatially varying fiber orientations.Computer Methods in Applied Mechanics and Engineering, 421:116809, 2024. doi: 10.1016/j.cma.2024. 116809

  19. [19]

    An efficient GPU solver for 3D topology optimization of continuous fiber-reinforced composite structures.Computer Methods in Applied Mechanics and Engineering, 435:117675, 2025

    Tianyuan Qi, Junpeng Zhao, and Chunjie Wang. An efficient GPU solver for 3D topology optimization of continuous fiber-reinforced composite structures.Computer Methods in Applied Mechanics and Engineering, 435:117675, 2025. doi: 10.1016/j.cma.2024.117675

  20. [20]

    Efficient large-scale 3D topology optimization with matrix-free MATLAB code.Structural and Multidisciplinary Optimization, 68(9):174, 2025

    Junpeng Wang, Niels Aage, Jun Wu, Ole Sigmund, and Rüdiger Westermann. Efficient large-scale 3D topology optimization with matrix-free MATLAB code.Structural and Multidisciplinary Optimization, 68(9):174, 2025. doi: 10.1007/s00158-025-04127-3

  21. [21]

    Topology optimization methods for 3d structural problems: A comparative study.Archives of Computational Methods in Engineering, 29(3):1525–1567, 2022

    Daniel Yago, Juan Cante, Oriol Lloberas-Valls, and Javier Oliver. Topology optimization methods for 3d structural problems: A comparative study.Archives of Computational Methods in Engineering, 29(3):1525–1567, 2022. doi: 10.1007/s11831-021-09626-2

  22. [22]

    A parallel geometric multigrid method for adaptive topology optimization.Structural and Multidisciplinary Optimization, 66(10):225, 2023

    David Herrero-Pérez and Sebastián Ginés Picó-Vicente. A parallel geometric multigrid method for adaptive topology optimization.Structural and Multidisciplinary Optimization, 66(10):225, 2023. doi: 10.1007/s00158-023-03675-w

  23. [23]

    Jan Christoph Krüger and Benedikt Kriegesmann. On the application of the multigrid method to topology optimization with orthotropic material with varying orientation.Struc- tural and Multidisciplinary Optimization, 68(8):158, 2025. doi: 10.1007/s00158-025-04102-y

  24. [24]

    Effi- cienthybridtopologyoptimizationusingGPUandhomogenization-basedmultigridapproach

    Arya Prakash Padhi, Souvik Chakraborty, Anupam Chakrabarti, and Rajib Chowdhury. Effi- cienthybridtopologyoptimizationusingGPUandhomogenization-basedmultigridapproach. Engineering with Computers, 39(5):3593–3615, 2023. doi: 10.1007/s00366-022-01771-x

  25. [25]

    Matrix-Free 3D SIMP Topology Optimization with Fused Gather-GEMM-Scatter Kernels

    Shaoliang Yang, Jun Wang, and Yunsheng Wang. Matrix-free 3D SIMP topology optimiza- tion with fused gather-GEMM-scatter kernels.arXiv preprint arXiv:2604.18020, 2026. doi: 10.48550/arXiv.2604.18020. URLhttps://arxiv.org/abs/2604.18020

  26. [26]

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

    Shaoliang Yang, Jun Wang, and Yunsheng Wang. A matrix-free Galerkin multigrid solver and failure-mode screen for single-GPU 3D SIMP linear systems.arXiv preprint arXiv:2604.26441, 2026. doi: 10.48550/arXiv.2604.26441. URLhttps://arxiv.org/abs/ 2604.26441

  27. [27]

    Erin Carson and Nicholas J. Higham. A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems.SIAM Journal on Scientific Computing, 39(6):A2834–A2856, 2017. doi: 10.1137/17M1122918

  28. [28]

    Azzam Haidar, Harun Bayraktar, Stanimire Tomov, Jack Dongarra, and Nicholas J. Higham. Mixed-precision iterative refinement using tensor cores on GPUs to accelerate solution of linear systems.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476(2243):20200110, 2020. doi: 10.1098/rspa.2020.0110

  29. [29]

    Higham and Theo Mary

    Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. Acta Numerica, 31:347–414, 2022. doi: 10.1017/S0962492922000022

  30. [30]

    McCormick, Joseph Benzaken, and Rasmus Tamstorf

    Stephen F. McCormick, Joseph Benzaken, and Rasmus Tamstorf. Algebraic error analysis for mixed-precision multigrid solvers.SIAM Journal on Scientific Computing, 43(5):S392–S419,

  31. [31]

    Three-precision algebraic multigrid on GPUs.Future Generation Computer Systems, 149:280–293, 2023

    Yu-Hsiang Mike Tsai, Natalie Beams, and Hartwig Anzt. Three-precision algebraic multigrid on GPUs.Future Generation Computer Systems, 149:280–293, 2023. doi: 10.1016/j.future. 2023.07.024

  32. [32]

    Element connectivity parameterization for topology optimization of geometrically nonlinear structures.International Journal of Solids and Structures, 42(7):1983–2009, 2005

    Gil Ho Yoon and Yoon Young Kim. Element connectivity parameterization for topology optimization of geometrically nonlinear structures.International Journal of Solids and Structures, 42(7):1983–2009, 2005. doi: 10.1016/j.ijsolstr.2004.09.005

  33. [33]

    Tortorelli

    Anna Dalklint, Mathias Wallin, and Daniel A. Tortorelli. Structural stability and artificial buckling modes in topology optimization.Structural and Multidisciplinary Optimization, 64 (4):1751–1763, 2021. doi: 10.1007/s00158-021-03012-z

  34. [34]

    Zhengfang Zhang, Yanqiang Dong, and Weifeng Chen. A non-ersatz material approach for the topology optimization of elastic structures based on piecewise constant level set method.Communications in Computational Physics, 30(5):1370–1389, 2021. doi: 10.4208/ cicp.OA-2020-0231

  35. [35]

    A level set method for structural topology optimization.Computer Methods in Applied Mechanics and Engineering, 192(1–2): 227–246, 2003

    Michael Yu Wang, Xiaoming Wang, and Dongming Guo. A level set method for structural topology optimization.Computer Methods in Applied Mechanics and Engineering, 192(1–2): 227–246, 2003. doi: 10.1016/S0045-7825(02)00559-5

  36. [36]

    Structural optimization using sensitivity analysis and a level-set method.Journal of Computational Physics, 194(1): 363–393, 2004

    Grégoire Allaire, François Jouve, and Anca-Maria Toader. Structural optimization using sensitivity analysis and a level-set method.Journal of Computational Physics, 194(1): 363–393, 2004. doi: 10.1016/j.jcp.2003.09.032

  37. [37]

    MFEM: A modular finite element methods library.Computers and Mathematics with Applications, 81:42–74, 2021

    Robert Anderson, Julian Andrej, Andrew Barker, Jamie Bramwell, Jean-Sylvain Camier, Jakub Cerveny, Veselin Dobrev, Yohann Dudouit, Aaron Fisher, Tzanio Kolev, Will Pazner, Mark Stowell, Vladimir Tomov, Ido Akkerman, Johann Dahm, David Medina, and Stefano Zampini. MFEM: A modular finite element methods library.Computers and Mathematics with Applications, 8...

  38. [38]

    A matrix-free approach for finite-strain hyperelastic problems using geometric multigrid

    Denis Davydov, Jean-Paul Pelteret, Daniel Arndt, Martin Kronbichler, and Paul Steinmann. A matrix-free approach for finite-strain hyperelastic problems using geometric multigrid. International Journal for Numerical Methods in Engineering, 121(13):2874–2895, 2020. doi: 10.1002/nme.6336

  39. [39]

    High-order matrix- free incompressible flow solvers with GPU acceleration and low-order refined preconditioners

    Michael Franco, Jean-Sylvain Camier, Julian Andrej, and Will Pazner. High-order matrix- free incompressible flow solvers with GPU acceleration and low-order refined preconditioners. Computers and Fluids, 203:104541, 2020. doi: 10.1016/j.compfluid.2020.104541

  40. [40]

    Will Pazner, Tzanio Kolev, and Jean-Sylvain Camier. End-to-end GPU acceleration of low-order-refined preconditioning for high-order finite element discretizations.The International Journal of High Performance Computing Applications, 37(5):578–599, 2023. doi: 10.1177/10943420231175462

  41. [41]

    Stitt, Kenneth Weiss, Vladimir Z

    Arturo Vargas, Thomas M. Stitt, Kenneth Weiss, Vladimir Z. Tomov, Jean-Sylvain Camier, Tzanio Kolev, and Robert N. Rieben. Matrix-free approaches for GPU acceleration of a high-order finite element hydrodynamics application using MFEM, Umpire, and RAJA. The International Journal of High Performance Computing Applications, 36(4):492–509,

  42. [42]

    Ham, and Paul H

    Tianjiao Sun, Lawrence Mitchell, Kaushik Kulkarni, Andreas Klöckner, David A. Ham, and Paul H. J. Kelly. A study of vectorization for matrix-free finite element methods.The International Journal of High Performance Computing Applications, 34(6):629–644, 2020. doi: 10.1177/1094342020945005. 30

  43. [43]

    Large-scale multigrid with adaptive Galerkin coarsening, 2025

    Fabian Böhm, Nils Kohl, Harald Köstler, and Ulrich Rüde. Large-scale multigrid with adaptive Galerkin coarsening, 2025. arXiv preprint arXiv:2511.13109

  44. [44]

    A flexible inner-outer preconditioned GMRES algorithm.SIAM Journal on Scientific Computing, 14(2):461–469, 1993

    Yousef Saad. A flexible inner-outer preconditioned GMRES algorithm.SIAM Journal on Scientific Computing, 14(2):461–469, 1993. doi: 10.1137/0914028

  45. [45]

    Society for Industrial and Applied Mathematics, Philadelphia, PA, 2 edition, 2003

    Yousef Saad.Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2 edition, 2003. doi: 10.1137/1.9780898718003

  46. [46]

    Briggs, Van Emden Henson, and Steve F

    William L. Briggs, Van Emden Henson, and Steve F. McCormick.A Multigrid Tutorial. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2 edition, 2000. doi: 10.1137/1.9780898719505

  47. [47]

    Oosterlee, and Anton Schüller.Multigrid

    Ulrich Trottenberg, Cornelius W. Oosterlee, and Anton Schüller.Multigrid. Academic Press, San Diego, CA, 2000

  48. [48]

    AmgX: A library for GPU accelerated algebraic multigrid and preconditioned iterative methods.SIAM Journal on Scientific Computing, 37(5): S602–S626, 2015

    Maxim Naumov, Marat Arsaev, Patrice Castonguay, Jonathan Cohen, Julien Demouth, Joe Eaton, Simon Layton, Nikolay Markovskiy, Istvan Reguly, Nikolai Sakharnykh, Vijay Sellappan, and Robert Strzodka. AmgX: A library for GPU accelerated algebraic multigrid and preconditioned iterative methods.SIAM Journal on Scientific Computing, 37(5): S602–S626, 2015. doi:...

  49. [49]

    Olson, Jacob Schroder, and Ben Southworth

    Nathan Bell, Luke N. Olson, Jacob Schroder, and Ben Southworth. PyAMG: Algebraic multigrid solvers in Python.Journal of Open Source Software, 8(87):5495, 2023. doi: 10.21105/joss.05495

  50. [50]

    Erin Carson and Nicholas J. Higham. Accelerating the solution of linear systems by iterative refinement in three precisions.SIAM Journal on Scientific Computing, 40(2):A817–A847,

  51. [51]

    Boman, Erin Carson, Terry Cojean, Jack Dongarra, Alyson Fox, Mark Gates, Nicholas J

    Ahmad Abdelfattah, Hartwig Anzt, Erik G. Boman, Erin Carson, Terry Cojean, Jack Dongarra, Alyson Fox, Mark Gates, Nicholas J. Higham, Xiaoye S. Li, Jennifer Loe, Piotr Luszczek, Srikara Pranesh, Siva Rajamanickam, Tobias Ribizel, Barry F. Smith, Kasia Swirydowicz, Stephen Thomas, Stanimire Tomov, Yaohung M. Tsai, and Ulrike Meier Yang. A survey of numeric...

  52. [52]

    Mixed-precision numerics in scientific applications: sur- vey and perspectives.The Journal of Supercomputing, 82(5):287, 2026

    Aditya Kashi, Hao Lu, Wesley Brewer, David Rogers, Michael Matheson, Mallikarjun Shankar, and Feiyi Wang. Mixed-precision numerics in scientific applications: sur- vey and perspectives.The Journal of Supercomputing, 82(5):287, 2026. doi: 10.1007/ s11227-026-08264-4. 31

  53. [2016]

    doi: 10.1109/TVCG.2015.2502588. 28

  54. [2018]

    doi: 10.1137/17M1140819

  55. [2021]

    doi: 10.1137/20M1348571. 29

  56. [2022]

    doi: 10.1177/10943420221100262