Search

Unique continuation principles are fundamental properties of elliptic partial differential equations, giving conditions that guarantee that the solution to an elliptic equation must be uniformly zero. Since finite-element discretizations are a natural tool to help gain understanding into elliptic equations, it is natural to ask if such principles also hold at the discrete level. In this work, we prove a version of the unique continuation principle for piecewise-linear and -bilinear finite-element discretizations of the Laplacian eigenvalue problem on polygonal domains in $\mathbb{R}^2$. Namely, we show that any solution to the discretized equation $-\Delta u = \lambda u$ with vanishing Dirichlet and Neumann traces must be identically zero under certain geometric and topological assumptions on the resulting triangulation. We also provide a counterexample, showing that a nonzero \emph{inner solution} exists when the topological assumptions are not satisfied. Finally, we give an application to an eigenvalue interlacing problem, where the space of inner solutions makes an explicit appearance.

The numerical analysis of higher-order mixed finite-element discretizations for saddle-point problems, such as the Stokes equations, has been well-studied in recent years. While the theory and practice of such discretizations is now well-understood, the same cannot be said for efficient preconditioners for solving the resulting linear (or linearized) systems of equations. In this work, we propose and study variants of the well-known Vanka relaxation scheme that lead to effective geometric multigrid preconditioners for both the conforming Taylor-Hood discretizations and non-conforming ${\bf H}(\text{div})$-$L^2$ discretizations of the Stokes equations. Numerical results demonstrate robust performance with respect to FGMRES iteration counts for increasing polynomial order for some of the considered discretizations, and expose open questions about stopping tolerances for effectively preconditioned iterations at high polynomial order.

Space-time finite-element discretizations are well-developed in many areas of science and engineering, but much work remains within the development of specialized solvers for the resulting linear and nonlinear systems. In this work, we consider the all-at-once solution of the discretized Navier-Stokes equations over a space-time domain using waveform relaxation multigrid methods. In particular, we show how to extend the efficient spatial multigrid relaxation methods from [37] to a waveform relaxation method, and demonstrate the efficiency of the resulting monolithic Newton-Krylov-multigrid solver. Numerical results demonstrate the scalability of the solver for varying discretization order and physical parameters.

This work introduces and assesses the efficiency of a monolithic $ph$MG multigrid framework designed for high-order discretizations of stationary Stokes systems using Taylor-Hood and Scott-Vogelius elements. The proposed approach integrates coarsening in both approximation order ($p$) and mesh resolution ($h$), to address the computational and memory efficiency challenges that are often encountered in conventional high-order numerical simulations. Our numerical results reveal that $ph$MG offers significant improvements over traditional spatial-coarsening-only multigrid ($h$MG) techniques for problems discretized with Taylor-Hood elements across a variety of problem sizes and discretization orders. In particular, the $ph$MG method exhibits superior performance in reducing setup and solve times, particularly when dealing with higher discretization orders and unstructured problem domains. For Scott-Vogelius discretizations, while monolithic $ph$MG delivers low iteration counts and competitive solve phase timings, it exhibits a discernibly slower setup phase when compared to a multilevel (non-monolithic) full-block-factorization (FBF) preconditioner where $ph$MG is employed only for the velocity unknowns. This is primarily due to the setup costs of the larger mixed-field relaxation patches with monolithic $ph$MG versus the patch setup costs with a single unknown type for FBF.

Irksome is a library based on the Unified Form Language (UFL) that enables automated generation of Runge--Kutta methods for time-stepping finite element spatial discretizations of partial differential equations (PDE). Allowing users to express semidiscrete forms of PDE, it generates UFL representations for the stage-coupled variational problems to be solved at each time step. The Firedrake package then generates efficient code for evaluating these variational problems and allows users a wide range of options to deploy efficient algebraic solvers in PETSc. In this paper, we describe several recent advances in Irksome. These include alternate formulations of the Runge--Kutta time-stepping methods and optimized support for diagonally implicit (DIRK) methods. Additionally, we present new and improved tools for building preconditioners for the resulting linear and linearized systems, demonstrating that these can lead to efficient approaches for solving fully implicit Runge-Kutta discretizations. The new features are demonstrated through a sequence of computational examples demonstrating the high-level interface and obtained solver performance.

