Use of algebraic dual representations in domain decomposition methods for Darcy flow in 3D domain

In this work we use algebraic dual representations in conjunction with domain decomposition methods for Darcy equations. We define the broken Sobolev spaces and their finite dimensional counterparts. In addition, a global trace space is defined that connects the solution between the broken spaces. Use of dual representations results in a sparse metric free representation of the constraint on divergence of velocity, the pressure gradient term and on the continuity constraint across the sub domains. To demonstrate this, we solve two test cases: i) manufactured solution case, and ii) industrial benchmark reservoir modelling problem SPE10. The results demonstrate that the domain decomposition scheme, although with more unknowns, requires less memory and simulation time as compared to the continuous Galerkin formulation.


Introduction
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 method [1], mortar method [2], FETI method [3,4], DPG method [5][6][7], Steklov-Poincaré method [8,9], etc. For a comprehensive discussion of these methods we refer the reader to the review papers [10,11]. 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 [12], 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 [13] 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 [2], 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 N = 1, 2, 3 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 [14] 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 [15] 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 [14], and ii) the benchmark test case SPE10 [15]. We draw conclusions and discuss the scope for future work in Section 6.

Broken Sobolev spaces
Let Ω ⊂ R 3 be a bounded domain with Lipschitz boundary ∂Ω. Let L 2 (Ω) be the space of square integrable functions and [L 2 (Ω)] 3 the space of square integrable vector fields in 3D, and H 1 (Ω), H (div; Ω) be the Hilbert spaces defined as .
The trace spaces of H 1 (Ω) is defined as and we denote by H −1/2 (∂Ω) its dual space.
Let Ω be broken into T non-overlapping open sub domains M i with Lipschitz boundary ∂M i , i = 1 , . . . , T , such that Let Ω T be the set of sub domains, and ∂Ω T be the set of boundaries of these sub domains defined as We define the broken Sobolev spaces for the set of sub domains Ω T as .
Let the set of boundary faces, ∂Ω e be defined as Then we define the trace space

Notation
We will denote the L 2 -inner product by (· , ·). A set of coordinates is denoted by x. We denote by N the highest polynomial degree of basis functions used in an element is N − 1.
We use two set of finite element spaces, the primal representation and the dual representation from [13]. The primal representation of finite element space, basis functions, and degrees of freedom are dentoed by X (·), Ψ (·), N respectively. The dual representation of space X (·), its associated basis functions and expansion coefficients are denoted with a tilde as X (·), Ψ (·) and N (·), respectively. The basis functions are always represented as row vectors where d k is the dimension of polynomial vector space of X (·). The expansion coefficients are always represented as column vectors We will Lemma 2 from [13], which states Lemma 1. If Ψ k (x) and Ψ k (x) are basis functions from primal and dual representations respectively, then these bases are bi-orthogonal with respect to each other where I is the identity matrix of dimension d k .
Corollary 1. The inner product between variables of primal and dual representation is vector dot product of expansion coefficients.
Proof. Let p, q ∈ X (Ω) × X (Ω), the the inner product is given by In general, if not explicitly mentioned otherwise, we use Gauss-Lobatto-Legendre points for numerical integration.

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 N > 1 each element also consists of a GLL mesh. We will first introduce the finite dimensional spaces for the domain Ω as defined in [13]. The basis functions for domain Ω are obtained from mapping of basis functions on a reference domain Ω = [−1, 1] 3 . For further details on construction we refer the reader to [13, §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) D (Ω) ⊂ H (div; Ω), ii) S (Ω) ⊂ L 2 (Ω), iii) D b (∂Ω) ⊂ H −1/2 (∂Ω), and dual representations to define the finite element spaces: i) S (Ω) ⊂ L 2 (Ω), ii) D b (∂Ω) ⊂ H 1/2 (∂Ω).
These spaces will then be used to define the broken finite dimensional spaces that will be used for DD formulation.
2 (x) . . . (2) Here, N 2 i (u) = f i u · n dA denotes the net flux u through the face f i . Then we can represent u as u (x) = Ψ 2 (x) N 2 (u) .
Using this notation, for two elements u, v ∈ D (Ω) the L 2 -inner product is given as where, M (2) := Ω Ψ 2 (x) Ψ 2 (x) dΩ is the mass matrix associated with the basis functions Ψ 2 (x). If there is a symmetric positive definite permeability tensor K, then the weighted inner product is given as dΩ is the mass matrix of the weighted inner product.
Let p be any element of the space S (Ω) and N v be the total number of volumes in the discretized domain Ω. If Ψ 3 (x) and N 3 (p) form the row vector of basis functions and the column vector of expansion coefficients, given as 2 (x) . . .
then we can represent p as where N 3 i (p) = V i p dV denotes the integral of p over mesh volume V i . Using this notation, the L 2 -inner product for two elements p, q ∈ S (Ω) is given as where, M (3) := Ω Ψ 3 (x) Ψ 3 (x) dΩ is the mass matrix associated to the basis functions Ψ 3 (x).
is defined as the restriction of vector fields in D (Ω) to the domain boundary ∂Ω. If N 2 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 N 2 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 tru ∈ D b (∂Ω) 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 [13,Ex. 3]. The inclusion matrix is a sparse metric-free matrix, that consists of +1, −1 and 0 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.

