Florian Haeberlein
email: florian.haeberlein@yahoo.fr

Laurence Halpern
email: halpern@math.univ-paris13.fr

Anthony Michel
email: anthony.michel@ifpen.fr

Schwarz Waveform Relaxation and Krylov Accelerators for Reactive Transport

Keywords: 65B99, 65F08, 65N55, 65Y05 multi-species reactive transport, global implicit approach, time-space domain decomposition, optimised Schwarz Waveform Relaxation, CO2 geological storage, high performance computing

come   L'archive ouverte pluridisciplinaire

Introduction

Scientists commonly agree to the indeed existing impact of CO 2 on the greenhouse effect. The transition from classical, (i.e. fossil-based), to renewable energy sources has already begun, nevertheless, they still did not exceed 3.5 % on a global level see http: //www.iea.org/publications/freepublications/publication/KeyWorld2014.pdf . As a consequence, in the next decades, our energy will mainly be supplied by fossil fuels. Carbon Capture and Storage (CCS) is seen as a promising way to ensure the transition from fossil based to renewable energy: on the one hand, it allows to keep existing infrastructures in and use the gained know-how on fossil power plants. On the other hand, it allows to provide necessary time to develop experience and build infrastructure in the relatively new field of renewable energies. The CCS approach separates CO 2 from other gases during the energy production process and transforms it into a supercritical or liquid state. It is then eventually transported for short distances by pipelines to the injection well, where it is injected in the subsurface. During and after injection of CO 2 into saline aquifers, several physical and chemical processes appear. The injected CO 2 dissolves partially in water and changes the pH, the water becomes acid and attacks the rock matrix. This changes the geophysical system, e. g. important changes in the porosity and permeability and hence in the way of how the aquifer moves. In order to ensure the reliability of the technical processes and the consequent changes in the subsurface system, a preprocessing numerical simulation has to be undertaken in order to predict it numerically.

Reactive transport models can be used for large field of applications, CO 2 geological storage simulation is only one of them. Those models can be represented by a system of time-dependent transport equations modelled by partial differential equations which are coupled by nonlinear functions that represent the source terms resulting from chemical reactions. In this paper, we concentrate on a two species reactive transport system which is a subsystem of the multispecies reactive transport system presented in [START_REF] Kräutle | General Mutli-Species Reactive Transport Problems in Porous Media: Efficient Numerical Approaches and Existence of Global Solutions[END_REF][START_REF] Häberlein | Time-Space Domain Decomposition Methods for Reactive Transport Applied to CO 2 Geological Storage[END_REF].

The numerical simulation of such reactive transport models in the context of CO 2 geological storage is a quite challenging task. This is mainly due to two different issues. On the one hand, the aim of the simulation is to provide information on a large area within a long time period in order to ensure the reliability of the storage process. The regarded scales are tremendously large, i. e. the interesting time scale is hundreds and even thousands of years and some hundreds of kilometres in space. On the other hand, the characteristics of the problem itself which contains highly different time scales and quite different levels of numerical complexity.

For the discretisation of the partial differential equations, finite volumes and close methods have established themselves for different reasons. On the one hand, they are easy to code, they work on very general meshes which are often based more on geological than on mathematical criteria, and they provide naturally interesting properties like flux conservation which is highly desirable in this context, see [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes, sushi: a scheme using stabilisation and hybrid interfaces[END_REF]. Nevertheless, all different classes of methods lead to a large system of algebraic equations which are nonlinear, due to chemical coupling terms. Efficient approaches to solve those problems have been developed and tested in different applications. Robust and fast methods are now available which combine nonlinear and linear solvers and which have only frugal need of memory capacity (e. g. Jacobian-free Newton-Krylov methods, see [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF]). With the new generation of multicore-multiprocessor computers, new parallel algorithms, based on domain decomposition in space, appeared, designed for linear and nonlinear problems (see the books [START_REF] Smith | Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations[END_REF][START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF][START_REF] Quarteroni | Domain Decomposition Methods for Partial Differential Equations[END_REF] and [START_REF] Ovtchinnikov | One level Newton-Krylov-Schwarz algorithm for unsteady non-linear radiation diffusion problem[END_REF] in this context. Nonetheless, besides the development in the previously mentioned fields, there are several issues, which have not been solved. One important field is the heterogeneity in space and time. As it has been mentioned, only a small part of the simulation domain is highly reactive while most of the domain is close to an equilibrium state. The chemical reactivity reflects in strong nonlinear coupling terms which have an impact on the nonlinear solver. Usually, there are two ways out of this problem. Either, one accepts a high number of iterations in the nonlinear solver with the risk that no convergence can be achieved. Or, one adapts the time step in order to ensure faster convergence with higher probability. In practice, the second choice is done since this is the a good compromise between higher costs and convergence probability. Cutting the time step is quite easy when it is done globally. But, as the reason for cutting the time step is localised, the time step should be cut only locally. Otherwise, the part of the simulation domain, which is not dominated by strong nonlinearities is solved with a much higher precision than needed.

In the context of CCS, it is also important to be able to use different mesh sizes and different time steps in some parts of the computational domains. The Schwarz waveform relaxation (SWR) methods are in their own conception able to deal with different time steps in different subdomains. The domain is split in several subdomains and the numerical approximation is done over a whole time window, iteratively on the subdomains which are coupled by transmission conditions. The choice of the transmission conditions is a crucial issue in order to achieve high-performing algorithms. So called Robin or Fourier transmission conditions, or even second order operators (Ventcel type) have proven to produce much faster algorithms than the "classical" Schwarz boundary operator, for which Dirichlet data exchange the informations between neighbouring subdomains. A procedure for the optimization of the coefficients has been set, theoretical and numerical results have been produced in the last decade, for the linear advection-diffusion-reaction equation, see [START_REF] Halpern | Optimized Schwarz waveform relaxation: roots, blossoms and fruits[END_REF][START_REF] Bennequin | Optimized schwarz waveform relaxation for advection reaction diffusion equations in two dimensions[END_REF] . An analysis of the optimized transmission conditions to be used for the two species reactive transport system considered here can be found in [START_REF] Häberlein | Optimized Schwarz waveform relaxation for nonlinear systems of parabolic type[END_REF]. Note that these algorithms are efficient for overlapping or nonoverlapping subdomains. They can very well handle the coupling between two different models with different space and time steps, see [START_REF] Halpern | Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF] for an ocean-atmosphere coupling.

The Schwarz waveform relaxation algorithm for nonlinear problems has be analysed in various frames, first for sublinear problems, see [START_REF] Gander | A waveform relaxation algorithm with overlapping splitting for reaction diffusion equations[END_REF]. Optimized Schwarz waveform relaxation algorihtm, with nonlinear transmission conditions were first introduced in [START_REF] Halpern | Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF], for the semilinear wave equation. In [START_REF] Caetano | Schwarz waveform relaxation algorithms for semilinear reaction-diffusion[END_REF], the semilinear advection-diffusion reaction equation in two dimensions was considered, ∂ t u -ν∆u + f (u) = 0, where f is constrained only to be in C 2 (R), with f (0) = 0. Nonoverlapping Robin-Schwarz and Ventcell-Schwarz were proposed and analyzed. The main difficulty in this case is that each iterate in a subdomain is solution of a nonlinear problem, whose solution has to be defined properly, and has its own time of existence T n j , where j is the number of the subdomain, and n the number of the iterate. The sequence (T n j ) n is decreasing, and it has been shown that there is a lower bound T * for these times. Then the convergence is achieved inside (0, T * ). From a numerical point of view, a nonlinear system has to be solved in each subdomain at every step, which has been implemented with P 1 finite elements in space, and a linearly implicit Euler scheme in time. It turns out that the requirement of small time interval given by the existence analysis is not compelling. Furthermore nonlinear transmission condition where the coefficients p and q depend on the iterates through the formulas of [START_REF] Bennequin | Optimized schwarz waveform relaxation for advection reaction diffusion equations in two dimensions[END_REF] were successfully implemented. For the nonlinear reactive transport system, with suitable assumptions on the coefficients, a similar analysis holds, this is the Classical Approach we introduce precisely in §3.

