{"id":"eab1c973-19d6-4709-9d48-d038d50a4818","arxiv_id":"2508.20255","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"Irksome now automates Runge-Kutta-Nyström time stepping for second-order PDEs in Firedrake, producing smaller stage-coupled systems and faster solves than first-order Runge-Kutta methods.","lead":"Researchers extended the Irksome finite element library to automatically build Runge-Kutta-Nyström time-stepping schemes for PDEs with second-order time derivatives. This avoids rewriting the equations as larger first-order systems, and makes fully implicit time stepping competitive with explicit methods for wave, plate, and poroelasticity problems.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Mixed-order RKN claim rests on an unproved assertion that zero non-physical initial values for w_t and p_t do not affect accuracy, while Section 6.3 concedes the DAE theory does not apply.","rationale":"The paper is a solid software-extension contribution: the RKN stage reduction in Section 2 (equations 2.13-2.18) is standard and mathematically correct, the reduction from two stage vectors to one for a 2-stage method on the wave equation is real, and the plate experiments (Figure 9) support the claimed roughly 2x run-time advantage of RKN over first-order RK. The wave-equation results (Figure 6) support the claim that fully implicit GL(2) can compete with explicit methods for this linear constant-coefficient problem, with the acknowledged caveat that matrices and preconditioners are reused across steps. Code is available via Zenodo, which is genuine independent support for reproducibility. The load-bearing weak point is the mixed-order extension that underlies the poroelasticity section and the abstract's 'mixed temporal order' claim. The paper's own Section 6.3 flags that the DAE convergence theory does not directly apply to the RKN schemes, and Section 3's assertion that non-physical zero initial values for w_t and p_t do not affect accuracy is asserted without proof. The RKN stage replacement was derived for pure second-order systems; applying it to first-order variables introduces an O(dt) perturbation at the initial step whose effect on order is tableau-dependent and unanalyzed. The numerical evidence is consistent with the authors' heuristic, but the product norm in Figures 10-11 has limited spatial order that can mask temporal-order defects, so the experiments do not fully settle the question. A targeted experiment using true initial values for w_t and p_t, plus temporal-error-dominated refinement, would settle whether the assertion holds. The abstract's '(essentially) required' phrasing is also overstated: the paper's own Figure 9b shows first-order RK works for the plate, merely about 2x slower, so RKN is not required in any strict sense; this is a rhetorical overreach rather than a technical failure. The dependence of the poroelasticity solver results on the unpublished Adler et al. 2025 solver is a reproducibility caveat, but it does not affect the RKN-versus-RK comparison, which uses the same solver for both. These considerations support the reader's CONDITIONAL verdict: the central construction is correct and the performance advantage is real, but the mixed-order accuracy claim and the abstract's overstatements should be addressed before the paper is treated as a definitive benchmark.","tokens_in":24750,"tokens_out":14158,"duration_ms":119198,"concrete_test":"Re-run the Section 6.3 manufactured-solution Biot experiments (Figures 10-11) with w_t(0) and p_t(0) set to the values obtained by differentiating the closed-form solution and evaluating (3.16) at t = 0, instead of zero, keeping all other settings identical. If product-norm errors and convergence orders are unchanged to solver tolerance, the Section 3 assertion holds for these tests; if errors shift by more than the observed grid-to-grid differences or orders drop, the non-physical initialization contaminates the results and the abstract's mixed-order claim must be softened.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's scope claim (Section 2) that RKN methods handle \"coupled systems mixing first- and second-order time derivatives\" and the abstract's poroelasticity claim both depend on the RKN treatment of the Biot system (3.16), where u has u_tt but w and p have only first-order time derivatives. Section 3 applies the RKN stage replacement to these first-order variables by introducing w_t and p_t as if they were velocities, requiring \"non-physical initial values for w_t and p_t, which we take to be zero,\" and asserting \"the values specified do not impact the resulting accuracy.\" This assertion is unproved, and Section 6.3 explicitly concedes that the cited DAE convergence theory \"does not directly apply to the Runge-Kutta-Nyström schemes used here.\" The stage approximations for w and p contain the term c_i * dt * (w_t,0 or p_t,0), so zero initialization perturbs the stage values by O(dt) at the first steps; whether stiff accuracy (Radau IIA), symplecticity (Gauss), or other tableau structure absorbs this perturbation, or whether it degrades the formal order, is exactly the missing analysis. The observed second-order (Radau) and first-order (Gauss) rates in Figures 10-11 are consistent with the authors' heuristic, but the product norm has limited spatial order (second order for the spaces used, with n_t = n_x/4), which can mask temporal-order defects; the agreement is therefore not decisive evidence. If the assertion fails, the poroelasticity accuracy results and the \"perform effectively even as we approach the incompressible limit\" claim are unreliable, and the general claim that the RKN formalism extends to mixed-order systems is unsupported.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper reports an extension of the Irksome/Firedrake time-stepping framework to Runge–Kutta–Nyström (RKN) methods for semidiscrete finite element discretizations of PDEs containing second-order time derivatives. Section 2 derives the RKN stage equations from an underlying RK tableau via A-bar = A^2 and b-bar = A^T b, yielding one set of stage unknowns rather than two. Section 3 presents wave, telegraph, dynamic plate, and fully dynamic Biot poroelasticity formulations; Section 4 surveys stage-coupled preconditioners; Section 5 describes UFL-based automated generation of stage-coupled variational problems and their PETSc solver integration. Numerical experiments compare RKN with explicit central differences and with RK applied to first-order forms, including energy conservation and a factor-of-two speedup for RKN over RK. The paper claims smaller algebraic systems, better runtime, competitiveness with explicit methods for wave equations, and essentiality of RKN for higher-order spatial derivatives.","tokens_in":24971,"tokens_out":10723,"duration_ms":96000,"significance":"The core derivation is standard, parameter-free, and cleanly presented, and the software contribution is potentially valuable: a user can write a UFL form with Dt(u,2) and obtain an implicit RKN discretization with automatically generated stage-coupled variational problems. The energy-conservation benchmark for Gauss-Legendre methods and the reproducible Zenodo artifact are concrete strengths. If the mixed-order poroelasticity extension is rigorously supported, the paper would demonstrate a genuinely broader applicability. At present, however, the central novelty claim for mixed first-/second-order systems rests on an unproved assertion, and the abstract overstates the necessity of RKN, so the contribution is promising but not yet fully established.","major_comments":[{"comment":"The extension of the RKN stage replacement to the mixed first-/second-order Biot system is asserted rather than proved. The text states that introducing non-physical initial values for w_t and p_t and taking them to be zero does not affect accuracy, but no order or stability analysis is given for the resulting mixed-order partitioned scheme, and Section 6.3 explicitly concedes that the cited DAE convergence theory \"does not directly apply to the Runge–Kutta–Nyström schemes used here.\" This is load-bearing because the abstract and Section 2 advertise coupled systems mixing first- and second-order time derivatives, and the poroelasticity experiments are a central validation. The authors should supply either a convergence analysis for this RKN treatment of first-order variables or numerical evidence that isolates temporal errors, for example temporal refinement at a fixed fine spatial mesh and a comparison of zero versus consistent non-zero initial values for w_t and p_t.","section":"Section 3 (Biot system, Eq. (3.16)) and Section 6.3"},{"comment":"The abstract's claim that RKN methods are \"(essentially) required to discretize wave-type equations with higher-order spatial derivatives\" is contradicted by the paper's own results: Figure 9 shows that the equivalent first-order system with RK-GL(2) also solves the plate problem, with RKN being roughly twice as fast but not uniquely capable. The unconditional stability of Gauss-Legendre is a property of the underlying RK method and is available in the first-order formulation as well. The wording should be softened to a comparative claim, such as \"substantially more efficient in this setting,\" rather than a claim of necessity.","section":"Abstract and Section 6.2"},{"comment":"The observed temporal convergence rates in the Biot experiments are not decisive, because the study refines time and space simultaneously with n_t = n_x/4, and the product norm has at most second-order spatial accuracy for the chosen elements, as the text itself acknowledges. A first-order temporal defect could be masked by a second-order spatial error floor. A temporal-only refinement study at a fixed fine spatial mesh would make the claimed RKN orders (second order for Radau, first order for Gauss) convincing and would also directly test the sensitivity to the non-physical initial values for w_t and p_t.","section":"Section 6.3, Figures 10 and 11"},{"comment":"The poroelasticity convergence rates and solver performance rest on the monolithic multigrid solver and DAE theory of [Adler et al. 2025], which is cited as \"in preparation.\" This work is used both to solve the linear systems and to justify the expected convergence rates, so the manuscript is not self-contained on a load-bearing point. The authors should either include sufficient detail on the solver and the applicable convergence theory or clearly mark these results as preliminary and provide alternative supporting evidence.","section":"Section 6.3"}],"minor_comments":[{"comment":"The clamped boundary conditions are written as u|_Omega = ∂u/∂n|_Omega = 0; the restriction should be to ∂Omega, the boundary of the domain, rather than to Omega itself.","section":"Equation (3.13)"},{"comment":"The reference to \"Clough and Toucher 1965\" should read \"Clough and Tocher 1965\" to match the standard name of the Hsieh-Clough-Tocher macroelement.","section":"References"},{"comment":"The product norm contains a term 10^{-6} h_t ‖w‖^2, where h_t is a time step; the mixed length/time scaling of this term should be clarified or the notation adjusted to avoid dimensional confusion.","section":"Section 6.3, product norm definition"},{"comment":"The sentence noting that GL(2) \"produces some inaccuracy with Q2\" on coarse meshes would benefit from a one-sentence explanation, since the phenomenon is not reflected in the asymptotic discussion that follows.","section":"Section 6.1, Figure 6b"}],"recommendation":"major_revision","confidential_remarks":"The main risks are the unproved mixed-order RKN claim in Section 3 and the dependence of Section 6.3 on the unpublished [Adler et al. 2025] solver and theory. I would ask the authors to add a temporal-refinement study, address the sensitivity to zero initial values for w_t and p_t, and soften the \"(essentially) required\" claim in the abstract before publication. The central RKN derivation and the software contribution are sound and valuable."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a solid, useful software paper. The RKN stage reduction (Abar = A^2, bbar = A^T b) is textbook, but the paper's contribution is the automation in Irksome, the dDAE boundary treatment, and the demonstration that RKN gives smaller stage-coupled systems with real runtime savings. The numerical experiments are honest and mostly well constructed. The paper deserves a serious referee.\n\nWhat's genuinely new: Irksome previously automated first-order RK; this extends it to Dt(u,2), handles mixed-order systems like Biot, and adds the differentiated-DAE boundary condition variant. The code is available via Zenodo. The plate experiments show RKN roughly halves runtime versus RK on the first-order system, which is a concrete, credible win.\n\nSoft spots, in order of real weight. First, the abstract claims RKN methods are \"(essentially) required\" for wave-type equations with higher-order spatial derivatives. The paper's own Figure 9 shows RK-GL(2) works fine on the plate problem, just slower. That claim is simply false as stated; it should be softened to \"advantageous\" or \"more efficient.\" Second, the Biot mixed-order treatment: Section 3 introduces non-physical initial values for w_t and p_t, taken to be zero, and asserts they don't affect accuracy. Section 6.3 then concedes the DAE convergence theory does not directly apply. That is a genuine gap. The observed rates in Figures 10-11 are consistent with the heuristic, but the product norm has limited spatial order, which can mask temporal-order defects. If the assertion fails, the poroelasticity accuracy claims would be in trouble. This needs either a proof, a more precise reference, or a much more hedged statement. Minor: the performance conclusions in Section 6.3 lean on an unpublished solver paper (Adler et al. 2025, in preparation); that's acceptable for a preprint but should be resolved before publication.\n\nThe math in Section 2 is clean. The boundary condition taxonomy (ODE, DAE, dDAE) is genuinely useful and the explicit-rank discussion is thoughtful. Citation pattern looks fine; self-cites are to prior Irksome work and are appropriate.\n\nBottom line: conditional acceptance is the right verdict. The core contribution is solid and the flaws are overstatement and one unproved heuristic, not a broken central argument. This should go to peer review. I'd cite it if I were building on Irksome or doing RKN preconditioning.","headline":"Solid Irksome extension automating RKN time stepping with real runtime savings, but the abstract oversells RKN as 'essentially required' and the mixed-order Biot extension rests on an unproved heuristic about zero initial derivatives.","tokens_in":25635,"tokens_out":1888,"would_cite":true,"duration_ms":16567,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65L06","65M20","65M60"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper shows that discretizing second-order time derivatives directly with Runge–Kutta–Nyström methods, rather than rewriting the equation as a first-order system, halves the algebraic unknowns and yields faster finite-element solves…","keywords":["Runge-Kutta-Nyström methods","second-order time derivatives","finite element methods","Irksome","Unified Form Language","wave equation","poroelasticity","implicit Runge-Kutta preconditioning"],"falsifier":"Run the Biot manufactured-solution experiment with non-zero smooth initial values for $\\boldsymbol{w}_t$ and $p_t$ (for example, $\\sin(t)$ evaluated at $t=0$), keeping everything else fixed, and compare the errors and convergence orders; if the errors change at the reported rates, the claim that these artificial initial values do not affect accuracy is false.","tokens_in":24426,"feed_emoji":"🌊","tokens_out":5651,"duration_ms":52524,"temperature":0.7,"pith_summary":"This paper extends the Irksome library so that users can write a finite-element semidiscrete form containing second-order time derivatives and have a Runge–Kutta–Nyström (RKN) method applied automatically. The central claim is that discretizing the second-order form directly, instead of rewriting it as a first-order system and applying an ordinary Runge–Kutta method, produces smaller stage-coupled algebraic systems and better run time. With effective preconditioning, the paper reports, fully implicit RKN methods can match or beat explicit methods for wave equations and are essentially required for wave-type equations with fourth-order spatial derivatives. The paper also reports that the approach handles fully dynamic poroelasticity, a mixed-order system, effectively even near the incompressible limit.","feed_headline":"Implicit wave solves, skipping the first-order rewrite","feed_subtitle":"Runge–Kutta–Nyström methods halve the stage system and beat explicit methods on wave equations.","key_machinery":"The key mechanism is the extended Runge-Kutta-Nyström tableau, where the second-derivative stage equation uses $\\bar A = A^2$ and the update uses $\\bar b = A^T b$; this reduces the stage system to a single set of unknowns. Around that identity, the paper builds an automated UFL preprocessing pass using a DAG traverser that commutes time differentiation with spatial operations, plus two algebraic solver strategies: the Clines block-triangular preconditioner based on an LDU approximation of the stage matrix, and monolithic multigrid with stage-coupled relaxation. The combination turns a user's high-level second-order weak form into a solved stage-coupled system.","core_discovery":"The paper's central discovery is that RKN stage equations can be written using only one set of stage variables approximating $u_{tt}$, with the solution and velocity updates built from the derived tableau coefficients $\\bar A = A^2$ and $\\bar b = A^T b$. This avoids auxiliary velocity variables: for a two-stage method on the wave equation, the stage system has $2\\times\\dim V_h$ unknowns instead of $4\\times\\dim V_h$. Numerically, a two-stage Gauss-Legendre RKN method takes far fewer time steps than explicit central differences or explicit Nyström, and with the Clines block-triangular preconditioner it is faster than both explicit methods on quadratic hexahedral meshes. On an elastodynamic plate problem with Hsieh-Clough-Tocher elements, RKN runs nearly twice as fast as applying the underlying Runge-Kutta method to the first-order system. For the Biot poroelasticity system, RKN with RadauIIA attains expected convergence orders as the Poisson ratio approaches $0.5$, a regime where the first-order Runge-Kutta formulation with the same solver fails to converge on finer meshes.","pith_inferences":["Extrapolating from the reported plate results, other wave-type PDEs currently integrated in first-order form could see similar runtime savings; a natural test would be nonlinear elasticity or acoustic-wave imaging problems.","Because the RKN advantage comes from shrinking the stage system rather than changing spatial discretization, the benefit should grow with the number of stages, so higher-order implicit RKN methods deserve a dedicated preconditioning study.","The poroelasticity conclusions rest on an unpublished monolithic multigrid solver; replacing it with a publicly available robust preconditioner for the first-order form would clarify whether the performance gap is intrinsic to RKN or depends on the solver."],"forward_implications":["Irksome users can write second-order semidiscrete forms directly and receive an RKN stage-coupled variational problem automatically, without hand-coding a time-stepping loop.","With effective preconditioning, fully implicit RKN methods can outperform explicit central differences and explicit Nyström for wave equations at the tested accuracies.","For equations with fourth-order spatial derivatives, implicit RKN avoids the explicit stability restriction $\\Delta t = O(h^2)$ and runs in roughly half the time of the equivalent first-order Runge-Kutta formulation.","For fully dynamic Biot poroelasticity, the RKN-RadauIIA scheme shows the expected convergence orders as $\\nu \\to 0.5$, where the first-order RK formulation with the same monolithic solver fails to converge on finer grids."],"supporting_citations":[{"why":"Introduces the classical explicit RKN method whose extended tableau is used as a benchmark and as an example of a method not derivable from a first-order RK tableau.","marker":"[Nyström 1925]"},{"why":"Supplies the order, stability, and DAE convergence theory that motivates the choice of Gauss-Legendre and RadauIIA tableaux and the expected convergence rates.","marker":"[Wanner and Hairer 1996]"},{"why":"The original Irksome framework for generating stage-coupled variational problems, extended here to second-order and RKN tableaux.","marker":"[Farrell et al. 2021a]"},{"why":"Previous Irksome extension whose preconditioner implementation and UFL manipulation patterns are reused for the RKN stage systems.","marker":"[Kirby and MacLachlan 2025]"},{"why":"Provides the LDU-based block-triangular preconditioner for RKN stage systems used in the wave-equation experiments.","marker":"[Clines et al. 2022]"},{"why":"Establishes monolithic multigrid with stage-coupled relaxation, the foundation for the plate and poroelasticity solvers.","marker":"[Van Lent and Vandewalle 2005]"},{"why":"Supplies the three-field Biot spatial discretization and the time-continuous Galerkin baseline that the RKN poroelasticity experiments build on.","marker":"[Kraus et al. 2025]"},{"why":"Provides the Hsieh-Clough-Tocher macroelement implementation used for the dynamic plate problem.","marker":"[Brubeck and Kirby 2025]"},{"why":"The monolithic multigrid solver with Vanka-style relaxation that carries the poroelasticity convergence and performance results.","marker":"[Adler et al. 2025]"}],"fun_headline_variants":["RKN time stepping halves wave-stage unknowns","Fewer unknowns per stage: RKN outruns explicit wave solvers","Direct second-order discretization shrinks stage systems","One stage set instead of two: RKN accelerates wave PDEs"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The paper assumes that applying the RKN stage replacement to variables whose true time derivative is only first-order, with artificially imposed zero initial values for those derivatives, preserves the expected accuracy; the paper explicitly notes that the available DAE convergence theory does not directly apply to these RKN schemes.","fun_headline_variants_meta":{"raw":{"variants":["RKN time stepping halves wave-stage unknowns","Fewer unknowns per stage: RKN outruns explicit wave solvers","Direct second-order discretization shrinks stage systems","One stage set instead of two: RKN accelerates wave PDEs"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000262,"raw_usage":{"total_tokens":1654,"prompt_tokens":1062,"completion_tokens":592,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":678,"completion_tokens_details":{"reasoning_tokens":524}},"tokens_in":678,"tokens_out":592,"duration_ms":6121,"temperature":1.0,"reasoning_tokens":524,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T16:48:22.331255+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the Biot manufactured-solution experiment with non-zero smooth initial values for $\\boldsymbol{w}_t$ and $p_t$ (for example, $\\sin(t)$ evaluated at $t=0$), keeping everything else fixed, and compare the errors and convergence orders; if the errors change at the reported rates, the claim that these artificial initial values do not affect accuracy is false.","supporting_citations":[],"review_version":1}