The divergence operator E 3,2
For any element u ∈ D (Ω), the divergence operation on u is defined as (see [13, §4.4] where E 3,2 is the discrete representation of the divergence operator that acts on the expansion coefficients of u. The divergence operation changes the degrees of freedom and the basis functions to those of space S (Ω).

Example 1. Divergence operator for 2D mesh
In Figure 1 on the left plot we show a square domain divided into 3 × 3 elements. The expansion coefficients of u ∈ D (Ω), are defined across the edges in the mesh. The divergence operation on any element K, 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 +1, −1, 0 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 E 2,1 remains unchanged.

Algebraic dual representations
In this section we will introduce the finite dimensional sub spaces for L 2 (Ω) and H 1/2 (∂Ω) using the algebraic dual representations as defined in [13].
For the finite dimensional space S (Ω) let S (Ω) be the corresponding algebraic dual representation . For any element q ∈ S (Ω), if Ψ 0 (x) and N 0 ( q) are the associated basis functions and the expansion coefficients, then we can represent q as Using Corollary 1, the L 2 -inner product between the elements p, q ∈ S (Ω) × S (Ω) is then given by We see that the L 2 -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.
It is known that the Sobolev spaces H 1/2 (∂Ω) and H −1/2 (∂Ω) are dual to each other. To replicate this duality also in the discrete setting, we choose D b (∂Ω) as finite dimensional sub space of H 1/2 (∂Ω) which is the algebraic dual and B 0 (λ) are the basis functions and the expansion coefficients then, we can represent λ as where The inner product between the elements λ , tr u ∈ D b (∂Ω) × D b (∂Ω), is given by

The gradient of dual representations
For a scalar field p ∈ S (Ω) , andp ∈ D b (∂Ω) its values on the domain boundary, the gradient operation for dual representations is defined as [13,Def. 19 The expansion coefficients of grad (p,p) are given as

Broken finite dimensional spaces
In this section we define the finite dimensional subset of Broken Sobolev spaces.
For any two elements u , v ∈ D (Ω T ), and symmetric positive definite permeability tensor K, the weighted L 2 -inner product is given by where, N 2 (u) and N 2 (v) are the column vector of assembled expansion coefficients of all the sub domains and the mass matrix M (2) K −1 is given as where Ψ 2 i (x) are the basis functions associated to the finite dimensional space D (M i ).
For elements p, q ∈ S (Ω T ), the L 2 -inner product is given by where N 3 (p), N 3 (q) are the column vectors of assembled expansion coefficients of all the sub domains, and the mass matrix M (3) is the block diagonal mass matrix given by where Ψ 3 i (x) are the basis associated with the finite dimensional space S (M i ).
Let p, q ∈ S (Ω T ) × S (Ω T ), then the inner product between the elements is given by where, N 0 (q) , N 3 (p) are the column vectors of assembled expansion coefficients. Then the finite dimensional sub space D b (∂Ω e ) ⊂ H 1/2 (∂Ω e ) is defined as For elements λ, tr u ∈ D b (∂Ω e ) × D b (∂Ω T ), we have the inner product given by where B 0 (λ), B 2 (µ) are the assembled coefficients over all sub domains. If µ = u · n| ∂M i for i = 1, . . . , T , i.e. flux component at sub domain boundaries, then using (5) we have that and we can write (24) as where, N 2 (u) is the column vector of assembled expansion coefficients and N 2 is the assembled trace matrix of all the sub domains.

The divergence operator for broken Sobolev spaces
Let q, u ∈ S (Ω T ) × D (Ω T ) be the elements of the broken finite dimensional spaces. Then the divergence operation on the vector field u is given by where, N 0 (q) , N 2 (u) are the column vectors of the assembled expansion coefficients, and E 3,2 is the assembled divergence operator. If E 3,2 i is the discrete representation of divergence operator on the domain M i , then we have If the topology of sub domain discretizations is also the same, see Example 1, then we have that E 3,2 T . In this paper we have only used the case with the same topology for all the sub domains, therefore where E 3,2 i is the topological divergence operator for any of the sub domains.

The gradient of dual representations
Let p be a scalar field represented by algebraic dual representation S (Ω T ), andp be the boundary value of the scalar field on the sub domain boundary faces γ ∈ ∂Ω e . The gradient operation for dual representations p ,p ∈ S (Ω T ) × D b (∂Ω e ) is then defined as, grad : For the pair of elements p ,p ∈ S (Ω T ) × D (∂Ω e ) the H 1 -norm is then defined as

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 u is the velocity, K is the symmetric positive definite permeability tensor, p is the pressure, f is the given right hand side term, n is the outward unit normal vector,û is the given velocity boundary condition imposed on Neumann boundary Γ N andp is the pressure boundary condition imposed on the Dirichlet boundary Γ D . The Lagrange functional for continuous formulation of Darcy equations is given by The algebraic system for (33), using algebraic dual representations of [13], is given by The above system can be solved for unknowns of p, and the unknowns of u as For the DD formulation of (32) we break the domain Ω into T 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 u, 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. Proof. In (37) if we take variations with respect to u, we get For any sub domain M i , we have If the solution is sufficiently smooth, then using integration by parts on the second term, we get Now, combining the first and the second term, and the third, fourth and fifth term we get This implies that K −1 u + grad p = 0 in the L 2 -sense, i.e. almost everywhere, and λ − p = 0 in the H 1/2 -sense, i.e. between the sub domains and p =p along Γ D . As this should hold for all v ∈ H (div; M i ), we have that λ = p, i.e. the Lagrange multipliers are the pressure boundary values on the sub domain boundaries, ∂M i . The optimality conditions at continuous level for the Lagrange functional (37) are given by: For given f ∈ L 2 (Ω), p ∈ H 1/2 (Γ D ),û ∈ H −1/2 (Γ N ), find u ∈ H (div; M i ), p ∈ L 2 (M i ), λ ∈ H 1/2 (∂Ω e ), such that The finite dimensional problem is then given by: For given f ∈ L 2 (Ω) In (39) we see that all terms, except one -the weighted inner product term, are the inner product between the primal and the dual representations. Consequently, these terms do not require evaluation of (dense) mass matrices. The only matrices associated with these terms are the sparse, metric-free divergence operator, or the inclusion matrix. See also (40).

The algebraic formulation
The inner product for the left hand side terms of (39) are evaluated as and the inner product for the right hand side terms of (39) are evaluated as Using above relations we can write the algebraic formulation for (39) as In (40), we see that the matrices E 3,2 , N 2 , are sparse and metric-free. By metric-free, we mean independent of the size of the elements, the shape of the elements (orthogonal or highly curved) or the order of the approximation. All the metric dependence is contained in the basis functions and therefore in the mass matrix M (2) K −1 . Using static condensation (40) can be solved efficiently for the trace variables B 0 (λ) only. We first solve for a global system of B 0 (λ) by where, The matrix A is block diagonal and can be constructed efficiently by evaluating for each sub domain separately as The inverse matrices in (42) are symmetric, positive definite and are evaluated using Cholesky decomposition. The left hand side of (41) is also symmetric, positive definite and the λ system is solved using Cholesky decomposition. The local degrees of freedom in sub domains M i , i = 1, . . . , T are then evaluated as The inverse terms in (43) and (44) are already evaluated in (42), and therefore evaluation of expansion coefficients of velocity and pressure field is simply a matrix multiplication step.

Test cases
In this section we present the computational results for two test cases using the DD formulation. All the simulations are executed on MATLAB release 2020b on a Macintosh machine with 2.6 GHz Intel Core i7 processor using a single processor.

Test case I: Manufactured solution
In this section we will solve a test problem from [14] and compare the results of DD formulation with the continuous formulation for hexahedral elements of order N = 1, 2, 3 with varying mesh refinements. We will show that i) the results from both the formulation are same up to machine precision, ii) that the DD formulation has optimal convergence rates, iii) the speed-up in simulation run times using DD formulation, and iv) comparison of condition number of the (35) and (41).
The domain for the problem is obtained by mapping the reference domain (ξ , η , ζ) ∈ Ω = [−1, 1] 3 by x =x + 0.03 cos (3πx) cos (3πŷ) cos (3πẑ) y =ŷ − 0.04 cos (3πx) cos (3πŷ) cos (3πẑ) z =ẑ + 0.05 cos (3πx) cos (3πŷ) cos (3πẑ) where, In Figure 2, in the left plot we show the reference domain Ω and on the right plot we show the domain of the problem which is obtained using (45). The permeability tensor, K, and the exact solution, p ex , are given by and p ex (x, y, x) = x + y + z − 1.5 , and the right hand side term is given by As in [14], we impose the Dirchlet boundary conditionsp atx = 0, 1 faces and Neumann boundary conditionsû at y = 0, 1,ẑ = 0, 1 faces, wherep For DD formulation of this test case we decompose the domain Ω with equal number of sub domains, K1, in each direction. For N = 1, 2 cases each of the sub domains is discretized into 2 × 2 × 2 elements. For N = 3 case each sub domain consists of a single element only. Therefore, total number of elements in one direction, K, are 2K1 for N = 1, 2 case, and K1 for N = 3 case. We define h = 2/K as the size of a non-deformed element of the reference domain Ω = [−1, 1] 3 .
In Figure 3 we compare the pressure and velocity profiles of the continuous formulation and the DD formulation for K = 3, N = 2 case. In the first row we plot the pressure, in the second row we plot the x-component of the velocity, in the third row we plot the y-component of the velocity, and in the fourth row we plot the z-component of the velocity. In the first column we plot the results from continuous element formulation, in the second column we plot the results from DD formulation and in the third column we plot the difference of results between both the formulations. In the third column we see that the maximum difference between the continuous and the DD formulation for any of the pressure or velocity profiles is 10 −13 .
In Figure 4 we show the error convergence results from the DD formulation. At the top-left we show the convergence of the L 2 -error for the constraint (div u − f ex ), at top-right we show the convergence of the error in the H (div; Ω) norm for the velocity field, and in the bottom-centre we show the convergence of the error in the H 1 (Ω) norm (using (31)), for the pressure field. On the x-axis we have the length of the non-deformed element h. All the error plots show optimal rate of convergence of order O (N).   In Table 1, Table 2, Table 3 we give the average simulation times to solve the continuous formulation, the DD formulation and the comparison between the two cases respectively. The simulation time (in seconds) is the average for five simulation runs.
In Table 1, we give the average simulation time to solve for continuous formulation for elements of order N = 1, 2, 3. In the first column we have the number of elements in one direction, in the second column we have the time taken to set-up the matrices and in the third column the time taken to solve (35) and (36). The simulation times that are less than 0.1 second, have not been measured with further accuracy. For N = 1 the system run out of memory for K = 64, for N = 2 at K = 32, and for N = 3 at K = 16.
In Table 2, we give the average simulation time to solve for DD formulation for elements of order N = 1, 2, 3. In the first column, we give the total numer of elements in one direction. In the second column we give the number of sub domains in one direction. In the third column we give the number of elements in each direction within each sub domain. In the fourth column we give the time to set-up the matrices. In the fifth column we give the time to solve (41). In the sixth column we give the time to solve (43) and (44). In the seventh column we give the total time, and in the eigth column we give the % amount of time spent to evaluate (41). For the cases where simulation times are less  We observe in the last column that for all the three cases i.e. N = 1, 2, 3, as we increase the total number of elements the percentage of time spent on solution of the Lagrange multiplier system increases. In the last case listed in the table for N = 1, 2, 3, the time spent in solution of (41) is above 80%. Given that DD formulation is for large simulations, we see that in this case the majority of the time is spent on solution of the Lagrange multiplier system (41).
In Table 3 we compare the total time taken to solve the continuous formulation and the DD formulation for elements of order N = 1, 2, 3. The comparison is made of cases with same total number of elements in one direction, such that topology of the mesh discretization remains the same. In the first column we have the total number of elements in one direction K, in the second column we have the total simulation time from continuous elements formulation, in the third column we have the total simulation time from DD formulation. In the fourth column we Firstly, we see that for all N, continuous formulation runs out of memory for lower values of K. That is, the DD formulation requires relatively less memory than the continuous formulation. Secondly, the simulation time required, for the same refinement, is less for DD method than for continuous formulation in all cases. Thirdly, the speed-up factor increases as the mesh is more refined, i.e. the total number of elements are increased. In Figure 5 we compare the condition number of the global system of (35) and (41). We observe that in both the cases, for N = 1, 2, 3, and same discretization, the condition number of the Lagrange multiplier system (41) is higher than the condition number of the pressure unknowns of (35). Moreover, the rate of growth of condition number for N = 1, 2, 3 is similar in both the cases.

