Since the early days of computational fluid dynamics there has been an immense increment in computational power. However, the need for development of fast numerical algorithms has persisted consistently. To that motive, the class of domain decomposition (DD) methods has played an important role in reducing the run times and memory requirements for numerical simulations. The general aim for these methods is that the domain of the problem is broken up into a set of many small independent sub domains. New degrees of freedom are introduced, the Lagrange multipliers, that impose the appropriate continuity constraints across these domains. It is then possible to solve for the coupled system of Lagrange multiplier equations only, which is a smaller system compared to the continuous unbroken formulation. The local solution in the sub domains is obtained independently of each other as a function of the solution of the Lagrange multipliers. This process reduces the computational burden of the problem in the sense that it is less demanding in terms of required memory and computational time, without compromising on the accuracy of the results. Over the last decades, many variants of DD methods have been developed, for example, hybrid method1977Raviart, mortar method 2000Wohlmuth, FETI method 1994farhat; 2001farhat, DPG method 2010Demkowicz; 2011Demkowicz; 2016Carstensen, Steklov-Poincaré method 1988Agoshkov; 1991Quarteroni, etc. For a comprehensive discussion of these methods we refer the reader to the review papers 2015Cockburn; 2016Pietro. The two main challenges for DD methods are: i) to find a subset of suitable finite dimensional spaces such that the system matrices are not singular, i.e. they do not produce spurious kinematic modes 2016Almeida, especially with respect to the trace spaces of the Lagrange multipliers; and ii) how to efficiently solve the global system of Lagrange multiplier equations which can become large for practical applications and is characterized by high condition numbers. The primary objective of this paper is to address the first challenge, i.e. to define the framework and extend the use of algebraic dual representations introduced in 2020jain for DD formulation of Darcy flow. The DD formulation used in this work is based on the hybrid method which is a form of discontinuous Galerking formulation. It is shown that the use of algebraic dual representation results in a sparse metric free representation for the divergence constraint on the velocity, the pressure gradient and the continuity constraint across the sub domain, even for high order spectral element method. We also see that the continuity constraints are local to the boundary face of adjoining sub domain elements. This construction is similar to the use of dual spaces in 2000Wohlmuth, where they have non-confirming sub domain boundaries and the continuity constraints are also local to the elements on the adjoining sub domain boundaries only.
In this paper we use spectral elements of order and demonstrate the advantage of using DD formulation with algebraic dual representations using two test cases. The first test case is a manufactured solution taken from 2012Wheeler on a randomly deformed smooth domain. Using this case we first show that using the DD formulation gives the same solution as the continuous formulation and have speed-up in simulation times. It is observed that on mesh refinement the speed-up in simulation time is larger in the case of high order elements (see Table 3). We also show optimal rates of convergence for varying mesh refinements. The second test case is a benchmark reservoir modelling test case SPE10 spe10 that demonstrates the use of this formulation on a practical application. The purpose of this work is fourfold: i) to demonstrate that the algebraic dual representations can be extended to the framework of DD methods without compromising on the accuracy of results, ii) to show that using DD method we can go to higher levels of mesh refinement with same memory hardware, iii) the use of DD formulation is more efficient in terms of simulation time, and iv) to show the applicability of the method for industry benchmark problems such as the SPE10 case.
This paper is structured as follows: In Section 2 we define the broken Sobolev spaces for DD formulation. In Section 3 we define the finite dimensional subset of the broken Sobolev spaces. In Section 4 we state the Darcy problem and the weak formulation of this problem for DD method. The algebraic formulation for this problem and the solution steps are also described in this section. In Section 5 we present the results for the two test cases, i) the manufactured test case from 2012Wheeler, and ii) the benchmark test case SPE10 spe10. We draw conclusions and discuss the scope for future work in Section 6.
2 Broken Sobolev spaces
Let be a bounded domain with Lipschitz boundary . Let be the space of square integrable functions and
the space of square integrable vector fields in 3D, and, be the Hilbert spaces defined as
The trace spaces of is defined as
and we denote by its dual space.
Let be broken into non-overlapping open sub domains with Lipschitz boundary , , such that
Let be the set of sub domains, and be the set of boundaries of these sub domains defined as
We define the broken Sobolev spaces for the set of sub domains as
Let the set of boundary faces, be defined as
Then we define the trace space
We will denote the -inner product by . A set of coordinates is denoted by . We denote by the highest polynomial degree of basis functions used in an element is .
We use two set of finite element spaces, the primal representation and the dual representation from 2020jain. The primal representation of finite element space, basis functions, and degrees of freedom are dentoed by , , respectively. The dual representation of space , its associated basis functions and expansion coefficients are denoted with a tilde as , and , respectively. The basis functions are always represented as row vectors
where is the dimension of polynomial vector space of . The expansion coefficients are always represented as column vectors
We will Lemma 2 from 2020jain, which states
If and are basis functions from primal and dual representations respectively, then these bases are bi-orthogonal with respect to each other
where is the identity matrix of dimension
is the identity matrix of dimension.
The inner product between variables of primal and dual representation is vector dot product of expansion coefficients.
Let , the the inner product is given by
In general, if not explicitly mentioned otherwise, we use Gauss-Lobatto-Legendre points for numerical integration.
3 Finite dimensional spaces
In this section we will define the finite dimensional spaces for hexahedral elements in 3D domains. The domain is discretized into a mesh that consists of points, edges, surfaces and volumes. For each element also consists of a GLL mesh. We will first introduce the finite dimensional spaces for the domain as defined in 2020jain. The basis functions for domain are obtained from mapping of basis functions on a reference domain . For further details on construction we refer the reader to (2020jain, §4.5). Here we directly use the mapped basis functions and the relevant properties of these spaces. We will use the primal representations to define finite element spaces: i) , ii) , iii) , and dual representations to define the finite element spaces: i) , ii) .
These spaces will then be used to define the broken finite dimensional spaces that will be used for DD formulation.
3.1 Primal representations
3.1.1 Finite element space
Let be the flux component of a vector field. For an element , let , form the row vector of basis functions and the column vector of expansion coefficients, given as
Here, denotes the net flux through the face .
Then we can represent as
Using this notation, for two elements the -inner product is given as
where, is the mass matrix associated with the basis functions .
If there is a symmetric positive definite permeability tensor, then the weighted inner product is given as
where is the mass matrix of the weighted inner product.
3.1.2 Finite dimensional space
Let be any element of the space and be the total number of volumes in the discretized domain . If and form the row vector of basis functions and the column vector of expansion coefficients, given as
then we can represent as
where denotes the integral of over mesh volume .
Using this notation, the -inner product for two elements is given as
where, is the mass matrix associated to the basis functions .
The trace space is defined as the restriction of vector fields in to the domain boundary . If is the discrete representation of the inclusion map, that maps degrees of freedom defined on the boundary to the global degrees of freedom of the element, and is the discrete representation of the trace operator thaat restricts global degrees of freedom to the boundary, then the basis functions on the boundary and degrees of freedom on the boundary are given by
If represents the restriction of flux component of vector field to the domain boundary, then the solution on the boundary can be represented by
where tr is the trace operator.
For construction of inclusion matrix, see (2020jain, Ex. 3). The inclusion matrix is a sparse metric-free matrix, that consists of and entries only, and is independent of shape and size of elements as long as the topology, numbering of the degrees of freedom and the orientation of the mesh remains the same.
3.1.4 The divergence operator
For any element , the divergence operation on is defined as (see (2020jain, §4.4)): , such that
where is the discrete representation of the divergence operator that acts on the expansion coefficients of . The divergence operation changes the degrees of freedom and the basis functions to those of space .
Divergence operator for 2D mesh
In Figure 1 on the left plot we show a square domain divided into elements. The expansion coefficients of , are defined across the edges in the mesh. The divergence operation on any element , is then defined as
If we assemble (7) for all the nine elements, with appropriate numbering, we get the discrete divergence operator as
It is a sparse metric free matrix that consists of entries only. If the domain is deformed, for example see the right plot of Figure 1, but the connection between the nodes, edges, surfaces, and volumes, remains the same, relation (7) remains the same and consequently the matrix remains unchanged.
3.2 Algebraic dual representations
In this section we will introduce the finite dimensional sub spaces for and using the algebraic dual representations as defined in 2020jain.
3.2.1 Finite dimensional space
For the finite dimensional space let be the corresponding algebraic dual representation . For any element , if and are the associated basis functions and the expansion coefficients, then we can represent as
Using Corollary 1, the -inner product between the elements is then given by
We see that the -inner product in (4) requires the evaluation of mass matrix, whereas the inner product in (8) requires only the vector product between the expansion coefficients which makes the discrete system more sparse and easier to set up.
3.2.2 Finite dimensional trace space
It is known that the Sobolev spaces and are dual to each other. To replicate this duality also in the discrete setting, we choose as finite dimensional sub space of which is the algebraic dual representation of .
For an element if and are the basis functions and the expansion coefficients then, we can represent as
The inner product between the elements , is given by
3.3 The gradient of dual representations
For a scalar field and its values on the domain boundary, the gradient operation for dual representations is defined as (2020jain, Def. 19), , such that
The expansion coefficients of are given as
expanded in the basis .
3.4 Broken finite dimensional spaces
In this section we define the finite dimensional subset of Broken Sobolev spaces.
3.4.1 Finite dimensional space
For the set of domains we define the finite element space as
For any two elements , and symmetric positive definite permeability tensor , the weighted -inner product is given by
where, and are the column vector of assembled expansion coefficients of all the sub domains and the mass matrix is given as
where are the basis functions associated to the finite dimensional space .
3.4.2 Finite dimensional space
For the set of domains we define the finite element space as
For elements , the - inner product is given by
where , are the column vectors of assembled expansion coefficients of all the sub domains, and the mass matrix is the block diagonal mass matrix given by
where are the basis associated with the finite dimensional space .
3.4.3 Finite dimensional space
For the set of boundaries we define the trace space as
3.4.4 Finite dimensional space
The algebraic dual representation of is denoted by and defined as
Let , then the inner product between the elements is given by
where, are the column vectors of assembled expansion coefficients.
3.4.5 Finite dimensional space
Let be any boundary face of the sub domains. Let be the space of restriction of basis of to and be its algebraic dual representation .
Then the finite dimensional sub space is defined as
For elements , we have the inner product given by
where , are the assembled coefficients over all sub domains.
3.4.6 The divergence operator for broken Sobolev spaces
Let be the elements of the broken finite dimensional spaces. Then the divergence operation on the vector field is given by
where, are the column vectors of the assembled expansion coefficients, and is the assembled divergence operator. If is the discrete representation of divergence operator on the domain , then we have
If the topology of sub domain discretizations is also the same, see Example 1, then we have that . In this paper we have only used the case with the same topology for all the sub domains, therefore
where is the topological divergence operator for any of the sub domains.
3.4.7 The gradient of dual representations
Let be a scalar field represented by algebraic dual representation , and be the boundary value of the scalar field on the sub domain boundary faces . The gradient operation for dual representations is then defined as, , such that
For the pair of elements the -norm is then defined as
4 Model anisotropic diffusion problem
In this section we will use the broken finite dimensional spaces defined in Section 3 to derive the algebraic formulation for DD formulation of Darcy problem.
The equations for Darcy problem in the domain are given by
where is the velocity, is the symmetric positive definite permeability tensor, is the pressure, is the given right hand side term, is the outward unit normal vector, is the given velocity boundary condition imposed on Neumann boundary and is the pressure boundary condition imposed on the Dirichlet boundary .
The Lagrange functional for continuous formulation of Darcy equations is given by
The algebraic system for (33), using algebraic dual representations of 2020jain, is given by
The above system can be solved for unknowns of , and the unknowns of as
For the DD formulation of (32) we break the domain into sub domains, see (1), and we use the Lagrange multipliers to enforce the required continuity across the sub domains. The weak formulation for (32) is then obtained using the Lagrange functional
where, on the left hand side, the first term is the kinetic energy term, the second term imposes the constraint on divergence of velocity field , the third term imposes continuity of flux across the sub domain faces, the fourth term imposes the Neumann boundary condition and the fifth term imposes the Dirichlet boundary conditions.
The Lagrange multipliers are pressure boundary values on sub domains.
In (37) if we take variations with respect to , we get
For any sub domain , we have
If the solution is sufficiently smooth, then using integration by parts on the second term, we get