Numerical and stochastic models of flow Pichot, B. Poirriez in 3D - - PowerPoint PPT Presentation

SLIDE 1

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Numerical and stochastic models of flow in 3D Discrete Fracture Networks

J-R. de Dreuzy 1, J. Erhel 2, G. Pichot2, B. Poirriez 3

1CNRS, G´

eosciences Rennes, 2INRIA, Rennes, 3INSA, IRISA, Rennes

Special Semester on Multiscale Simulation and Analysis in Energy and the Environment Linz, October 2011

SLIDE 2

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Context of the SAGE team

common team to INRIA and IRISA in Rennes, Brittany http ://www.irisa.fr/sage MICAS project (funded by ANR) collaboration with Geosciences Rennes, University of Lyon, University of Poitiers, University of Rennes. study of flow in fractured media and fractured porous media study of flow and transport in heterogeneous porous media uncertainty quantification : random transmissivty fields and random networks

• f fractures

software platform H2OLab, with libraries MPFRAC and PARADIS Some other related projects GEOFRAC, with Estime team at Inria Rocquencourt, on fractured porous media GRT3D, global reactive transport model, with Andra and Momas group (nuclear waste disposal) H2OGuilde, a user interface for H2OLab Hydromed, with Morocco and Tunisia, on uncertainty quantification and inverse problems Arphymat, with Archeosciences Rennes and IPR, on prehistoric fires

SLIDE 3

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Outline

1

Problem description Geometry Flow equations

2

Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method

3

Linear solvers

4

Domain Decomposition method Domain decomposition of a DFN PCG applied to the Schur complement Results with Schur

5

Conclusion

SLIDE 4

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Geometry Flow equations Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Stochastic Generation of DFN

Following a Discrete Fracture Network approach, fractures are planes with the following random properties : Parameter Random distribution length power law shape disks / ellipses position uniform

• rientation

uniform Conductivity homogeneous / correlated log-normal

Example of DFN with 217 fractures

SLIDE 5

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Geometry Flow equations Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Stochastic Generation of DFN

The broad natural fracture length distribution is modeled by a power law distribution (Bour et al, 2002) : p(l)dl = 1 a − 1 l−a l−a+1

min

dl, where p(l)dl is the probability of observing a fracture with a length in the interval [l, l + dl], lmin is the smallest fracture length, and a is a characteristic exponent.

SLIDE 6

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Geometry Flow equations Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Flow model

Current assumptions : The rock matrix is impervious : flow is only simulated in the fractures, Study of steady state one phase flow, There is no longitudinal flux in the intersections of fractures. Flow equations within each fracture Ωf : ∇ · u(x) = f(x), for x ∈ Ωf , u(x) = −T (x)∇p(x), for x ∈ Ωf , p(x) = pD(x),

• n ΓD ∩ Γf ,

u(x).ν = qN(x),

• n ΓN ∩ Γf ,

u(x).µ = 0,

• n Γf \{(Γf ∩ ΓD) ∪ (Γf ∩ ΓN)},

ν (resp. µ ) outward normal unit vectors T (x) a given transmissivity field (unit [m2.s−1]), f(x) ∈ L2(Ωf ) sources/sinks. Continuity conditions in each intersection : pk,f = pk,

• n Σk, ∀f ∈ Fk,
• f ∈Fk

uk,f .nk,f = 0,

• n Σk,

with Fk the set of fractures with Σk (the k-th intersections) on the boundary,

SLIDE 7

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Conforming mesh at the intersections

Mixed-Hybrid Finite Element Method (MHFEM) for DFNs

• Ref. J. Erhel et al., SISC, Vol. 31, No. 4, pp. 2688-2705, 2009

Makes it easy to deal with complex geometry ; Conforming mesh at the fracture intersections ; A linear system with only trace of pressure unknowns : AΛ = b, with A a symmetric positive definite matrix, the flux at the edges and the mean pressure are then easily derived locally on each triangle. Specific mesh generation :

1

A first discretization of boundaries and intersections is done in 3D by using elementary cubes.

2

The discretization of the boundaries and intersections within the fracture f is obtained by a projection of the previous voxel discretization within the fracture plane.

3

Some local corrections are made to ensure topological properties.

4

Once the borders and intersections are discretized, a 2D mesh of each fracture is generated, using triangular elements.

SLIDE 8

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Channelling in large fractures

In DFN, flow is highly channelled = an opportunity to reduce the number of unknowns and the computational cost, by using a non conforming mesh at intersections. How to apply a mortar method [Arbogast et al., SINUM 37(4) (2000)] to DFN ?

SLIDE 9

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Mixed-Hybrid Mortar Finite Element Method applied to DFN

Mixed-Hybrid Mortar Finite Element Method (MHMFEM) for DFNs

