2.2.2. PETSc-Based Solver

The solution of the linear systems arising from the FEM discretization of the coupled velocity-pressure system is carried out taking advantage of PETSc, the Portable, Extensible Toolkit for Scientific Computation [63]. This step is key to the overall efficiency and viability of the method as an appropriate candidate to tackle Newtonian and non-Newtonian multiphase flows in which large density and viscosity ratios are involved, as those taking place in experiments dealing with drop and bubble deformation. As the ratio increases, so does the condition number of the resulting velocity matrix, which hampers the convergence of traditional iterative approaches to the solution of the coupled velocity-pressure system such as preconditioned conjugate gradient techniques [64] or traditional "splitting" techniques, usually problematic at low (viscosity-dominant) Reynolds numbers [65]. Instead, we use a physics-based, block-preconditioning approach presented by Elman and collaborators [66] and implemented in the PETSc-FieldSplit preconditioner: we define separate fields for the discrete velocity and pressure and combine them using a Schur-complement factorization [67] of the total system matrix, which is efficiently stored as a MatNest structure K=[A B; B<sup>T</sup> 0], in which the corresponding sub-matrices for the discrete velocity, discrete gradient, and discrete divergence are the four entries K\_ij, for i,j=1,2. The following snippet collects a sample function to build the block matrix.

```
PetscErrorCode fun_petsc_construct_block_matrix_PETSC_Stokes_system (...)
{
  ...
  /* Build first sub-matrix "K00" */
  ierr = MatDuplicate(K00,MAT_COPY_VALUES,& (Stokes_Problem_system->K_ij[0]));CHKERRQ(ierr);ierr = MatAssemblyBegin(Stokes_Problem_system->K_ij[0],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);ierr = MatAssemblyEnd(Stokes_Problem_system->K_ij[0],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ...
  /* Build complete, block matrix "K" from sub-matrices "K_ij" as a nested matrix */
  ierr = MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, \
               Stokes_Problem_system->K_ij, & (Stokes_Problem_system->K)); CHKERRQ(ierr);
  ierr = MatAssemblyBegin(Stokes_Problem_system->K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Stokes_Problem_system->K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ...
}
```
The previous approach can now be used to tackle the saddle-point problem resulting from our FEM-based discretization of the Navier–Stokes equations using the PETSc-FieldSplit preconditioner for the system matrix, collecting all the information in a suitable structure such as the PETSc\_Saddle\_point\_system sample structure included in the following code snippet.

```
/* Structure with all the information required to solve a saddle-point problem */
typedef struct
{
 Mat K; /*MatNest matrix with 4 sub-matrices: K = [K00 K01; K10 K11]*/
 Mat K_ij[4]; /*each of the four sub-matrices of block matrix "K"*/
 Mat Pmat; /*MatNest precond-matrix: Pmat = [Pmat00 Pmat01; Pmat10 Pmat11]; */
 Mat Pmat_ij[4]; /*each of the four sub-matrices of precond-matrix "Pmat"*/
 Vec x; /*solution of system*/
 Vec b; /*right hand side vector*/
 KSP ksp; /*solver context*/
 MatNullSpace nullsp; /*Nullspace for pressure (in Stokes solved by PCFieldSplit)*/
 PC pc; /*preconditioner context*/
 IS isg[2]; /* index sets (rows) of splits "0" and "1" (e.g. velocity, pressure)*/
}
PETSC_Saddle_point_system;
```
The CHOLMOD package [68] is used to perform the sparse direct Cholesky factorization of the positive-definite block matrix K00=A, while the Least-Squares Commutator (LSC) preconditioner [67] takes care of the second field (pressure) involved in the upper factorization of the Schur complement, using a combination of the conjugate gradient and additive-Schwarz methods [69] for the multigrid levels of the ML (Multi-Level) preconditioner package [70] used for LSC. The following options passed to PETSc were found to offer excellent performance in demanding multiphase flow problems with high density ratios such as those found in Section 3.2.

```
/* KSP00=preonly (cholesky,CHOLMOD); KSP11:preonly (lsc); ksp_lsc: cg (ml(+asm)). */
char options_pc_stokes[] = "-stokes_ksp_rtol 5.e-9 -stokes_ksp_diagonal_scale \
```
These are some of the general configuration options and data structures that are most responsible for the performance boost we find for the kind of complex, polymer flows studied in this work. Given the platform-independent model of C and PETSc, the snippets could be included with little modification in any existing C code used by any researcher with a similar underlying mathematical problem, independently of the particular field of knowledge.
