{"id":"147d5cd9-0945-4683-8104-9f3a4df0583e","arxiv_id":"1908.07071","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A low-order refined finite element mesh with ordered ILU-smoothed geometric multigrid and overlapping additive Schwarz yields matrix-free preconditioners whose iteration counts stay bounded in h, polynomial degree p, and DG penalty parameter eta.","lead":"This paper designs preconditioners that let high-order finite element and discontinuous Galerkin solvers run without ever forming the dense system matrix, using a coarser low-order mesh plus tailored smoothers. A smart generalist should read it because it removes a major memory bottleneck for high-order methods, which promise accuracy but are often too costly to solve in practice.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Corollary 1's p-robust bound depends on assumption (17) that the ILU(0)-smoothed local multigrid V-cycle is a uniform preconditioner in h and p, which is asserted rather than proved; a single p or aspect-ratio regime where the ordered-ILU smoother degrades would invalidate the headline claim.","rationale":"The reader's verdict is CONDITIONAL, with the weakest assumption identified as the unproved uniform local-solver bound (17). I agree that this is the single most load-bearing concern. The central claim, as stated in Corollary 1, is an unconditional bound kappa(B A_p) <= C independent of h and p. Its proof via Theorem 1 explicitly assumes (17), and the paper's only support for (17) is numerical evidence on a limited parameter range. The paper itself flags the risk: the refined mesh is anisotropic with aspect ratio growing like p (Remark 2), point smoothers degrade under Gauss-Lobatto refinement (Figure 1), and the ILU ordering is critical since AMD-ordered ILU fails while MDF and RCM work (Figure 2, Section 3.4). No analytical estimate is given for the contraction of the ILU-smoothed V-cycle on these anisotropic patches. Therefore, the uniform bound could fail in a regime not covered by the tables. Other aspects of the paper, including the Schwarz decomposition Lemma 2 and the DG norm-equivalence argument, appear internally consistent and are either proved or properly cited to established results, so they are not the weakest point. A concrete high-p test of the local solver would settle whether the concern lands. Since the reader already marked the result CONDITIONAL rather than ACCEPT and the identified risk is unchanged, no verdict adjustment is needed.","tokens_in":25487,"tokens_out":7030,"duration_ms":75708,"concrete_test":"Measure the local solver condition number kappa(R_j A_j) or the V-cycle contraction factor on a single Gauss-Lobatto-refined patch (one element with p-by-p subelements, and one vertex patch), for p = 20, 40, 80, 160 in 2D and 3D, using the same MDF-ordered ILU(0) pre/post-smoothing and coarse solver as in the paper. If kappa(R_j A_j) grows systematically with p, or if the full preconditioner CG iteration counts in Tables 2 and 3 grow beyond a bounded plateau past p=20, assumption (17) is violated and Corollary 1 is not supported. A growth rate like O(p) would refute the uniform bound; a flat curve would strengthen, though not prove, uniformity.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 3.2 derives Theorem 1 and Corollary 1 conditional on (17): there exist C0,C1 independent of h,p such that C0 a(uj,uj) <= a(Rj Aj uj,uj) <= C1 a(uj,uj) for every local patch solver. The paper states \"we will assume that Rj is a uniform preconditioner for Aj\" and never proves this for the specific element-structured V-cycle with MDF-ordered ILU(0) smoothing introduced in Section 3.4. The Schwarz decomposition analysis (Lemma 2, Theorem 1) is sound for the chosen patches, but it only transfers the burden to the unproved local-solver uniformity. Figures 1-2 and Tables 2-3 provide numerical evidence up to p=20 on Cartesian and unstructured meshes; however, the local mesh aspect ratio grows like p (Remark 2), and Section 3.4 itself shows that the choice of ordering matters (AMD-ordered ILU degrades while MDF/RCM do not), so the uniformity of MDF-ordered ILU(0) as p grows is a genuine open question rather than a consequence of the analysis. If some p or aspect ratio defeats the ILU smoother, as AMD ordering already does at moderate p, then the local solvers cease to be uniform in p, and the constant C in Corollary 1 either grows with p or does not exist. The stated corollary would then be false for the implemented preconditioner, even though the Schwarz framework itself remains stable.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper designs and analyzes a matrix-free preconditioner for high-order continuous and discontinuous Galerkin discretizations of elliptic problems. The method is based on the FEM-SEM equivalence: a high-order operator is preconditioned by a low-order (p=1) finite element operator on a refined mesh of Gauss-Lobatto points. The resulting anisotropic low-order system is solved by an overlapping additive Schwarz method with vertex patches, using an element-structured geometric multigrid V-cycle with ordered ILU(0) smoothing as local solvers. The analysis proves stability of the Schwarz decomposition on the p-anisotropic mesh (Lemma 2) and derives a condition-number bound (Theorem 1) that is uniform in h and p conditional on assumption (17), namely that each local solver R_j is a uniform preconditioner for its subspace operator. The method is extended to interior-penalty and BR2 DG discretizations, where a norm-equivalence argument yields a bound independent of the penalty parameter (Proposition 3 and Corollary 2). Extensive numerical experiments on Cartesian, unstructured, variable-coefficient, and anisotropic meshes, as well as DG penalty studies and timing scalings, support the practical claims.","tokens_in":25862,"tokens_out":12577,"duration_ms":122131,"significance":"If the local-solver uniformity in (17) holds, the paper delivers a practical, matrix-free solver achieving bounded iteration counts under h- and p-refinement and under DG penalty variation, with near-optimal setup and application costs. The analysis of the Schwarz decomposition on p-anisotropic meshes in Lemma 2 is a genuine technical contribution, and the DG norm-equivalence argument in Proposition 3 is clean. The numerical study is unusually thorough for a solver paper, covering many of the regimes that matter in practice. The main caveat is that the key uniformity of the ILU-smoothed multigrid local solver is assumed rather than proved; this is the single load-bearing assumption on which the headline robustness claims rest.","major_comments":[{"comment":"The p-robust bound κ(BA_p) ≤ C in Corollary 1 is obtained from Theorem 1 together with the FEM-SEM equivalence, but Theorem 1 is conditional on assumption (17): each local solver R_j must satisfy C0 a(u_j,u_j) ≤ a(R_j A_j u_j, u_j) ≤ C1 a(u_j,u_j) with C0, C1 independent of h and p. The manuscript explicitly says in Section 3.2 that 'we will assume that R_j is a uniform preconditioner for A_j' and provides no proof of this for the specific element-structured V-cycle with MDF-ordered ILU(0) smoothing developed in Sections 3.3 and 3.4. This gap is load-bearing because Section 3.4 itself shows that uniformity is not generic: AMD-ordered ILU degrades while MDF/RCM orderings do not, so the property is not a consequence of the general Schwarz framework. Consequently, Corollary 1 as stated overclaims an unconditional bound for the implemented preconditioner; it should either be accompanied by a proof of (17) or be reworded to present the bound as conditional on a numerically verified property.","section":"Section 3.2, Eq. (17); Theorem 1; Corollary 1"},{"comment":"The DG result κ(B_DG A_BR2) ≤ C in Corollary 2 inherits the same unproved assumption. Theorem 2 requires R_p to be a uniform approximation to A_p^{-1}, and the paper establishes this hypothesis only through Theorem 1, which itself rests on (17). Therefore the claimed robustness with respect to h, p, and the penalty parameter η is conditional on the same local-solver uniformity. If (17) fails in some parameter regime, the DG corollary and the associated robustness claims would also fail. The author should make this conditional structure explicit in the statements of the DG results.","section":"Section 3.5; Theorem 2; Corollary 2"},{"comment":"The numerical evidence offered in support of (17) is supportive but not conclusive. In the vertex-patch tests on the 32×32 mesh, the CG iteration count increases from 26 at p=2 to 38 at p=20, a growth of about 46%. Moreover, at p=2 the iteration count grows steadily under h-refinement (4, 10, 14, 20, 26 for nx=2 through 32). The paper characterizes these trends as pre-asymptotic, but the data are also consistent with slow, unbounded growth of the condition number. If the author intends the headline claims to be unconditional, then a proof of (17) or at least a much more extensive empirical study (larger p, direct condition-number estimates, and a study of the dependence on aspect ratio) is needed. Otherwise, the wording should state explicitly that the observed robustness is an empirical finding rather than a proven uniform bound.","section":"Section 4.1; Table 3"}],"minor_comments":[{"comment":"The abstract and conclusions state that the preconditioner is 'robust in h and p' and 'robust in the penalty parameter' without qualification. These statements should be aligned with the conditional nature of Theorem 1 and Corollaries 1-2, for example by adding a clause such as 'assuming the local solvers are uniform preconditioners, which is verified numerically.'","section":"Abstract and Section 5"},{"comment":"In the vertex-patch column of Table 3, the p=2 row shows a monotonic increase in iteration counts with h-refinement (4, 10, 14, 20, 26 for nx=2 to 32). Even if this is a pre-asymptotic effect, the author should comment on whether the growth persists at larger nx and how it interacts with the claimed h-robustness of the parallel version.","section":"Table 3 and Section 4.1"},{"comment":"The symbol h is used both for the coarse-mesh spacing (in Lemma 2 and the product-rule bound) and for the refined-mesh element size (in Remark 2). A brief sentence distinguishing these two uses, e.g., h for the coarse mesh and ~h/p for the refined mesh, would improve readability.","section":"Lemma 2 and Remark 2"},{"comment":"The step in which Ih is said to be bounded on Q2([0,1]^d) in the H^1_i seminorms is stated as 'it is clear'; a short justification with the explicit constant c_I (for instance by a direct calculation on the reference cube) would make the proof more self-contained.","section":"Proof of Lemma 2"}],"recommendation":"major_revision","confidential_remarks":"The paper is a solid contribution to matrix-free preconditioning of high-order methods, and the conditional analysis plus numerical study make a case for publication. The single most important issue is assumption (17): the uniformity of the ILU-smoothed multigrid local solver is asserted without proof, and it is the key to the advertised p-robustness. I would suggest the editor ask the author either to prove (17) for the proposed smoother (possibly under structural assumptions on Gauss-Lobatto meshes and the MDF ordering) or to clearly downgrade the unconditional wording of Corollaries 1-2 and the abstract. A related, milder concern is the apparent h-growth in Table 3 at low p, which should be addressed explicitly."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The short version: this is a well-engineered, honest paper on matrix-free preconditioners for high-order CG and DG, but the uniform p-robustness claim is conditional on an unproved assumption about the local solver. I think it deserves a serious referee.\n\nWhat is actually new: the paper combines known ingredients (FEM-SEM equivalence, additive Schwarz, multigrid, ILU) in a non-routine way. The element-structured geometric multigrid with MDF-ordered ILU(0) smoothing works as an algebraic substitute for line smoothing on unstructured meshes, and the unified treatment of CG, IP, and BR2 with penalty robustness is clean. Lemma 2 is the real theoretical contribution: it proves a stable decomposition for the p-anisotropic refined mesh, which is nontrivial because the aspect ratio grows like p. The numerical study is broad: h- and p-refinement, variable coefficients, curved meshes, DG penalties up to 10^4, and a direct comparison with BoomerAMG showing AMG degrades with p.\n\nThe soft spots, in proportion. The main one is exactly what the stress-test note flags: Equation (17) is an assumption, not a theorem. The paper explicitly says “we will assume that Rj is a uniform preconditioner for Aj,” and the proof of Theorem 1 transfers the burden to that assumption. The numerical evidence up to p=20 is consistent, and the paper is honest enough to show that ordering matters (AMD fails), so the positive evidence is not trivial. Still, Corollaries 1 and 2 are conditional on that local-solver uniformity, and a single p or aspect-ratio regime where MDF-ordered ILU degrades would break the headline claim. A referee should push for a proof in a simplified model setting or at least a sharper empirical characterization. Two smaller things: the theory covers affine meshes while curved elements are empirically validated only, and no code or data are provided, which makes independent reproduction harder.\n\nNone of this is fatal. The paper does not overclaim: it labels the local-solver assumption as such, and the Schwarz framework plus Lemma 2 stand on their own. This is a useful contribution for anyone working on high-order matrix-free solvers, especially in memory-constrained environments.\n\nRecommendation: send it to peer review. Ask the authors to clarify the status of (17), add a remark on what would be needed to prove or disprove it, and consider releasing a small reproducible example. I would also bring it to a reading group; the discussion of ILU orderings in the context of anisotropy is worth seeing.","headline":"Solid, honest paper on matrix-free high-order preconditioning, but the headline p-robustness rests on an explicit unproved assumption about the local ILU-smoothed V-cycle.","tokens_in":26314,"tokens_out":2778,"would_cite":true,"duration_ms":30220,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N30","65N55","65F08","65N35"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper constructs a matrix-free preconditioner whose condition number stays bounded as mesh size, polynomial degree, and DG penalty vary.","keywords":["matrix-free preconditioning","high-order finite elements","discontinuous Galerkin","additive Schwarz","FEM-SEM equivalence","geometric multigrid","ILU smoothing","interior penalty"],"falsifier":"Fix a $2\\times2$ coarse mesh with constant coefficient and a single patch covering the whole domain, then run the proposed preconditioner on the Poisson problem at $p=40, 80, 160$ and count conjugate gradient iterations to a fixed residual reduction, say $10^{-8}$. If the count grows without bound as $p$ increases, the uniform condition number bound of Corollary 1 is false, since the only unproved ingredient is the uniform local-solver assumption (17). Equivalently, compute the extreme eigenvalues of $R_j A_j$ on one patch and check whether they stay within fixed bounds as $p$ grows.","tokens_in":25292,"feed_emoji":"🧮","tokens_out":9962,"duration_ms":87620,"temperature":0.7,"pith_summary":"Solving high-order finite element and discontinuous Galerkin systems usually costs too much memory if the system matrix is assembled, since the number of nonzeros grows like $O(p^{2d})$ in dimension $d$ and polynomial degree $p$. This paper constructs a preconditioner that avoids forming the high-order matrix: it replaces the degree-$p$ operator by a low-order ($p=1$) operator on a mesh refined at Gauss-Lobatto points, which is equivalent in energy with constants independent of $h$ and $p$, and solves that sparse problem with an overlapping additive Schwarz method whose local solvers are geometric multigrid V-cycles with ordered ILU smoothing. The central claim is that the condition number of the preconditioned system stays bounded as the mesh size $h$, the polynomial degree $p$, and (for interior penalty and BR2 discontinuous Galerkin methods) the penalty parameter $\\eta$ vary. A sympathetic reader would care because bounded condition numbers mean a bounded number of conjugate gradient iterations, so the solve time scales like the matrix-free operator application and the memory cost is linear in the number of degrees of freedom.","feed_headline":"Bounded solver iterations for matrix-free high-order finite elements","feed_subtitle":"Low-order refined preconditioners keep conjugate gradient steps flat as h, p, and penalty vary.","key_machinery":"The load-bearing object is the low-order refined operator $K_h$: a $p=1$ finite element stiffness matrix on the mesh obtained by subdividing each high-order element into parallelepipeds whose vertices are adjacent Gauss-Lobatto nodes. FEM-SEM equivalence (Proposition 2) makes $K_h$ spectrally equivalent to the high-order $K_p$, so a uniform solver for $K_h$ is a uniform preconditioner for $K_p$. The solver is the additive Schwarz preconditioner $B=\\sum_j R_j Q_j$ from (13), with one coarse space and one subspace per overlapping vertex patch; the local $R_j$ are geometric multigrid V-cycles on the structured Gauss-Lobatto hierarchy using minimum-discarded-fill ordered ILU(0) smoothing, which plays the role of line relaxation on unstructured meshes. Lemma 2 is the key technical step: it bounds the energy of the interpolated partition-of-unity pieces by $h^{-2}\\|w\\|_{L^2}^2 + |w|_{H^1}^2$ with constants independent of element aspect ratio, which is what lets the Schwarz theory survive the $O(p)$ anisotropy. For DG, the machinery is the decomposition $V_{DG} = V_B + V_p$, with a point Jacobi solve on the bubble space $V_B$ and the conforming additive Schwarz solver on $V_p$; norm equivalence of IP and BR2 (Proposition 3) transfers the bound.","core_discovery":"On its own terms, the paper establishes two uniform bounds. For continuous Galerkin, Corollary 1 asserts that the additive Schwarz preconditioner $B$ defined by (13) satisfies $\\kappa(B A_p) \\le C$ with $C$ independent of $h$ and $p$, where $A_p$ is the high-order stiffness operator. For discontinuous Galerkin, Corollary 2 asserts $\\kappa(B_{DG} A_{BR2}) \\le C$ independent of $h$, $p$, and the penalty $\\eta$, and Theorem 2 gives the same for the interior penalty form. The argument chains three ingredients: FEM-SEM spectral equivalence between $K_p$ and the low-order refined $K_h$; a Schwarz space decomposition whose stability is proved in Lemma 2 despite the non-shape-regular refined mesh; and uniform local solvers whose uniformity is assumed in (17) and verified numerically.","pith_inferences":["Editorial extension: the unproved uniformity assumption on the ILU-smoothed multigrid V-cycle, equation (17), is the only step separating the theorem from a fully rigorous $h$-$p$ bound; replacing ILU with exact line solves on the structured patches, or proving a smoothing property for MDF-ordered ILU, would close the gap.","Editorial extension: Lemma 2 needs only isotropic overlap, not shape regularity, so the same additive Schwarz construction should apply to hp-adaptive meshes, boundary-layer meshes, and other anisotropic refinements beyond Gauss-Lobatto grids.","Editorial extension: because the preconditioner never sees the high-order matrix entries, it is natural to test it on high-contrast variable coefficients and on vector-valued systems such as elasticity or Maxwell equations, where the scalar Schwarz theory does not directly apply.","Editorial extension: the paper's cost model assumes one geometric multigrid V-cycle per patch; on massively parallel hardware, the overlap factor and coarse-solve communication could change the practical balance, so a distributed-memory implementation study would test whether the bounded iteration count translates to bounded wall-clock time."],"forward_implications":["For continuous Galerkin, a bounded $\\kappa(BA_p)$ means the number of conjugate gradient iterations is bounded as $h\\to 0$ and $p\\to\\infty$, so total solve time is $O(p\\,n_{\\rm dof})$ with $O(n_{\\rm dof})$ memory.","For interior penalty and BR2 discontinuous Galerkin methods, the same boundedness holds as the penalty $\\eta$ grows, so the stabilization parameter can be chosen for accuracy rather than solver convenience.","The low-order refined matrices are sparse with a constant number of nonzeros per row, independent of $p$, which is what makes forming them cheap in time and memory.","The vertex-patch additive Schwarz version parallelizes the local solves, with iteration counts roughly 1.5 to 2 times the single-patch case in the reported tables.","Numerical tests on variable coefficients, curved unstructured meshes, and anisotropic element strips show the same bounded iteration counts, provided subdomain overlap is made isotropic on anisotropic meshes as described in Section 4.5."],"supporting_citations":[{"why":"Supplies the FEM-SEM spectral equivalence between the high-order and low-order refined operators (Proposition 2), the foundation for replacing $K_p$ by $K_h$.","marker":"[23]"},{"why":"Provides the one-dimensional $L^2$ and $H^1$ interpolation equivalences at Gauss-Lobatto points used to prove the spectral equivalence.","marker":"[18]"},{"why":"Supplies the additive Schwarz framework and stable decomposition estimates used in Theorem 1.","marker":"[29]"},{"why":"Domain decomposition theory used for the subspace correction analysis and condition number bounds.","marker":"[59]"},{"why":"Parallel subspace correction framework justifying the additive Schwarz preconditioner and the $c_1/c_0$ bound.","marker":"[63]"},{"why":"Lions' lemma, which reduces the lower bound to constructing a stable decomposition of the refined space.","marker":"[43]"},{"why":"Supplies the DG space decomposition $V_{DG}=V_B+V_p$ and the penalty-independent IP preconditioning result stated as Theorem 2.","marker":"[6]"},{"why":"Provides the trace and inverse-trace inequalities used to prove equivalence of the IP and BR2 norms.","marker":"[17]"},{"why":"Supplies boundedness and coercivity of the BR2 form with lifting operators, used in Lemma 3 and Corollary 2.","marker":"[13]"}],"fun_headline_variants":["Preconditioner keeps CG iterations flat for high-order FEM","Robust preconditioner for h, p, and penalty in DG and CG","Uniform bounds for high-order FEM via p=1 preconditioning","Two corollaries: uniform preconditioner for CG and DG"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"Everything rests on the assumption that the local approximate solve on each overlapping patch is uniformly accurate as the mesh is refined and the polynomial degree grows; the paper demonstrates this by experiment, but does not prove it.","fun_headline_variants_meta":{"raw":{"variants":["Preconditioner keeps CG iterations flat for high-order FEM","Robust preconditioner for h, p, and penalty in DG and CG","Uniform bounds for high-order FEM via p=1 preconditioning","Two corollaries: uniform preconditioner for CG and DG"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000676,"raw_usage":{"total_tokens":3059,"prompt_tokens":915,"completion_tokens":2144,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":531,"completion_tokens_details":{"reasoning_tokens":2071}},"tokens_in":531,"tokens_out":2144,"duration_ms":13680,"temperature":1.0,"reasoning_tokens":2071,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T12:27:12.319896+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Fix a $2\\times2$ coarse mesh with constant coefficient and a single patch covering the whole domain, then run the proposed preconditioner on the Poisson problem at $p=40, 80, 160$ and count conjugate gradient iterations to a fixed residual reduction, say $10^{-8}$. If the count grows without bound as $p$ increases, the uniform condition number bound of Corollary 1 is false, since the only unproved ingredient is the uniform local-solver assumption (17). Equivalently, compute the extreme eigenvalues of $R_j A_j$ on one patch and check whether they stay within fixed bounds as $p$ grows.","supporting_citations":[{"cited_title":"Canuto , Stabilization of spectral methods by finite element bubble functions , Computer Meth- ods in Applied Mechanics and Engineering, 116 (1994), pp","cited_arxiv_id":null,"evidence_quote":"Provides the one-dimensional $L^2$ and $H^1$ interpolation equivalences at Gauss-Lobatto points used to prove the spectral equivalence."},{"cited_title":"Dryja and O","cited_arxiv_id":null,"evidence_quote":"Supplies the additive Schwarz framework and stable decomposition estimates used in Theorem 1."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Lions' lemma, which reduces the lower bound to constructing a stable decomposition of the refined space."},{"cited_title":"Brezzi, G","cited_arxiv_id":null,"evidence_quote":"Supplies boundedness and coercivity of the BR2 form with lifting operators, used in Lemma 3 and Corollary 2."}],"review_version":1}