However, there is still space for acceleration. Before optimised transmission conditions have been in the spotlight, either different overlap sizes of the subdomains (see [START_REF] Bjørstad | To overlap or not to overlap: A note on a domain decomposition method for elliptic problems[END_REF] for a detailled study of the influence of the overlap size) and/or Krylov-type methods like BiCGStab ( [START_REF] Van Der | Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems[END_REF]) or GMRES ( [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]) have been applied in order to accelerate the convergence of Schwarz-type methods independently of the chose transmission conditions. Brakkee and Wilders ( [START_REF] Brakkee | The Influence of Interface Conditions on Convergence of Krylov-Schwarz Domain Decomposition for the Advection-Diffusion Equation[END_REF]) have studied the influence of interface conditions on the convergence of Krylov-Schwarz domain decomposition methods and showed that an application of a Krylov-type method on the interface problem has no overhead compared to a standard approach but accelerates significantly the convergence speed of the algorithm for all considered transmission condition types. In this work, we are interested in applying this well-established technique to nonlinear time-dependent problems in order to benefit from its accelerating properties known in the linear context. Optimized Schwarz Waveform relaxation for the advection-diffusion equation have been implemented as a preconditioner, for a Krylov method for the interface problem presented in §4, in [START_REF] Bennequin | Optimized schwarz waveform relaxation for advection reaction diffusion equations in two dimensions[END_REF]. The present paper proposes a new space-time parallel algorithm for the Nonlinear Reactive Transport. It relies on the interpretation of the Schwarz algorithm as a Jacobi algorithm for the space-time interface problem. This interface problem can in turn be solved by Newton-Krylov or Krylov-Newton space-time methods. We present two of those methods. For sake of simplicity, the presentation is made in the case of two subdomains. The paper is organised as follows: In section 2 we set up the problem to solve. In section 3 we describe the Schwarz waveform relaxation algorithm. In section 4 the reduction to the interface variables is presented. The two new approaches are described in section 5 . Numerical issues as implementation details and the used framework as well as numerical results in 2D and 3D are given in section 6. Finally in section 7, we present a test-case developed in the industrial platform Arcane, involving space-time adaptation of the classical approach.

Problem Description

We consider the model problem of a coupled two-species reactive transport system

∂ t (φu) + div(-a∇u + bu) -R(u, v) = f u on Ω × [0, T ], ∂ t (φv) +R(u, v) = f v on Ω × [0, T ], (2.1)
The spatial domain is Ω ⊂ R d , d = 1, 2, 3 and the simulation period is [0, T ]. φ(x) > 0 is the porosity of the medium. The mobile species u is subject to a linear transport, defined by the operator Lu := div(-a∇u + bu), including diffusion described by a positive diffusion coefficient a > 0 and advection described by a Darcy field b ∈ R d . The fixed species v is coupled to the mobile species u by a nonlinear coupling term R(u, v). Both species are subject to a force, defined by the right hand side terms f u , f v . Initial conditions are given for the mobile and fixed species

u(x, 0) = u 0 (x), v(x, 0) = v 0 (x) on Ω, (2.2)
and a boundary condition for the mobile species u

Gu = g(x, t), on ∂Ω × (0, T ], (2.3)
where the linear boundary operator G can be of different types, e. g. Dirichlet, Neumann.

This problem arises as a subproblem of multispecies reactive transport systems and appears to be the most challenging subsystem since mobile and fixed species are coupled by nonlinear reaction terms resulting of kinetic reactions. In [START_REF] Kräutle | General Mutli-Species Reactive Transport Problems in Porous Media: Efficient Numerical Approaches and Existence of Global Solutions[END_REF], S. Kräutle shows that, after reducing a general multispecies reactive transport systems with kinetic and equilibrium reactions, the resulting nonlinearly coupled equations are of the same type as system (2.1). He also established the well-posedness of the problem.

For the sake of readability of the forthcoming analysis, we will use the following notation to concentrate the problem defined by equations (2.1, 2.2, 2.3):

∂ t (φw) + Lw + F(w) = f w on Ω × [0, T ], w(x, 0) = w 0 (x)
on Ω,

J w = g(x, t) on ∂Ω × (0, T ].
(2.4)

By setting w = (u, v), F(•) = (-R(•), R(•)), f w = (f u , f v ) w 0 = (u 0 , v 0 ), J =

The Schwarz Waveform Relaxation Algorithm

We decompose the domain Ω into two non-overlapping domains Ω 1 and Ω 2 and call the common boundary Γ = Ω 1 ∩ Ω 2 the interface. We introduce the following Schwarz waveform relaxation algorithm to approximate the solution of (2.4): Suppose the iterates w k-1 1 and w k-1 2 to be given, then one step of the algorithm consists in solving for both subdomains i = 1, 2

∂ t (φw k i ) + Lw k i + F(w k i ) = f w on Ω i × [0, T ], w k i (x, 0) = w 0 (x) on Ω i , J w k i = g(x, t) on (∂Ω i \ Γ) × (0, T ], B i w k i = B i w k-1 3-i on Γ × [0, T ]. (3.1)
Note that for the sake of readability we do not introduce additional indices for the right hand side source term f w of the transport equation and the initial and boundary conditions w 0 and g in Ω i . To finish the definition of the algorithm, we provide an initial guess w 0 i . This initial guess could be the result of a computation on a coarser grid. The interface transmission operators B i , i = 1, 2 are linear operators. Different choices are conceivable, we limit our developments here to linear differential operators of Robin type, i. e.

B i = ∂ ni + p,
with n i the unit outward normal of Ω i on Γ and p ∈ R, p > 0 a constant. The convergence of the algorithm for the linear problem was proved in [START_REF] Häberlein | Time-Space Domain Decomposition Methods for Reactive Transport Applied to CO 2 Geological Storage[END_REF]. As for the nonlinear problem, the convergence can be proved using the tools in [START_REF] Caetano | Schwarz waveform relaxation algorithms for semilinear reaction-diffusion[END_REF]. Other transmission conditions have been considered in those two works, using second order Ventcell operators for instance, but not in the nonlinear context yet. The parameter p can be defined

for two subdomains such as to minimize the convergence factor of the linearized problem.A full theory was set in [START_REF] Häberlein | Optimized Schwarz waveform relaxation for nonlinear systems of parabolic type[END_REF] in the case of semilinear term

R(u, v) = k(v-Ψ(u)) see (6.1), linearized into R(u, v) = k(v -cu) .
In this work, we use the same p across the subdomains, but they also can be chosen different in different subdomains, improving the convergence of the algorithm, see [START_REF] Dubois | Optimized Schwarz Methods for the Advection-Diffusion Equation and for Problems with Discontinuous Coefficients[END_REF].

4 The Classical Approach: A Fixed Point Algorithm for the Interface Problem

It is possible to reduce algorithm (3.1) to the so-called interface variables. Define the operator

M i : (λ, f ) → w i solution of          ∂ t (φw i ) + Lw i + F(w i ) = f w on Ω i × [0, T ], w i (x, 0) = w 0 (x) on Ω i , J w i = g(x, t) on (∂Ω i \ Γ) × (0, T ], B i w i = λ on Γ × [0, T ],
where f = (f w , w 0 , g) represents all source terms excepting the ones on the interface that are represented separately by λ. Define the interface variable at iteration k by

λ k i := (∂ ni + p)w k i .
Owing to the transmission conditions of algorithm (3.1) on the interface and the relation

n 1 = -n 2 ,
we obtain for i = j,

λ k+1 i = (∂ ni + p)w k+1 i = (∂ ni + p)w k j = -(∂ nj + p)w k j + 2pw k j = -λ k j + 2p M j (λ k j , f
). Algorithm (3.1) can therefore be rewritten as

λ k 1 = -λ k-1 2 + 2p M 2 (λ k-1 2 , f ), λ k 2 = -λ k-1 1 + 2p M 1 (λ k-1 1 , f ). (4.1)
The Schwarz algorithm (3.1) is therefore a fixed point algorithm for the interface problem

λ 1 λ 2 = -λ 2 + 2p M 2 (λ 2 , f ) -λ 1 + 2p M 1 (λ 1 , f ) . (4.2)
In a linear context, suppose the operator F(.) to be linear, the operator M i is linear in both arguments and the algorithm (3.1) is the Jacobi method applied to the linear system

