{"id":"ac66d732-8aa2-45a2-8ef8-c4f19f24ae6c","arxiv_id":"2506.04157","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":8,"one_line_summary":"A matrix-free, multigrid-preconditioned solver with a V-cycle BFBT Schur complement handles compressible, high-contrast-viscosity mantle convection Stokes problems at billion-unknown scale.","lead":"Geoscientists and applied mathematicians present a matrix-free solver for the compressible viscous-flow equations that drive mantle convection, replacing direct matrix assembly and solution by a combination of multigrid, Uzawa preconditioning, and a new Schur-complement approximation. They demonstrate the approach on a 2.7-billion-unknown spherical-shell simulation of Earth's mantle over 400 million years, showing roughly 69 percent of compute time goes to the Stokes solve.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (21) approximates B A^{-1} B^T, not the actual Schur complement (B+C) A^{-1} B^T; no bound on the C mismatch is given, so the claimed compressible robustness is untested outside the single density profile.","rationale":"The reader's weakest-assumption analysis points to the unproved V-cycle BFBT approximation for the non-symmetric compressible saddle point system. My stress test sharpens this into a specific algebraic defect: with C != 0, the exact Schur complement is \\bar B A^{-1} B^T, while Eq. (21) is built from B only, so the approximation is not even algebraically consistent at C = 0 in the limiting sense unless C is small. This is more concrete than a generic concern about missing spectral analysis, and it identifies a specific term that the numerical study never isolates. The proposed test directly measures sensitivity to the C term by varying \\lambda \\nabla \\ln \\rho; if the method is robust, the preconditioner's neglect of C is benign for this problem class, and the claim stands. If not, the conclusion in Sec. 8 overstates robustness. The paper has genuine independent support: reproducible source code is provided, the solver structure is clearly documented, and the strong/weak scaling data are internally consistent. Those strengths keep the result plausible, so I do not escalate beyond CONDITIONAL. My specific concern is a refinement, not a rejection, of the reader's verdict, so the verdict remains UNCHANGED.","tokens_in":28529,"tokens_out":3879,"duration_ms":36510,"concrete_test":"Using the open-source HyTeG code linked in Appendix A, run the Sec. 7.4 2D Stokes test with the V-cycle BFBT parameters mV=1, degV=1, tolV-BFBT=0.1 fixed, and scale the compressibility term by replacing \\nabla \\ln \\rho with \\lambda \\nabla \\ln \\rho for \\lambda = 0, 1, and 2, keeping the viscosity profile, mesh, and tolerances unchanged. If FGMRES iterations and wall time to reach tol(u,p)=1e-5 remain within about 20% across \\lambda, the C-mismatch concern does not land. If they degrade substantially, Eq. (21) should either incorporate \\bar B = B + C in the appropriate factor or the robustness claim must be qualified.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim that the solver is robust and scalable for the compressible TALA system with C != 0 rests on the V-cycle BFBT approximation in Eq. (21), which defines \\hat S_V^{-1} = (B \\hat A_C^{-1} B^T)^{-1} (B \\hat A_C^{-1} A \\hat A_C^{-1} B^T) (B \\hat A_C^{-1} B^T)^{-1}. For the generalised saddle point system (15) with off-diagonal block \\bar B = B + C, the exact Schur complement is S = \\bar B A^{-1} B^T, not B A^{-1} B^T. Thus Eq. (21) targets the symmetric incompressible Schur complement, and correctness for C != 0 requires the extra term C A^{-1} B^T to be negligible in a norm compatible with the Krylov iteration. The paper explicitly notes in Sec. 5.2.1 that B can be replaced by B + C in the Uzawa updates, but Eq. (21) never does so. No spectral estimate or numerical ablation isolating the compressibility term is provided. The numerical evidence uses one fixed density profile with Di/Gamma_0 about 0.378 and one viscosity model (24), so the demonstrated robustness to C is in-sample and unsupported outside this parameter point. This is a concrete algorithmic gap rather than a disagreement with consensus: the BFBT identity is exact only when C = 0, and the manuscript supplies neither analysis nor an experiment showing that the C term is harmless at the tested, let alone more extreme, compressibility.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a matrix-free iterative solver for a compressible (TALA) Stokes system coupled to a temperature advection-diffusion equation. Spatial discretisation uses P2-P1 Taylor-Hood elements on blended annulus and spherical-shell meshes; temporal discretisation combines BDF2 for diffusion with a particle-based modified method of characteristics for advection. The Stokes subproblem is solved by FGMRES with Uzawa-type block preconditioners, Chebyshev-smoothed geometric multigrid for the velocity block, and a proposed 'V-cycle BFBT' Schur approximation in which the diagonal of A is replaced by a single GMG V-cycle. The paper reports second-order temporal convergence, a comparison of Schur approximations and Uzawa variants, strong and weak scaling studies, and a 400 Myr mantle simulation with 2.11e9 velocity-pressure unknowns on 15,360 cores.","tokens_in":29044,"tokens_out":8340,"duration_ms":83314,"significance":"If the solver behaviour is as robust as claimed, the paper is a useful contribution to geodynamic modelling: matrix-free, scalable Stokes solvers that tolerate high viscosity contrast and a compressible mass balance are of genuine practical importance. The manuscript is careful in its internal solver comparisons, and the source code is made available with a DOI, which is a concrete reproducibility strength. The central robustness claim, however, is supported only by a single viscosity and density configuration, and the proposed V-cycle BFBT approximation for the non-symmetric compressible system is not backed by spectral analysis or an ablation isolating the compressibility term. The scalability evidence is reported per FGMRES iteration rather than as total time to solution, with no iteration counts for the weak scaling runs. These gaps limit confidence in the general robustness and scalability claims as currently stated.","major_comments":[{"comment":"Eq. (21) defines the inverse Schur operator using B \\hat A_C^{-1} B^T and B \\hat A_C^{-1} A \\hat A_C^{-1} B^T, i.e., it targets the incompressible Schur complement B A^{-1} B^T. For the generalised saddle point system (15) with off-diagonal block \\bar B = B + C, the exact Schur complement is \\bar B A^{-1} B^T. No bound or numerical experiment is given showing that the C A^{-1} B^T correction is negligible, and the note in Sec. 5.2.1 that B can be replaced by B + C in the Uzawa updates does not repair Eq. (21), because the inner products in the BFBT formula still use B only. The Sec. 8 claim that the solver is 'robust with respect to this term' is therefore not established outside the tested parameter point.","section":"5.2.3, Eq. (21)"},{"comment":"The V-cycle BFBT approximation rests on two unproved algorithmic assumptions: that a single Chebyshev-smoothed GMG V-cycle with mV = degV = 1 is an accurate enough surrogate for A^{-1} inside the BFBT identity, and that the resulting \\hat S_V is a sufficiently accurate Schur complement for the non-symmetric compressible system. No spectral analysis or convergence theorem is given, and the numerical support is limited to the Frank-Kamenetskii viscosity (24) with the radial base profile of Fig. 6 and one density profile with Di/Gamma_0 approximately 0.378. If either assumption degrades under other high-contrast viscosity profiles or stronger compressibility, the central robustness claim fails; a parameter study varying Di/Gamma_0 and viscosity contrast, or an ablation isolating the C-block contribution, is needed.","section":"5.2.3 and Sec. 7.2"},{"comment":"The weak scaling efficiency in Fig. 15a is computed from the average FGMRES time per iteration, not from the total saddle point solve time, and the paper does not report iteration counts for the weak scaling runs. If the number of FGMRES iterations grows with the number of unknowns, the per-iteration metric overstates the true time-to-solution scalability. The strong scaling plot (Fig. 13) likewise reports time per FGMRES iteration for the saddle point part. Total solve times and iteration counts should be reported for both scaling studies to support the scalability claim.","section":"7.5"}],"minor_comments":[{"comment":"The text contains the typo 'FMGRES' where the method is FGMRES; please correct both occurrences.","section":"7.2 and 7.5"},{"comment":"The left subcaption 'This is to make the scaling of the text similar.' appears to be a leftover sentence and should be removed or rewritten.","section":"Fig. 13"},{"comment":"The unknown count is given as '2 .106· 10^9' in Sec. 7.2, while the abstract and Sec. 8 state 2.11e9; please harmonise the notation.","section":"7.2"},{"comment":"The legend abbreviations SM, SV and S_w are not defined in the caption; please add a sentence identifying them with Eqs. (19), (21) and (20).","section":"Fig. 12"},{"comment":"The statement that omega = 0.0125 'worked best' for the weighted BFBT suggests parameter tuning for that variant; please state explicitly which parameters were used for each variant in Fig. 12 and whether the tuned value was used only for weighted BFBT.","section":"7.4"},{"comment":"The manufactured-solution test sets rho = 1 and uses a prescribed velocity, so it validates the temporal splitting of the temperature equation but not the coupled compressible Stokes solver; this limitation should be stated explicitly.","section":"7.1"},{"comment":"The claim that the degree-6 polynomial approximation of the exponential 'does not negatively impact the convergence rate' is not backed by a comparison; please add an iteration-count or residual-history comparison.","section":"5.2.2"}],"recommendation":"major_revision","confidential_remarks":"The paper fits the scope of the journal and the code availability is a clear strength. The main risk is that the central robustness claim for the V-cycle BFBT preconditioner in the compressible case rests on a single configuration, with no analysis of the C-block mismatch and no ablation. If the revision adds such an ablation and reports iteration counts and total solve times for the scaling experiments, I would be comfortable moving toward acceptance. I do not see citation or novelty issues requiring editorial action."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Worth reading if you work on Stokes preconditioners or large-scale mantle convection. The genuinely new piece is Eq. (21), the V-cycle BFBT approximation that replaces diag(A) with a cheap GMG V-cycle; combined with symmetric Uzawa and FGMRES it delivers a 2.11e9-unknown compressible TALA Stokes solve in about 10.5 days on 15,360 cores. That is a real capability data point.\n\nThe paper is unusually transparent about solver parameters (Table B1), gives a full source-code DOI, reports strong and weak scaling, and includes an internal comparison against mass-matrix and weighted-BFBT Schur approximations. The temporal convergence test is clean second order. The authors also do not oversell: they explicitly note that the compressible off-diagonal block makes the saddle-point system non-symmetric and that B can be replaced by B+C in the Uzawa updates.\n\nThe stress-test concern is legitimate. Eq. (21) is written for B A^{-1} B^T, while the actual generalized saddle-point Schur complement is (B+C) A^{-1} B^T. The paper supplies no bound, no spectral estimate, and no experiment that isolates the effect of C by varying compressibility. The numerical evidence is one density profile with Di/Gamma_0 about 0.378 and one viscosity model. That does not make the claim false, but it makes “robust for compressible TALA” an in-sample statement rather than a demonstrated property. A second, smaller issue: solver parameters such as omega, tolerances, l_eta, and l_min are tuned on the test problems, and the performance comparisons appear to be single-run. That tempers the strong-scaling numbers, though it is not a disqualifier. There is also no external benchmark against ASPECT or a community compressible-convection test; again a limitation, not a fatal one.\n\nWho should read this: numerical linear algebra people looking for practical preconditioners for high-contrast Stokes problems, and geodynamics code developers. The paper deserves serious peer review: the contribution is concrete, reproducible, and the central algorithmic claim is worth testing. A referee should push for a compressibility ablation and a clear statement about why omitting C in Eq. (21) is harmless, or a corrected formula that includes it.","headline":"Solid engineering paper with a genuinely new Schur complement variant and an honest large-scale demonstration; the main caveat is that the V-cycle BFBT approximation ignores the compressibility term C, and robustness beyond the single tested density profile is not established.","tokens_in":29509,"tokens_out":2725,"would_cite":true,"duration_ms":26697,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N30","65N55","65Y05","65F10","65Z05"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper proposes a matrix-free iterative solver for compressible Stokes flow with extreme viscosity contrasts and demonstrates its robustness on a 2.11-billion-unknown mantle convection simulation.","keywords":["Stokes problem","multigrid methods","high contrast","matrix-free methods","Schur complement","block preconditioner","geophysics","mantle convection"],"falsifier":"Compute the eigenvalues of $\\hat S_V^{-1} S$ or, equivalently, run the FGMRES method on a spherical-shell Stokes problem with a thin, stiff slab of viscosity contrast $10^6$ and with the compressibility term strengthened by increasing the dissipation number or density gradient. If the iteration count per time step grows linearly with the contrast or with mesh refinement, or if the V-cycle BFBT approximation degrades relative to the diag($A$) version, the central robustness claim is refuted.","tokens_in":28314,"feed_emoji":"🌋","tokens_out":6409,"duration_ms":57040,"temperature":0.7,"pith_summary":"The paper aims to show that the compressible Stokes system at the heart of mantle convection simulations, with viscosity contrasts spanning many orders of magnitude and a non-symmetric compressibility term, can be solved robustly and scalably by a fully matrix-free iterative solver. The key step is a new Schur complement approximation, the V-cycle BFBT, which replaces the diagonal of the velocity block in the standard BFBT formula with a single cheap geometric multigrid V-cycle. The authors support the claim with a 2.11e9-unknown Stokes solve inside a 400-million-year mantle simulation, completed in 10.52 days on 15,360 cores, together with strong and weak scaling tests. A sympathetic reader would care because faster forward solvers make higher-resolution mantle convection runs and iterative inverse problems feasible.","feed_headline":"Billion-scale mantle flow solver handles extreme viscosity","feed_subtitle":"A matrix-free Stokes solve manages 2.11 billion unknowns in 10.52 days, enabling 400-Myr mantle convection runs.","key_machinery":"The load-bearing object is the V-cycle BFBT Schur complement approximation of Eq. (21): $\\hat S_V^{-1} = (B \\hat A_C^{-1} B^T)^{-1} (B \\hat A_C^{-1} A \\hat A_C^{-1} B^T) (B \\hat A_C^{-1} B^T)^{-1}$, where $\\hat A_C$ is one geometric multigrid V-cycle with Chebyshev smoothing, using one pre- and post-smoothing step, in place of the diagonal diag($A$). It approximates the inverse Schur complement inside an FGMRES-accelerated symmetric Uzawa block preconditioner; the outer application also solves $B \\hat A_C^{-1} B^T$ by a CG solver preconditioned by the inverse-viscosity mass matrix to a relative tolerance of 0.1. The idea is that the V-cycle captures the action of $A^{-1}$ more faithfully than the diagonal, at a cost low enough for a matrix-free, massively parallel setting.","core_discovery":"The central claim is that Eq. (21), called the V-cycle BFBT approximation, provides a sufficiently accurate inverse Schur complement for the generalized saddle point system $(A \\quad B^T; \\, B+C \\quad 0)$ of compressible TALA mantle flow, even when viscosity varies over several orders of magnitude. In this approximation, the diagonal of $A$ used in the classical diag($A$)-BFBT formula is replaced by $\\hat A_C$, a single Chebyshev-smoothed geometric multigrid V-cycle with one pre- and post-smoothing step. Wrapped in an FGMRES outer iteration preconditioned by a symmetric Uzawa block preconditioner, the resulting solver is matrix-free and is demonstrated to be robust for the temperature- and space-dependent Frank-Kamenetskii viscosity model considered, reaching a Stokes solve with 2.11e9 unknowns in a 400-myr mantle convection simulation.","pith_inferences":["The robustness claim is algorithmic rather than proven: Eq. (21) is not backed by a spectral analysis of $\\hat S_V^{-1} S$, so a natural next test is to compute the eigenvalue distribution of $\\hat S_V^{-1} S$ for the given viscosity model and for stiffer ones.","The same 'replace diag(A) by one V-cycle' trick could plausibly transfer to other saddle point problems with heterogeneous coefficients, such as ice-sheet flow or magma dynamics, where mass-based Schur complements are known to degrade.","The comparison data suggest a practical plateau: for loose tolerances of $10^{-3}$ to $10^{-1}$, cheaper Schur approximations are competitive, so production codes that only need modest accuracy per time step may not need the extra machinery.","A concrete stress test would be to push the compressibility parameter, for example by increasing the dissipation number or the density gradient, and check whether the iteration count per FGMRES step stays bounded; the paper only demonstrates one parameter set."],"forward_implications":["The Stokes subproblem with 2.11e9 unknowns can be solved inside a 400-myr mantle simulation in 10.52 days on 15,360 cores, with 7.27 days spent on the saddle point FGMRES solve.","When only a small residual reduction is required per time step, all considered Schur complement approximations perform comparably; the V-cycle BFBT is clearly fastest when a larger reduction is needed.","Strong and weak scaling tests show that both the saddle point solver and the advection-diffusion solver scale nearly ideally to 15,360 cores, with some super-ideal behaviour attributed to smoother viscosity on finer grids.","The BDF2 plus particle-advection operator splitting retains second-order temporal convergence with a bounded error constant even for small diffusion coefficients, so advection-dominated regimes do not spoil the time accuracy.","Because the whole method is matrix-free, the limiting resource becomes memory rather than operator assembly, which is the property that makes higher-resolution simulations with 2.7e9 spatial degrees of freedom feasible."],"supporting_citations":[{"why":"Supplies the flexible GMRES outer iteration that permits variable preconditioning inside the saddle point solver.","marker":"[20]"},{"why":"Provides the Uzawa type block preconditioners and the convergence conditions used to build the solver.","marker":"[17]"},{"why":"Shows that BFBT-type Schur complement approximations are suitable for high-viscosity-contrast Stokes flow, motivating the new variant.","marker":"[62]"},{"why":"Introduces the weighted BFBT approximation and its asymmetric viscous rescaling, the main algorithmic baseline compared against the V-cycle BFBT.","marker":"[63]"},{"why":"Uses a similar saddle point solver structure for mantle convection and serves as a comparison for solver performance.","marker":"[6]"},{"why":"Documents Stokes solvers at extreme scale on hybrid grids, providing context for the scalability claims.","marker":"[42]"},{"why":"Supports matrix-free monolithic multigrid methods for Stokes and generalized Stokes problems, which the present solver extends to the compressible setting.","marker":"[18]"},{"why":"Shows a matrix-free multilevel preconditioner for the generalized Stokes problem with discontinuous viscosity, a direct antecedent for high-contrast robustness.","marker":"[19]"},{"why":"Develops a Schur complement approach robust to large viscosity jumps, another precursor of the BFBT strategy.","marker":"[56]"}],"fun_headline_variants":["Matrix-free Stokes solver scales to 2.11B unknowns","Billion-scale mantle flow with a matrix-free Uzawa solver","Solver for extreme viscosity contrasts reaches 2B+ unknowns","No matrix assembly: mantle convection at billion-scale"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that replacing diag($A$) in the BFBT formula with one cheap Chebyshev-smoothed multigrid V-cycle still gives a Schur complement approximation accurate enough for the non-symmetric compressible saddle point system; no spectral or convergence proof is provided, and the numerical evidence covers a single viscosity model and parameter set.","fun_headline_variants_meta":{"raw":{"variants":["Matrix-free Stokes solver scales to 2.11B unknowns","Billion-scale mantle flow with a matrix-free Uzawa solver","Solver for extreme viscosity contrasts reaches 2B+ unknowns","No matrix assembly: mantle convection at billion-scale"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000793,"raw_usage":{"total_tokens":3478,"prompt_tokens":916,"completion_tokens":2562,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":532,"completion_tokens_details":{"reasoning_tokens":2493}},"tokens_in":532,"tokens_out":2562,"duration_ms":19210,"temperature":1.0,"reasoning_tokens":2493,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T10:47:24.552954+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute the eigenvalues of $\\hat S_V^{-1} S$ or, equivalently, run the FGMRES method on a spherical-shell Stokes problem with a thin, stiff slab of viscosity contrast $10^6$ and with the compressibility term strengthened by increasing the dissipation number or density gradient. If the iteration count per time step grows linearly with the contrast or with mesh refinement, or if the V-cycle BFBT approximation degrades relative to the diag($A$) version, the central robustness claim is refuted.","supporting_citations":[{"cited_title":"SIAM Journal on Matrix Analysis and Applications 39(2), 932–960 (2018) https://doi.org/10.1137/16M1106304","cited_arxiv_id":null,"evidence_quote":"Provides the Uzawa type block preconditioners and the convergence conditions used to build the solver."},{"cited_title":"SIAM Journal on Scientific Computing 39(5), 272–297 (2017) https://doi.org/10.1137/16M108450X","cited_arxiv_id":null,"evidence_quote":"Introduces the weighted BFBT approximation and its asymmetric viscous rescaling, the main algorithmic baseline compared against the V-cycle BFBT."},{"cited_title":"Journal of Computational Science 17(3), 509–521 (2016) https://doi.org/10.1016/j.jocs.2016.06.006","cited_arxiv_id":null,"evidence_quote":"Documents Stokes solvers at extreme scale on hybrid grids, providing context for the scalability claims."},{"cited_title":"Journal of Computational Physics 230(24), 8835–8851 (2011) https: //doi.org/10.1016/j.jcp.2011.09.007","cited_arxiv_id":null,"evidence_quote":"Develops a Schur complement approach robust to large viscosity jumps, another precursor of the BFBT strategy."}],"review_version":1}