In recent years, solvers for finite-element discretizations of linear or linearized saddle-point problems, like the Stokes and Oseen equations, have become well established. There are two main classes of preconditioners for such systems: those based on block-factorization approach and those based on monolithic multigrid. Both classes of preconditioners have several critical choices to be made in their composition, such as the selection of a suitable relaxation scheme for monolithic multigrid. From existing studies, some insight can be gained as to what options are preferable in low-performance computing settings, but there are very few fair comparisons of these approaches in the literature, particularly for modern architectures, such as GPUs. In this paper, we perform a comparison between a block-triangular preconditioner and a monolithic multigrid method with the three most common choices of relaxation scheme - Braess-Sarazin, Vanka, and Schur-Uzawa. We develop a performant Vanka relaxation algorithm for structured-grid discretizations, which takes advantage of memory efficiencies in this setting. We detail the behavior of the various CUDA kernels for the multigrid relaxation schemes and evaluate their individual arithmetic intensity, performance, and runtime. Running a preconditioned FGMRES solver for the Stokes equations with these preconditioners allows us to compare their efficiency in a practical setting. We show monolithic multigrid can outperform block-triangular preconditioning, and that using Vanka or Braess-Sarazin relaxation is most efficient. Even though multigrid with Vanka relaxation exhibits reduced performance on the CPU (up to $100\%$ slower than Braess-Sarazin), it is able to outperform Braess-Sarazin by more than $20\%$ on the GPU, making it a competitive algorithm, especially given the high amount of algorithmic tuning needed for effective Braess-Sarazin relaxation., Comment: submitted to IJHPCA

In this paper, we are concerned with efficiently solving the sequences of regularized linear least squares problems associated with employing Tikhonov-type regularization with regularization operators designed to enforce edge recovery. An optimal regularization parameter, which balances the fidelity to the data with the edge-enforcing constraint term, is typically not known a priori. This adds to the total number of regularized linear least squares problems that must be solved before the final image can be recovered. Therefore, in this paper, we determine effective multigrid preconditioners for these sequences of systems. We focus our approach on the sequences that arise as a result of the edge-preserving method introduced in [6], where we can exploit an interpretation of the regularization term as a diffusion operator; however, our methods are also applicable in other edge-preserving settings, such as iteratively reweighted least squares problems. Particular attention is paid to the selection of components of the multigrid preconditioner in order to achieve robustness for different ranges of the regularization parameter value. In addition, we present a parameter culling approach that, when used with the L-curve heuristic, reduces the total number of solves required. We demonstrate our preconditioning and parameter culling routines on examples in computed tomography and image deblurring.

We investigate a novel monolithic algebraic multigrid (AMG) preconditioner for the Taylor-Hood ($\pmb{\mathbb{P}}_2/\mathbb{P}_1$) and Scott-Vogelius ($\pmb{\mathbb{P}}_2/\mathbb{P}_1^{disc}$) discretizations of the Stokes equations. The algorithm is based on the use of the lower-order $\pmb{\mathbb{P}}_1\text{iso}\kern1pt\pmb{\mathbb{P}}_2/\mathbb{P}_1$ operator within a defect-correction setting, in combination with AMG construction of interpolation operators for velocities and pressures. The preconditioning framework is primarily algebraic, though the $\pmb{\mathbb{P}}_1\text{iso}\kern1pt\pmb{\mathbb{P}}_2/\mathbb{P}_1$ operator must be provided. We investigate two relaxation strategies in this setting. Specifically, a novel block factorization approach is devised for Vanka patch systems, which significantly reduces storage requirements and computational overhead, and a Chebyshev adaptation of the LSC-DGS relaxation is developed to improve parallelism. The preconditioner demonstrates robust performance across a variety of 2D and 3D Stokes problems, often matching or exceeding the effectiveness of an inexact block-triangular (or Uzawa) preconditioner, especially in challenging scenarios such as elongated-domain problems.

While constraints arise naturally in many physical models, their treatment in mathematical and numerical models varies widely, depending on the nature of the constraint and the availability of simulation tools to enforce it. In this paper, we consider the solution of discretized PDE models that have a natural constraint on the positivity (or non-negativity) of the solution. While discretizations of such models often offer analogous positivity properties on their exact solutions, the use of approximate solution algorithms (and the unavoidable effects of floating -- point arithmetic) often destroy any guarantees that the computed approximate solution will satisfy the (discretized form of the) physical constraints, unless the discrete model is solved to much higher precision than discretization error would dictate. Here, we introduce a class of iterative solution algorithms, based on the unigrid variant of multigrid methods, where such positivity constraints can be preserved throughout the approximate solution process. Numerical results for one- and two-dimensional model problems show both the effectiveness of the approach and the trade-off required to ensure positivity of approximate solutions throughout the solution process.