Id Id -2p M 2 (•, 0) Id -2p M 1 (•, 0) Id • λ 1 λ 2 = 2p M 2 (0, f ) 2p M 1 (0, f ) , (4.3)
where Id denotes the identity operator. Krylov subspace acceleration can be done by solving the linear interface problem by a Krylov type method like GMRES or CG instead of using a splitting method like Jacobi or Gauß-Seidel (cf. [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF]).

Note finally, that the interface variables are time-space variables, i. e. they represent the transmission values on the entire interface for both subdomains globally in time.

The New Approach

In the nonlinear context, i. e. F(.) is nonlinear, the algorithm described on the interface variables by (4.1) is equivalent to the application of a "bloc-wise" fixed point iteration method for the nonlinear system (4.3). Note that the operators M i are nonlinear. We use an iterative method (Newton's method for instance) in order to realise them.

Schwarz-Newton-Krylov or Nested Iteration Approach

The first approach consists in treating system (4.2) by a Newton-Krylov approach. We seek the zeros of the nonlinear function

Ψ(λ) := -λ 2 + 2p M 2 (λ 2 , f ) -λ 1 + 2p M 1 (λ 1 , f ) - λ 1 λ 2 .
One step n → n + 1 of Newton's method consists in solving the linear system

Ψ (λ n ) • (λ n+1 -λ n ) = -Ψ(λ n ).
The derivative of Ψ is given by

Ψ (•) = -Id -Id +2p∂ λ M 2 (•, f ) -Id +2p∂ λ M 1 (•, f ) -Id .
Owing to the definition of a linearised operator

M lin i : (A, h, f ) → w i solution of          ∂ t φw i + Lw i + Aw i = f w on Ω i × [0, T ], w i (x, 0) = w 0 (x) on Ω i , J w i = g(x, t) on (∂Ω i \ Γ) × (0, T ], B i w k+1 i = h on Γ × [0, T ], (5.1)
one can state the derivative

Ψ (•) = -Id -Id +2p M lin 2 (F (M 2 (λ n 2 , f )), •, 0) -Id +2p M lin 1 (F (M 1 (λ n 1 , f )), •, 0) -Id .
The entire procedure of the approach is then described by the following algorithm:

Input: λ 0 (initial guess), ε (precision), n max (maximum iterations) Return: λ (solution) n = 0 repeat // Set up RHS if n = 0 then
No previous iterate is available, use the previous time step as initial guess for the actual time step during evaluation of operators M i else Use the previous iterate globally in time as initial guess during the evaluation of operators

M i end if -Ψ(λ n ) = λ n 1 + λ n 2 -2p M 2 (λ n 2 , f ) λ n 1 + λ n 2 -2p M 1 (λ n 1 , f )
// Solve the linear system Ψ δ λ n = -Ψ(λ n ) with a Krylov-type method for every Krylov-iteration k do realise a Matrix-vector multiplication by

Ψ δ k λ n = -δ k λ n 1 -δ k λ n 2 + 2p M lin 2 (F (M 2 (λ 2 , f )), δ k λ n 2 , 0) -δ k λ n 1 -δ k λ n 2 + 2p M lin 1 (F (M 1 (λ 1 , f )), δ k λ n 1 , 0) end for // Update variables λ n+1 = λ n + δ λ n n = n + 1 until n = n max or δ λ n < ε or b < λ = λ n
The approach requires in every iteration of the outer loop (indices in n) to set up a right hand side-vector that demands to solve two nonlinear problems in the subdomains. Therefore, a nested iterative procedure is necessary (Newton for instance), for this reason, we call this approach "Nested Iteration Approach" (NIA) due to the split iterative approaches of the nonlinear interface problem and the nonlinear subproblems. The name "Schwarz-Newton-Krylov" can be used in order to explain the order of application of the different tools: The global problem is first attacked by a Schwarz-type domain-decomposition method. The resulting nonlinear interface problem is attacked by a Newton-type method where, at every iteration, the resulting linear system is solved by a Krylov-type method. Unfortunately, the name "Newton-Krylov-Schwarz" has already been widely used for another type of methods and therefore "Schwarz-Newton-Krylov" may be confusing. On these methods for stationary problems, see [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF].

Newton-Schwarz-Krylov or Common Iteration Approach

The second approach is not based on the nonlinear interface problem but attacks the global problem (2.4) up from the beginning. We apply Newton's method on this system and solve in every iteration n → n + 1 the linear system

(∂ t φ + L + F (w n ))(w n+1 -w n ) = -(∂ t w n + Lw n + F(w n ) -f w ) on Ω × [0, T ], w n+1 (x, 0) = w 0 (x)
on Ω,

J (w n+1 -w n ) = -J w n + g(x, t)) on ∂Ω × (0, T ],
which can be reformulated to

(∂ t φ + L + F (w n ))w n+1 = F (w n )w n -F(w n ) + f w on Ω × [0, T ], w n+1 (x, 0) = w 0 (x) on Ω, J (w n+1 ) = g(x, t) on ∂Ω × (0, T ],
owing to the linearity of the operators L and J .

We apply then the domain decomposition method on this linear system as we have done in section 3 and reduce the problem to the interface variables. Note that no further iteration index has to be introduced for the domain decomposition algorithm since no numerical method has yet been assigned to the resulting linear interface problem. The resulting linear problem for every iteration is then given by

Id Id -2p M lin 2 (F (w n 2 ), •, 0) Id -2p M lin 1 (F (w n 1 ), •, 0) Id • λ n+1 1 λ n+1 2 = = 2p M lin 2 (F (w n 2 ), 0, (F (w n 2 )w n 2 -F(w n 2 ) + f w , w 0 , g)) 2p M lin 1 (F (w n 1 ), 0, (F (w n 1 )w n 1 -F(w n 1 ) + f w , w 0 , g))
.

Logically, the values for (u n 1 , u n 2 ) that are needed to evaluate the RHS and the matrixvector multiplication at every iteration have to be provided by a nonlinear solution with operators M i (λ n i , f ). By giving an initial guess also for the first iterate (w -1 1 , w -1 2 ), this procedure can be replaced by calculating the values only by the linearised operators

w n i = M lin i (F (w n-1 i ), λ n i , (F (w n-1 i )w n-1 i -F(w n-1 i ) + f w , w 0 , g)),
because, suppose the solution has converged, the linearised operator gives the same solution as the nonlinear operator.

The entire procedure of the approach is then described by the following algorithm:

Input: λ 0 , (w -1 1 , w -1 2 ) (initial guess), ε (precision), n max (maximum iterations) Return: λ (solution) n = 0 repeat // Update subdomain solutions w n 1 = M lin 1 (F (w n-1 1 ), 0, (F (w n-1 1 )w n-1 1 -F(w n-1 1 ) + f w , w 0 , g) w n 2 = M lin 2 (F (w n-1 2 ), 0, (F (w n-1 2 )w n-1 2 -F(w n-1 2 ) + f w , w 0 , g) // Set up RHS b = 2p M lin 2 (F (w n 2 ), 0, (F (w n 2 )w n 2 -F(w n 2 ) + f w , w 0 , g)) 2p M lin 1 (F (w n 1 ), 0, (F (w n 1 )w n 1 -F(w n 1 ) + f w , w 0 , g))
// Solve the linear system Aλ n+1 = b with a Krylov-type method for every Krylov-iteration k do realise a Matrix-vector multiplication by

Aλ n+1,k = λ n+1,k 1 + λ n+1,k 2 -2p M lin 2 (F (w n 2 ), λ n+1,k 2 , 0) λ n+1,k 1 + λ n+1,k 2 -2p M lin 1 (F (w n 1 ), λ n+1,k 1 , 0) end for // Update variables δ λ n = λ n+1 -λ n n = n + 1 until n = n max or δ λ n < ε λ = λ n
The approach requires in every iteration of the outer loop (indices in n) to set up a right hand side-vector that demands to solve two linear problems in the subdomains. Moreover, in the matrix-vector multiplication inside the Krylov-method, only linear problems in the subdomains are evaluated. No nested nonlinear iterative method is needed. For this reason and in contrast to the first approach, we call this approach "Common Iteration Approach" (CIA) due to the common iterative approach of the nonlinear character of the monodomain problem. The name "Newton-Schwarz-Krylov" can be used in order to explain the order of application of the different numerical tools: The global problem is first attacked by a Newton-type method. At every iteration, the resulting linear problem is decomposed by a Schwarz-type algorithm where the problem is reduced to the interface variables. The resulting linear system is then solved by a Krylov-type method. As in the first case, we recommend not to use this name as the name "Newton-Krylov-Schwarz" has already been used for another type of methods.

