Automated Galerkin time stepping in Irksome

Problem Statement
The discretization of time-dependent partial differential equations (PDEs) has traditionally relied on Runge–Kutta (RK) methods. While effective, these methods often struggle to preserve non-quadratic invariants and dissipation laws inherent in geometric numerical integration, particularly for Hamiltonian and dissipative systems. Although Galerkin-in-time (or variational-in-time) discretizations, such as Discontinuous Galerkin (DG) and Continuous Petrov–Galerkin (CPG) schemes, offer superior structure-preserving properties, their adoption has been hindered by a lack of general-purpose software. Implementations of these schemes have historically been bespoke, cumbersome, and limited to specific applications, often requiring the "repurposing" of spatial discretization tools not designed for time integration. Furthermore, the algebraic systems resulting from these methods are monolithic, coupling all temporal degrees of freedom within a time step, which complicates the derivation of residuals, Jacobians, and efficient solvers.

Methodology
This work extends the Irksome library, a time-stepping tool for the Firedrake finite element framework, to automate both DG and CPG time stepping for semidiscrete variational problems. The core methodology relies on syntactic manipulation of the Unified Form Language (UFL).

1. Variational Formulations:

   • DG Schemes: The method discretizes the time derivative and jump terms over a time interval $I_n$ using a polynomial space of degree $s$. It employs a "strong" form derived via integration by parts, handling discontinuities at time steps through jump terms.
   • CPG Schemes: These schemes enforce continuity of the solution across time steps and strongly enforce initial data. The test functions are sought in a lower-degree discontinuous space ($s-1$), allowing the method to be applied to general semidiscrete forms, including those without explicit time derivatives on all variables.
   • Auxiliary Variables: To handle differential-algebraic equations (DAEs) and mixed formulations (e.g., incompressible Navier–Stokes), the framework introduces auxiliary variables (such as pressure) in discontinuous-in-time spaces of degree $s-1$, distinct from the primal variables in continuous spaces of degree $s$.

2. Automation Strategy:

   • The user provides a semidiscrete variational form in UFL. Irksome automatically constructs the coupled stage equations by expanding the solution and test functions in temporal bases.
   • Quadrature weights are absorbed into the test functions via linearity to simplify the UFL manipulation.
   • The system is assembled into a block-structured algebraic system where the temporal stages are coupled, similar to fully implicit RK methods.
   • Strong boundary conditions are handled via $L^2$ projections onto the temporal basis.

3. Algebraic Solvers:

   • The paper leverages Firedrake's solver infrastructure to apply monolithic multigrid methods to the resulting stage-coupled systems.
   • It utilizes additive Schwarz relaxation schemes (specifically Vanka-type decompositions) adapted for the block structures of DG and CPG systems, allowing for scalable solutions of large-scale problems.

Key Contributions

• First General-Purpose Package: This work presents the first general-purpose package to natively support both DG and CPG time stepping with a unified interface.
• Seamless Integration: Users can switch between Runge–Kutta, DG, and CPG schemes with minimal code changes (often a single line), maintaining the same high-level UFL interface.
• Support for DAEs and Mixed Systems: The implementation explicitly supports auxiliary variables, enabling the structure-preserving discretization of complex systems like incompressible and compressible flows where variables have different continuity requirements in time.
• Flexible Quadrature: The framework allows for flexible temporal quadrature, enabling the user to choose different rules for different terms to optimize cost or preserve specific conservation laws.

Numerical Results
The paper validates the implementation through four representative PDE systems:

1. Gross–Pitaevskii Equation (Hamiltonian System): Demonstrates that CPG schemes exactly conserve energy (up to solver tolerance) but not wave function normalization, whereas implicit Gauss–Legendre RK methods conserve normalization but not energy.
2. Allen–Cahn Equation (Gradient Descent System): Shows that CPG schemes preserve unconditional energy dissipation. DG schemes also exhibit monotonic energy decay, though the theoretical guarantee is more complex.
3. Incompressible Navier–Stokes (Flow Past a Cylinder):
   • Validates the implementation against benchmark drag and lift data.
   • Highlights that treating pressure as a discontinuous auxiliary variable is crucial for CPG convergence in DAE systems.
   • Identifies that standard CPG can accumulate divergence errors; a stabilization technique using right-Radau quadrature for the incompressibility constraint is proposed and shown to restore convergence.
   • Comparisons show DG generally yields lower errors for this specific problem, but CPG remains competitive with appropriate stabilization.
4. Compressible Navier–Stokes–Fourier: Demonstrates the preservation of mass, momentum, energy, and entropy generation in a supersonic perturbation test, properties that traditional RK methods fail to maintain.

Significance and Claims
The authors claim that this work removes the barrier to entry for using high-order, structure-preserving Galerkin-in-time methods. By automating the derivation of complex, stage-coupled variational forms and integrating them with scalable monolithic solvers, Irksome enables researchers to explore the theoretical benefits of these schemes (such as conservation of non-quadratic invariants) in practical, large-scale PDE simulations. The paper concludes that while the software successfully realizes these structure-preserving properties in practice, challenges remain in developing even more efficient linear and nonlinear solvers specifically tailored to the unique block structures of these time-stepping schemes.