• Ref. G. Pichot et al., Applicable Analysis, Volume 89, 1629-1643, 2010
• Ref. G. Pichot et al., SISC, In revision

A method to mesh the fractures independently and to refine the chosen fractures using a posteriori estimators. Same advantages as MHFEM A simple mesh generation A reduced number of unknowns while keeping a solution of good quality a complex numerical method with Mortar conditions A new specific mesh generation : For each fracture f, choose a mesh step, then

1

A first discretization of boundaries and intersections is done in 2D, by using elementary squares, it leads to a stair-case like discretizations.

2

Some local corrections are made to ensure some topological properties.

3

Once the borders and intersections are discretized, a 2D mesh of each fracture is generated, using triangular elements.

SLIDE 10

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Meshing procedure : Example

The discretization of intersections is non matching ⇒ Mortar conditions are required to ensure the continuity of heads and fluxes.

SLIDE 11

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Mortar principle

Mortar method principle : It consists in choosing arbitrarily for each intersection a master fracture (m) and a slave fracture (s). Classical Mortar approach : each edge is either master or slave Specific Mortar method : some edges have several master or slave properties

SLIDE 12

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Master and slave unknowns

Hydraulic head Notation Cell P Interior edge Λin Intersection edge ΛΣ Master edge Λm Slave edge Λs Jump of flux Notation Intersection edge QΣ Master edge Qm Slave edge Qs Conforming mesh : Λm = Λs = ΛΣ Classical Mortar conditions : relations between master and slave unknowns Specific Mortar conditions : new relations with intersection unknowns

SLIDE 13

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Mortar global conditions

Trace of hydraulic head Jump of flux Λs = CΛm Qm + CT Qs = 0 ΛΣ = AmΛm + AsΛs AmT QΣ = Qm AsT QΣ = Qs with C a block matrix of dimension NsxNm, with blocks (Ck) of dimension Nk,sxNk,m for which each block represents the L2-projection from the master side to the slave side with coefficients Cln, l ∈ {1, ..., Nk,s}, n ∈ {1, ..., Nk,m} : Cln =

• |E m

n ∩ E s l |

|E s

l |

• ,

where the notation |E| stands for the length of the edge E. As and Am are ponderation matrices that gives ΛΣ as the mean of Λm and Λs.

SLIDE 14

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Equations

                   D P − Rin RΣ(Am + AsC) Λin Λm

• = f,

Min MΣ(Am + AsC) (MΣ(Am + AsC))T (Am + AsC)T BΣ(Am + AsC) Λin Λm

• −

Rin RΣ(Am + AsC) T P − v = 0.

• btained by inverting locally Poiseuille’s law on each triangle and by expressing the

fluxes in terms of traces of hydraulic head and mean heads in the following equations : First set of equations : mass conservation Second set of equations : continuity of flux through inner edges Third set of equations : continuity of flux through intersection edges

SLIDE 15

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Linear system

For all cases, with conforming or non conforming mesh, we get a linear system of the form   D −R −RT M  

• P

Λ

• =
• f

v

• Then the system reduces to :

AΛ = c, with A = M − RT D−1 R, Λ =

• Λin

Λm

• and c = v + RT D−1f.

Assuming the transmissivity is locally symmetric positive definite, and the network is connected, the matrix J =   D −R −RT M   is symmetric and, with the presence of Dirichlet boundary conditions within at least

• ne fracture, it is positive definite.

Then A is also symmetric positive definite.

SLIDE 16

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Consistency of the results

Criteria checked for all simulations : Null sum of the fluxes over all the system Null sum of the fluxes over all intersections between fractures Boundary conditions satisfied Continuity of the flux on inner edges (that is egdes that are not intersection).

50 fractures, 315 intersections, 41967 edges, Mean head computation

SLIDE 17

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Convergence criterium

Numerical convergence is estimated via a discrete relative L2 error :

1

A computation is performed on a fine mesh Tη that gives a reference mean pressure pη

2

Simulations are performed on coarsened grids Th of mesh step h > η. The mean head obtained on coarse meshes, ph, are then compared with pη [Martin et al. 2005] : ||ph − pη||2

L2(Ω) =

• Tη∈Tη(Πηph − pη)2|Tη|
• Tη∈Tη(pη)2|Tη|

, |Tη| the area of the triangle Tη ∈ Tη Πηph the projection of ph onto the fine mesh Tη.

SLIDE 18

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Mixed Hybrid Finite Element Method Mortar method Linear solvers Domain Decomposition method Conclusion

Convergence analysis

On 25 Monte-Carlo simulations : Parameter Random distribution length power law shape disks position uniform

• rientation

uniform Parameter Value a 3.5 L/lmin 2 NMC 25 Mesh step from 0.05 to 0.09 Density 2

Log10 of the Convergence criterium vs mesh scale

SLIDE 19

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Solving the linear system

