GridapSolvers.jl: Scalable multiphysics finite element solvers in Julia

Summary and statement of need

The ever-increasing demand for resolution and accuracy in mathematical models of physical processes governed by systems of Partial Differential Equations (PDEs) can only be addressed using fully-parallel advanced numerical discretization methods and scalable solvers, able to exploit the vast amount of computational resources in state-of-the-art supercomputers.

One of the biggest scalability bottlenecks within Finite Element (FE) parallel codes is the solution of linear systems arising from the discretization of PDEs. The implementation of exact factorization-based solvers in parallel environments is an extremely challenging task, and even state-of-the-art libraries such as MUMPS (Amestoy et al., 2001, 2019) and PARDISO (Schenk & Gärtner, 2011) have severe limitations in terms of scalability and memory consumption above a certain number of CPU cores. Hence, the use of iterative methods is crucial to maintain scalability of FE codes. Unfortunately, the convergence of iterative methods is not guaranteed and rapidly deteriorates as the size of the linear system increases. To retain performance, the use of highly scalable preconditioners is mandatory. For simple problems, algebraic solvers and preconditioners (i.e., those based solely on the algebraic system) are enough to obtain robust convergence. Many well-known libraries providing algebraic solvers already exist, such as PETSc (Balay et al., 2021), Trilinos (Trilinos Project Team, 2020), and Hypre (Hypre, n.d.). However, algebraic solvers are not always suited to deal with more challenging problems.

In these cases, solvers that exploit the physics and mathematical discretization of the particular problem are required. This is the case of many multiphysics problems involving differential operators with a large kernel such as the divergence (Arnold et al., 1997) and the curl (Arnold et al., 2000). Examples can be found amongst highly relevant problems such as Navier-Stokes, Maxwell, and Darcy. Scalable solvers for this type of multiphysics problems rely on exploiting the block structure of such systems to find a spectrally equivalent block-preconditioner, and are often tied to a specific discretization of the underlying equations.

As a consequence, high-quality open-source parallel finite element packages like FEniCS (Logg et al., 2012) and deal.II (Arndt et al., 2021) already provide implementations of several state-of-the-art physics-informed solvers (Cui et al., 2024; Farrell et al., 2021). The Gridap ecosystem (Badia & Verdugo, 2020) aims to provide a similar level of functionality within the Julia programming language (Bezanson et al., 2017).

To this end, GridapSolvers is a registered Julia software package which provides highly scalable physics-informed solvers tailored for the FE numerical solution of PDEs on parallel computers within the Gridap ecosystem of packages. Emphasis is put on the modular design of the library, which easily allows new preconditioners to be designed from the user’s specific problem.