Test case II: SPE 10
In this section we show the results for solution of SPE10 benchmark problem [15]. It is often used for validation of numerical schemes for reservoir modelling applications because of its challenging permeability field. The domain of the problem and the natural logarithm of the permeability field is shown in Figure 6. The size of the domain is 1200 f t×2200 f t×170 f t. The domain is divided into equal blocks of size 20 f t×10 f t×2 f t, i.e. 60×220×85 = 1122000 blocks. Each block has a constant isotropic permeability tensor K i = k i I. The top 70 f t of the domain represents the Tarbert formation that has smooth changes in the permeability field component, and the bottom 100 f t represents the Upper Ness formation that has sharp changes in the permeability field component.
The numerical solution of this problem using the continuous formulation was not possible because the system ran out of memory and therefore we only present the results with DD method. This shows that DD method can also be advantageous in cases to reduce the memory requirements for large simulations. For numerical solution using the DD method, we use two different approaches. In Case 1, we divide the domain into 15 × 55 × 17 sub domains, and each sub domain is further divided into 4 × 4 × 5 elements. In Case 2, we divide the domain into 12 × 44 × 17 sub domains, and each sub domain is further divided into 5 × 5 × 5 elements. We also considered a third case with 6 × 22 × 17 sub domains, and each sub domain divided into 10 × 10 × 5 elements, but this case exceeded the memory bounds of the system. The total number and topological configuration of the elements remain the same in all the cases. We use only the lowest order elements, i.e. N = 1, in both the cases. The results of the pressure field from both the cases are shown in Figure 7. In the left plot we see the results from Case 1, in the centre plot we see the results from Case 2, and in the right plot we see the difference between the two results. We observe that the maximum difference in results is of order 10 −12 .
In Table 4 we give the average run time of five simulation runs to solve the SPE10 case using DD formulation. In the second column we have the total number of elements, in the third column we have the total number of sub domains, in the fourth column we have the total number of elements within each sub domain, in the fifth column we have the set-up time of matrices, in the sixth column we have the time taken to solve (41), in the seventh column we have the time taken to solve p & u. (44) & (43), in the eighth column we have the total time, and in the ninth column we have the % of time taken to solve the Lagrange multiplier system (41) with respect to the total time. The results show that for the same total number and topology of elements, the choice of sub domain decomposition also affects the simulation times.
This test case demonstrates that use of algebraic dual representations with DD method can be applied to practical applications.

Conclusions
In this paper we have presented the use of algebraic dual spaces for DD method for Darcy flow. We have defined the broken Sobolev spaces and their finite dimensional counterparts. We have also defined the global finite dimensional trace space that connects the broken spaces. These spaces are used to solve the DD formulation of Darcy equations. It is shown that using algebraic dual representations the matrix representation of continuity constraint across the sub domains becomes sparse and metric-free. The first test case is a manufactured solution where it is shown that i) the results from continuous formulation and DD formulation are precisely the same and that DD formulation also has optimal rate of convergence of errors, ii) the DD formulation is more memory efficient, iii) the DD formulation is more efficient in terms of simulation run times. In the second test case we show that DD scheme can be used to reduce the memory requirements and solve for large practical applications.
We have demonstrated that algebraic dual spaces can be used with DD schemes. In future this work will be extended to broken H 1 (Ω) and H (curl; Ω) spaces, to address problems such as vector Laplacian, Stokes flow, and linear elasticity.