Linear system to solve : In the previous section, we have seen that using MHFEM or MHFEM with Mortar leads to a linear system in terms of trace of hydraulic head unknowns : AΛ = c with A a symmetric positive definite matrix. Solvers tested Direct solver :Cholesky factorization : it is robust but memory and CPU expensive Iterative solver : multigrid method : it is efficient in most cases but it fails in some cases Iterative solver : Preconditioned Conjugate Gradient (PCG) method : it is robust with a multigrid preconditioner

SLIDE 20

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Experiments with stochastic DFN

• Ref. B. Poirriez, Ph-D thesis, University of Rennes, December 2011

Generation of 417 random connected Discrete Fracture Networks Experiments using MPFRAC software and scientific platform H2OLab Sequential computations done with a multiprocessor with 16 GB memory, windows server 10 random simulations for each set of parameters Increasing the density of fractures Parameter value dimension L/lmin 3 power law exponent 2.5 ;3.5 ;4.5 density 2 ;3 ;4 ;5 ;6 ;7 mesh step 0.04 Reducing the mesh step Parameter value dimension L/lmin 3 power law exponent 2.5 ;3.5 ;4.5 density 4 mesh step 0.01 ;0.02 ;0.03 ;0.04 ;0.05 ;0.06 ;0.07 ;0.08 ;0.09

SLIDE 21

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Results with a direct solver (UMFPACK)

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

10 10

1

10

2

CPU Time system size n UMFPACK

Increasing the density of fractures

6x10

4 8x10 4 10 5

2x10

5

4x10

5

6x10

5 8x10 5 10 6

2x10

6

10 10

1

10

2

10

3

UMFPACK

CPU Time mesh nb

Reducing the mesh step 387 systems can be solved with this computer CPU and memory requirements increase rapidly with the size

SLIDE 22

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Results with Algebraic MultiGrid (Boomer-AMG from Hypre)

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

10 10

1

10

2

AMG CPU Time system size n

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

100 200 300 400 500 600

Nb of cycles system size n AMG

Increasing the density

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

10 10

1

10

2

10

3

AMG CPU Time system size n

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

100 200 300 400 500 600

AMG Number of cycles system size n

Reducing the mesh step 375 systems can be solved with less than 1000 cycles convergence rate highly depends on the network

SLIDE 23

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Results with PCG + AMG (from Hypre)

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

10 10

1

10

2

PCG CPU Time system size n

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

1 2 3 4 5 6

Nb of iterations system size n PCG

Increasing the density

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

10 10

1

10

2

10

3

PCG CPU Time system size n

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

1 2 3 4 5 6 7 8 9 10

PCG Number of iterations system size n

Reducing the mesh step all systems can be solved convergence rate depends on the network

SLIDE 24

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

A geometric partition of the network

The network : Ω = ∪iΩi, i = 1,..., Nf , Nf total number of fractures. A specific matrix shape : block angular form            A1,1 . . . . . . . . . A1,Nf +1 . . . ... . . . . . . Ai,i Ai,Nf +1 . . . ... . . . AT

1,Nf +1

. . . AT

i,Nf +1

. . . ANf +1,Nf +1            One block for each fracture Same memory complexity as a 2D problem Many fractures : they can be grouped into fracture sets

SLIDE 25

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Domain decomposition of a Fracture Network

Partition into fracture sets The network is decomposed into Nk fracture sets Fk, k = 1, ...Nk, Fk = ∪iΩi, i ∈ Ik Each subdomain is a connected fracture set to ensure that the submatrix is irreducible Some subdomains are floating, with no Dirichlet boundary condition block angular form of the matrix The matrix is permuted and partitioned into blocks B = PT AP =

Nk

• k=1

Bk, with P a permutation matrix.

SLIDE 26

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Schur complement with subdomains

Local block matrix : For a block k, k = 1, ...Nk : Bk =         . . . . . . . . . . . . Bk,k Bk,Nk +1 . . . ... . . . . . . B

T k,Nk +1

. . . B

(k) Nk +1,Nk +1

        , B =

• k

Bk The block Bk,k is symmetric positive definite (SPD). Schur Complement and local Schur complement : S =

Nk

• k=1

RT

k SkRk, with Sk = B(k) Nk +1,Nk +1 − BT k,Nk +1B−1 k,kBk,Nk +1

The Schur complement matrix S is SPD but Sk can be singular. Equivalent system : Sx = b, with b = bNk +1 −

• k

BT

k,Nk +1B−1 k,kbk

SLIDE 27

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Preconditioned Conjugate Gradient

M spd matrix Initialisation : Choose x0 r0 = b − S x0 z0 = Mr0 p0 = z0 Iterations : For j = 0, 1, ... until convergence qj = Spj ⇐ Computation of matrix-vector product αj = rT

