Skip to main content
Have a personal or library account? Click to login
FEniCSx-pctools: Tools for PETSc Block Linear Algebra Preconditioning in FEniCSx Cover

FEniCSx-pctools: Tools for PETSc Block Linear Algebra Preconditioning in FEniCSx

By: Martin Řehoř and  Jack S. Hale  
Open Access
|Sep 2025

Figures & Tables

Figure 1

Diagram summarising the workflow to implement a custom fieldsplit preconditioner using the FEniCSx-pctools library. (Top left block, UFL) The symbolic UFL bilinear form defining a block operator is first defined as a 2D list of UFL forms a_ufl. (Top right block, DOLFINx) The list of UFL forms is just-in-time compiled to a list of compiled finite element forms a_dolfinx, which is then assembled into a block matrix K in PETSc’s “mpiaij” format (bottom right block, PETSc/petsc4py). We now enter the domain of FEniCSx-pctools, bottom left block. Instead of directly associating the matrix K with the linear solver, we first pass the original UFL form a_ufl and K to create_splittable_matrix_block to create a splittable matrix K_splittable. This is a necessary step to allow block extraction based on arbitrarily combined index sets provided by the fieldsplit preconditioner fs_wrapped. The extraction itself is implemented in the createSubMatrix method of the splittable matrix’ Python-context K_ctx, while the index sets are configurable via setUp method of the preconditioner’s Python-context fs_ctx. (Bottom right block, PETSc/petsc4py) The methods in fs_ctx and K_ctx are transparently executed by the linear solver upon calling its solve method with the right-hand side vector b and a vector x to store the solution.

Figure 2

Standard DOLFINx code for assembling the block linear system Kx = b for the mixed Poisson problem. We define a mesh consisting of 1024 × 1024 quadrilateral cells. For the flux space Qh we choose Brezzi-Douglas-Marini elements of first-order, and for the pressure space Ph discontinuous Lagrange elements of zeroth-order. The right-hand side forcing term f is drawn from a uniform distribution and then the mixed Poisson variational problem is defined using UFL.

Figure 3

Continuation of Figure 2. The fundamental FEniCSx-pctools factory function create_splittable_matrix_block takes the DOLFINx assembled matrix along with the array of bilinear forms that defines it, and returns a PETSc Mat of type “python” with the necessary functionality to apply block preconditioning strategies.

Figure 4

Continuation of Figure 3. We first specify that PETSc should use the custom class fenicsx_pctools.pc.WrappedPC to wrap the fieldsplit preconditioner, and GMRES as an outer solver for Pupper1Kx=Pupper1b. At the next level, we ask for upper Schur complement preconditioning with the structure of the preconditioner specified by the user.

Figure 5

Continuation of Figure 4. We tell PETSc to approximate A˜1 using one application of block Jacobi preconditioning with the matrix A, the already-assembled upper-left block of K.

Figure 6

Continuation of Figure 5. We first setup a class SchurInv with a method apply that will apply the approximate inverse of the discontinuous Galerkin Laplacian matrix S to the vector x and place the result in y. We then tell PETSc to use this method when it needs the action of S˜1.

Figure 7

Continuation of Figure 6. Finally, we set all of the options on the PETSc objects and solve. This solver setup gives a nearly mesh independent number of GMRES iterations (17) tested up to a mesh size of 1024 × 1024.

Figure 8

Visualisation of the solution to the Rayleigh-Bénard problem [18]. The fluid velocity is shown as streamlines coloured by the velocity magnitude, the fluid pressure as an isocontour field, and the fluid temperature on the back surface only. Fluid is held at a higher temperature on the left and a lower temperature on the right of the cube. Less dense fluid rises on the left, and denser fluid sinks on the right (the setup takes the gravitational acceleration vector pointing up). At steady-state this forms a distinctive circulation pattern which is clearly shown in the fluid velocity streamlines.

Table 1

Performance metrics for the Rayleigh-Bénard problem [18] with customised PCD-AMG preconditioning. Weak scaling at 100k DOF per process. Aion Cluster, 50% core utilisation per node. The number in brackets is the average iterations per outer linear solve. Wall time for SNES solve is the time taken to execute SNES.solve, which includes PETSc linear algebra operations, DOLFINx assembly operations and the coupling implemented by FEniCSx-pctools.

DOF (×106)MPI PROCESSESNONLINEAR ITERATIONSLINEAR ITERATIONSNAVIER–STOKES ITERATIONSTEMPERATURE ITERATIONSWALL TIME FOR SNES SOLVE (S)
6.3596428115 (14.4)49 (6.1)26.7
12.612828117 (14.6)49 (6.1)27.2
25.6425629133 (14.8)56 (6.2)31.2
101.7102427103 (14.7)43 (6.1)25.9
203.5204827102 (14.6)44 (6.3)26.1
408.940962582 (16.4)31 (6.2)22.4
816.8819226102 (17)41 (6.8)27.1
DOI: https://doi.org/10.5334/jors.494 | Journal eISSN: 2049-9647
Language: English
Submitted on: Nov 24, 2023
Accepted on: Sep 5, 2025
Published on: Sep 22, 2025
Published by: Ubiquity Press
In partnership with: Paradigm Publishing Services
Publication frequency: 1 issue per year

© 2025 Martin Řehoř, Jack S. Hale, published by Ubiquity Press
This work is licensed under the Creative Commons Attribution 4.0 License.