{'fraction_non_alphanumeric': 0.036293866700966834, 'fraction_numerical': 0.016249705232919586, 'mean_word_length': 4.7405857740585775, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Composable block solvers for the four-field double porosity/permeability model

The objective of this paper is twofold. First, we propose two composable block solver methodologies to solve the discrete systems that arise from finite element discretizations of the double porosity/permeability (DPP) model. The DPP model, which is a four-field mathematical model, describes the flow of a single-phase incompressible fluid in a porous medium with two distinct pore-networks and with a possibility of mass transfer between them. Using the composable solvers feature available in PETSc and the finite element libraries available under the Firedrake Project, we illustrate two different ways by which one can effectively precondition these large systems of equations. Second, we employ the recently developed performance model called the Time-Accuracy-Size (TAS) spectrum to demonstrate that the proposed composable block solvers are scalable in both the parallel and algorithmic sense. Moreover, we utilize this spectrum analysis to compare the performance of three different finite element discretizations (classical mixed formulation with H(div) elements, stabilized continuous Galerkin mixed formulation, and stabilized discontinuous Galerkin mixed formulation) for the DPP model. Our performance spectrum analysis demonstrates that the composable block solvers are fine choices for any of these three finite element discretizations. Sample computer codes are provided to illustrate how one can easily implement the proposed block solver methodologies through PETSc command line options.


M. G. Knepley
Associate Professor, State University of New York at Buffalo One of the main contributions of this paper is performing Time-Accuracy-Size spectrum analysis of various finite element discretizations for the DPP mathematical model. One of the main findings under this analysis is that the continuous stabilized Galerkin (CG-VMS) formulation has a higher Digits of Efficacy (DoE), which measures the combined effect of numerical accuracy and computational time, than the H(div) formulation. The graphical abstract illustrates this finding for the case of four-node tetrahedron element and for the macro-pressure under the DPP model.

INTRODUCTION
Due to recent growth in the exploration of hydrocarbons from unconventional sources (i.e., oil and gas from tight shale), there has been a growing interest to understand and to model flows in porous media with complex pore-networks [Straughan, 2017]. This interest has been amplified due to recent advances in additive manufacturing, which allow for creating materials with complex porenetworks for various applications ranging from water purification filters to composite manufacturing. The porous materials in the aforementioned applications typically exhibit two or more dominant pore-networks with each pore-network displaying distinctive hydro-mechanical properties [Delage et al., 1996;Schmidt and Mcdonald, 1979]. Moreover, due to presence of fissures, there could be mass transfer across the pore-networks [Barenblatt et al., 1960].
To address flows in these types of porous materials, several mathematical models have been proposed in the literature (see [Straughan, 2017] and references therein). A class of models, which is commonly referred to as double porosity/permeability (DPP) models, have been found to be particularly attractive in modeling flows in porous media with two pore-networks (e.g., see [Barenblatt et al., 1960;Boutin and Royer, 2015;Choo et al., 2016;Dykhuizen, 1990;Warren and Root, 1963]). Recently, a DPP mathematical model with strong continuum thermomechanics underpinning has been derived in . This model, which will be central to this paper and will be referred to as the DPP model from here on, describes the flow of a single-phase incompressible fluid in a rigid porous medium with two distinct pore-networks, with possible mass transfer across the pore-networks. The governing equations form a boundary value problem in terms of four-fields and the nature of the PDE is elliptic under steady-state responses.
Except for some academic problems, it is not possible to obtain analytical solutions for the governing equations under the DPP model. Hence, there is a need to resort to numerical solutions. Recently, several numerical formulations have been developed to solve the governing equations under the DPP model; which include [Choo and Borja, 2015;Joshaghani et al., 2018]. However, these works addressed small-scale problems.
The problems that arise in subsurface modeling and other applications involving flow through porous media are typically large-scale in nature. These large-scale problems cannot be solved on a standard desktop or by employing direct solvers; as such a computation will be prohibitively expensive. But large-scale problems from subsurface modeling are routinely tackled using parallel computing tools and by employing iterative linear solvers. The current iterative solver methodologies have been developed and successfully employed for either single-field problems (e.g., Poisson's equation, linear elasticity) or for two-field problems (e.g., Darcy equations, Stokes equations) using two-field composable solvers [Brown et al., 2012;Rathgeber et al., 2016;Saad, 2003]. However, there is a gap in knowledge when one wants to solve large-scale problems under the DPP model, which involves four independent field variables. Unlike Darcy equations, the governing equations under the DPP model cannot be written as a single-field Poisson equation solely in terms of pressures , or even as a two-field problem.
To facilitate solving large-scale problems under the DPP model, we present two four-field composable block solver methodologies. Appealing to PETSc's composable solver features [Balay et al., 2017;Brown et al., 2012] and Firedrake Project's finite element libraries , we will show that the proposed composable block solvers can be effectively implemented in a parallel setting. The two salient features of the proposed block solvers are: they are scalable in both the algorithmic and parallel senses. They can be employed under a wide variety of finite element discretizations. Both these features will be illustrated in this paper using representative two-and three-dimensional problems.
In order to illustrate that the proposed composable solvers can be used under a wide variety of finite element discretizations, we will employ three popular finite element discretizations -the classical mixed formulation (which is based on the Galerkin formalism) using H(div) elements, the CG-VMS stabilized formulation  and the DG-VMS stabilized formulation [Joshaghani et al., 2018]. We will consider H(div) discretizations for simplicial elements (triangle [TRI] and tetrahedron [TET]) and non-simplicial elements (quadrilateral [QUAD] and hexahedron [HEX]). In particular, we employ the lowest-order Raviart-Thomas spaces for simplicial elements [Boffi et al., 2013;Raviart and Thomas, 1977] 1 . For non-simplicial elements, the velocity spaces for QUAD and HEX elements are, respectively, RCTF 1 and NCF 1 [Arnold and Logg, 2014; 2 . The CG-VMS formulation is based on the variational multi-scale (VMS) formalism [Hughes, 1995] and is stable under any arbitrary interpolation order for velocity and pressure fields. The DG-VMS formulation is a discontinuous version of the CG-VMS formulation and is built by combining the VMS formalism and discontinuous Galerkin techniques. The DG-VMS formulation has been shown to accurately capture physical jumps in flow profiles across highly heterogeneous porous media [Joshaghani et al., 2018].
Recently, [Chang et al., 2018a] have proposed the Time-Accuracy-Size (TAS) performance spectrum model, which is an enhanced version of the original spectrum model proposed in [Chang et al., 2018b] obtained by incorporating accuracy into the spectrum model. The TAS spectrum model can be used to study performance of numerical formulations in a parallel setting. Herein, we will utilize the TAS model specifically to achieve the following: (i) We show that the proposed composable solvers are algorithmically scalable. (ii) We compare the performance of the two proposed composable solvers on a particular hardware. (iii) We discuss how the choice of finite element mesh type could affect the solver performance. (iv) We compare the performance of the chosen three finite element discretizations (the classical mixed formulation with H(div) elements, the CG-VMS stabilized formulation and the DG-VMS stabilized formulation) for solving the governing equations under the DPP model.
The work reported in this paper will be valuable to subsurface modelers on three fronts. First and the obvious one is that the proposed composable block solver methodologies facilitate solving large-scale problems involving flow through porous media with multiple pore-networks. Second, our work can guide an application scientist to choose a finite element discretization among several choices. Third, our work illustrates on how to utilize performance metrics other than the commonly used metric -the total time to solution -in subsurface modeling. A couple of these other metrics include Digits of Efficacy (DoE) and the total Degrees-of-Freedom (DoF) processed per second (DoF/s).
The rest of this paper is organized as follows. The governing equations under the DPP model and convenient grouping of the field variables are presented in Section 2. The weak forms of the three finite element formulations that are employed in this paper are presented in Section 3. The proposed two block solver methodologies are discussed in detail in Section 4. The framework of performance spectrum model along with the guidelines on how to interpret the resulting diagrams are presented in Section 5. The performance of the proposed block solvers is illustrated using numerical simulations in Section 6. In the same section, we also compare the performance of the chosen finite element discretizations using the TAS performance spectrum model, which provides guidance to an application scientist with respect to several metrics (e.g., time-to-solution, digitsof-efficacy). Finally, conclusions are drawn in Section 7.

GOVERNING EQUATIONS AND GROUPING OF FIELD VARIABLES
We now document the most important equations under the DPP mathematical model. In this paper, we restrict our treatment of the model to a steady-state response; however, the proposed composable block solvers are equally applicable in a transient setting. We refer the two dominant pore-networks as macro-pore and micro-pore, and the quantities associated with the two porenetworks are, respectively, identified by subscripts 1 and 2. The porous medium is denoted by Ω. Mathematically, Ω ⊂ R nd is assumed to be open and bounded, where nd denotes the number of spatial dimensions. In this paper, nd = 2 or 3. The gradient and divergence operators with respect to a spatial point x ∈ Ω are, respectively, denoted by grad [·] and div [·]. The pressure and the discharge (or Darcy) velocity fields in the macro-pore network are, respectively, denoted by p 1 (x) and u 1 (x), and the corresponding fields in the micro-pore network are denoted by p 2 (x) and u 2 (x). The viscosity and true density of the fluid are denoted by µ and γ, respectively.
The governing equations for a steady-response under the DPP mathematical model take the following form: where k 1 (x) and k 2 (x), respectively, denote the (isotropic) permeabilities of the macro-pore and micro-pore networks, β is a dimensionless characteristic of the porous medium, and b(x) denotes the specific body force. In the above equations, the mass transfer per unit volume, χ(x), from the macro-pore network to the micro-pore network is modeled as follows: As one can see from equations (2.1a)-(2.1d), the flow in one pore-network is coupled with its counterpart in the other pore-network through the inter-pore mass transfer. Unlike Darcy equations, it is not possible to rewrite the governing equations under the DPP model in form of a single-field Poisson's equation in terms of pressures. One has to deal with the governing equations in their mixed form.
2.1. Grouping of field variables in continuum setting. We now discuss two ways of grouping the field variables, which form the basis for the proposed composable block solvers. Under the first approach, the field variables are grouped based on the scale of the pore-network. That is, all the field variables (i.e., velocity and pressure) pertaining to the macro-pore network are placed in one group, and the field variables of the micro-pore network are placed into another. We refer to this splitting of field variables as the scale-split and the associated grouping takes the following form: The governing equations of the DPP model under the scale-split can be compactly written as follows: In the above equation, the differential operator takes the following form: where I denotes the identity tensor, O denotes the zero tensor, and the forcing function takes the following form: Under the second approach, the field variables are grouped based on the nature of the fields. That is, field variables of a similar kind are placed in the same group. We refer to this splitting of field variables as the field-split and the associated grouping takes the following form: The governing equations of the DPP model under the field-split can be compactly written as follows: where the differential operator takes the following form: (2.9) and the forcing function can be written as follows: 2.2. Boundary conditions. The boundary of the domain will be denoted by ∂Ω. The unit outward normal to the boundary at x ∈ ∂Ω is denoted by n(x). The velocity boundary condition (i.e., no penetration boundary condition) for each pore-network takes the following form: where Γ u i denotes that part of the boundary on which the normal component of the velocity is prescribed in the macro-pore (i = 1) and micro-pore (i = 2) networks, and u n1 (x) and u n2 (x) denote the prescribed normal components of the velocities on Γ u 1 and Γ u 2 , respectively. The pressure boundary condition for each pore-network reads: in which Γ p i is that part of the boundary on which the pressure is prescribed in the macro-pore (i = 1) and micro-pore (i = 2) networks, and p 01 (x) and p 02 (x) denote the prescribed pressures on Γ p 1 and Γ p 2 , respectively. For mathematical well-posedness, we assume that

CLASSICAL AND STABILIZED MIXED WEAK FORMULATIONS
As mentioned earlier, we will employ three different mixed formulations -the classical mixed formulation using H(div) discretizations, a continuous stabilized mixed formulation and a discontinuous stabilized mixed formulation. Under a mixed formulation, velocities and pressures are taken to be the primary variables. However, for numerical stability, a mixed formulation should either satisfy or circumvent the Ladyzhenskaya-Babuška-Brezzi (LBB) inf-sup stability condition [Brezzi and Fortin, 2012]. This naturally places all the mixed formulations into either of two categories. A mixed formulation in the first category is built on the classical mixed formulation (which is based on the Galerkin formalism) but places restrictions on the interpolation functions for the independent field variables to satisfy the LBB condition. To put it differently, not all combinations of interpolation functions for the field variables satisfy the LBB condition under the classical mixed formulation. A mixed formulation in the second category augments the classical mixed formulation with stabilization terms so as to circumvent the LBB condition and to render a stable formulation. In this paper, we consider one mixed formulation from the first category and two from the second category.
We denote the set of all square-integrable functions on Ω by L 2 (Ω). The standard L 2 innerproduct over a set K is denoted as follows: For convenience, the subscript K will be dropped if K = Ω in the case of classical mixed formulation and the stabilized mixed continuous Galerkin formulation. In the case of the stabilized mixed discontinuous Galerkin formulation, the subscript K will be dropped if K = Ω, which will be defined later in equation (3.10). We now provide details of the three chosen mixed finite element discretizations. Let us define the following function spaces for the velocities and pressures fields as follows: where H 1 (Ω) is a standard Sobolev space, and H −1/2 (·) is the dual space corresponding to H 1/2 (·). Rigorous discussion of Sobolev spaces are accessible in [Krylov, 2008]; and further discussion of function spaces are provided by [Brezzi and Fortin, 2012].
The classical mixed formulation based on discretizations from the lowest-order Raviart-Thomas spaces is commonly referred to as the RT0 formulation; which is frequently used in subsurface modeling [Chen et al., 2006]. The unknowns under the RT0 formulation on a triangle are fluxes at the midpoints of edges of the element and element-wise constant pressures. The finite dimensional subspaces for each velocity and pressure fields under the lowest-order Raviart-Thomas discretization on a triangle, which are collectively denoted by RTF 1 ⊕ DP 0 , take the following form: where T h is a triangulation on Ω. These subspaces on a tetrahedron, which are denoted by N1F 1 ⊕ DP 0 , take the following form: where T h , in this case, is a tetrahedralization on Ω.
In addition to H(div) discretizations on simplicial meshes, we also consider the corresponding discretizations on non-simplicial element -QUAD and HEX. The velocity spaces for QUAD and HEX elements are, respectively, RCTF 1 and NCF 1 [Arnold and Logg, 2014;. The (macro-and micro-) pressures are element-wise constants, and DG 0 is commonly used to denote element-wise constant discretization on non-simplicial elements. See Figure 1 and Table  1 for a description of these discretizations. The finite dimensional subspaces for the RCTF 1 and NCF 1 discretizations can be written precisely using the language of finite element exterior calculus.
But such a description needs introduction of additional jargon and notation, which is beyond the scope of this paper. We, therefore, refer the reader to [Arnold and Logg, 2014;Arnold et al., 2006Arnold et al., , 2010. However, to guide the reader, the degrees-of-freedom for these discretizations are shown in Figure 1 and Table 1.
3.2. Stabilized mixed continuous Galerkin formulation (CG-VMS). The weak form of the CG-VMS formulation can be written as follows: Find (u 1 (x), u 2 (x)) ∈ U 1 ×U 2 and (p 1 (x), p 2 (x)) ∈ Q such that we have where the bilinear form and the linear functional are defined, respectively, as follows: An attractive feature of the CG-VMS formulation is that nodal-based equal-order interpolation for all the field variables (micro-and macro-velocities and pressures) is stable, which is not the case with the classical mixed formulation. The stability is achieved by the addition of stabilization terms, which circumvent the LBB condition.

Stabilized mixed discontinuous Galerkin formulation (DG-VMS).
Formulations under the discontinuous Galerkin (DG) method inherit attractive features of both finite element and finite volume methods by allowing discontinuous basis functions (e.g., in the form of piecewise polynomials) [Hesthaven and Warburton, 2007]. The DG method supports non-matching grids and hanging nodes, and hence ideal for hp adaptivity [Cockburn et al., 2000]. Moreover, the method can naturally handle jumps in the profiles of the solution variables [Hughes et al., 2006;Joshaghani et al., 2018]. Since the DG method offers several attractive features, we choose the DG formulation that is recently proposed by [Joshaghani et al., 2018] for the DPP model, which will be referred to as the DG-VMS formulation, as one of the three representative formulations to illustrate the performance of the proposed composable block solvers.
We now document the weak form under the DG-VMS formulation. To this end, we decompose the domain into N ele open subdomains such that where ω i denotes the i-th subdomain. The union of all open subdomains will be denoted by We denote the two adjacent subdomains sharing a given interior edge by ω + and ω − . The unit normal vectors on the shared interface Γ ± pointing outwards to ω + and ω − are, respectively, denoted by n + and n − . The jump and average operators on an interior facet for a scalar field ϕ(x) are, respectively, defined as follows: For a vector field τ (x) these operators are defined as follows: τ := τ + · n + + τ − · n − and {τ } : where τ + and τ − are defined similar to equation (3.12). We denote the set of all square-integrable functions on ω i by L 2 (ω i ). We denote the set of all functions that belong to L 2 (ω i ) and are continuously differentiable by H 1 (ω i ). We then introduce the following broken Sobolev spaces (which are piece-wise discontinuous spaces): The weak form under the DG-VMS formulation reads as follows: Find (u 1 (x), u 2 (x)) ∈ U dg × U dg , (p 1 (x), p 2 (x)) ∈ Q dg such that we have where the bilinear form and the linear functional are defined, respectively, as follows: where η u and η p are non-negative, non-dimensional numbers. The DG-VMS formulation also circumvents the LBB condition and is stable under arbitrary combinations of interpolation functions for the field variables.

PROPOSED FOUR-FIELD SOLVERS
The fully discrete formulations for the DPP model can be assembled into the following linear problem: where K is the stiffness matrix, u is the vector of unknown velocities and pressure, and f is the corresponding forcing or RHS vector. Solving the system of equations (4.1) in a fast and scalable way requires careful composition and manipulation of the four different physical fields. In this section, we demonstrate how this can be done through PETSc [Balay et al., 2017[Balay et al., , 2018Dalcin et al., 2011] and its composable solver capabilities [Brown et al., 2012]. The individual block components of the stiffness matrix K for the mixed Galerkin formulation using H(div) elements can be categorized into the following: For the CG-VMS formulation, the individual block components of the stiffness matrix can be categorized into the following: Likewise, the block components of the stiffness matrix for the DG-VMS formulation read: The components of the corresponding RHS vector f for equations (4.2), (4.3) and (4.4) are respectively. Specifically, we employ PETSc's block solver capabilities, in the PCFIEDLSPLIT class, taking two fields at a time. However, the global DPP model is a four field problem so we subdivide our problem recursively such that we end up with 2×2 blocks. Conceptually, PETSc can employ a wide variety of block solver methodologies on a 2×2 matrix: where A, B, C, and D are individual block matrices which also consist of 2×2 blocks. Although equation (4.8) is conceptually a 4×4 block matrix, PETSc's field-splitting capabilities enables us to break the system down dynamically at runtime into two levels of 2×2 blocks. We now propose two different ways one can compose scalable and efficient solvers and preconditioners for blocks A, B, C, and D with the individual components shown in equations (4.2), (4.3), and (4.4).
4.1. Method 1: splitting by scales. One option is to split the global problem by scales. That is, each macro-or micro-scale 2×2 block will contain its corresponding velocity and pressure fields. Under this solver strategy, equation (4.1) is then rewritten as: where 0 is a zero matrix, u 1 and p 1 are the respective macro-scale velocity and pressure vectors, v 2 and p 2 are the respective micro-scale velocity and pressure vectors. The individual 2×2 blocks from equation (4.8) would be Although the off diagonal blocks B and C contain the inter-scale pressure coupling terms, they are very sparse so we will ignore these blocks for now. The composition of the A and D blocks are similar to the classical mixed Poisson problem so the Schur complement approach outlined in Mapakshi et al., 2018] and the references within can be applied.
The task is to individually precondition the decoupled A and D blocks. We note that they admit factorizations of where I is the identity matrix and are the Schur complements for the A and D blocks, respectively. The inverses can therefore be written as The task at hand is to approximate the inverses of the K 1 vv , K 2 uu , S 1 , and S 2 blocks. The first two blocks are simply mass matrices so we can invert them using the ILU(0) (incomplete lower upper) solver. For the Schur complement blocks, we employ a diagonal mass-lumping of K 1 uu and K 2 uu to estimate K 1 uu −1 and K 2 uu −1 because they are spectrally equivalent to the identity. That is, to precondition the inner solvers responsible for inverting S 1 and S 2 . For these blocks we employ the multigrid V-cycle on S 1 p and S 2 p from the HYPRE BoomerAMG package [Falgout and Yang, 2002]. We expect these to work because the S blocks are spectrally equivalent to the Laplacian, modulo the penalty terms. In [Mapakshi et al., 2018] it turns out the presence of the VMS stabilization terms in the K 1 pp and K 2 pp blocks do not drastically affect the performance or scalability of this solver strategy.
Instead of completely solving for the K −1 uu and S p of both scales, we apply only a single sweep of ILU(0)/block Jacobi and V-cycle, respectively, and rely on GMRES [Saad and Schultz, 1986] to solve the entire 4×4 block system. Thus this outer GMRES is able to pick up the inter-scale pressure coupling blocks B and C. The PETSc command-line options for this solver methodology is given in listing 1.

Method 2: splitting by fields.
Another option is to group the velocities and pressures of both scales into two different blocks. If this approach is taken, equation (4.1) is then rewritten as: (4.19) and the individual blocks in equation (4.8) would now look like Unlike the previous methodology, we can work directly with the above stiffness matrix, which admits a factorization of where the Schur complement S is The inverse can therefore be written as Although A is a 2×2 block containing velocities spanning across two different scales, we can still approximate A −1 by inverting the entire A block using ILU(0) because the off-diagonal blocks are zero and the diagonal blocks consist of only mass matrices. Approximating S −1 is a little trickier because equation (4.22) is a dense 2×2 block with off-diagonal terms. However, we can still employ a diagonal mass-lumping of A to estimate A −1 because it is again spectrally equivalent to the identity. The preconditioner needed for S −1 is: The off-diagonal blocks only consist of mass-matrix terms but the decoupled diagonal blocks are identical to equations (4.17) and (4.18). Thus, we individually employ multigrid V-cycle on each of the diagonal blocks. As in the previous solver methodology, only a single sweep of ILU(0) and the two multigrid V-cycles are needed for the A −1 matrix and the two diagonal terms within the S p matrix, respectively, and the GMRES method is employed to solve the entire block system. The PETSc implementation is shown in listing 2.

Computer implementation.
The finite element capabilities are provided by the Firedrake Project package Homolya and Ham, 2016;Homolya et al., 2017Homolya et al., , 2018Luporini et al., 2015Luporini et al., , 2017Rathgeber et al., 2016] with GNU compilers. This sophisticate finite element simulation package and its software dependencies can be found at [Zenodo/COFFEE, 2017;Zenodo/FIAT, 2018;Zenodo/FInAT, 2018;ZENODO/firedrake, 2018;Zenodo/PETSc, 2018;Zenodo/petsc4py, 2018;Zenodo/PyOP2, 2018;Zenodo/TSFC, 2018;Zenodo/UFL, 2018]. The computational meshes are built on top of the DMPlex unstructured grid format [Knepley and Karpeev, 2009;Lange et al., 2015 and partitioned through the Chaco package [Hendrickson and Leland, 1995]. Pictorial descriptions of the specific elements represented by this mesh format and utilized in this paper are illustrated in Figure 1. The DMPlex data structure interfaces very nicely with PETSc's suite of parallel solvers and provide excellent scalability across thousands of MPI processes [Chang et al., , 2018b. Sample Firedrake codes for some of these benchmark problems can be found in Appendix A.
In our PETSc implementation, the same global matrix will be assembled for both solvers. The preconditioners differ by the subblocks which are extracted. The different sparsity pattern of the subblocks contributes to the performance differences seen in the solvers, but the overall assembly time remains unchanged for either solvers.

PERFORMANCE SPECTRUM MODELING
To understand the parallel performance and algorithmic scalability of the proposed DPP composable block solver methodologies for the three finite element formulations, a performance model is needed. The performance model based on the Time-Accuracy-Size (TAS) spectrum analysis outlined in [Chang et al., 2018a] shall be used as the basis for understanding the quality of these finite element formulations with the proposed block solvers. We now briefly highlight the performance metrics used in this section and why they are each important in each of their own ways.

Mesh convergence.
This criterion uses the convergence notion to account for numerical accuracy of a solution in the performance spectrum. In this paper, we are adopting L 2 norm of the error defined as: where u is the exact solution, u h is the the finite element solution, and h is measure of element size. Based on theory, most finite element discretizations will have an upper-bound for L 2 error norm as follows: where α is known as convergence rate and C is some constant. When reporting and comparing how much accuracy is attained for each discretization, we use the notation of Digits of Accuracy (DoA) defined as: and plot DoA against Digits of Size (DoS), which is defined as: Noting that for most formulations DoF = Dh −d , where d is the spatial dimension and D is some constant, the slope of DoA vs DoS plot is in the order of α d . Any tailing off from the line plot is an indicator of incorrect implementation or solver convergence tolerances being too relaxed. Furthermore, the ratio DoA/DoS can be a good indicator of how much accuracy is achieved per DoF.

Strong-scaling.
In this basic parallel scaling while the size of problem remains unchanged, the number of processes increases. In general, this metric comments on the marginal efficiency of each additional processes assigned to a problem. It is conventional to plot number of processes against the parallel efficiency defined as: where proc is the number of MPI processes, T 1 is the total wall-clock time needed on a single MPI processes, and T p is the total wall-clock time needed with proc MPI processes. However, this metric must be interpreted carefully for the following reasons: (1) Solver iteration counts: The number of solver iterations may fluctuate as the number of MPI processes changes. This can happen for a number of reasons, whether it is algorithmic implementation or relaxed convergence criterion. It is necessary to also report the number of KSP iterations required as the number of processes changes. (2) Problems too small: If the DoF count is too small for a particular MPI concurrency, communication time will swamp the computation time, thus reducing the parallel efficiency. This issue may arise when making comparative studies between different finite element discretizations, as different formulations have different DoF counts for a given h-size. Furthermore, for Python-based simulation packages like Firedrake, overheads from just-in-time compilation and instantiation of objects can also affect the strong-scaling . (3) Problems too large: If the DoF count is too large for a particular MPI concurrency, the problems not only drop out of the various levels of cache in the memory hierarchy but also invoke several expensive cache misses which can slow down the overall performance. This may result in superlinear speedups, like the BLMVM bound-constrained optimization solver in .
Lastly, the global problem size for this scaling analysis is fixed, so we need an additional parallel scaling metric which explains whether the performance of our solvers might degrade due to increased KSP iteration counts or memory contention as the problem size increases.
5.3. Static-scaling. As described in [Chang et al., 2018b], static-scaling is a scaling analysis where the MPI concurrency is fixed but the problem size is increased. The essential metric for this analysis is the computation rate (DoF over Time). In this paper, we run a series of problem sizes at a fixed parallelism and plot the computation rate against the wall-clock time. Note that the time need not be the total time to solution, instead one could look at various phases like the finite element assembly or solver computation rates.
Static-scaling returns information on performance and scalability of software and solvers across different hardware architectures. This scaling analysis also captures both strong-scaling and weakscaling effects. Assuming that the block solvers are of O(N ) scalability, where N = DoF , optimal scaling is indicated by a horizontal curve. Any tail offs at small problem sizes suggests strongscaling effects whereas tail offs at large problem sizes indicate suboptimal algorithmic or memory effects. The exact reasoning for the tail offs towards the right can be verified through arithmetic intensity, which is the measure of the total work over the total bytes transferred (see [Chang et al., 2018b] and the references within).

Digits-of-Efficacy (DoE).
The final metric needed for our performance spectrum study is the Digits of Efficacy (DoE). This metric measures the accuracy production by a particular scheme in a given amount of time. The DoE could be defined as: Assuming that straight lines are captured in both the mesh convergence and static-scaling diagrams, the DoE has a linear dependence on problem size and returns a slope of d − α (see [Chang et al., 2018a] for details on the exact derivation). This efficacy measure is analogous to the action of a mechanical system, that is the product of energy and time. In the TAS spectrum analysis, the DoE represents an analogous action for computation, and we speculate that an optimal algorithm minimizes this product over its runtime. Since the DoE takes the negative logarithm of action, a higher DoE is desirable.

REPRESENTATIVE NUMERICAL RESULTS
In this section, after clarifying the terminology and framework adopted for the performance spectrum model, we solve the four-field DPP model in two-and three-dimensional settings in order to demonstrate the implementation of the proposed composable block solvers and gauge their performances. The two-dimensional problem will be conducted in serial (one MPI process) on a dual socket Intel Xeon E5-2609v3 server node. The three-dimensional problems will be conducted on a dual socket Intel Xeon E5-2698v3 server node and will utilize up to 16 MPI processes (8 MPI processes per socket). On different performance metrics, H(div), CG-VMS, and DG VMS formulations are compared for both simplicial (TRI, TET) and non-simplicial (QUAD, HEX) meshes. Both two-dimensional and three-dimensional problems were adopted by  for the convergence analysis of continuous stabilized mixed formulation (CG-VMS) and by [Joshaghani et al., 2018] for the convergence analysis of discontinuous stabilized mixed formulation (DG-VMS) for the DPP model. We are generating three series of outputs for first-order CG-VMS, DG-VMS with η p = η u = 10, and H(div) formulations.
6.1. Two-dimensional study. For this first problem, let us consider a two-dimensional DPP boundary value problem with governing equations stated in equations (2.1)-(2.2) and (2.11)-(2.12). The homogeneous (i.e., constant macro and micro-permeabilities) bi-unit square computational   640  20  2646  20  2646  20  14400  20  9600  20  4080  20  2480  40  10086  40  10086  40  57600  40  38400  40  16160  40  9760  80  39366  80  39366  80  230400  80  153600  80  64320  80  38720  160  155526  160  155526  160  921600  160  614400  160  256640  160  154240 domain and boundary conditions for this study are shown in Figure 2(a), and the corresponding parameters are described in Table 2. The analytical solution for the pressure and velocity fields takes the following form: u 1 (x, y) = −k 1 e πx sin(πy) e πx cos(πy) − η βk 1 e ηy (6.1a) p 1 (x, y) = µ π e πx sin(πy) − µ βk 1 e ηy (6.1b) u 2 (x, y) = −k 2 e πx sin(πy) e πx cos(πy) + η βk 2 e ηy (6.1c) p 2 (x, y) = µ π e πx sin(πy) + µ βk 2 e ηy (6.1d) where η is defined as: For two-dimensional performance spectrum analysis, all three finite element formulations will start off with the same h-sizes and will be refined up to 6 times. The initial TRI and QUAD coarse meshes are shown in Figures 2(b) and 2(c) and the corresponding DoF counts for each formulation is shown in Table 3. The mesh convergence results with respect to DoA and DoS are performed under field-splitting solver and are shown in Figures 3 and 4 respectively for TRI and QUAD meshes. It should be noted that by applying scale-splitting solver, very same results could be obtained and for brevity, we decided not to plot them in figures. It can be seen in these diagrams that the CG-VMS and DG-VMS lines exhibit a slope α d = 1, which verifies that our Firedrake implementation of these discretizations is correct. The H(div) lines exhibit a slope of 0.5 for TRI meshes but appear to have superlinear convergence for the QUAD meshes, which has also been observed in other Firedrake endeavors [Gibson et al., 2018]. It can also be seen that if the solver tolerances are not strict enough, the mesh convergence lines will tail off. Nonetheless, the CG-VMS and DG-VMS have the highest ratios of DoA over DoS in most of these diagrams which suggests that each DoF in VMS formulations has a greater level of contribution to the overall numerical accuracy than their H(div) counterparts.
Static-scaling results for both block solver strategies are shown in Figure 5, and we see that the total wall clock time is almost equally distributed among the assemble and solve phases. The field-splitting methodologies are slightly worse than their scale-splitting counterparts for the VMS formulations. However, the difference in performance is almost negligible when we look at the total time. The DoF counts are too small as the line curves for both the assembly and solve phases flatten out when all three formulations have roughly 10K DoF or more. No matter which mesh is utilized, the H(div) formulation processes its DoF count faster than either VMS formulations. Figures 6 and 7 contain DoE diagrams for TRI and QUAD meshes, respectively. Although H(div) appears to have the highest computation rates, it has a lower DoA than its VMS counterparts which results in a much smaller DoE. The QUAD mesh on the other hand has a very high DoA and it beats out its VMS counterparts for all the fields.
First, we consider when all discretizations have an h-size = 1/16. Strong-scaling results for both field-splitting and scale-splitting block solver methodologies for H(div), CG-VMS, and DG-VMS can be found in Tables 6, 7, and 8, respectively. All three Tables indicate that the KSP iteration counts between field-splitting and scale-splitting are identical whereas the wall-clock time for scale-splitting is slightly smaller. The KSP counts for H(div) and CG-VMS do not change much when the number of MPI processes increases, whereas DG-VMS's KSP counts increase drastically. This increase in KSP iteration counts will affect the parallel efficiency so one has to be careful when interpreting these results. Nonetheless, we see that the DG-VMS parallel efficiency is the highest, even with its proliferated KSP counts. This is attributed to the fact that the DoF count for DG-VMS is larger than CG-VMS and H(div). All three tables indicate that higher DoF counts bring in more efficiency in the parallel sense.
Second, we consider the case when all discretizations contain approximately 200K degrees-offreedom. Strong-scaling results for both field-splitting and scale-splitting block solver methodologies for H(div), CG-VMS, and DG-VMS can be found in Tables 9, 10, and 11, respectively. Like with the same h-size case, the scale-splitting method appears to be more efficient in terms of wall-clock time needed despite having the same KSP counts as the field-splitting method. It can also be seen that the H(div) and CG-VMS KSP counts do not fluctuate much with MPI processes and that DG-VMS KSP counts still increase dramatically. However, tuning the mesh sizes such that all finite element discretizations have the same DoF count enables us to have better understanding of the parallel performance, especially for three-dimensional problems. It can be seen that H(div) requires the least amount of wall-clock time resulting in the lowest parallel efficiency, but that does not mean this is a bad formulation. In order to understand the quality of the H(div) discretizations, we need to take into consideration the numerical accuracy and perform a TAS spectrum analysis.
6.2.2. Test 2: TAS Spectrum Analysis. For the TAS spectrum analysis, we consider a range of problems, shown in Table 12, such that all finite element formulations in each refinement step have roughly the same DoF count. The mesh convergence results with respect to DoA and DoS, for both  TET and HEX meshes, are shown in Figure 9. CG-VMS and DG-VMS lines indicate a slope of 2 3 , which again corroborates that our Firedrake implementation of these formulations are correct. The H(div) lines exhibit a slope of 1 3 for TET mesh. However, similar to the two-dimensional problem for non-simplicial element QUAD, H(div) exhibits super linear convergence for the HEX meshes. We are not observing any tail-offs in these results as the solver relative convergence tolerance of 1e − 7 was strict enough. The observation that both CG-and DG-VMS have the highest DoA over DoS ratio for almost all velocity and pressure fields implies that they have greater levels of contribution to the overall numerical accuracy than the H(div) schemes. Static-scaling results for both block solver strategies are presented in Figure 10. Flat lines appear in all six subfigures, indicating that the proposed block-solver methodologies are scalable under the chosen h-sizes and hardware environment. It is a common belief among application scientists that a solver exhibits worse scaling than an assembly procedure, since assembly is almost entirely local. However, the results show that for all the chosen discretizations-no matter what solver methodology is employed-time to assemble stiffness matrix is higher than the solver time. This infers that we have successfully optimized solvers to such an extent that the assembly procedure is more dominant.
Analogous to the two-dimensional problem, the scale-splitting methodologies are slightly better than their field-splitting counterparts for the all formulations. Evidently, this disparity is more clear for VMS formulations at the solve time level. However, the difference in performance is almost inconsequential when we look at the total time. It can be seen that the DoF counts are sufficient to level out the curves when all three formulations have roughly 20K DoFs or more. Regardless of the mesh type, the H(div) formulation processes its DoF count faster than either VMS formulations. Figures 11 and 12 contain DoE diagrams for TET and HEX meshes, respectively. For the case of TET mesh type, in spite of H(div) having the fastest computation rates, it has a lower DoA than its VMS counterparts which in turns lead to a much smaller DoE with steep declining curve. On the contrary, for the case of the HEX mesh, H(div) surpasses its VMS counterparts due to its high DoA values. These diagrams demonstrate how numerical accuracy can have a drastic effect on the overall computational performance of these various finite element formulations.

CLOSURE
We have developed two block solver methodologies which are capable of solving large-scale problems under the four-field DPP mathematical model. We have also presented a systematic performance analysis of various finite element discretizations for the DPP model using the recently  proposed Time-Accuracy-Size (TAS) spectrum model, which takes into consideration important metrics such as mesh convergence, static-scaling, and Digits of Efficacy (DoE). We have also identified strong-scaling issues one needs to be cognizant of when the block solvers are applied to various finite elements. In our numerical studies, two-and three-dimensional problems had analogous performance trends, despite their marked discrepancy in time to solution. Some of the salient features of the proposed composable block solver methodologies are as follows:  (1) Both composable solvers are compatible with different kinds of mixed finite element formulations: H(div) and non-H(div) elements, simplicial and non-simplicial elements, nodeand edge-based discretizations, and continuous and discontinuous approximations. (2) Both composable solvers are scalable in both parallel and algorithmic senses.
(3) The solvers can be implemented seamlessly using the existing PETSc's composable solver options. Hence, one can leverage on the existing parallel computing tools to implement these composable solvers into existing simulators.
Some of the main conclusions from the performance analysis based on the TAS spectrum model are as follows: (1) Scale-split vs. field-split. For a fixed problem size, the scale-splitting methodology tends to be slightly more efficient in terms of wall-clock time needed despite having the same KSP counts as the field-splitting method. However, selecting either solver methodologies will be left to the programmer's convenience and limitations as switching from one strategy to another exerts negligible overall effects on performance metrics.
(2) H(div) vs. VMS formulations. (a) No matter what mesh type is chosen, DoFs are processed the fastest under the H(div) formulation compared to the CG-VMS or DG-VMS formulations. (b) The VMS formulations yield much higher overall numerical accuracy for all velocity and pressure fields than their H(div) counterparts. The exception is for nonsimplicial meshes, where the H(div) formulation exhibits super linear convergence.