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.

• 3 publications
• 1 publication
• 2 publications
05/30/2020

Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows

Proper orthogonal decomposition (POD) stabilized methods for the Navier-...
10/18/2019

An exactly mass conserving space-time embedded-hybridized discontinuous Galerkin method for the Navier-Stokes equations on moving domains

This paper presents a space-time embedded-hybridized discontinuous Galer...
12/17/2021

HHO methods for the incompressible Navier-Stokes and the incompressible Euler equations

We propose two Hybrid High-Order (HHO) methods for the incompressible Na...
09/13/2019

Divergence-free tangential finite element methods for incompressible flows on surfaces

In this work we consider the numerical solution of incompressible flows ...
03/15/2020

A novel staggered semi-implicit space-time discontinuous Galerkin method for the incompressible Navier-Stokes equations

A new high order accurate staggered semi-implicit space-time discontinuo...
12/22/2021

Domain Decomposition in space-time of 4D-VAR Data Assimilation problem: a case study on the ROMS software

Domain Decomposition of 4D-VAR Data Assimilation (DD-4DVAR) is made up o...
02/04/2021

Aitken-Schwarz heterogeneous Domain Decomposition for EMT-TS Simulation

In this paper, a Schwarz heterogeneous domain decomposition method (DDM)...

1 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

1977Raviart, 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

 H1/2(∂Ω):={λ∈L2(∂Ω):∃ p∈H1(Ω)s.t.λ=p|∂Ω},

and we denote by its dual space.

Let be broken into non-overlapping open sub domains with Lipschitz boundary , , such that

 (1)

Let be the set of sub domains, and be the set of boundaries of these sub domains defined as

 ΩT={Mi}i=1,…,T∂ΩT={∂Mi}i=1,…,T.

We define the broken Sobolev spaces for the set of sub domains as

 L2(ΩT):= ∏M∈ΩTL2(M) H(div;ΩT):= ∏M∈ΩTH(div;M). H−1/2(∂ΩT):= ∏Γ∈∂ΩTH−1/2(Γ)