Boundary value problems based on the convection-diffusion equation arise naturally in models of fluid flow across a variety of engineering applications and design feasibility studies. Naturally, their efficient numerical solution has continued to be an interesting and active topic of research for decades. In the context of finite-element discretization of these boundary value problems, the Streamline Upwind Petrov-Galerkin (SUPG) technique yields accurate discretization in the singularly perturbed regime. In this paper, we propose efficient multigrid iterative solution methods for the resulting linear systems. In particular, we show that techniques from standard multigrid for anisotropic problems can be adapted to these discretizations on both tensor-product as well as semi-structured meshes. The resulting methods are demonstrated to be robust preconditioners for several standard flow benchmarks.

Clustering is a commonplace problem in many areas of data science, with applications in biology and bioinformatics, understanding chemical structure, image segmentation, building recommender systems, and many more fields. While there are many different clustering variants (based on given distance or graph structure, probability distributions, or data density), we consider here the problem of clustering nodes in a graph, motivated by the problem of aggregating discrete degrees of freedom in multigrid and domain decomposition methods for solving sparse linear systems. Specifically, we consider the challenge of forming balanced clusters in the graph of a sparse matrix for use in algebraic multigrid, although the algorithm has general applicability. Based on an extension of the Bellman-Ford algorithm, we generalize Lloyd's algorithm for partitioning subsets of Rn to balance the number of nodes in each cluster; this is accompanied by a rebalancing algorithm that reduces the overall energy in the system. The algorithm provides control over the number of clusters and leads to "well centered" partitions of the graph. Theoretical results are provided to establish linear complexity and numerical results in the context of algebraic multigrid highlight the benefits of improved clustering., Comment: 29 pages, 14 figures

Domain decomposition methods (DDMs) are popular solvers for discretized systems of partial differential equations (PDEs), with one-level and multilevel variants. These solvers rely on several algorithmic and mathematical parameters, prescribing overlap, subdomain boundary conditions, and other properties of the DDM. While some work has been done on optimizing these parameters, it has mostly focused on the one-level setting or special cases such as structured-grid discretizations with regular subdomain construction. In this paper, we propose multigrid graph neural networks (MG-GNN), a novel GNN architecture for learning optimized parameters in two-level DDMs\@. We train MG-GNN using a new unsupervised loss function, enabling effective training on small problems that yields robust performance on unstructured grids that are orders of magnitude larger than those in the training set. We show that MG-GNN outperforms popular hierarchical graph network architectures for this optimization and that our proposed loss function is critical to achieving this improved performance.

Algebraic Multigrid (AMG) methods are often robust and effective solvers for solving the large and sparse linear systems that arise from discretized PDEs and other problems, relying on heuristic graph algorithms to achieve their performance. Reduction-based AMG (AMGr) algorithms attempt to formalize these heuristics by providing two-level convergence bounds that depend concretely on properties of the partitioning of the given matrix into its fine- and coarse-grid degrees of freedom. MacLachlan and Saad (SISC 2007) proved that the AMGr method yields provably robust two-level convergence for symmetric and positive-definite matrices that are diagonally dominant, with a convergence factor bounded as a function of a coarsening parameter. However, when applying AMGr algorithms to matrices that are not diagonally dominant, not only do the convergence factor bounds not hold, but measured performance is notably degraded. Here, we present modifications to the classical AMGr algorithm that improve its performance on matrices that are not diagonally dominant, making use of strength of connection, sparse approximate inverse (SPAI) techniques, and interpolation truncation and rescaling, to improve robustness while maintaining control of the algorithmic costs. We present numerical results demonstrating the robustness of this approach for both classical isotropic diffusion problems and for non-diagonally dominant systems coming from anisotropic diffusion.

Sparse matrix representations are ubiquitous in computational science and machine learning, leading to significant reductions in compute time, in comparison to dense representation, for problems that have local connectivity. The adoption of sparse representation in leading ML frameworks such as PyTorch is incomplete, however, with support for both automatic differentiation and GPU acceleration missing. In this work, we present an implementation of a CSR-based sparse matrix wrapper for PyTorch with CUDA acceleration for basic matrix operations, as well as automatic differentiability. We also present several applications of the resulting sparse kernels to optimization problems, demonstrating ease of implementation and performance measurements versus their dense counterparts.