Numerical Approach

Krylov-subspace accelerators for the linear case have become very popular because they are easy to implement and have no significant overhead compared to a standard approach for the Schwarz domain decomposition algorithm. In the nonlinear case, they are as easy to implement as in the linear case: Once the nonlinear solver for a problem in a subdomain is available, the effort for adding the linearised solver and the implementation of the matrix-vector multiplication and the right hand side vector is negligible. Moreover, in many environments, sophisticated and ready-to-use Krylovtype methods like GMRES or BiCGStab are available.

Implementation

In our prototype code for the two species nonlinear coupled reactive transport system (2.4), the equations are discretised by an implicit Euler scheme in time and a hybrid finite volume scheme (cf. [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes, sushi: a scheme using stabilisation and hybrid interfaces[END_REF]) in space including an upwind advection discretisation and a two-point diffusion discretisation on simple rectangular meshes. Nonlinear problems in the subdomain solver are treated with a global implicit approach by means of Newton's method. The arising linear system is solved by an exact LU-decomposition and convergence of Newton's method is controlled by the residual and the step size with a precision of ε = 10 -8 .

The Classical Approach has been implemented as presented in (4.1). We have implemented both methods, the Nested Iteration Approach and the Common Iteration Approach. Both methods make use of a Krylov-type method, GMRES for instance. In order both approaches to be competitive, we apply a precision strategy in the mood of inexact Newton methods, i. e. in the first iterations of the Newton method, we will not oversolve the linear system and can limit therefore the number of costly subdomain evaluations within a matrix-vector multiplication. The more we advance in the outer Newton iteration, the more precise the linear system has to be solved. The appearing linear system at iteration n is solved up to a precision of

max min 1 1 + n , Ψ(λ n ) , 10 -8 ,
for the Nested Iteration Approach and

max 10 -1 • min 10 -n , δ λ n , 10 -8 ,
for the Common Iteration Approach.

Note that the Common Iteration Approach needs to store the discretised values of a solution in both subdomains. This can be viewed as a huge drawback in high performance codes. The Nested Iteration Approach a priori does not suffer from this drawback. Anyway, it may be highly desirable to store the solution and use it as initial guess for the evaluation of the right hand side. In practice, using the solution of the previous iterate as initial guess reduces significantly the number of nested Newton iterations in the nonlinear subdomain solver.

Finally, concerning the stopping criterion of the outer Newton iteration, the Nested Iteration Approach can be classically controlled by both the residual and the step size norm. The Common Iteration Approach can though no longer be controlled by the residual norm since we have eliminated that term up from the beginning. Recalculating the residual term afterwards is possible but the cost for this may not be negligible.

The classical approach

We study now the numerical behaviour of the Schwarz waveform relaxation algorithm 3.1. In the first part, we concentrate on issues concerning the influence of the nonlinearity on the performance of the classical approach. In the second part, we compare the classical approach with the two new approaches presented in section 6.3. For the numerical results in this section we fix the time period t ∈ [0, 1] and the global domain

Ω = [0, 1] × [0, 1] ∈ R 2 . Discrete steps are ∆t = ∆x = ∆y = 2 • 10 -2 . The physical pa- rameters are porosity φ = 1, diffusion a = 1, advection is (b x , b y ) = (1 • 10 -2 , 5 • 10 -2 ).
Defining the function f (x, y, t) = (sin(πx) cos(πy) cos(πt)+cos(πx) sin(πy) cos(πt)+cos(πx) cos(πy) sin(πt)+1)/2, we can set the initial values to u 0 = f (x, y, 0), v 0 = f (x, y, 0)/c for (x, y) ∈ Ω and we impose Dirichlet boundary conditions with values set to u b (x, y, t) = f (x, y, t) for (x, y) ∈ ∂Ω. The function f provides a heterogeneity in space and time so that we can ensure that the exact solution that we reconstruct numerically does not degenerate to a stationary problem. We decompose Ω into non-overlapping subdomains Ω 1 = [0, 0.5] × [0, 1] and Ω 2 = [0.5, 1] × [0, 1]. We impose a random initial guess on the interface Γ 1 in order to ensure the presence of a wide range of possible frequencies in the error.

Therefore we study, as in the linear case, the convergence behaviour of the Schwarz waveform relaxation algorithm with Robin transmission conditions using different parameters p for the transmission condition. We proceed ten iterations of the algorithm and focus on the resulting error on the interface values which indicate us the numerical performance of the transmission condition with respect to the parameter p. We study different nonlinear coupling function. First, we consider an adsorption process that is modelled by a BET isotherm law:

Ψ(u) = Q s K L u (1 + K L u -K S u)(1 -K S u) . (6.1)
BET theory is a rule for the physical adsorption of gas molecules on a solid surface and serves as the basis for an important analysis technique for the measurement of the specific surface area of a material (cf. Brunauer et al. [START_REF] Brunauer | Adsorption of Gases in Multimolecular Layers[END_REF]). This law is insofar mathematically interesting as it is neither convex nor concave (cf. figure 6.1). The coupling term is given by R(u, v) = 100(v -Ψ(u)) with Q S = 2, K S = 0.7, K L = 100. We study another nonlinear function that is given by an exponential equilibrium model

R(u, v) = exp(10(2v -3u)) -1,
which models the reaction rate of a kinetic reaction of the form 2v ⇐⇒ 3u where the forward reaction 2v -→ 3u is much faster than the backward reaction 2v ←-3u. In this case, we have no asymptotic analysis at hand, and we use the optimised parameter of the single equation of advection-diffusion and the therefore obtained transmission condition behaves well. In figure 6.2 we plot the error at iteration 10 varying the parameter p of the Robin transmission condition.The square locates the optimised parameter of the linear system established in [START_REF] Häberlein | Optimized Schwarz waveform relaxation for nonlinear systems of parabolic type[END_REF]. In both cases, we observe that the theoretically optimised parameter is close to the best performance for the nonlinear coupled advection-diffusion-reaction system. Tests with other nonlinear functions confirmed this behaviour. The best results are obtained when we have the exact parameter from the linear analysis, that is when the nonlinearity is not too strong.

In the linear case, domain decomposition methods are often used with many subdomains in order to distribute data and computational efforts to several processors. In the nonlinear context, domain decomposition methods can be used in order to localise time step constraints resulting from heterogeneity. In the context of CO 2 geological storage modelling, highly reactive moving fronts appear in the geochemical system. Those regions of strong chemical disequilibrium induce highly localised constraints on the time step using a global implicit approach. If the time step is chosen large, the number of Newton iterations for one time step increases drastically and if the time step is chosen too large, the standard Newton approach fails to converge. The Schwarz waveform relaxation approach provides the possibility to use different discretisations and numerical approaches in the subdomains. We exploit the possibility to choose different time grids in the subdomains. By this way, we can select a "reactive domain" with small time steps in order to keep the number of Newton iterations for the time steps small. In a "non reactive domain" we can chose much larger time steps and the number of Newton iterations stays acceptable. By this way we can limit the time step restrictions to the area where they appear instead of letting them influence the global time step of the whole simulation domain.

We exemplify this feature with the following test case: the time period is [0, 1] and the global domain is

Ω = [0, 1]×[0, 1] ∈ R 2 . Discrete steps are ∆x = ∆y = 2•10 -2 . The diffusion parameter is ν = 5•10 -2 , advection is (b x , b y ) = (1.5, 1.0).
The reaction term is realised by use of the BET isotherm. The initial values are set to (u 0 , v 0 ) = (0, 0) which is an equilibrium state. We model the entry of a reactive front by imposing a Dirichlet boundary condition on x = 0 with values g(x = 0, y, t) = sin(yπ). All other boundaries are set to be of no diffusion type, i. e. we impose no concentration gradient on the boundary. By the incoming reactive front, the chemical system is subject to a strong disequilibrium perturbation. In order the number of Newton iterations not to be too excessive (less than ten), one has to chose a time step of 10 -1 , i. e. ten time steps. The global monodomain approach has a linear system of 5200 discrete unknowns to solve in every time step and every Newton iteration. The first time step needs 9 Newton iterations to reach convergence, the following time steps (2nd to 10th) need each 7 Newton iterations to reach convergence where we suppose the solution of the previous time step as initial guess of the Newton iteration. We measure the effort of the global approach by the effort of the inversion for one matrix multiplied by the number of matrix inversions since this is the most costly operation. The global effort is hence

(1 • 9 + 9 • 7) • (5200) 3 = 1.01 • 10 13 .
In a domain decomposition approach, we can chose the reactive subdomain to be Ω 1 = [0, 0.4 + ∆x] × [0, 1] (2242 discrete unknowns) and the non reactive subdomain to be Ω 2 = [0.4, 1] × [0, 1] (3160 unknowns discrete). Both subdomains have an overlap of one layer of cells. For the reactive subdomain we chose the same time step as in the global monodomain approach, i. e. ten time steps, while for the non reactive subdomain, we can use only one time step without the number of Newton iterations to become important. We impose the initial state as initial guess for the interface values in the Schwarz waveform relaxation iteration and impose the solution of the previous Schwarz waveform relaxation iteration in the Newton iterations. In the first iterate we proceed as in the global approach, i. e. we impose the solution of the previous time step as initial guess for the Newton iterates. In total, we need three Schwarz waveform relaxation iterations in order to reach convergence. In the first iteration, we need for the reactive subdomain 9 iterations in the first time step and 7 iterations for each following time step. In the non reactive subdomain we need 4 iterations. In the second Schwarz iteration, we need 2 iterations for the first 7 time steps and 3 iterations for the last 3 time steps in the reactive subdomain and 3 iterations in the non reactive subdomain. In the third iteration we need 2 iterations for all time steps in the reactive and non reactive subdomains. The total effort is hence

(1 • 9 + 9 • 7 + 17 • 2 + 3 • 3) • (2242 3 ) + (4 + 3 + 2) • (3160) 3 = 1.58 • 10 12 .
The effort for the domain decomposition solution is hence by a factor of 10 smaller than the effort for a global monodomain solution. Note that this estimation is quite optimistic, since using more performing linear solvers for the linear problems during the Newton iterations may reduce the gain of a domain decomposition approach compared to the global approach.

In figures 6.3 and 6.4 we plot the concentration of u and v at t = 0.5 and t = 1.0 at the third iteration of the domain decomposition algorithm.

The new approach

We use a similar testing approach as before and set k = 100, Q S = 2, K S = 0.7 and

K L = 100. Initial values are set to (u 0 , v 0 ) = ( 1 2 , 1 3
). By defining the function g(x, y, t) = (sin(πx) cos(πy) cos(2πt)+cos(πx) sin(πy) cos(2πt)+cos(πx) cos(πy) sin(2πt)+1)/2, we impose Dirichlet boundary conditions with values set to u b (x, y, t) = g(x, y, t) for (x, y) ∈ ∂Ω.

In a first time, we are interested in the sensibility of the three approaches with respect to the parameter p of the Robin transmission condition. We discretise the numerical domain with ∆x = ∆y = 1/40 and ∆t = 1/10 and impose a random initial guess on the interface for the first iteration. As both subdomains have the same size, the number of overall matrix inversions in the linear and nonlinear subdomain solvers for three approaches is a meaningful criterion to measure the numerical performance. We proceed the three approaches for different parameters p of the Robin transmission condition and plot in figure 6.5 the number of matrix inversions. One observers first that the performance of the classical approach, i. e. a fixed point method on the nonlinear interface problem depends highly on the parameter p that one chooses for the Robin transmission condition. The best parameter is p * ≈ 40. The two new approaches, NIA and CIA also show the best performance at p * but are much less sensitive to the choice of the parameter. Especially if, in realistic test cases where one has no idea of the best parameter, one underestimates the unknown parameter p * , the new approaches do not suffer from the exponential loss of performance. Second, one observes that the CIA is always more performing than the classical and the NIA approach. The NIA is, in a wide range of parameters, more performing than the classical approach but the last one stays more performing in an important range of parameters around the optimal It turned out that the two new methods, NIA and CIA, have a cost overhead that becomes non negligible if space discretisations are chosen too coarse. For this reason, we study the asymptotic behaviour of the three approaches using always the optimal parameter of the classical approach. We refine the problem in space using always ∆x = ∆y. Note that we keep the time step constant at ∆t = 0.1. Refining the discretisation also in time would lead to a problem that is quasi linear at every time step since we use a global implicit approach. The negligible nonlinearity would result in a minimal number of nested Newton iterations and the overhead cost would become more important. We measure again the overall number of matrix inversions in the three approaches. One observes that the overhead cost of the two new approaches compared to the classical approach becomes negligible up from a discretisation with about 150 grid points per dimension for the NIA and about 20 grid points per dimension for the SIA. For problems finer than the respective thresholds, the new approaches are always more performing compared to the classical approach with the best parameter for the transmission condition. Moreover, the finer the discretisation, the larger the problem, the more important the accelerating property of the two new approaches. Note that both new approaches have the same slope of O(N 1/7 ) in the asymptotic behaviour which is considerably less than the slope of the classical approach which behaves like O(N 1/2.75 ). The vertical translation of the curves for the two new approaches indicates their overhead cost. Again here, a theoretical anlysis will be interesting.

In order to exemplify the accelerating property of the two new approaches, we perform a simulation with N x = N y = 200 points in each dimension keeping the number of time steps constant and compare the convergence behaviour of the stopping criteria of the three methods. In figure 6.7 we plot the convergence criterion versus the number of matrix inversions. One observes the quadratic convergence of the new approaches since they are Newton-based, the quadratic convergence is observed late in the history since the initial guess (randomly chosen) is far from the exact solution. The classical approach shows only a superlinear convergence, also in this case, the superlinear character appears late in the convergence history.

Finally, we want to apply the two new approaches to a benchmark test case in the context of CO 2 geological storage. The 3D test case is based on the benchmark for the SHPCO2 project (Simulation haute performance pour le stockage géologique du CO 2 ) which is described in [START_REF] Michel | Cas Tests SHPCO2 n • 3 -Test synthétique pour le transport réactif dans le cadre du stockage géologique de CO2[END_REF]. The global domain is set to Ω = [0, 4750] × [0, 3000] × [-1100, -1000] and is decomposed into the two nonoverlapping subdomains Ω 1 = [1000, 2500] × [0, 3000] × [-1050, -1000] and Ω 2 = Ω \ Ω 1 . We call Ω 1 the reactive subdomain since in this subdomain an injection of the mobile species u is modelled by a source term. The initial state is zero for the mobile and immobile species. We consider again the BET isotherm law as nonlinear coupling term. The injected mobile species is partially absorbed by the reaction and partially transported by mainly advection. Simulation time is [0, 100]. As the waveform relaxation approach allows to use different discretisations in the subdomains, we chose to use ten time steps in the reactive subdomain Ω 1 and only five time steps in the subdomain Ω 2 . This choice is insofar justified since the rapid injection in the reactive subdomain restricts the time step size by imposing a maximum number of Newton iterates of ten. As in the subdomain Ω 2 , the mobile species appears only by transport processes on a slower time scale than the injection, one can chose a larger time step in order to respect the maximum number of Newton iterations. Concerning the parameter of the Robin transmission condition, we use the same as before. The initial guess on the interface is zero for both subdomain interfaces. In figure 6.8 we plot the convergence histogram, i. e. the stopping criterion in a logarithmic scale versus the CPU time (normalised to the CPU time of the classical approach). Note that both subdomains have a different size of unknowns and therefore the number of matrix inversions, as used in the previous examples, is no longer a valid tool to measure the effort. One observes clearly that the classical approach converges only linearly while the two new approaches show a quadratic convergence. Moreover, the two new approaches need only about 20 percent of the CPU time of the classical approach. Note that the stopping criterion of the Common Iteration Approach seems to crash down in the last iteration as one observes in figure 6.8. This behaviour is due to the fact that the GMRES solver provided in the last iteration the same solution as in the last but one iteration due to the precision strategy. As a consequence, the norm of the variation of the iterates is zero and this indicates that the algorithm has converged. To indicate this behaviour in the plot, we set the convergence criterion of the last iterate to the overall precision 10 -8 .

In the first test case using the BET isotherm we observed that the Common Iteration Approach is always more performing than the Nested Iteration Approach. However, in the SHPCO2 test case the opposite holds. Note that there is no a priori more performing method, all depends on the problem to treat. We performed several tests with different problems and we can state several points: First, both new approaches, the Nested Iteration Approach and the Common Iteration Approach show always much less sensibility to the parameter of the Robin transmis-sion condition than the classical approach. Second, both new approaches show always a more favourable slope in the case of asymptotic mesh refinement in space on the one hand and on the other hand they show always the same slope up to a vertical translation which indicates the different overhead cost. Third, the question which of both new approaches is more performing on a problem is a challenge between the complexity to solve the problem with a domain decomposition approach and the difficulty to solve it with a nonlinear approach. Problems which are easy to solve with a domain decomposition approach since they are for example quite advective in normal interface direction and less diffusive have concentrated their main challenge in the nonlinearity. This is the case of the SHPCO2 test where the main difficulty lies in the massive injection of mobile species in the reactive domain which induces an elevated number of Newton iterations. The Nested Iteration approach showed to be more performing on this cases where the nonlinearity is the main challenge. On the other hand, for test cases which are difficult to solve with a domain decomposition approach, for example highly diffusive cases like the synthetic case in the previous part, the Common Iteration Approach seems to be more performing. Altogether, we can summarise that both new approaches are more performing than the classical approach whenever the overhead cost due to coarse discretisations become negligible. Moreover, which one of the new approaches is more performing lies on the challenging character of the problem and can hardly be predicted in advance. Nevertheless, passing from one approach to another is very simple since they only inverse two loops.

First numerical tests in industrial environment Arcane

The work presented in this article is part of an international project called SHPCO2 (high performance simulation of CO 2 geological storage), initiated by the French National Research Agency (ANR). One part of the project is to study different domain decomposition approaches for multi-species reactive transport problems in the context of high performance computing for CO 2 geological storage simulation. Arcane is the underlying industrial simulation platform for HPC applications which is under co-development of the two national French institutes IFP Energies nouvelles and CEA Energies alternatives, see [START_REF] Häberlein | Time-Space Domain Decomposition Methods for Reactive Transport Applied to CO 2 Geological Storage[END_REF] for specificities.

A first approach of Schwarz Waveform Relaxation methods has been implemented into the Arcane development platform to show the ability to treat these kind of complex problems in a HPC environment. In the context of a domain specific language for parallel finite volumes (cf. [START_REF] Prud'homme | Basic concepts to design a DSL for Parallel Finite Volume Applications[END_REF]), the classical fixed-point iteration type approach as described in equation (4.2) has been implemented. Spatial grids on the interface are conformal but time-stepping can be different in the sub-domains. This allows for individual approaches of combined time-stepping and nonlinearity control in the different sub-domains. Currently, two sub-domains are treated, one which contains most of the non-linear problems due to chemical reactions. The other subdomain is essentially characterised by the absence of strong nonlinearities.

We have tested the approach with the SHPCO2 test case. This synthetic test case contains all major physical and chemical effects which are characteristic for CO 2 geological storage modelling. The full definition of the test case can be found in CO 2 (aq)

Cl - H + N a + Ca ++ SiO 2 (aq)             , q = q 1 q 2 = Calcite Quartz ,
and secondary species

x = x 1 x 2 = HCO - 3 OH -, z = z 1 = CO 2 (solid) .
The first group being an independent set of species while the second group forms species whose concentration can be determined by solving the associated equilibrium problem. The equilibrium reactions are represented in a Morel tableau indicating the stoechiometric indices of the associated equilibrium reactions:

c 1 c 2 c 3 c 4 c 5 c 6 c 7 c 8 q 1 q 2 x 1 1 1 -1 x 2 1 -1 z 1 1
.

Besides the equilibrium reactions, we include four reactions whose reaction speeds are modelled by different kinetic laws:

Dissolution of Calcite: Calcite + H + -→ Ca ++ + HCO - We simulate the test case on a 2D mesh with 76 cells in x-direction and 48 cells in y-direction. The initial subdomains are chosen such as the reactive domain covers the initial CO 2 -plume and is slightly blown-up. The non-reactive domain contains the part where the system is in chemical equilibrium. Both domains do overlap in the initial state (see figure 7.1). For the interface transmission conditions, we impose Robin conditions with a localised one-dimensional approximation of 0th order of the optimal parameter. The Schwarz waveform relaxation algorithm is applied to a timewindow and iteration is proceeded between the sub-domains until the change of the concentration at the end of a time window is smaller than 10 -4 . We perform a simulation of the first 95.13 years with three equally sized time windows and ten time steps per time window (i. e. ∆t = 3.17 years). In order to show the localisation of numerical difficulties related to kinetic reactions appearing mainly in the reactive subdomain, we chose no adaptive time strategy and concentrate on the number of nonlinear steps needed to solve the subdomains. In figure 7.2, we plot the pH-value at the end of the simulation. The initially solid CO 2 has partially dissolved into the liquid phase where it becomes acid resulting in low pH-values. The associated regions are chemically highly reactive, both from a kinetic and an equilibrium point of view. As a consequence, non-linearities originating from chemistry are concentrated in regions with low pH-value. Note that we have chosen the initial reactive sub-domain such that at the end of the simulation, the reactive domain covers the whole region with low pH-value.

The Schwarz waveform relaxation algorithm needs 2 iterations per time window to converge. The overall number of global nonlinear steps are 421 in the reactive and 306 in the non reactive subdomains. The linear systems the during the nonlinear iteration are more difficult to solve in the reactive domain than in the non reactive domain. The iterative solver of GMRES type preconditioned with Hypre's Euclid preconditioner needed overall 10234 iterations in the reactive and and only 3806 iterations in the non reactive subdomain. This means, that in the reactive subdomain, the linear solver needed 24.3 iterations in average and in the non reactive domain it needed only 12.4 iterations.

Conclusion

Basing on a nonlinear coupled reactive transport system we have developed two new approaches for solving the interface problem in the nonlinear case. The Nested Iteration Approach applies a Newton-Krylov method on the nonlinear interface problem of the Schwarz algorithm. The Common Iteration Approach applies a Newton method on the global problem and solves than in every iteration the appearing linear interface problem of the Schwarz algorithm with a Krylov method. We have implemented and tested both methods. Comparative results with the classical approach that consists in the application of a fixed point method on the nonlinear interface problem have been provided. The numerical tests showed that, besides an overhead cost for coarse space discretisations, both methods have an accelerating property and show much less sensibility with respect to the choice of the parameter of the Robin condition. The quadratic convergence behaviour of the new approaches outperform the superlinear convergence behaviour of the classical approach. An implementation of the Schwarz waveform relaxation algorithm in a complex industrial case has shown the usefulness of the method. An extension to the new approaches, together with a study of how to reduce the overhead cost will certainly be fruitful.

Figure 6 . 1 :

 61 Figure 6.1: BET Isotherm with Q S = 2, K S = 0.7, K L = 100

Figure 6 . 2 :

 62 Figure 6.2: Variation of the error of the 10th iterates versus the parameter p of the Robin transmission condition. Nonlinear function. Left, BET isotherm law, Right, exponential equilibrium model.

Figure 6 . 3 :Figure 6 . 4 :

 6364 Figure 6.3: Concentration u and v at t = 0.5 using the BET isotherm with an incoming reactive front. Solution of the third domain decomposition iteration.

Figure 6 . 5 :

 65 Figure 6.5: Number of matrix inversions versus parameter p of the Robin transmission conditions for the classical approach (fixed point on the nonlinear interface problem), Nested Iteration Approach and Common Iteration Approach

Figure 6 . 6 :Figure 6 . 7 :

 6667 Figure 6.6: Number of matrix inversions versus number of discrete points per dimension (N x = N y ) for the classical approach (fixed point on the nonlinear interface problem), Nested Iteration Approach and Common Iteration Approach

Figure 6 . 8 :

 68 Figure 6.8: Convergence history with 200 points per space dimension for the Classical Approach, Nested Iteration Approach and Common Iteration Approach

3 3

 3 Precipitation of Calcite: Calcite + H + ←-Ca ++ + HCO - Dissolution of Quartz:Quartz -→ SiO 2 (aq) Precipitation of Quartz: Quartz ←-SiO 2 (aq)

Figure 7 . 1 :

 71 Figure 7.1: Choice of subdomains in the SHPCO2 2D test case. Orange: position of the initially present CO 2 . Green: Areas with low porosity. Shaded grey: reactive subdomain. Violet solid line: Interface of the reactive subdomain. Hatched annulus: overlapping region between reactive and non reactive subdomain.

Figure 7 . 2 :

 72 Figure 7.2: pH of the SHPCO2 test case at t = 95.13 years (end of the first time-window for the domain decomposition simulation). The low pH-value can be seen as an indicator of highly chemically reactive zones.

this work was supported by the research french agency ANR (project SHP-CO 2 ).


Juliette Rambourg

Stéphane Conversy
email: stephane.conversy@enac.frjrambourg@houstonmethodist.org

Garbey Marc
email: mgarbey2@houstonmethodist.orghelene.gaspard-boulinc@enac.fr

Hélène Gaspard-Boulinc

Collaboration within the Surgical Suite: BoardProbe Design For and With the Surgical Team

Keywords: Computerization, surgical team, collaboration, ATC, decisionaids. ACM Classification H.5.m. Information interfaces and presentation (e.g., HCI): Miscellaneous

The expectations of computerizing hospitals are high, in the belief that it will improve health care quality, reduce costs, and increase administrative efficiency. However, effective computerization of hospitals is a real challenge. The explanations of this failure are typical of HCI misdesign. In this paper, we focus on the computer support of the activity within the surgical suite. Our research objectives are to understand the local activity at the DUNN Operating Rooms surgical suite of the Houston Methodist Hospital, propose a technological solution to instrument the staff activity, and deliver a design method and a toolkit to adapt our solutions to other hospitals.

Introduction

The expectations of computerizing hospitals are high, in the belief that it will improve health care quality, reduce costs, and increase administrative efficiency [START_REF] Himmelstein | Hospital Computing and the Costs and Quality of Care: A National Study[END_REF]. Computerization is the act of using a computer to do something that was done by people or other machines before (Cambridge Dictionary). Electronic Patient Records are an example that has been implemented in order to improve cost-effectiveness, distribution, reliability, and exchange [START_REF] Kim | A solution to the distribution and standardization of multimedia medical data in E-Health[END_REF].

However, effective computerization of hospitals is a real challenge: Himmelstein et al. [START_REF] Himmelstein | Hospital Computing and the Costs and Quality of Care: A National Study[END_REF] estimate that so far computerization did not bring any evidence of improvement in administrative costs neither overall costs, or brought modest improvements in quality of patient outcomes. Black et al. [START_REF] Black | The Impact of eHealth on the Quality and Safety of Health Care: A Systematic Overview[END_REF] highlight the gap between the postulated and empirically demonstrated benefits of eHealth technologies, as well as the lack of best practice guidelines in effective development and deployment strategies. The explanations of this failure are typical of HCI misdesign. Multiple works highlight the fact that computerization brings benefits to a specific context if the system has been especially customized for this given context. For example, Himmelstein et al. [START_REF] Himmelstein | Hospital Computing and the Costs and Quality of Care: A National Study[END_REF] states that computerization cannot decrease costs because what is saved is compensated by the purchase and the maintenance of the computer system (presumably evolution maintenance of the system to adapt to context particularities). The failure also remains in the lack of pertinent evaluation of the computer system and should include socio-technical factors to adapt and fit a local activity as best as possible [START_REF] Black | The Impact of eHealth on the Quality and Safety of Health Care: A Systematic Overview[END_REF]. Evidence of the beneficial impact of such systems is limited to a few academic clinical centers of excellence where the systems were developed in house, undergoing extensive evaluation with continual improvement, supported by a strong sense of local ownership by their clinical users [START_REF] Black | The Impact of eHealth on the Quality and Safety of Health Care: A Systematic Overview[END_REF]. The contrast between the success of these systems and the relative failure of much of the wider body of evidence is "striking" [START_REF] Black | The Impact of eHealth on the Quality and Safety of Health Care: A Systematic Overview[END_REF]. Himmelstein et al. [START_REF] Himmelstein | Hospital Computing and the Costs and Quality of Care: A National Study[END_REF] mentions that the few custom-built systems that have been deployed have improved quality of patient outcomes. However, the extent to which the results of these primary studies on custom-built systems can be generalized beyond their local environment is questionable [START_REF] Black | The Impact of eHealth on the Quality and Safety of Health Care: A Systematic Overview[END_REF].

In this paper, we focus on the computer support of the activity within the "surgical suite". The surgical suite is a group of one or more operating rooms and adjunct facilities, such as a sterile storage area, scrub room, and recovery (The free dictionary, medical dictionary). There is an important demand of course for improving patient outcomes in surgery but also reducing its high costs. Therefore, we target the management of the surgical suite, in particular the management of patients, operations, operating rooms and hospital personnel around a shared whiteboard.

Our research objectives are: 1. Understand the local activity at the DUNN Operating Rooms surgical suite of the Houston Methodist Hospital: how does the staff exchange, shift, evolve and adapt to the patient flow and staff and equipment allocation? What are their habits and how do they deal with unforeseen events? What are their current organizational and collaborative tools and how do they deal with them? What is the team synergy? 2. Propose a technological solution to instrument the staff activity the best possible. What do they need to improve their efficiency, effectiveness and satisfaction? What do they need to improve patient outcomes? What do they need to be more cost-effective? How can we measure the efficiency of these solutions? 3. Deliver a design method and a toolkit to adapt our solutions to other hospitals and other surgical suites, to help generalize effective computerization of hospitals beyond a particular one.

In the next section we describe the related work. We then present the three methods that we follow: the Technology Probes, the Air Traffic Control analogy and the ethnographic study of the activity. We finally expose our first prototypes and our perspectives for future work.

Related work

Surgical suite studies Several ethnographic studied collaboration in the surgical suite [START_REF] Bardram | Activity-based computing for medical work in hospitals[END_REF], [START_REF] Bardram | A web of coordinative artifacts: collaborative work at a hospital ward[END_REF], [START_REF] Scupelli | Supporting coordination in surgical suites: physical aspects of common information spaces[END_REF], [START_REF] Tang | Supporting transitions in work: informing large display application design by understanding whiteboard use[END_REF]. Bardram [START_REF] Bardram | Activity-based computing for medical work in hospitals[END_REF] mentions that the problem of computerization at hospitals is mainly because "the systems are "typically single-user oriented", whereas teamwork is omnipresent. These studies describe activities within the surgical suite of other hospitals, sometimes around the same kind of whiteboards as in DUNN Operating Rooms (OR). They also provide recommendations for designing collaborative tools for the surgical suite. However, these studies did not develop collaborative tools themselves or test their design, which is missing in CSCW for hospitals [START_REF] Fitzpatrick | A Review of 25 Years of CSCW Research in Healthcare: Contributions, Challenges and Future Agendas[END_REF]. Fitzpatrick et al. [START_REF] Fitzpatrick | A Review of 25 Years of CSCW Research in Healthcare: Contributions, Challenges and Future Agendas[END_REF] also mention that, for the past 20 years, commercial offerings of electronic whiteboards have been a failure. They claim that the mistake comes mainly from the lack of usability, the weakness of the quality of the display and hardware.

As recommended by Bossen et al. [START_REF] Bossen | Collaboration in-between: The Care Hotel and Designing for Flexible Use[END_REF], we should consider "moving from design-for-use to design-for-future-use", in order to deploy pliant systems. The design choice of Bossen et al. [START_REF] Bossen | Collaboration in-between: The Care Hotel and Designing for Flexible Use[END_REF] to improve collaboration in the care hotel was an electronic whiteboard, a large device with good visibility. The whiteboard of the surgical suite is what Bossen et al. [START_REF] Bossen | Collaboration in-between: The Care Hotel and Designing for Flexible Use[END_REF] call the "common ground". This tool is essential to coordination and understanding between teammates. The situation awareness of the staff in the surgical suite relies partly on what is written on the whiteboard. According to Branham et al. [START_REF] Branham | Let's Go From the Whiteboard: supporting transitions in work through whiteboard capture and reuse[END_REF] whiteboards are poor tools for archival and reuse, but they are familiar, pliant and allow for quick creation and modification.

Air Traffic Control studies

Air Traffic Control and the surgical suite activity are both critical environments. They share specific requirements related to collaboration. Air Traffic Control research projects involved the development of new concepts that imitated and improved the current way of working of the staff by minimizing the difference with work habits [START_REF] Chatty | Pushing the limits of ATC user interface design to avoid S&M interaction: the DigiStrips Experience International conference on Air Traffic Management[END_REF], [START_REF] Conversy | Supporting Air Traffic Control Collaboration with a TableTop System[END_REF], [START_REF] Hurter | Strip'TIC: Exploring Automatic Paper Strip for Air Traffic Controllers[END_REF]. We will inspire from these studies to design our collaborative tools for the surgical suite.

Methods

Inquiries and Observations

To capture the local activity at DUNN OR, 1. We realized contextual inquiries and shadowing within the surgical suite, with the main actors. 2. And we set up a camera that can capture a whole whiteboard.

We have already started conducting contextual inquiries and direct observations within the surgical suite. We interviewed two Board Runners, two surgeons, two nurses, two Anesthesiologists and a Member of the DUNN OR surgical suite council. Because the staff is very busy, the Board Runners, the Surgeons and the Member of the council were interviewed out of the context of work. In a second part, Board Runners, Anesthesiologists and Nurses were interviewed and observed in their work environment.

With the camera, we can capture 1000 photos per day, taken with 1min delay between each photo or 30s. On the photos, we can see the symbols that are written/put on the whiteboard, the modifications/evolutions (updates, adds, writings, erasings, see Figure 4) brought to the whiteboard along the day, as well as the user-whiteboard (Figures 2 and Figure 5) and the user-user interactions (Figures 3 and4).

From these studies, we extracted 20 identified scenarios (an example Figure 6), we started to identify 25 first user requirements and we identified more than 30 different "symbols" used on the whiteboard (which is a local language).

Technology Probe

Our approach is to equip the surgical suite information hall with BoardProbe, a whiteboard-like multi-touch screen application, by following the Technology Probe design approach. Such a probe has 3 dimensions [START_REF] Boehner | How HCI interprets the probes[END_REF], [START_REF] Hutchinson | Technology probes: inspiring design for and with families[END_REF]:

1. The social science goal of collecting information about the use and the users of the technology in operational, track how users deal with the probe over time: we must understand the local activity of the staff at the surgical suite DUNN OR, Houston Methodist Hospital. 2. The engineering goal of testing the technology in operational: the technology has to be functional and allow them to perform their daily tasks such as managing the schedule, adding and communicating information through the day. How can we measure their efficiency through this probe? 3. The design goal of inspiring the staff and the design team to create adapted features answering the needs of the surgical team.

In order to capture the users' interaction with BoradProbe, we will use the webcam on the screen, automatic screenshots and logs of the input.

Air Traffic Control Analogy

An analogy can be drawn between the surgical suite computerization and Air Traffic Control computerization issues linked to critical environments (Figures 5, Figure 6 and Figure 7). In terms of social science, we will compare the collaborative activity of the surgical staff and their whiteboard with the collaborative activity of the controllers and their strips board. In terms of technology/engineering issues, we encounter the same problems due to the difficulty to set experiments and new tools for activities in operational context. This difficulty is one of the reasons why we chose to use a Technology Probe methodology: to allow the staff to actually use a very simple but functional device to gather information, test the technology and evolve the design.

Results of first observations

The surgical suite The central organizational and collaborative tool consists of two large whiteboards placed in the surgical suite. On Figure 8, "X" represents where the large whiteboards are. "A" represents the Post Operating Rooms (PostOp). There are also other rooms where the patient scheduled for surgery can be before or after the surgery: Holding Room -Holding, Pre Operating room -PreOp and PostOp. "B" represents the Coordinative Area. It is physically separated in 2 parts: the Information Hall with the whiteboards and the Control Room where necessary information about cases is dispatched. Information within the Control Room is available for individual consultation. Information that is displayed in the Information Hall has for purpose to be shown publicly -to the whole staff of the surgical suite. "C" represents the operating rooms, where the surgeries take place. DUNN OR contains about 23 operating rooms, which requires a very strict organization in terms of management of the resources -human and equipment.

The surgical team

The cases are updated throughout the day: new scheduling, shifts, patient location, etc. All necessary information about the organization of the case is contained in the whiteboard, which makes it a large source of information, allowing synchronous as well as asynchronous collaboration among the team, made of Board Runners, Nurses, Surgeons, Anesthesiologists and the Members of the DUNN OR surgical suite council (some examples Figure 9). The case schedule for the next day is edited late in the evening on the whiteboard, and starts to be used by the staff early in the morning (Figure 4). [Board Runner]The "name alert" magnet is used when two patients have the same name. Blue is usually for the anesthesia information and black is for the first edition of the schedule in the evening.

[Nurse] As soon as we have time: when a patient arrives in PreOp, we cross out the time of his case line so that everybody can know he is in here. We cross it then a second and third time, as a star, when the surgery has begun.

[Anesthesiologist] I do not really input information myself on the board, I read it to know what I have to do, this is more the other anesthesiologist Dr.

A. that manages whiteboard inputs.

[Surgeon] When I arrive the morning in the surgical suite, I like to see this big whiteboard with all the information I need for the day from where I am. I check that my patients and their surgery procedures are well written by the nurses on the whiteboard.

[Council member] We need objective data to understand the staff and know how efficient we are to improve patient outcomes and reduce costs of surgery.

Early Prototypes

Our guidelines So far, we determined the necessary "minimum" features to set up the BoardProbe within the surgical suite and observe the use (see Current features). The purpose is that the staff could be able to use the prototype as a replacement of the whiteboard. This set was progressively tested by the DUNN OR staff and people from the Computational Surgery laboratory, Houston Methodist Hospital.

Current features

We develop an application (Figures 1 and Figure 14), running on a 70 inches MultiTouch screen. Our purpose is to reach a maximum of flexibility and similarity compare to the original whiteboard. The global organization is exactly the same as the original whiteboard (Figures 5 and Figure 14). We first developed the basic feature "write-erase" that is mandatory for flexibility purpose and technology transition (Figures 12 and Figure 13). MAGNETS Magnets are items that can be dragged and dropped freely on a case -Name conflict alert, or allergy to latex, (Figure 10). ADDS-ON AND REORGANIZATION Add-ons (Figure 11) are additional cases that will be added during the day, as an emergency or not. Their insertion in the schedule is always a difficult decision and a "physically" complicated task to perform, due to the lack of space on the board.

Future work

The project started in January 2016; we began with ethnologic studies and design of prototypes, with users' feedback, observations and tests constantly in the loop. We are planning to set up the experiment in operational at DUNN OR within a couple months. Before that we need to make sure the application is robust enough. The application will be improved according to the user satisfaction, efficiency and effectiveness during the user experiments.

Figure 3 .

 3 Figure 3. Two Nurses working together on preparing the whiteboard the night before the next day.

Figure 2 .

 2 Figure 2. A Nurse, taking notes from the whiteboard on personal support.

Figure 4 .

 4 Figure 4. Staff walking in front of the whiteboard, sending a text to coworkers.

Figure 6 .

 6 Figure 6. Example of a scenario of Add-on arrival.

Figure 5 .

 5 Figure 5. Comparison Surgical Suite activity and Air Traffic Controllers. activity.

Figure 7 .

 7 Figure 7. Extract of the analogy tabs between Air Traffic Control and the Surgical Suite activity.

Figure 8 .Figure 9 .

 89 Figure 8. Reconstruction of a part of the surgical suite of DUNN OR, Houston Methodist Hospital. Juliette Rambourg, Vid Fikfak

Figure 13 .

 13 Figure 13. Erasing function with two fingers.

Figure 14 .

 14 Figure 14. Prototype overview: same columns and visual architecture as the whiteboard.

Figure 10 .

 10 Figure 10. Magnets feature.

Figure 12 .

 12 Figure 12. Writing feature.

Figure 11 .

 11 Figure 11. Adds-on feature.

IHM'16 Travaux en Cours

Acknowledgments

Preparation of this paper was supported by National Science Foundation under the I/UCRC for Cyber-Physical Systems for the Hospital Operating Room Grant No. 106022. We also thank the Center for Computational Surgery members and the DUNN OR surgical suite staff of the Houston Methodist Hospital.