Let the set of boundary faces, be defined as

 ∂Ωe={γij=∂¯¯¯¯¯¯¯Mi⋂∂¯¯¯¯¯¯¯Mjfori
 ∂ΩΓ={γi=∂Mi⋂∂Ωfori=1,…,T}.

Then we define the trace space

 H1/2(∂Ωe)={λ|γ∈L2(γ),γ∈∂Ωe}=∏Γ∈∂ΩeH1/2(Γ).

Notation

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

 Ψk(x)=(ϵ(k)1(x)ϵ(k)2(x)…ϵ(k)dk(x)),

where is the dimension of polynomial vector space of . The expansion coefficients are always represented as column vectors

 Nk(u)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Nk1(u)Nk2(u)⋮Nkdk(u)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

We will Lemma 2 from 2020jain, which states

Lemma 1.

If and are basis functions from primal and dual representations respectively, then these bases are bi-orthogonal with respect to each other

 ∫Ω˜Ψ3−k(x)Ψk(x)dΩ=I, (2)

where

is the identity matrix of dimension

.

Corollary 1.

The inner product between variables of primal and dual representation is vector dot product of expansion coefficients.

Proof.

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 D(Ω)⊂H(div;Ω)

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

 Ψ2(x)=(ϵ(2)1(x)ϵ(2)2(x)…ϵ(2)d2(x)),N2(u)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝N21(u)N22(u)⋮N2d2(u)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,respectively. (3)

Here, denotes the net flux through the face .

Then we can represent as

 u(x)=Ψ2(x)N2(u).

Using this notation, for two elements the -inner product is given as

 (u,v)=∫Ωu⊺vdΩ=N2(u)⊺(∫ΩΨ2(x)⊺Ψ2(x)dΩ) N2(v)=N2(u)⊺M(2)N2(v),

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

 (u,v)K−1=∫Ωu⊺K−1vdΩ=N2(u)⊺(∫ΩΨ2(x)⊺K−1(x)Ψ2(x)dΩ) N2(v)=N2(u)⊺M(2)K−1N2(v),

where is the mass matrix of the weighted inner product.

3.1.2 Finite dimensional space S(Ω)⊂L2(Ω)

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

 Ψ3(x)=(ϵ(3)1(x)ϵ(3)2(x)…ϵ(3)d3(x)),N3(p)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝N31(p)N32(p)⋮N3d3(p)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,% respectively,

then we can represent as

 p(x)=Ψ3(x)N3(p),

where denotes the integral of over mesh volume .

Using this notation, the -inner product for two elements is given as

 (p,q)=∫Ωp⊺q dΩ=N3(p)⊺(∫ΩΨ3(x)⊺Ψ3(x)dΩ) N3(q)=N3(p)⊺M(3) N3(q), (4)

where, is the mass matrix associated to the basis functions .

3.1.3 Db(∂Ω)⊂H−1/2(∂Ω)

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

 Ψ2b(x):=Ψ2(x)N2andB2(u⋅n):=N⊺2N2(u). (5)

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 E3,2

For any element , the divergence operation on is defined as (see (2020jain, §4.4)): , such that

 div u=Ψ3(x)E3,2N2(u), (6)

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 .

Example 1.

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

 ∫Kdiv u dK=∫∂Ku⋅n dγ=−N1(u)left+N1(u)right−N1(u)bottom+N1(u)top. (7)

If we assemble (7) for all the nine elements, with appropriate numbering, we get the discrete divergence operator as

 E2,1=\resizebox436.87215pt$⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−100100000000−100100000000000−1001000000−10010000000000000−10010000−10010000000−10010000000000−1001000000000−100100000000−100100000000000−1001000000−100100000−1001000000000000−10010000000−10010000000000−1001000000000−100100000000−1001⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.$

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 ˜S(Ω)⊂L2(Ω)

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

 \definecolor[named]pgfstrokecolorrgb0,0,0\pgfsys@color@gray@stroke0\pgfsys@color@gray@fill0q(x)=˜Ψ0(x)˜N0(\definecolor[named]pgfstrokecolorrgb0,0,0\pgfsys@color@gray@stroke0\pgfsys@color@gray@fill0q),

where

 ˜Ψ0(x)=Ψ3(x)(M(3))−1,and,˜N0(\definecolor[named]pgfstrokecolorrgb0,0,0\pgfsys@color@gray@stroke0\pgfsys@color@gray@fill0q)=M(3)N3(q).

Using Corollary 1, the -inner product between the elements is then given by

 (q,p)=∫Ωq⊺p dΩ=˜N0(\definecolor[named]pgfstrokecolorrgb0,0,0\pgfsys@color@gray@stroke0\pgfsys@color@gray@fill0q)⊺N3(p). (8)

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 ˜Db(∂Ω)⊂H1/2(∂Ω)

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

 λ(x)=˜Ψ0b(x)˜B0(λ). (9)

where

 ˜Ψ0b(x)=Ψ2b(x)(M(2)b)−1˜B0(λ)=M(2)bB2(λ),with,M(2)b=∫Ψ2b(x)⊺Ψ2b(x)dΓ. (10)

The inner product between the elements , is given by

 (11)

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

 ∫Ω˜grad(p,^p)u dΩ=−∫Ωp(div u) dΩ+∫∂Ω^p(u⋅n) dΓ∀u∈D(Ω). (12)

The expansion coefficients of are given as

 (13)

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 D(ΩT)⊂H(div;ΩT)

For the set of domains we define the finite element space as

 D(ΩT):=∏M∈ΩTD(M). (14)

For any two elements , and symmetric positive definite permeability tensor , the weighted -inner product is given by

 (15)

where, and are the column vector of assembled expansion coefficients of all the sub domains and the mass matrix is given as

 (16)

where are the basis functions associated to the finite dimensional space .

3.4.2 Finite dimensional space S(ΩT)⊂L2(ΩT)

For the set of domains we define the finite element space as

 S(ΩT):=∏M∈ΩTS(M). (17)

For elements , the - inner product is given by

 (p,q)=∫Ωp⊺q dΩ=T∑i=1∫Mip⊺q dM=N3(p)⊺M(3)N3(q), (18)

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

 M(3)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣M(3)1M(3)2⋱M(3)T⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,with,M(3)i=∫MiΨ3i(x)⊺Ψ3i(x)dM, (19)

where are the basis associated with the finite dimensional space .

3.4.3 Finite dimensional space Db(∂ΩT)⊂H−1/2(∂ΩT)

For the set of boundaries we define the trace space as

 Db(∂ΩT):=∏M∈ΩTDb(∂M). (20)

3.4.4 Finite dimensional space ˜S(ΩT)⊂L2(ΩT)

The algebraic dual representation of is denoted by and defined as

 ˜S(ΩT):=∏M∈ΩT˜S(M). (21)

Let , then the inner product between the elements is given by

 (p,q)=T∑i=1∫Miq⊺p dM=˜N0(q)⊺N3(p), (22)

where, are the column vectors of assembled expansion coefficients.

3.4.5 Finite dimensional space ˜Db(∂Ωe)⊂H1/2(∂Ωe)

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

 ˜Db(∂Ωe)=∏γ∈∂Ωe˜Db(γ). (23)

For elements , we have the inner product given by

 (λ,tr u)=T∑i=1∫∂Miλ⊺tr u dΓ=T∑i=1˜B0i(λ)⊺B2i(μ), (24)

where , are the assembled coefficients over all sub domains.

If for , i.e. flux component at sub domain boundaries, then using (5) we have that

 B2(μ)=B2(u⋅n)=N⊺2N2(u), (25)

and we can write (24) as

 (λ,u⋅n)=˜B0(λ)⊺N⊺2N2(u), (26)

where, is the column vector of assembled expansion coefficients and is the assembled trace matrix of all the 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

 (q,div u)=T∑i=1∫Miq div u dM=˜N0(q)⊺E3,2N2(u)∀q∈˜S(ΩT), (27)

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

 E3,2=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣E3,21E3,22⋱E3,2T⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (28)

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

 E3,2=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣E3,2iE3,2i⋱E3,2i⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (29)

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

 (31)

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

 (32)

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

 L(u,p;f,^p,^u)=∫Ω12u⊺K−1u dΩ−∫Ωp(div u−f)dΩ+∫ΓD^p(u⋅n)dΓ+∫ΓN(u⋅n−^u)dΓ. (33)

The algebraic system for (33), using algebraic dual representations of 2020jain, is given by

 (34)

The above system can be solved for unknowns of , and the unknowns of as

 ˜N0(p) = (E3,2M(2)K−1−1E3,2⊺)−1(N3(f)+E3,2M(2)K−1−1N2˜B0(^p)) (35) N2(u) = M(2)K−1−1(E3,2⊺˜N0(p)−N2˜B0(^p)) (36)

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

 L(u,p,λ;f,^p,^u) = (37) +∫∂Mi∩ΓNλ(u⋅n−^u)dΓ+∫∂Mi∩ΓD^p(u⋅n)dΓ},

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.

Lemma 2.

The Lagrange multipliers are pressure boundary values on sub domains.

Proof.

In (37) if we take variations with respect to , we get

 T∑i=1{∫Miv⊺K−1u dM−∫Mi(div v)p dM+∫∂Mi∖ΓD(v⋅n)λ dΓ+∫∂Mi∩ΓD^p(v⋅n)dΓ}=0∀v∈H(div;Ωt).

For any sub domain , we have

 ∫Miv⊺K−1u dM−∫Mi(div v)p dM+∫∂Mi∖ΓD(v⋅n)λ dΓ+∫∂Mi∩ΓD(v⋅n)^pdΓ=0∀v∈H(div;Mi).

If the solution is sufficiently smooth, then using integration by parts on the second term, we get