j rj

pT

j qj

xj+1 = xj + αjpj rj+1 = rj − αjqj zj+1 = Mrj+1 ⇐ Preconditioning βj = rT

j+1zj+1

rT

j zj

pj+1 = zj+1 + βjpj

SLIDE 28

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Neumann Neumann preconditioning

Mn = ∆

Nk

• k=1

RT

k ˜

S−1

k Rk∆

If the subdomain is not floating, the matrix Sk is non singular and ˜ Sk=Sk If the subdomain is floating, Sk singular ⇒ Non singular approximation ˜ Sk ∆ is a diagonal matrix which accounts for the number of subdomains involved Mn is not given explicitly Solving z = Mnr z = Mnr = ∆

Nk

• k=1

RT

k zk, with zk = ˜

S−1

k

rk, with rk = Rk∆r Solving ˜ Skzk = rk is equivalent to solve Bk

• xk

zk

• =
• rk
• with

Bk =

• Bk,k

˜ Bk,Nk +1 ˜ BT

k,Nk +1

˜ B(k)

Nk +1,Nk +1

• The matrix Bk is symmetric positive definite.

SLIDE 29

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Computing the local products

Matrix-vector computation q = Sp =

Nk

• k=1

qk, with qk = RT

k SkRkp

qk = RT

k

• B(k)

Nk +1,Nk +1Rkp − BT k,Nk +1B−1 k,kBk,Nk +1Rkp

• Cholesky factorization : before iterations

Bk = LkLT

k with Lk =

• Lk,k

Lk,Nk +1 L(k)

Nk +1,Nk +1

• This gives the Cholesky factorization Bk,k = Lk,kLT

k,k

and the Cholesky factorization ˜ Sk = L(k)

Nk +1,Nk +1L(k)T Nk +1,Nk +1

SLIDE 30

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Deflation preconditioning

• Ref. J. Erhel et al., SIMAX, Volume 21, no. 4, pp. 1279-1299, 2000

M is any preconditioner, for example M = Mn Projection matrix Z a basis of a coarse subspace restriction onto the subspace : Sd = ZT SZ Projection P = I − SZS−1

d ZT

Balancing preconditioning Mb = PT MP + ZS−1

d ZT

Deflation preconditioning Md = PT MP with the initial choice x0 = ZS−1

d ZT b

Deflation and balancing are equivalent for this initial choice [R. Nabben and C. Vuik, SISC, 27(5), 2006] Sketch of the proof ZT r0 = 0 and ZT rk = 0 Thus Mdrk = PT Mrk And, during PCG iterations, Md = Mb = PT M

SLIDE 31

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

SCHUR-Solve implementation(in MPFRAC)

as many subdomains as fractures Nf mutual factorization deflation with the signature of subdomains : Z of size Nf implemented as a library and interfaced with MPFRAC in H2OLab experiments with a conforming mesh more details in Baptiste’s Poirriez Ph-D manuscript to be defended in December 2011

SLIDE 32

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Results with SCHUR-Solve

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

10 10

1

10

2

Schur CPU Time system size n

10

5

2x10

5

3x10

5

4x10

5

5x10

5

6x10

5 7x10 5

20 40 60 80 100

Schur Nb of iterations system size n

Increasing the density

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

10 10

1

10

2

10

3

Schur CPU Time system size n

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

20 40 60 80 100

Schur Nb of iterations system size n

Reducing the mesh step

SLIDE 33

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Domain decomposition

• f a DFN

PCG applied to the Schur complement Results with Schur Conclusion

Comparison of solvers

6x10

4

8x10

4

10

5

2x10

5

4x10

5 6x10 5

8x10

5

10

6

2x10

6

4x10

6 6x10 6

8x10

6

0,1 1

PCG AMG Umf Schur relative CPU time system size n

t < 1E5 t < 1.55E5 t < 2.3E5 t < 3.6E5 t < 5.5E5 t < 8E5 t < 6E5

10 20 30 40 50 60 70 80 90 100

Fastest(%) system size n umf amg pcg schur

Schur is reliable and fast

SLIDE 34

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

On-going and future work

Build a posteriori estimators to optimize mesh generation Develop a parallel software Run large scale DFNs simulations to study scalability Run Monte Carlo simulations to derive upscaling rules

SLIDE 35

Numerical and stochastic models of flow in 3D Discrete Fracture Networks J-R. de Dreuzy , J. Erhel , G. Pichot, B. Poirriez Outline Introduction Problem description Derivation of the linear system Linear solvers Domain Decomposition method Conclusion

Important announcement

DD21 will be in Rennes, in June 25-29, 2012. International conference on Domain Decomposition methods co-organized by Inria and University of Caen with an unforgettable excursion to Mont Saint-Michel in Normandy Save the date ! Sponsors are welcome...