A key consideration in the development of numerical schemes for time-dependent partial differential equations (PDEs) is the ability to preserve certain properties of the continuum solution, such as associated conservation laws or other geometric structures of the solution. There is a long history of the development and analysis of such structure-preserving discretisation schemes, including both proofs that standard schemes have structure-preserving properties and proposals for novel schemes that achieve both high-order accuracy and exact preservation of certain properties of the continuum differential equation. When coupled with implicit time-stepping methods, a major downside to these schemes is that their structure-preserving properties generally rely on exact solution of the (possibly nonlinear) systems of equations defining each time-step in the discrete scheme. For small systems, this is often possible (up to the accuracy of floating-point arithmetic), but it becomes impractical for the large linear systems that arise when considering typical discretisation of space-time PDEs. In this paper, we propose a modification to the standard flexible generalised minimum residual (FGMRES) iteration that enforces selected constraints on approximate numerical solutions. We demonstrate its application to both systems of conservation laws and dissipative systems.

A long-standing issue in the parallel-in-time community is the poor convergence of standard iterative parallel-in-time methods for hyperbolic partial differential equations (PDEs), and for advection-dominated PDEs more broadly. Here, a local Fourier analysis (LFA) convergence theory is derived for the two-level variant of the iterative parallel-in-time method of multigrid reduction-in-time (MGRIT). This closed-form theory allows for new insights into the poor convergence of MGRIT for advection-dominated PDEs when using the standard approach of rediscretizing the fine-grid problem on the coarse grid. Specifically, we show that this poor convergence arises, at least in part, from inadequate coarse-grid correction of certain smooth Fourier modes known as characteristic components, which was previously identified as causing poor convergence of classical spatial multigrid on steady-state advection-dominated PDEs. We apply this convergence theory to show that, for certain semi-Lagrangian discretizations of advection problems, MGRIT convergence using rediscretized coarse-grid operators cannot be robust with respect to CFL number or coarsening factor. A consequence of this analysis is that techniques developed for improving convergence in the spatial multigrid context can be re-purposed in the MGRIT context to develop more robust parallel-in-time solvers. This strategy has been used in recent work to great effect; here, we provide further theoretical evidence supporting the effectiveness of this approach., Comment: Updates v1 with a change of LaTeX template, some minor re-wording, and small bibliographic changes

The fourth-order PDE that models the density variation of smectic A liquid crystals presents unique challenges in its (numerical) analysis beyond more common fourth-order operators, such as the classical biharmonic. While the operator is positive definite, the equation has a "wrong-sign" shift, making it somewhat more akin to an indefinite Helmholtz operator, with lowest-energy modes consisting of plane waves. As a result, for large shifts, the natural continuity, coercivity, and inf-sup constants degrade considerably, impacting standard error estimates. In this paper, we analyze and compare three finite-element formulations for such PDEs, based on $H^2$-conforming elements, the $C^0$ interior penalty method, and a mixed finite-element formulation that explicitly introduces approximations to the gradient of the solution and a Lagrange multiplier. The conforming method is simple but is impractical to apply in three dimensions; the interior-penalty method works well in two and three dimensions but has lower-order convergence and (in preliminary experiments) seems difficult to precondition; the mixed method uses more degrees of freedom, but works well in both two and three dimensions, and is amenable to monolithic multigrid preconditioning. Our analysis reveals different behaviours of the error bounds with the shift parameter and mesh size for the different schemes. Numerical results verify the finite-element convergence for all discretizations, and illustrate the trade-offs between the three schemes.

Domain decomposition methods are widely used and effective in the approximation of solutions to partial differential equations. Yet the optimal construction of these methods requires tedious analysis and is often available only in simplified, structured-grid settings, limiting their use for more complex problems. In this work, we generalize optimized Schwarz domain decomposition methods to unstructured-grid problems, using Graph Convolutional Neural Networks (GCNNs) and unsupervised learning to learn optimal modifications at subdomain interfaces. A key ingredient in our approach is an improved loss function, enabling effective training on relatively small problems, but robust performance on arbitrarily large problems, with computational cost linear in problem size. The performance of the learned linear solvers is compared with both classical and optimized domain decomposition algorithms, for both structured- and unstructured-grid problems.