Ahmed Kayad Moussa
email: ahmed.kayad_moussa@mines-paristech.fr

Dominique Bruel
email: dominique.bruel@mines-paristech.fr

Role of temperature change in micro seismic activity during fluid injections in faulted and fractured zones. Part 1: Updating the thermal modelling in a DFN model using a double media approach

Keywords: induced seismic, heat transfer, fracture system, geothermal system, hot dry rock,

Introduction

According to the general agreement that carbon emissions use should be strongly controlled, governments are now committed to diversify the primary fossil energy resource and particularly use more renewable resources. This is even more vital for countries that depend on imported oil products. Developing electricity production obtained from geothermal fields is therefore a priority at many places in the world. A step forward in this domain has been achieved with the Enhanced Geothermal System (EGS) concept, which was demonstrated at a pilot plant scale at the Soultz sous Forêts site [START_REF] Genter | Overview of the current activity of the European EGS SOultz project: from exploration to electricity production[END_REF]. In such reservoirs, developed in deep hard rocks, the overall flow pattern is controlled by the pre-existing network of geological discontinuities. Access to the reservoir can be improved thanks to hydraulic stimulation techniques resulting in better connections between boreholes and these natural pathways. Controlling the efficiency of stimulation techniques upon the natural structures forming the bulk part of the reservoir is still a matter of research, as results obtained at various sites are highly variable. In all cases, raising the pore pressure, thereby reducing the shear strength of fractures zones and promoting shear failure and some increase in fracture permeability, is accompanied by micro seismic activity and acoustic emissions (AE). These signals are recorded and analysed. The spatial localization of rupture sources often delineates planar structures where pressure perturbations were propagated. It is therefore confirmed that flow occurs in a 3D network of 2D structures and that specific quantitative models are required to capture the coupled hydraulic-mechanical processes at issue. During the past decade, attention was first given to the engineering of the reservoir using various combinations of hydraulic tests, and to the evaluation of the hydraulic improvement of the pathways developed the fractured system, using tracer tests [START_REF] Sanjuan | Tracer testing of the geothermal heat exchanger at Soultz-sous-Forêts (France) between 2000 and 2005[END_REF]. Recent works are re-focussing at physical processes, such as the determination of fracture/matrix transfer area, (K. [START_REF] Pruess | Thermal Single-Well Injection-Withdrawal Tracer Tests for Determining Fracture-Matrix Heat Transfer Area[END_REF] or at comparisons of tracer and thermal transport in fractured reservoirs (Juliusson et Horne, 2010) to investigate the extent to which tracer return can be used to predict thermal breakthrough. Very few new programmes were dedicated to long term thermal behaviour, assuming that a reservoir with sufficient size would last long enough for industrial purposes. The phenomenon of large induced seismic events (LISE) then received considerable interest after the generation in 2003 of a Magnitude 2.9 earthquake at Soultz sous Forêts (France). In 2006, a M 3.4 earthquake was generated during the hydraulic development of a new EGS reservoir, near the Swiss city of Basel, that stopped the project. This issue of large seismic event is also raised at other places where natural steam reservoirs exist and are exploited for decades. At the Geysers site (California) it is clearly observed that earthquakes with magnitudes larger than 4 become more frequent, in correspondence with continuous re-injection of fluid performed to maintain the reservoir pressure. In this case, it was argued and demonstrated (Rudqvist et al., 2002) ° that the thermal poro-elastic cooling below re-injection zone is responsible for mechanical changes, sufficient to activate dislocation mechanisms along pre-existing fractures. However the working hypothesis of equivalent porous media did not offer the possibility to the authors of going further ahead with fracture failure analysis in a realistic fractured reservoir. Other pioneer 3D modelling works can be mentioned. Conceptual models based on deterministic regular fracture patterns or simplified geometries (Swensson 1990[START_REF] Hicks | A hydrothermo-mechanical numerical model for HDR geothermal reservoir evaluation[END_REF][START_REF] Kolditz | Quasi-steady-state strategy for numerical simulation of geothermal circulation in hot dry rock fractures[END_REF][START_REF] Kolditz | Numerical simulation of flow and heat transfer in fractured crystalline rocks: Application to the hot dry rock site in rosemanowes (U.K.)[END_REF] were first developed. Since fractures are unevenly distributed in space, forming clusters, with widely variable size (see [START_REF] Genter | Over ten years of geological investigations within the HDR soultz project[END_REF] at Soultz sous Forêts site ), and since the heat to mine out is stored in the matrix surrounding the fractures, continuum models based on equivalent porous materials are not adequate. Randomly distributed networks [START_REF] Robinson | A A discrete fracture model for a hot dry rock geothermal reservoir[END_REF][START_REF] Bruel | Heat extraction modelling from forced fluid flow through stimulated fractured rock masses. Application to the Rosemanoves Hot Dry Rock reservoir[END_REF][START_REF] Bruel | Impact of induced thermal stresses during circulation tests in an engineered fractured geothermal reservoir. Example of the Soultz sous Forêts, European hot fractured rock geothermal project[END_REF]Willis-Richard et al., 1996) embedded in the rock mass have been developed. Most of the recent efforts are attempts to mix both representations, a discrete randomly distributed network embedded in a porous system. A recent example is given by (E. Juliusson and R.N. Horne, 2010) using Finite Element type numerical methods, and working on sophisticated unstructured computational grids. One of their conclusions is however that numeric for a realistic (2D) network was at the margin of being feasible, and the authors claim for the need of new appropriate upscaling algorithms and new non meshing techniques.

The purpose of this paper is to go back to the understanding of the long term potential impacts of coupled Thermo-Hydro-Mechanical processes (THM) in EGS geological settlements, using a discrete fracture network (DFN) specific approach. This discrete approach is able to handle the specific geometrical characteristics of fracture networks and faulted zone. Such an approach was already introduced for EGS modelling by [START_REF] Bruel | Heat extraction modelling from forced fluid flow through stimulated fractured rock masses. Application to the Rosemanoves Hot Dry Rock reservoir[END_REF][START_REF] Bruel | Impact of induced thermal stresses during circulation tests in an engineered fractured geothermal reservoir. Example of the Soultz sous Forêts, European hot fractured rock geothermal project[END_REF] for the study of a 3 year long circulation test at the Rosemanoves site (UK) and for the simulation of a four months long circulation tests at 3.5 km depth performed in 1997 in the upper reservoir of the Soultz sous Forets site (France). For the latter case, no long term fluid circulation was planned for this EGS sites, so that only local thermo-elastic couplings were considered, with no interactions in between temperature changes and fluid properties.

The work presented here after is an update of the FRACAS DFN code [START_REF] Baujard | Numerical study of the impact of fluid density on the pressure distribution and stimulated volume in the Soultz HDR reservoir Geothermics[END_REF] where a double media is introduced to solve for thermal conductive heat transfer at the reservoir scale coupled with local heat exchange along the fractures embedded in the rockmass and heat convection along the flowing pathways. These efforts are clearly required to tackle new research aspects concerning long term seismic risk analysis. This should be also particularly useful for new sites as those planned in faulted rift zones of eastern Africa, where the prevailing tectonic regime makes fractures very close to instability. As an example, C. Doubre et G. Peltzer (2007) analyse natural micro-seismic activity and show that present day fracture propagation at Asal rift (Djibouti) is probably linked with natural flow along fractures. The Asal hydrothermal region has been identified as a favorable area to develop the geothermal resources and projects are already underway. But the first evaluation of the production potential considers that re-injection of (cold) fluids will be necessary to compensate for reservoir drawdown [START_REF] Daher | Geothermal Resource Assessment of Asal Field[END_REF].

Theoritical background of the FRACAS Discrete Fracture Network flow code

The FRACAS modelling approach is based on the assumption that fluid moves through a rock mass within a system of interconnected fractures and that flow in the rock matrix is negligible by comparison. The three-dimensional hydraulically conductive network of planar, disc-shaped fractures generation is inherited from [START_REF] Cacas | Modeling fracture flow with a discrete fracture network ; calibration and validation. 1. the flow model[END_REF] modelling work. Series of individual fractures are generated within a rectangular block of rock, based on stochastic descriptions of fracture density, fracture orientation, and fracture size. Improvements were achieved to capture more realistic geological systems. Specific fracture sets can be defined as models of faults zones. Relay structures or simple planar fault segments can be generated as well as combinations of adjacent blocks with different fracture network properties.

Flow rules

Hydraulic conductivities

The overall flow model is based on the premise that flow in granitic basements primarily occurs on channels within fractures. The 3D structure resulting of the superimposition of the connected portions of planes is thus treated as a 3D network of 1D pipes. Equivalence of this channel flow approach with others solving true 2D flow in planes has been discussed (Dershowitz and Fidelibus, 1999). Time-dependant analysis requires assumptions to be made concerning the form of fluid flow within the fracture network. The general form of fluid flow assumed in each fracture is based on an analytical solution, known as the « cubic law », for fluid flow between approximately parallel surfaces. Modified forms of the cubic law are used to account for the effects of changes in the morphology of contact between fracture surfaces. As the effective stress across a fracture increases flow through the fracture becomes confined to a limited number of channels. Thus, in FRACAS, the single phase volumetric flux (m 3 s -1 ) in the x-direction through a length l (m) of a fracture has the form :

dx dh glF a Q ν 12 3 0 - = ( 1)
where a 0 is the hydraulic aperture (m) of the fracture at zero effective stress, g is acceleration due to gravity (ms -2 ), ν is the kinematic viscosity (m 2 s -1 ), dh/dx is the hydraulic head gradient driving flow through the fracture, and F is a dimensionless function dependent on effective stress. The effective stress is defined as σ'=σ π -p, where σ π (Pa) is the rock stress normal to the fracture surface and p (Pa) is the fluid pressure in the fracture. Empirical expressions for the closure law F, in which F decreases as the effective stress increases (F=1 at zero effective stress) are presented by [START_REF] Jeong | Numerical experiments of flow and transport in a variable-aperture fracture subject to effective normal stresses[END_REF]. This coupled approach was tested against various in situ experiments performed in hard rocks in the frame work of a nuclear waste storage research program, and predictive results were compared with those of other numerical approaches [START_REF] Rejeb | Hydromechanical effects of shaft sinking at the Sellafield site[END_REF]. The approach was extended [START_REF] Bruel | Using the migration of the induced seismicity as a constraint for fractured hot dry rock reservoir modelling[END_REF] to account for shear loading and Coulomb frictional failure, in an attempt to simulate reservoir hydraulic stimulation processes and micro-seismic predictions.

However, in this paper the objective is to incorporate the potential effects of a short term thermal perturbation to evaluate how significantly it may alter the reservoir performances. Therefore F will be re-derived later, according a new procedure described in the section devoted to the thermo-hydro-mechanical interactions, after the basics for heat exchange modeling in the FRACAS discrete fracture network model have been introduced.

Specific storage of the fractures

The specific storage S (m -1 ) of a fracture with R as radius (m) is given by:

0 0 ' 1 S d da a g S + = σ ρ ( 2)
where ρ is the fluid density (kgm -3 ) and g the acceleration due to gravity (ms -2 ). S 0 denotes the potential specific storage of the adjacent weathered rock that can be identified on cores and in situ, using sonic logs. These zones exist on both fracture sides and can extend up to a metric scale [START_REF] Genter | Over ten years of geological investigations within the HDR soultz project[END_REF]. In deriving the specific storage term, the rock stress has been assumed constant, such that changes in effective stress are caused only by fluid pressure changes. Also, fluid compressibility is generally much smaller than fracture compressibility and, thus, has been neglected. The form the of the term da i /dσ' calculated at the center of the fracture disc, depends on the assumed form of the relationship between the hydraulic fracture aperture and the effective stress, as discussed above.

Thermal modelling

Convective and conductive heat transfer along individual fractures

Heat exchange along the circulated fractures can be evaluated assuming that heat conduction develops perpendicular to each fracture plane. No energy is retained by the volume of fluid within the fractures. The conservation of energy is written at each fracture centre as a balance between energy transferred by convection from or to neighbouring fractures and energy dissipated by conduction in the current cell. At any time, continuity of the temperature at the fracture wall is assumed, that is :

θ f = θ m (y=0) ( 3)
Subscript f and m respectively stand for fluid and matrix. If q ij and q ik [m 3 s -1 ] denotes fluxes entering and leaving a given disc i at time t, with respective temperatures θ fi (t) and θ fi (t+dt), the energy exchanged during dt is dQ i [J], given by :

dQ i = Σ j ρ f C f q ij θ fj dt -Σ k ρ f C f q ik θ fi dt ( 4)
This energy is related to the heat flux at the fracture wall, Φ i , across the exchange area S i [m 2 ], that has to dissipate in an adjacent rock volume according to a diffusive equation. The shape of this adjacent block is cylindrical, with a radius equal to that of the fracture-disc and a length l i chosen so that the total volume cumulating all the cylinders is equal to the reservoir volume in order to fulfil the global heat capacity conservation.

The temperature at the face opposite to the circulated fracture is a new variable. θ ml = θ m (y=l).

This initial value of this local variable is set to the temperature estimated at depth for the corresponding fracture, given the natural measured temperature profile.. The developments insuring block to block thermal interactions for predicting long term thermal depletions are proposed in the next section. We will assume that interactions are controlled by the thermal diffusivity coefficient, which is in the order of 1. m 2 s -1 for rocks. Temperature changes may propagate at significant distances greater than 10 m within time periods no longer than a year. The equations presented here below are solved sequentially.

dQ i = Φ i S i dt,
(5)

Φ i = K m (dθ m /dy) y=0 (6) t C K m m m m m ∂ ∂ = ∆ θ θ ρ (7)
At each disc i in the network, the diffusivity equation is written in a discrete form following a standard finite difference scheme, and the set of equations is solved for θ fi and θ mi at time t+dt knowing similar quantities and boundary θ mli values at time t. The temperature of the injected fluid is continuously prescribed. Parameters are matrix heat conductivity K m [W m -1 K -1 ], heat capacity C f , C m [J kg -1 K -1 ] for both fluid and matrix, density ρ f , ρ m [kg m -3 ] for both fluid and matrix.

Conductive heat transfer at the global scale and derivation of internal boundary values θ θ θ θ ml

The basic idea suggested in this section is that the heat extracted during a given time step from a sub volume containing a number of connected fractures by a fluid forced to circulate along these fractures can be approximated by an equivalent density term, source or sink, depending on the sign of the calculated heat balance. Therefore, the method consists at any time step, when fluxes, fluid temperature and local matrix temperature profiles are known at each fracture disk i,, in solving a succession of tasks as follows:

• Define a regular 3D grid, with cells cubic in shape, overlapping the fractured reservoir, each cell being large enough to contain a sub-network of fractures • Identify the discs belonging to each cell • Calculate the net heat balance per cell, from the differences in temperature profiles at t and t+dt at each disc within the cell • Derive a source term from the temperature change for each cell, accounting for cell volume and time step duration. • Solve for a heat conduction problem over the 3D grid, using prescribed temperature at the top of the reservoir, regional geothermal heat flux at the base of the reservoir, heat conduction and capacity coefficients for the rock mass, and the discrete volumic heat source at each cell.

• Use the temperature obtained at each cell at date t as the new θ θ θ θ ml value for all the discs belonging to that cell.

The two step algorithm is sequential, but the time stepping for the conductive behaviour at large scale can be larger than the time stepping for solving the local convective-conductive heat exchange at the local scale. This new algorithm has been iùplemented in our FRACAS code. Numemerical verifications are now focussing on the hydro-thermal global behaviour of the double media, using benchmarking against other numerical codes that can handle some fracture patterns (GEOFRACK, TOUGH2). An outcome of the FRACAS code development is that experimental testing programmes running over years but combining successions of injection/production phases separated by month long rest periods can be modelled in a single run.

Preliminary conclusion

Apart from standard numerical checking, a series of verification examples against true data sets is under construction. The first one will be based on results obtained on one of the first Hot Dry Rock research programme run at the Rosemanoves quarry site (Cornwall, UK), extensively documented by RH Parker (1988).

A main aspect in this data base is that it contains a long term circulating test that lasted about 3 years and that shows a thermal breakthrough partly due to thermo mechanical interactions. Fracture network is well described, as a dense regular network made of 3 directional sets (blocky pattern). A large number of hydraulic tests in between the different wells have been performed as well as tracer tests to discuss heterogeneity and role of main flow paths. A main interest in running a model in such a system is to revisit the potential impact of fracture stimulation at the network scale and rediscuss some of the remarks concerning single well injection-withdrawal tracer tests done by [START_REF] Pruess | Thermal Single-Well Injection-Withdrawal Tracer Tests for Determining Fracture-Matrix Heat Transfer Area[END_REF], stating that temperature return is insensitive to fracture aperture changes. As thermal processes may combine with mechanical processes, we also intend to derive rock thermal shrinkage at the reservoir scale, and obtain the dramatic short-circuit effect observed at Rosemanoves site. Given shear rupture modelled along thermo-mechanical perturbed pathways, we will suggest, as done in [START_REF] Bruel | Using the migration of the induced seismicity as a constraint for fractured hot dry rock reservoir modelling[END_REF], a catalogue for the thermally induced micro-seismicity.


Optimized Sparse Matrix Formats on GT200 and Fermi GPUs

List of authors: M. Hugues 1 and S. Petiton CNRS/LIFL, Laboratoire d'Informatique Fondamentale de Lille, France maxime.hugues@lifl.fr, serge.petiton@lifl.fr 2 The rise of computing power now goes through the multiplication of cores on a single-chip because of the frequencies limitation. Multicore processors will compose future peta and exascale supercomputers which will allow applications to sustain many petaflops and thus to make breakthrough on scientific challenge problems. One massively multicore processor already available on the market is the GPU and has been making the rebirth of data parallel programming which is focused on data mapping and movements. GPU is a high peak performance architecture which gives the opportunity to accelerate scientific applications. Dense linear algebra has already shown that it could benefit from this architecture [START_REF] Volkov | Benchmarking GPUs to Tune Dense Linear Algebra[END_REF] and some libraries such as CUFFT or CUBLAS are available to easily develop high performance applications. In opposite, sparse matrices have irregular data structures which generate a memory bottleneck by accessing data in a non-contiguous fashion. The main optimization strategy is to get the most efficient memory accesses by arranging data in a well-suited manner depending on the matrix structure. Sparse formats define this data arrangement in memory. The sparse matrix vector product (SpMV or Ax = y) is one of the most important kernel in scientific computing. This kernel is intensively used and is the main performance key of iterative methods, such as Conjugate Gradient or GMRES, for solving partial differential equations, large linear systems or eigenvalues problems. Thereby, many researches have been conducted to optimize this kernel implementation through various sparse formats on many-core processors. On multicore processors, the Block Compressed Sparse Row (BCSR) format has been evaluated and optimized on different platforms like Cell BE, Intel, AMD and Sun processors [START_REF] Williams | Optimization of Sparse Matrix-Vector Multiplication on Emerging Multicore Platforms[END_REF]. It clearly achieves the best performances by reusing sparse matrix elements in the different cache levels thanks to the register and cache blocking techniques. On massively parallel architectures, some implementations of SpMV kernel have been evaluated by using various sparse formats on NVIDIA GPUs [START_REF] Bell | Implementing Sparse Matrix-Vector Multiplication on Throughput-Oriented Processors[END_REF].

They have shown that their HYB Format, a mix between the COO and Ellpack/ITpack format, achieves the best results on a set of matrices previously used on multicore processors and the DIA and Ellpack/ITpack is the most efficient format for a finite differences discretization of the Laplacian. In [START_REF] Baskaran | Optimizing Sparse Matrix-Vector Multiplication on GPUs[END_REF], authors have also proposed an optimized implementation of the Compressed Sparse Row (CSR) format by using vectorization and padding which increase the efficiency and outperform the implementation of Nvidia. A sparse linear solver [START_REF] Buatois | Concurrent Number Cruncher: a GPU Implementation of a General Sparse Linear Solver[END_REF] method has been developed on GPU and has proposed the implementation of the SpMV by using the BCSR format for register blocking. In spite the fact, GPU is called data parallel architecture, like Connection Machine, the past results obtained on this type of architecture are not fully reusable as is.

1 2 1
We propose in this talk to study some implementations on GPU of sparse format which has been successful on data parallel machine (CSC, SGP [START_REF] Ferng | Basic Sparse Computations on Data Parallel Computers[END_REF]) and present a performance comparison of common sparse formats for the SpMV kernel with our optimizations. This work also compares the performance achievement on the GT200 and the new Nvidia architecture, Fermi, through two evaluation methodologies The first method is based on a set of matrices from the University of Florida Sparse Matrix Collection and a finite difference discretization of the Laplace operator in 2D with different stencil in order to compare our results to previous of [START_REF] Bell | Implementing Sparse Matrix-Vector Multiplication on Throughput-Oriented Processors[END_REF]. The second one is focused on the sparse matrix format performances following the distribution of the non zeros values in the matrix by using a matrix pattern, named C-Diagonal q-perturbed. These patterns was introduced to evaluate sparse format performances on data parallel supercomputers following communications which depends on the distribution of nnz [START_REF] Petiton | Massively Parallel Preconditioners for the Sparse Conjugate Gradient Method[END_REF]. The C-Diagonal matrix is defined as a matrix with nc diagonals from the main diagonal to the right. The C-Diagonal q-perturbed matrix is built from a C-Diagonal matrix by modifying the data structure with a random perturbation. Each element of the matrix has a perturbation called q which defines whether the value has to be set at a random distance from the main diagonal, thus the element column index is modified.

On the set of matrices and the GT200, we demonstrate Ellpack/ITpack format with a columnmajor order in memory is the best format for common matrices by achieving 18.7Gflop/s in single precision and 11Gflop/s in double precision. Our variant of Ellpack/ITpack with a row-major order performs also well, depending on the matrix structure. For matrices with a nnz per row lower than 16, this implementation and all vectorized formats, such as the CSR implementation, perform poorly because they can not realize a coalescing with a half-warp. The proposed vectorized version of BCSR format achieves also good performance on matrices with a sufficient nnz to have an optimized blocking. Results on compressed column formats such as CSC, Ellpack/ITpack with column compression and SGP show that type of format perform less better than compressed row format. This comes from the summing of all elements belonging to the same row which makes irregular data access and degrades performances. The performance study of sparse formats following the distribution of nnz also shows interesting results about the choice of format for different matrix structure. Ellpack/ITpack implementation with columnmajor order stays the most efficient sparse format by performing 21Gflop/s in single precision and 13.8Gflop/s in double precision with q = 0. Nevertheless, we have seen for sparse matrices with a strong distribution of nnz (q > 0.5) that the various formats achieve approximately the same performance by reaching 2.0 Gflop/s in single and double precision. This can be interesting for iterative methods which require high precision to keep stability and converge faster without loosing performance. On Fermi, the performance behavior observed is the same with an improvement for some sparse formats such as CSR and BCSR. The first one benefits of a 1.7 SpeedUp in comparison to the GT200 by achieving 17Gflop/s and the second one realizes 25Gflop/s with a SpeedUp of 2. BCSR outperform Ellpack/ITpack format on Fermi in single and double precision by using the new cache levels on the GPU. The C-diagonal q-perturbed benchmark shows that Ellpack/ITpack is still the best because as the distribution of nnz grows the number of Block for BCSR also grows and creates extra zeros to compute. Moreover, we can notice the performance are quiet the same but they decline faster when q is growing than

on the GT200 architecture and finally reach less figures for q = 1.0 in single precision. In double precision, we observe a SpeedUp of 1.5 when the distribution of nnz increases and demonstrates the effort of Nvidia on Fermi to optimize double precision computations is clearly significant.


Stephane Attal

Francesco Petruccione

Christophe Sabot

Ilya Sinayskiy

Open Quantum Random Walks

come

Introduction

Nowadays quantum random walks, such as the Hadamard quantum random walk, are quite a successful domain of research, with important applications in Quantum Information Theory (see [Kem] for a survey). These quantum random walks are particular discrete-time quantum dynamics on a state space of the form H ⊗ C Z d . The space C Z d stands for a state space labelled by a net Z d , while the space H stands for the degrees of freedom given on each point of the net. The quantum evolution concerns pure states of the system which are of the form

|Ψ = i∈Z d |ϕ i ⊗ |i .
After one step of the dynamics, this state is transformed into another pure state,

|Ψ = i∈Z d |ϕ i ⊗ |i .
Each of these two states gives rise to a probability distribution on Z d , the one we would obtain by measuring the position on C Z d :

Prob({i}) = ϕ i 2 .
The point is that the probability distribution associated to |Ψ cannot be deduced from the distribution associated to |Ψ by "classical rules", that is, there is no classical probabilistic model (such as a Markov transition kernel, or else) which gives the distribution of |Ψ in terms of the one of |Ψ . One needs to know the whole state |Ψ in order to compute the distribution of |Ψ . These quantum random walks, that we shall call Unitary Quantum Random Walks (for a reason which will appear clear in Section 10) have been successful for they give rise to strange behaviors of the probability distribution as time goes to infinity. In particular one can prove that they satisfy a rather surprising Central Limit Theorem whose speed is n, instead of √ n as usually, and the limit distribution is not Gaussian, but more like functions of the form (see [Kon])

x → √ 1 -a 2 (1 -λx) π (1 -x 2 ) √ a 2 -x 2 ,
where a is a constant.

The purpose of this article is to introduce a new type of quantum random walks, that we suggest to call Open Quantum Random Walks (O.Q.R.W.).

These quantum random walks also deal with a quantum dynamics on a state space H ⊗ C Z d , but they consider the evolution of density matrices

ρ = i∈Z d ρ i ⊗ |i i| .
More or less, the principle is the same as above, and the dynamics leads to a new density matrix ρ =

i∈Z d ρ i ⊗ |i i| .
To each of them is associated the probability distribution obtained when measuring the position

Prob({i}) = Tr (ρ i ) , i ∈ Z d .
This new type of quantum random walks is very different from the Unitary Quantum Random Walks. It seems that there is no inclusion whatsoever, though we prove in Section 10 a very strong link between the two walks, in the way they can be physically implemented.

Actually, the limit behavior of Open Quantum Random Walks shows up a dissipative character, it tends to converge to a classical behavior, that is, it seems to give rise to classical Central Limit Theorems: one can see the distribution converging to Gaussian limits, or to mixtures of Gaussian limits.

The point to be stressed is the generality of our setup. It allows to consider a very wide class of quantum random walks on nets as well as on graphs. Our setup is the exact quantum generalization of the construction of a classical Markov chain on a net, or on a graph. By the way, we shall show that Open Quantum Random Walks contain all the classical Markov chains as particular cases.

Our conviction is that this type of quantum random walks gives rise to a vast field of exploration for the behavior of open quantum systems. It may be as rich as the one of classical Markov chains and it shall give rise to the same type of questions: existence of invariant states, ergodic behavior, Central Limit Theorems, Large Deviation Principle, recurrence and transience, etc.

Many of the examples that we have explored lead us to think that these quantum random random walks may apply in many realistic physical situations. Their dissipative behavior makes them physically more realistic, while keeping a very quantum behavior. For example, some of the examples that we shall explore in this article make us think of possible applications, such as heat conduction and quantum Fourier's law for a one dimensional model (such as the quantum version of the "Simple Exclusion Process", see [Bod]) and realistic model for excitation transport on a chain of quantum systems.

Note that the main physical implications of this article have already been announced and summarized in a letter [APSS].

Also, it has to be said that the idea of considering matrices of completely positive maps such that the lines (or columns, depending on the point of view) form a so-called quantum operation, appeared earlier in [Gud]. This approach is presented as a "quantum Markov chain". These objects present clearly several common points in their structure with our Open Quantum Random Walks, but they are not studied as giving rise to quantum random walks. Except at the end of the article where an incorrect parallel with Unitary Quantum Random Walk is claimed.

General Setup

We now introduce the general mathematical and physical setup of the Open Quantum Random Walks. For sake of completeness we recall in this section several technical lemmas which ensure that our definitions are consistent. The proofs of these lemmas are postponed to Section 13. We are given a set V of vertices, which might be finite or countable infinite. We consider all the oriented edges {(i, j) ; i, j ∈ V}. We wish to give a quantum analogue of a random walk on the associated graph (or lattice).

We consider the space K = C V , that is, the state space of a quantum system with as many degrees of freedom as the number of vertices; when V is infinite countable we put K to be any separable Hilbert space with an orthonormal basis indexed by V. We fix an orthonormal basis of K which we shall denote by (|i ) i∈V .

Let H be a separable Hilbert space; it stands for the space of degrees of freedom (or chirality as they call it in Quantum Information Theory) given at each point of V. Consider the space H ⊗ K.

For each edge (i, j) we are given a bounded operator B i j on H. This operator stands for the effect of passing from j to i. We assume that, for each

j i B i j * B i j = I , (1)
where the above series is strongly convergent (if infinite). This constraint has to be understood as follows: "the sum of all the effects leaving the site j is I ". It is the same idea as the one for transition matrices associated to Markov chains: "the sum of the probabilities leaving a site j is 1". By Lemma 2.1 which follows, to each j ∈ V is associated a completely positive map on the density matrices of H:

M j (ρ) = i B i j ρB i j * .
Lemma 2.1 Let (B i ) be a sequence of bounded operators on a separable Hilbert space H such that the series i B * i B i converges strongly to a bounded operator T . If ρ is a positive trace-class operator on H then the series

i B i ρB * i is trace-norm convergent and Tr i B i ρB * i = Tr (ρT ) .
The operators B i j act on H only, we dilate them as operators on H ⊗ K by putting M i j = B i j ⊗ |i j| . The operator M i j encodes exactly the idea that while passing from |j to |i on the lattice, the effect is the operator B i j on H. By Lemma 2.2, which follows, the series i,j M i j * M i j converges strongly to the operator I.

Lemma 2.2 Let K and H be separable Hilbert spaces. Consider an orthonormal basis (|i ) of K. Assume that B i j are bounded operators on H such that, for all j, the series i B i j * B i j is strongly convergent to I. Define the bounded operators

M i j = B i j ⊗ |i j| on H ⊗ K. Then the series i,j M i j * M i j converges strongly to I.
As a consequence we can apply Lemma 2.1 to the set of operators (M i j ) i,j and the mapping

M(ρ) = i j M i j ρ M i j * (2)
defines a completely positive map on H ⊗ K.

We shall especially be interested in density matrices on H ⊗ K with the particular form

ρ = i ρ i ⊗ |i i| , (3)
where each ρ i is not exactly a density matrix on H: it is a positive and trace-class operator but its trace is not 1. Indeed the condition that ρ is a state aims to

i Tr (ρ i ) = 1 . ( 4
)
The importance of those density matrices is justified by the following.

Proposition 2.3 Whatever is the initial state ρ on H⊗K, the density matrix M(ρ) is of the form (3).

Before proving this proposition, let us recall a basic result on partial traces. We can now come back to the proof of the proposition.

Proof [of Proposition 2.3] We have M(ρ) = i,j B i j ⊗ |i j| ρ B i j * ⊗ |j i| = i,j (B i j ⊗ I)(I ⊗ |i j|) ρ (I ⊗ |j i|)(B i j * ⊗ I) .
If we put ρ j = j| ρ |j (as in Lemma 2.4), we get

M(ρ) = i,j (B i j ⊗ I)(ρ j ⊗ |i i|)(B i j * ⊗ I) = i,j B i j ρ j B i j * ⊗ |i i| .
Each of the operators B i j ρ j B i j * is positive and trace-class, hence so is the

operator j≤M B i j ρ j B i j * . But we have Tr j≤M B i j ρ j B i j * = j≤M Tr (ρ j B i j * B i j ) .
As i B i j * B i j = I, each of the operators B i j * B i j is smaller than I (in the sense that I -B i j * B i j is a positive operator). Hence, Tr (ρ j B i j * B i j ) ≤ Tr (ρ j ), as can be easily checked. This shows that

j Tr B i j ρ j B i j * < ∞
and that j≤M B i j ρ j B i j * converges in trace-norm to a positive trace-class

operator j B i j ρ j B i j * which satisfies Tr j B i j ρ j B i j * = j Tr (ρ j B i j * B i j ) .
In particular by Lemma 2.1 we have,

i Tr j B i j ρ j B i j * = i j Tr (ρ j B i j * B i j ) = j Tr ρ j i B i j * B i j = j Tr (ρ j ) = 1 .
This means that the series (in the variable i)

i j B i j ρ j B i j * ⊗ |i i|
is trace-norm convergent. We now immediately have the relation

i,j B i j ρ j B i j * ⊗ |i i| = i j B i j ρ j B i j * ⊗ |i i| .
This proves that M(ρ) is of the form (3).

The states of the form (3) are mixtures of initial states ρ i on each site i, but they express no mixing between the sites. An immediate consequence of the proof of Proposition 2.3 is the following important formula.

Corollary 2.5 If ρ is a state on H ⊗ K of the form ρ = i ρ i ⊗ |i i| , then M(ρ) = i j B i j ρ j B i j * ⊗ |i i| . (5)
This is exactly the quantum analogue of a usual random walk: after one step, on the site i we have all the contributions from those pieces of the state which have travelled from j to i.

Open Quantum Random Walks

If the state of the system H ⊗ K is of the form

ρ = i ρ i ⊗ |i i| ,
then a measurement of the "position" in K, that is, a measurement along the orthonormal basis (|i ) i∈V , gives the value |i with probability Tr (ρ i ) .

As proved in Corollary 2.5, after applying the completely positive map M the state of the system H ⊗ K is

M(ρ) = i j B i j ρ j B i j * ⊗ |i i| .
Hence a measurement of the position in K would give that each site i is occupied with probability

j Tr B i j ρ j B i j * . (6)
For each fixed choice of the initial states ρ i , the probability to find the particle on the site i after one step is the usual one for the usual random walk where the probabilities to go from j to i are Tr B i j ρ j B i j * . Now, let us see what would have happened if the measurement happens after two steps only. The state of the system is

M 2 (ρ) = i j k B i j B j k ρ k B j k * B i j * ⊗ |i i| .
Hence measuring the position, we get the site |i with probability

j k Tr B i j B j k ρ k B j k * B i j * . ( 7
)
Had we repeated the rule for the first step, as for a classical random walk, we would have found

j k Tr B i j ρ j B i j * Tr B j k ρ k B j k * ,
which is clearly different from (7), in general. Actually there is no way to understand the probability measure given in (7) for two steps with the help of only the probability measure on one step (6).

The random walk which is described this way by the iteration of the completely positive map M is not a classical random walk, it is a quantum random walk. The rules for jumping from a site to another are dictated by the sites, but also by the chirality. This is what we call an "Open Quantum Random Walk".

Let us resume these remarks in the following proposition, which follows easily from the previous results and remarks.

Proposition 3.1 Given any initial state ρ (0) on H ⊗ K, then for all n ≥ 1 the states ρ

(n) = M n (ρ (0) ) are all of the form ρ (n) = i ρ (n) i ⊗ |i i| .
They are given inductively by the following relation:

ρ (n+1) i = j B i j ρ (n) j B i j * .
For each n ≥ 1, the quantities

p (n) i = Tr (ρ (n) i ) , i ∈ V define a probability distribution p (n) on V, it is called the " probability distri-
bution of the open quantum random walk at time n".

Examples

Before going ahead with the properties of the Open Quantum Random Walks, let us introduce a few examples that we shall discuss in more detail in Sections 11 and 12.

An Example on Z

First of all, it is very easy to define a stationary open quantum random walk on Z. Let H be any Hilbert space and B, C be two bounded operators on H such that

B * B + C * C = I .
Then we can define an open quantum random walk on Z by saying that one can only jump to nearest neighbors: a jump to the left is given by B and a jump to the right is given by C.

In other words, we put

B i-1 i = B and B i+1 i = C
for all i ∈ Z, all the others B i j being equal to 0. Starting with an initial state ρ (0) = ρ 0 ⊗ |0 0|, after one step we have the state

ρ (1) = Bρ 0 B * ⊗ |-1 -1| + Cρ 0 C * ⊗ |1 1| .
The probability of presence in |-1 is Tr (Bρ 0 B * ) and the probability of presence in |1 is Tr (Cρ 0 C * ).

After the second step, the state of the system is

ρ (2) = B 2 ρ 0 B 2 * ⊗ |-2 -2| + C 2 ρ 0 C 2 * ⊗ |2 2|+ + (CBρ 0 B * C * + BCρ 0 C * B * ) ⊗ |0 0| .
The associated probabilities for the presence in |-2 , |0 , |2 are then

Tr (B 2 ρ 0 B 2 * ), Tr (CBρ 0 B * C * + BCρ 0 C * B * ) and Tr (C 2 ρ 0 C 2 * ) ,
respectively.

One can iterate the above procedure and generate our open quantum random walk on Z.

As further example, that we shall discuss more later, take

B = 1 √ 3 1 1 0 1 and C = 1 √ 3 1 0 -1 1 .
The

ρ (0) = 1 0 0 0 ⊗ |0 0| ,
we find the following probabilities for the 4 first steps:

| -4 | -3 | -2 | -1 |0 | + 1 | + 2 | + 3 | + 4 n = 0 1 n = 1 1 3 2 3 n = 2

An Example on a Graph

In order to give an example on finite graphs it is useful to fix a notation. We shall denote the operators involved in the dissipative quantum random walk in a way similar to the notation of stochastic matrices for Markov chains. If the set of vertices is V = {1, . . . , V }, we shall denote the operators B i j inside a V × V -matrix as follows:

     B 1 1 B 2 1 . . . B V 1 B 1 2 B 2 2 . . . B V 2 . . . . . . . . . . . . B 1 V B 2 V . . . B V V      .
That is, on line j are all the operators for the contributions B i j which start from j and go to another site i. The usual property for stochastic matrices that the sum of each line is 1, is replaced by

i B i j * B i j = I for each line.
With this notation one can easily describe examples. On the graph with two vertices we consider the transition operators of the form

D 1 D 2 B C
where D 1 and D 2 are any diagonal matrices such that D * 1 D 1 + D * 2 D 2 = I and where

B = 0 √ p 0 0 and C = 1 0 0 √ 1 -p ,
for some p ∈ (0, 1). For example, with

D 1 = 1 √ 2 0 0 1 √ 3 and D 2 = 1 √ 2 0 0 2 3
and p = 1/4, with the initial state

1 2 I ⊗ |1 1| ,
we find the following probabilities for the first steps site:

|1 |2 n = 0 1 0 n = 1 0.42 0.58 n = 2 0.22 0.78 n = 3 0.16 0.84 n = 4
0.13 0.87

In Section 12 we shall study this example in more details.

Inhomogeneous Case

Note that our presentation of Open Quantum Random Walks is based on a behavior which is "homogeneous" in time. This is to say that our assumption is that the operators B i j for passing from site |j to site |i are always the same, step after step. This is the analogue of an homogenous Markov chain.

Of course one may be interested in even more complicated situations, such as non-homogeneous Markov chains. This would correspond here to choosing time-dependent operators B i j (n) for the transitions. For each n the constraint

i B i j (n) * B i j (n) = I
has to be respected for all j. This allows to define the time-dependent completely positive maps M(n) by

M(n)(ρ) = i,j B i j (n) ⊗ |i j| ρ B i j (n) ⊗ |i j| * .
The quantum random walk is then the one given by the products of these completely positive maps L(n):

ρ (n+1) = M(n + 1)(ρ (n) ) .

Recovering Classical Markov Chains

It is very interesting to notice that all the classical Markov chains can be recovered as particular cases of Open Quantum Random Walks. We only treat here the homogenous case, the discussion would be similar in the nonhomogeneous case.

Consider P = (P (j, i)) a stochastic matrix, that is P (j, i) are classical probability transitions on V. They express the transition probabilities of a Markov chain (X n ) on V, that is,

P (j, i) = P(X n+1 = i | X n = j) .
In particular, recall that i∈V P (j, i) = 1 for all j. Proposition 6.1 Put H = K = C V and consider any family of unitary operators U i j on C N , i, j ∈ V. Consider the operators

B i j = P (j, i) U i j .
They satisfy

i B i j * B i j = I
for all j. Furthermore, given any initial state ρ (0) , the associated open quantum random walk (M n ) has the same probability distributions (p (n) ) as the classical Markov chain (X n ) with transition probability matrix P and initial measure p

(0) i = Tr ( i| ρ (0) |i ) .
Proof The relation on the operators B i j is obvious for

i B i j * B i j = i P (j, i) U i j * U i j = i P (j, i)I = I .
Whatever is the initial state ρ, if we put ρ i = i| ρ |i , we get by Proposition 2.3 and its proof

M(ρ) = i k P (k, i)U i k ρ k U i k * ⊗ |i i| .
The probability to be located on site i is then

k P (k, i) Tr U i k ρ k U i k * = k P (k, i) Tr (ρ k ) .
That is, we get the classical transition probabilities for a classical Markov chain on the set V, driven by the transition probabilities P (i, j) and with initial measure p (0) i = Tr (ρ i ). After two steps, the probability to be located at site i is

k l P (l, k)P (k, i) Tr U i k U k l ρ l U k l * U i k * = k l P (l, k)P (k, i) Tr (ρ l ) .
That is, once again the usual transition probabilities for two steps of the above Markov chain. It is not difficult to get convinced, by induction, that this works for any number of steps.

Quantum Trajectories

Coming back to general Open Quantum Random Walks, we shall now describe a very interesting way to simulate them by means of Quantum Trajectories. This property seems very important when one wants to study the limit behavior of these quantum random walks. The principle of the quantum trajectories associated to an open quantum random walk is the following. Starting from any initial state ρ on H ⊗ K we apply the mapping M and then a measurement of the position in K. We end up with a random result for the measurement and a reduction of the wave-packet gives rise to a random state on H ⊗ K of the form

ρ i ⊗ |i i| .
We then apply the procedure again: an action of the mapping M and a measurement of the position in K.

Theorem 7.1 By repeatedly applying the completely positive map M and a measurement of the position on K, one obtains a sequence of random states on H ⊗ K. This sequence is a non-homogenous Markov chain with law being described as follows. If the state of the chain at time n is ρ ⊗ |j j|, then at time n + 1 it jumps to one of the values

1 p(i) B i j ρB i j * ⊗ |i i| , i ∈ V, with probability p(i) = Tr B i j ρB i j * .
This Markov chain (ρ (n) ) is a simulation of the master equation driven by

M, that is, E ρ (n+1) | ρ (n) = ρ = M(ρ) .
Furthermore, if the initial state is a pure state, then the quantum trajectory stays valued in pure states and the Markov chain is described as follows. If the state of the chain at time n is the pure state |ϕ ⊗ |j , then at time n + 1 it jumps to one of the values

1 p(i) B i j |ϕ ⊗ |i , i ∈ V, with probability p(i) = B i j |ϕ 2 .
Proof Let ρ ⊗ |j j| be the initial state. After acting by M the state is

i (B j i ρB j i * ) ⊗ |i i| .
Measuring the vertices, gives the site i with probability

p(i) = Tr (B j i ρB j i * ) .
By the usual wave-packet reduction postulate, the state after having been measured with this value is

1 p(i) (B j i ρB j i * ) ⊗ |i i| .
This state being given, if we repeat the procedure, then clearly the next step depends only on the new state of the system. We end up with a (nonhomogenous) Markov chain structure.

On average, the values of this Markov chain after one step is

E ρ (n+1) | ρ (n) = ρ ⊗ |j j| = i p(i) 1 p(i) (B j i ρB j i * ) ⊗ |i i| = i (B j i ρB j i * ) ⊗ |i i| = M(ρ (n) ) .
If ρ is a pure state |φ φ| ⊗ |i i|, then it stays a pure state at each step. Indeed, any initial pure state |φ φ| ⊗ |i i| will jump randomly to one of the states 1

p j i B j i |φ φ|B j i * ⊗ |j j| with probability p(i) = Tr (B j i |φ φ|B j i * ) .
In other words, it jumps from the pure state |φ ⊗|i to any of the pure states 1

p(i) B j i |φ ⊗ |j with probability p(i) = B j i |φ 2 .
We have a classical Markov chain valued in the space of wave functions of the form |φ ⊗ |i . On average, this random walk simulates the master equation driven by M.

Physical Implementation

It is natural to wonder how such Open Quantum Random Walks can actually be realized physically. We shall here discuss a way to achieve it.

For the sake of a simple discussion, we restrict ourselves in this section to the case where either V is finite, or the number of non-vanishing B i j 's is finite for every fixed j. This is the case in all our examples and makes all the sums finite in the following.

Consider an open quantum random walk on V with chirality space H and with associated transition operators B i j . Recall that we have supposed that

i∈V B i j * B i j = I
for all j ∈ V. Hence, for all j ∈ V there exists a unitary operator U (j) on H ⊗ K whose first column (we choose |1 to be the first vector) is given by

U i 1 (j) = B i j .
This unitary operator is a unitary operator that dilates the completely positive map

M j (ρ) = i B i j ρ B i j * .
In other words, the completely positive map M j on H is the partial trace of some unitary interaction between H and some environment E. It is wellknown that the dimension of the environment can be chosen to be the same as the number of Krauss operators appearing in the decomposition of M j , that is, in our case they are indexed by V. Hence the environment can be chosen to be E = K.

The state space on which one performs the physical implementation is H ⊗ K 1 ⊗ K 2 where K 1 and K 2 are two copies of K. Let us present the main ingredients which shall appear in the physical implementation.

Each unitary operator U (j) defined above acts on H ⊗ K 1 . We construct the unitary operator

U = j U (j) ⊗ |j j| which acts now on H ⊗ K 1 ⊗ K 2 .
We shall also need the so-called swap operator S on K 1 ⊗ K 2 defined by

S(|j ⊗ |k ) = |k ⊗ |j .
It is a unitary operator on K 1 ⊗ K 2 which simply expresses the fact of exchanging the two systems K 1 and K 2 . We shall also use a decoherence procedure on the space K 1 , along the basis (|i ). By this we mean the following: if the system is in a superposition of pure states

|ϕ = i λ i |i ,
then this system is coupled to an environment in such a way and in a sufficiently long time, for the state of K 1 to become

i |λ i | 2 |i i| .
This is to say that we have chosen a coupling of K 1 with some environment which makes the off-diagonal terms of the density matrix |ϕ ϕ| converge exponentially fast to 0. This kind of decoherence is now well-known in physics.

It is rather easy to describe an environment and an explicit Hamiltonian which will produce such a result. For example, by performing a coupling to an environment which gives rise on K to a Lindblad semigroup evolution with generator

L(ρ) = - 1 2 [Q , [Q , ρ] ]
where Q is the operator

Q = i i |i i| ,
for example. In that case it is easy to see that any initial state ρ 0 has its off-diagonal terms decaying like exp(-(i -j) 2 t/2) and the diagonal terms remaining constant. One can give an Hamiltonian description of such a decoherence procedure by using the continuous-time limit of repeated quantum interactions as developed in [A-P].

Finally, we shall need a refreshing procedure, that is, if K 1 is in any state ρ then we put it back to the state |1 1|.

Proposition 8.1 Consider the quantum system H ⊗ K 1 ⊗ K 2 , together with some initial state

ρ (0) = k ρ k ⊗ |1 1| ⊗ |k k| .
If we perform successively 1) an action of the unitary operator U , 2) a decoherence on the basis (|i ) of the system K 1 , 3) an action of the swap operator I ⊗ S 4) a refreshing of the system K 1 to the state |1 1| then the state of the system becomes

k l B k l ρ l B k l * ⊗ |1 1| ⊗ |k k| .
That is, one reads the first step of the dissipative quantum random walk on H ⊗ K 2 . By iterating this whole procedure one produces the dissipative quantum random walk on H ⊗ K 2 .

Proof The unitary operator U (k) admits a decomposition

U (k) = i,j U i j (k) ⊗ |j i| .
In particular we have

j U i j (k) * U i j (k) = δ i,i I .
On the space H⊗K 1 ⊗K 2 the operator U as defined above is then decomposed into

U = i,j,k U i j (k) ⊗ |j i| ⊗ |k k| .
Now starting in a pure state |φ ⊗ |1 ⊗ |k , we get

U (|φ ⊗ |1 ⊗ |k ) = j U j 1 (k)|φ ⊗ |j ⊗ |k = j B j k |φ ⊗ |j ⊗ |k .
This is the first step of the procedure. The second step consists in performing a decoherence on the first space K. The pure state

On the space H and the second space K one can now read the first step of our quantum random walk. Finally, refresh the first space K into the state |1 , we then end up with the state

j B j k |φ φ|B j k * ⊗ |1 1| ⊗ |j j| ,
on which one can apply our procedure again.

If the initial state is not a pure state but a density matrix, a mixture of pure states, it is not difficult to see that the procedure described above gives the right combination and the right final state.

To summarize, the quantum random walk is obtained in the following way. Dilate each of the maps L k into a unitary operator U (k) on H ⊗ K 1 , start in the desired initial state on H and the second space K 2 , with the first space K 1 being in the state |1 , then iterate the following procedure on H ⊗ K 1 ⊗ K 2 : -apply the unitary operator k U (k) ⊗ |k k|, -perform a decoherence on K 1 -apply the unitary shift I ⊗ S -refresh the first space K 1 into the state |1 . The dissipative quantum random walk now appears on H ⊗ K 2 .

Coming Back to Examples

Let us illustrate this physical implementation on the two physical examples developed in Section 4.

The Example on Z

In the case of stationary walks on Z the procedure can be considerably simplified, as follows. The procedure we describe below is slightly different from the one presented in Proposition 8.1, but it is actually the same one, presented in a different way, taking into account several simplifications offered by the model.

Consider an open quantum random walk on Z driven by two operators B and C on H. Consider a unitary operator U on H ⊗ C 2 of the form

U = B X C Y ,
that is, a dilation of the completely positive map driven by B and C. Let K = C Z and consider the space H ⊗ C 2 ⊗K. On the space C 2 ⊗K we consider the shift operator given by Applying the shift operator, the state becomes

S(|0 0| ⊗ |k k|) = |0 0| ⊗ |k -1 k -1| and S(|1 1| ⊗ |k k|) = |0 0| ⊗ |k + 1 k + 1| .
B|ϕ ϕ|B * ⊗ |0 0| ⊗ |k -1 k -1| + C|ϕ ϕ|C * ⊗ |1 1| ⊗ |k + 1 k + 1| .
Refreshing the space C 2 we end up with

B|ϕ ϕ|B * ⊗ |0 0| ⊗ |k -1 k -1| + C|ϕ ϕ|C * ⊗ |0 0| ⊗ |k + 1 k + 1| .
One can read the first step of the dissipative quantum random walk on H⊗K:

B|ϕ ϕ|B * ⊗ |k -1 k -1| + C|ϕ ϕ|C * ⊗ |k + 1 k + 1| .

The Example on a Graph

Let us now detail the second example of Section 4, the open quantum random walk on the 2-vertices graph.

Consider a two-level quantum system H coupled to another two-level quantum system K 1 via the Hamiltonian

H = iγ(a ⊗ a * -a * ⊗ a)
where a = 0 1 0 0 .

Then the unitary evolution associated to this Hamiltonian, for a time length t = 1 is given by

U = e -iH =     1 0 0 0 0 cos(γ) -sin(γ) 0 0 sin(γ) cos(γ) 0 0 0 0 1     .
Hence, for a good choice of γ, that is, for sin(γ) = √ p we have

U =     1 0 0 0 0 √ 1 -p - √ p 0 0 √ p √ 1 -p 0 0 0 0 1     .
In other words U is of the form

U = C X B Y
as a block matrix on K 1 , where B and C are those matrices associated to our example. This is to say that we have given here an explicit dilation of the completely positive map associated to the matrices B and C.

If D 1 and D 2 are two diagonal matrices satisfying D * 1 D 1 + D * 2 D 2 = I then assume, for simplicity only, that they have real entries

D 1 = a 0 0 α , D 2 = b 0 0 β , with a 2 + b 2 = α 2 + β 2 = 1.
Then, one can write a = cos(λ) and α = cos(µ).

Considering the Hamiltonian

K = λ 0 0 µ ⊗ 0 -i i 0 ,
we get that e -iK is of the form

V = D 1 X D 2 Y .
We have realized a concrete physical dilation of the completely positive map associated to D 1 and D 2 . Following Proposition 8.1, consider on

H ⊗ K 1 ⊗ K 2 = C 2 ⊗ C 2 ⊗ C 2 the unitary evolution V 0 0 U ,
written as a block matrix on K 2 . This is to say that H is coupled to K 1 with the Hamiltonian K when K 2 is in the state |1 1| and H is coupled to K 1 with the Hamiltonian H when K 2 is in the state |2 2|. In this context, the swap operator S takes the following simple form on

K 1 ⊗ K 2 S =     1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1     .
Then, following the four steps of Proposition 8.1 gives a realization of the associated quantum random walk on H ⊗ K 2 .

Unitary Quantum Random Walks

The Open Quantum Random Walks we have been describing up to now are actually very different from the well-known Unitary Quantum Random Walks, such as the Hadamard random walk (see Introduction for some references). This is to say that they produce probability distributions which are not of the same type as the ones usually observed with the Hadamard quantum random walks.

It seems to us that there is no way to produce limit distributions such as the one observed in the Hadamard quantum random walk central limit theorem, with open quantum random walks. The limit behaviors of Open Quantum Random Walks seems to be all Gaussian or mixtures of Gaussians. The dissipative character of our quantum random walks makes them very different from the unitary evolution describing the usual type of quantum random walks. In fact there is quite a surprising and strong link between the two types of quantum random walks. Let us develop it here.

Let V be a set of vertices, let H be a Hilbert space representing the chirality. For each pair (i, j) in V 2 we have a bounded operator B i j on H.

Instead of the usual condition

i B i j * B i j = I
for all j, we shall ask here a much stronger condition, namely for all j, j ∈ V

i B i j * B i j = δ jj I . ( 10
)
In other words, being given two starting points j and j , the sum of the "contributions" which go to the same points i ∈ V vanish, unless j = j in which case we recover the usual condition. Note that there is no classical analogue of this condition for classical Markov matrices.

Let us illustrate this condition with an example. For a stationary quantum random walk on Z we are given two operators B and C on H which represent the effect of making one step to the left or one step to the right. The usual condition, obtained by taking j = j gives

B * B + C * C = I .
Now, taking j = j + 1, we get a supplementary condition:

C * B = 0 .
This is the only new condition added to the usual one in that case. Note that these two conditions together imply in particular that B + C is unitary.

These two conditions are typically satisfied by the following family of examples. Let This is typically the case with Hadamard random walk where

U = 1 √ 2 1 1 1 -1 .
Let us see what happens, in the general context, with this additional condition. The point is the following, if we are given a pure state on H ⊗ K of the form

|ψ = i |ϕ i ⊗ |i with the condition ψ 2 = i ϕ i 2 = 1 then the state |ψ = i j B i j |ϕ j ⊗ |i
is of the same form and satisfies

ψ 2 = i j,j ϕ j , B i j * B i j ϕ j = j,j ϕ j , δ jj I ϕ j = j ϕ j 2 = 1 .
Hence, at each step we get a state of the form

|ψ = i |ϕ i ⊗ |i with the condition ψ 2 = i ϕ i 2 = 1 .
In particular it determines, at each step, a probability distribution on V by putting

P (i) = ϕ i 2 .
This is exactly the picture for the Unitary Quantum Random Walks, such as the Hadamard quantum random walk. Now the interesting point is the way one can physically realize such Unitary Quantum Random Walks and the way this construction is similar to the one of Open Quantum Random Walks.

Proposition 10.1 If the transition operators B i j satisfy the more restrictive condition (10), then applying the same physical procedure as in Proposition 8.1 without the decoherence step (step 2) gives rise to a unitary quantum random walk.

Proof Let us follow again the steps of the construction in Proposition 8.1. Starting in a pure state |φ ⊗ |1 ⊗ |k , we get

U (|φ ⊗ |1 ⊗ |k ) = j U j 1 (k)|φ ⊗ |j ⊗ |k = j B j k |φ ⊗ |j ⊗ |k .
This is the first step of the procedure. We now skip the decoherence part and apply I ⊗ S to the state. We get the state

j B j k |φ ⊗ |k ⊗ |j . (11)
On the space H and the second space K one can now read the first step of the quantum random walk. Finally, refresh the first space K into the state |1 , we then end up with the state

j B j k |φ ⊗ |1 ⊗ |j ,
on which one can apply our procedure again. We recognize the action of the type of quantum random walks we announced.

Examples on Z

We are now back to Open Quantum Random Walks. In this section and the following one, we shall review several concrete examples. We shall show up their limit behaviors, based on simulations.

A Walk With Only One Step to the Left

Recall the way one constructs dissipative quantum random walks on Z with the help of two operators B and C on H (Section 4). We start with a rather simple example. Take

B = 0 √ p 0 0 and C = 1 0 0 √ 1 -p ,
for some p ∈ [0, 1]. Put q = 1 -p. Consider the initial state

ρ 0 = α z z β ⊗ |0 0| .
After the first step, the state on the site | -1 is pβ 0 0 0 and the state on the site | + 1 is

α 0 0 qβ .
This means that the value | -1 is reached with probability pβ and the value | + 1 with probability α + qβ.

Let us compute the second step. As B 2 = 0 then | -2 is not reached. Only the sites |0 and | + 2 are reached, with respective states

(p + pq)β 0 0 0 = (1 -q 2 )β 0 0 0 and α 0 0 q 2 β .
The probability to reach | + 2 is α + q 2 β and the probability to reach |0 is (1 -q 2 )β.

By an easy induction we get that the state at time n is supported by |n -2 and |n only, with respective states (1 -q n )β 0 0 0 and α 0 0 q n β , with probability (1 -q n )β and α + q n β, respectively. Hence, computing for example the quantum trajectories associated to this open quantum random walk, we have the behavior of a random walk which goes straight to the right, with only one possible jump to the left. Before the first step to the left, the probability to go to the right, at step n, is α + q n+1 β α + q n β and the probability to jump to the left is (q n -q n+1 )β α + q n β .

After the first step to the left, the walk will only go to the right with probability one.

As n goes to +∞ the states on |n -2 and |n converge respectively to β 0 0 0 and α 0 0 0 .

A More Quantum Example

Take the same structure of quantum random walk on Z, but with

B = 1 √ 3 1 1 0 1 and C = 1 √ 3 1 0 -1 1 .
We have given in Section 4 the first few probability distributions associated to this quantum random walk. It appears immediately rather wild, asymmetric and uncentered.

The point is that, on numerical simulations, one can see these quantum distributions tend to a centered Gaussian. We start with the density matrix

ρ (0) = 1 0 0 0 ⊗ |0 0| .

Examples with Several Gaussians

One can produce examples where several Gaussians are appearing, including the case where Gaussians are reduced to Dirac masses. For example, taking

B = 0 0 0 2 √ 5 and C = 1 0 0 1 √ 5 ,
gives the following shapes (Figure 2), when starting with the initial state When n tends to +∞ then both the "Dirac soliton" on the right hand side and the Gaussian on the left hand side survive.

ρ (0) = 1 2 1 0 0 1 ⊗ |0 0| .
One can produce a perturbation of the above model, it will then converge to a different limit. Put

B =    0 β 0 4 5 -β 2 2    and C =    1 0 0 1 5 -β 2 2    .
We get the following shapes (Figures 3,4 and 5), where the Gaussian does not survive anymore at +∞.

An Example in Dimension 5

Another interesting family of behaviors is obtained with the following choice of B and C. For the sake of a compact notation, we put C 2 = cos(2t), C 4 = cos(4t), S 2 = sin(2t) and S 4 = sin(4t). Consider the matrices

B = 1 4          0 -2S 2 -S 4 0 2S 2 -S 4 0 -2S 2 -S 4 0 -2 3 2 S 4 0 2S 2 -S 4 0 -2 3 2 S 4 0 -2 3 2 S 4 0 2S 2 -S 4 0 -2 3 2 S 4 0 -2S 2 -S 4 0 2S 2 -S 4 0 -2S 2 -S 4 0          and C = 1 8       L 0 C 0 L 0 4(C 2 + C 4 ) 0 4(-C 2 + C 4 ) 0 C 0 2(1 + 3C 4 ) 0 C 0 4(-C 2 + C 4 ) 0 4(C 2 + C 4 ) 0 L 0 C 0 L      
, where

L = 3 + 4C 2 + C 4 , L = 3 -4C 2 + C 4 , C = - √ 6(1 -C 4 ) .
Simulations of this open quantum random walk indicates that the limit behavior exhibits two Gaussians plus a Dirac soliton. The two Gaussians get slowly constructed, point by point, as the soliton loses its mass. In Figure 6 we show the time evolution when the parameter t is equal to t = π/40, the initial state being ρ (0) = 1 5 I ⊗ |0 0| . Changing the parameter t makes the Gaussians moving at different speeds and even change their direction (Figures 7 and8). It is easy to produce Open Quantum Random Walks on Z 2 by specifying 4 matrices N, W, S, E on H which satisfy

N * N + W * W + S * S + E * E = I . (12)
Then, we ask the random walk to jump from any site to the four nearest neighbors, following N , W , S or E, respectively. One can for example combine two 1-dimensional Open Quantum Random Walks by asking them to act on the different coordinate axis. For example, take

N = √ λ 1 √ 3 1 1 0 1 and S = √ λ 1 √ 3 1 0 -1 1 together with W = (1 -λ) 0 cos(θ) cos(φ) 0 sin(θ) and E = (1 -λ) 1 0 0 cos(θ) sin(φ) , for some λ ∈ [0, 1].
One can obtain behaviors with a single Gaussian, as below (Figure 9), with λ = 1/6, θ = π/3, φ = π/7. 12). But it does not seem obvious to us that the properties of the two-dimensional walk constructed this way can be directly related to the properties of the two one-dimensional walks which helped its construction.

As an example, take

B = 1 √ 3 1 1 0 1 and C = 1 √ 3 1 0 -1 1 together with A =    0 β 0 4 5 -β 2 2    and D =    1 0 0 1 5 -β 2 2    ,
which are two Open Quantum Random Walks on Z that we have already considered above. Even when β = 0 one does not see anymore any sign of the Dirac soliton behavior in the corresponding 2-dimensional random walk.

In the following picture (Figure 11), the initial state is

1 2 I ⊗ |(0, 0) (0, 0)| .
D 1 D 2 B C , with B = 0 √ p 0 0 and C = 1 0 0 √ q ,
for some p ∈ (0, 1), q = 1 -p and

D 1 = a 0 0 b and D 2 = c 0 0 d
where we assume for simplicity that a, b, c, d are reals such that a

2 + c 2 = 1, b 2 +d 2 = 1, 0 < a 2 < 1, 0 < b 2 < 1, a = b, c = d, ab = √ q, a 2 = q, b 2 = q.
All the cases can be considered and easily solved, but it is not worth developing them here. The cases that we consider here are the nontrivial ones.

Proposition 12.1 Whatever is the initial state, the state of the open quantum random walk defined above converges to

1 0 0 0 ⊗ |2 2| . The speed of convergence is v n where v = min{a 2 , b 2 , |ab| , √ q}. Proof Let ρ (n) = ρ (n) 1 ⊗ |1 1| + ρ (n)
2 ⊗ |2 2| be the state at time n, where we put

ρ (n) i = α (n) i z (n) i z(n) i β (n) i .
The induction formula is then

ρ (n+1) 1 = D 1 ρ (n) 1 D * 1 + Bρ (n) 2 B * ρ (n+1) 2 = D 2 ρ (n) 1 D * 2 + Cρ (n) 2 C * ,
Then any initial state, for example any state of the form ρ 0 ⊗ |1 1|, will converge to the state 1 0 0 0 ⊗ |N N | .

Though this example is rather classical in its behavior, it is interesting for it gives a model of a sort of "excitation transport": giving any initial state on the site 1 only, the state will then be, more or less quickly, transported along the chain and will end up into the excited state on the site |N . In this case the probability distribution can be easily seen (and proved) to be a soliton of probability, supported by three sites only, traveling along the chain at speed 1 and ending up into the excited state on |N . For any initial state, for example

An Example on a 4-Graph

1 2 0 0 0 1 ⊗ |3 3| + 1 2 0 0 0 1 ⊗ |4 4| ,
all the "mass" of the states will get transferred to the sites 1 and 2 and end up in the ground state.

13 Appendix: Proofs of the Lemmas

As promised we end up with the complete proofs of the lemmas given in Section 2. The last term converges to 0 with N , for the sum j I ⊗ |j j| converges strongly to I. Hence choose a N such that this term is smaller than ε/2. The first term of the right hand side converges to 0 with M , for the operators T M ⊗ |j j| converge strongly to I ⊗ |j j|. Choose a M such that this term is smaller than ε/2. We have proved the announced strong convergence.

Proof

Proof [of Lemma 2.4] Consider the following sesquilinear form on H 2 (φ, ψ) → φ ⊗ j , ρ (ψ ⊗ j) . As ρ is trace-class, the right hand side above is finite. This means that the left hand side is finite for every orthonormal basis (e n ), hence ρ j is trace-class.

In particular we have We have proved all the assertions of the lemma.

Lemma 2. 4

 4 Let ρ be a trace-class operator on H ⊗ K and (|j ) be an orthonormal basis of K. The operator (I ⊗ |i j|) ρ (I ⊗ |j i|) can be written as ρ j ⊗ |i i| for some trace-class operator ρ j on H, which we shall denote by j| ρ |j . Furthermore we have Tr ( j| ρ |j ) = Tr ρ (I ⊗ |j j|) .

  operators B and C do satisfy B * B + C * C = I. Let us consider the associated open quantum random walk on Z. Starting with the state

  |j j| ⊗ |k k| . (8) Applying I ⊗ S to the state (8) we get the state j B j k |φ φ|B j k * ⊗ |k k| ⊗ |j j| .

  Now, let us detail the procedure. Starting with a state |ϕ ⊗ |0 ⊗ |k we apply the operator U ⊗ I and end up with the state B|ϕ ⊗ |0 ⊗ |k + C|ϕ ⊗ |1 ⊗ |k . Applying the decoherence on C 2 we get the state B|ϕ ϕ|B * ⊗ |0 0| ⊗ |k k| + C|ϕ ϕ|C * ⊗ |1 1| ⊗ |k k| .

  and C satisfy B * B + C * C = I and C * B = 0 .

Figure 1 Figure 1 :

 11 Figure 1 then shows the distribution obtained at times n = 4, n = 8 and n = 20.

Figure 2 :

 2 Figure 2: An O.Q.R.W. on Z which gives rise to a Gaussian travelling to the left and a Dirac mass travelling to the right (at times n = 5, n = 10, n = 20)

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: With a small perturbation (β = 1/5) of the model of Fig. 2, one can see the walk starting with almost the same shape, but the Gaussian moving towards left does not survive (at times n = 5, n = 10, n = 20)

Figure 6 :

 6 Figure 6: Two Gaussians moving towards right are constructed point by point while the soliton loses its mass (parameter t = π/40, times n = 50, n = 100, n = 200)

Figure 7 :Figure 8 :

 78 Figure 7: With the parameter t = 3π/40, one Gaussian is now moving towards left (at times n = 50, n = 100, n = 200)

Figure 9 :

 9 Figure 9: An O.Q.R.W. on Z 2 which exhibits a single Gaussian asymptotically (at time n = 50)

Figure 10 :

 10 Figure 10: With different parameters, the O.Q.R.W. on Z 2 gives rise to two Gaussian pics (at time n = 30)

Figure 11 :

 11 Figure 11: In this combination of two one-dimensional O.Q.R.W. the soliton travelling to the right of Fig. 2 and Fig. 5 has now disappeared (at time n = 50, for β = 0 and β = 1/2)

  The model can be very much sped up by replacing the transition matrix with

1 =ρ 2 =ρ 3 =ρ 4 =

 1234 Take a graph made of 4 vertices |1 , |2 , |3 and |4 . They are connected as follows. B and C are as in the first example and where D 1 , D 2 , D 3 , D 4 are diagonal matrices such thatD * 1 D 1 + D * 2 D 2 = D * 3 D 3 + D * 4 D 4 = I .The invariant state of this quantum random walk isCρ 2 C * + D 1 ρ 3 D * 1 Cρ 1 C * + D 3 ρ 4 D * 3 Bρ 1 B * + D 4 ρ 4 D * 4 Bρ 2 B * + D 2 ρ 3 D * 2 .It is not very difficult to check that the unique invariant state for this random

  [of Lemma 2.1] Each of the operators B i ρB * i is positive and trace-class. Put Y N = n≤N B i ρB * i , for all N ∈ N. For all N < M the operator Y M -Y N is positive and trace-class. Put T N = n≤N B * i B i . Then, for all N < M we haveY M -Y N 1 = Tr (Y M -Y N ) = N <n≤M Tr (B i ρB * i ) = Tr ρ(T M -T N ) .As ρ can be decomposed inton λ n |e n e n | ,with n λ n < ∞, the last term above is equal ton λ n e n , (T M -T N )e n .As the operators T N converge strongly to T , each of the terms e n , (T M -T N )e n converge to 0 as N and M go to +∞. Furthermore, they are all bounded independently of N and M , for the sequence (T N e n ) N ∈N is bounded. By Lebesgue's theorem Y M -Y N 1 tends to 0, that is, (Y N ) N ∈N isa Cauchy sequence in trace-norm. Hence it converges to a trace-class operator Y . The identity Tr (Y ) = Tr (ρT ) is now obvious. Proof [of Lemma 2.2] Put T N = n≤N B * i B i . We have |j j| . Then for all φ ∈ H ⊗ K we have j≤N T M ⊗ |j j| φ -φ ≤ j≤N T M ⊗ |j j| φ -j≤N I ⊗ |j j| φ + + j≤N I ⊗ |j j| φ -φ ≤ j≤N T M ⊗ |j j| φ -I ⊗ |j j| φ + j≤N I ⊗ |j j| φ -φ .

  Since| φ ⊗ j , ρ (ψ ⊗ j) | ≤ φ ψ ρ this sesquilinear form is obviously of the form φ , ρ j ψ for some bounded operator ρ j on H. For any orthonormal basis (e n ) of H we have n | e n , ρ j e n | = n | e n ⊗ j , ρ (e n ⊗ j) | .

  Tr (ρ j ) = n e n ⊗ j , ρ (e n ⊗ j) = n i e n ⊗ i , ρ (I ⊗ |j j|) (e n ⊗ i) = Tr ρ (I ⊗ |j j|) .Now, for any φ, ψ ∈ H and any x, y ∈ K we haveφ ⊗ x , (I ⊗ |i j|) ρ (I ⊗ |j i|) (ψ ⊗ y) = x , i i , y φ ⊗ j , ρ (ψ ⊗ j) = x , i i , y φ , ρ j ψ = φ ⊗ x , (ρ ⊗ |i i|) (ψ ⊗ y) .This proves that (I ⊗ |i j|) ρ (I ⊗ |j i|) = ρ j ⊗ |i i| .

* Work supported by ANR project "HAM-MARK", N • ANR-09-BLAN-0098-01, by South African Research Chair Initiative of the Department of Science and Technology and National Research Foundation

1

which gives

2 . This can be easily solved explicitly, giving

One can see easily that all the terms are then converging to 0, with speed at least v n , except the term α

(after a few lines of computations), that is, it tends to 1.

Excitation Transport

The above idea can be generalized to a chain of N sites connected as follows


B Rey

K Mocellin

L Fourment

A node-nested Galerkin multigrid method for metal forging simulation

Keywords: algebraic multigrid, mixed formulation, finite element method, 3D metal forming, mesh coarsening, node-nested meshes, automatic remeshing

A node-nested Galerkin multigrid method is developed to solve systems provided by mixed formulations of 3D metal forming problems. An algebraic approach is used where operators are built on node-nested coarse meshes obtained by an automatic coarsening algorithm. This blackbox multigrid preconditioner is developed within the PETSc library. It is plugged to the FORGE3® finite element software. The linear rate of convergence and the very high efficiency of the resulting multigrid solver are evaluated for large scale problems with nonlinear behaviour.

Introduction

This work is motivated by the success of the finite element methods for solving three dimensional metal forming problems at industrial scale. However, the constant need for more accuracy and reliability of the computations result in smarter spatial discretization that cannot be only made available by more powerful computers. In fact, reducing the mesh size by a factor of 2 in each space direction in 3D, which is not much, results into an increase of the number of nodes by a factor 8. As at least 75% of the computational time of a finite element simulation is spent for the resolution of linear systems, using an iterative solver which computational cost is proportional to N 3/2 (N being the number of degrees of freedom), this reduction will result into an increase of computational time by about a factor of 20. This shows the need to reduce the computational cost of systems resolution and more particularly its dependency to the number of degrees of freedom, independently on the increase of the computers power. Therefore, solver with a linear rate of convergence is quite an attractive solution. Multigrid method is known to be efficient for solving systems derived from finite element problems [START_REF] Okasanya | Algebraic multigrid for stabilized finite element discretizations of the Navier-Stokes equations[END_REF]. As discussed in [START_REF] Brandt | The Gauss Center research in multiscale scientific computation[END_REF], multi-scale strategies have the unique potential of solving most of the mathematical problems with N unknowns in O(N) works. This paper discusses the implementation of a multigrid solver into the Forge3® F.E. software that is dedicated to 3D thermo-mechanical forging simulation.

Forging simulation problem

More details on Forge3® model are available in [START_REF] Chenot | Three dimentional finite element modeling of forging process[END_REF]. For simplicity, the problem equations are presented under hot forging conditions when the elastic deformations can be neglected and the material assumed to be homogenous, isotropic and incompressible. The mechanical problem consists in the equilibrium equation, where inertia and gravity contributions can be neglected, and the incompressibility equation:

⎩ ⎨ ⎧ Ω = Ω = in 0 div in 0 div v σ (1)
where Ω is a time dependant domain, σ the stress tensor and v the velocity. This couple of equations has to be completed with boundary conditions where the boundary of domain Ω is decomposed into:

F C Ω ∂ ∪ Ω ∂ = Ω ∂ ; C Ω ∂
is the contact surface with forming tools and

F Ω ∂
is the free surface :

( ) ( ) [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ Ω ∂ = - Ω ∂ ≤ Ω ∂ ≤ - Ω ∂ = Ω ∂ ∆ ∆ - = - C n die C n C die F c s q s f n v v n v v n v v K on 0 . on 0 on 0 . on 0 on 1 σ σ σ α τ (2)
The friction power law gives the tangent shear stress τ at the workpiece/tool interface as a function of the relative tangent velocity s v ∆ defined by:

( ) [ ]n n v v v v v die die s . - - - = ∆ ( 3
)
f α is a friction coefficient, K the material consistency, and q a sensitivity coefficient. die v denotes the velocity of the die and n the outside normal to the workpiece. The second equality of ( 2) is the free surface condition. The last three equations of (2) are the Signorini conditions, which denote the unilateral contact condition and the non-penetration condition of the workpiece into the tools. The contact pressure n σ is defined by:

n n n . σ σ = (4)
The constitutive law of the material is regarded as viscoplastic. The Norton-Hoff law is considered here:

( ) ε ε σ & & 1 3 K 2 p - = + m I ( 5
)
where p is the hydrostatic pressure, I the identity tensor, ε & the equivalent strain rate, m the viscoplasticity coefficient and ε& the strain rate tensor.

Invoking the virtual work principle, a weak formulation of the mechanical problem formulated in terms of velocity and pressure is obtained. Tetrahedra and a P1+/P1 mini element mixed interpolation are selected to discretize the problem. First introduced in [START_REF] Arnold | A stable finite element for Stokes equations[END_REF] this interpolation ensures the stability condition (Brezzi Babuska condition) by enriching the linear velocity field with a bubble function at the center of the element. The unilateral contact condition is enforced to finite element nodes by a penalty method, so the additional penalty terms that appear in the system matrix are located on its block diagonal. The resulting nonlinear set of equations is solved by a Newton-Raphson algorithm. In the linear case (

1 = = q m
), the mini-element interpolation provides a discrete linear system to solve that has the following matrix form:

⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ∆ ∆ ∆ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ p b v b T b T b b R R R P V V B B B H B H 0 0 0 (6)
where V ∆ and P ∆ respectively are the nodal vector corrections to the initial values of the velocities and pressure, and b V ∆ is the velocity bubble contribution at the center of the elements. The internal velocity bubble terms are eliminated from the equations by static condensation, which provides the arrival of a new diagonal block matrix C defined by:

b b T b B H B C 1 - = (7)
The remaining linear parts of the correction of the velocity and pressure are then solution of the following system:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∆ ∆ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - b b T b p v T R H B R R P V C B B H 1 (8)
where H and C are both symmetric positive definite, so that the global matrix of the system is not definite.

The system unknowns are the three velocity components and the pressure for each node of the finite element mesh: the number of degrees of freedom (dof) of the system equals four times the number of nodes. This linear system is usually solved by a preconditioned conjugate residual method (PCR), where preconditioning is performed by an incomplete Cholesky factorization. This reference preconditioner is replaced by the multigrid method in the present work and more precisely by a single iteration of a multigrid V-cycle, as in [START_REF] Iwamura | An efficient algebraic multigrid preconditioned conjugate gradient solver[END_REF].

Multigrid

Notations

The algebraic multigrid approach [START_REF] Strüben | A review of algebraic multigrid[END_REF] is used to compute the coarse matrix, so it is defined by the Galerkin strategy:

P RA A fine coarse = (9)
The Ritz-Galerkin strategy is also selected:

T R P = (10)
The coarsest grid correction is compute by a direct solver that uses LU factorization.

Figure 1 shows the utilised n-grids V-cycle with a smoother ( ) b A S , , a restriction operator i R from the fine grid 1 i to the coarser grid i , and a prolongation operator T i R from the coarse grid i to the finer grid 1 i . ( ) b A D , is a direct solver used on the coarsest level. The system to solve is so denoted ---------------------------------------------------------------------------------------------------------------------------------------fig. 1: Multigrid V-cycle algorithm

0 0 0 b X A = where 0 X is the unknown. ---------------------------------------------------------------------------------------------------------------------------------------- function ( ) i i r A MG , if n i < then ( ) i i i b A S X , ← // pre-smoothing step i i i i X A b r - = // residual computation i i i r R r = +1 // residual restriction ) , ( 1 1 1 + + + = i i i r A MG X δ // correction computation 1 + = i T i i X R X δ δ // correction prolongation i i i X X X δ + ← // correction ( ) i i i b A S X , ← // post-smoothing step else ( ) i i i b A D X , ← // coarsest grid direct resolution return i X end -

Automatic coarse grid creation

According to equation [START_REF] Chan | Domain decomposition and multigrid algorithms for elliptic problems unstructured meshes[END_REF], only the restriction operators need to be evaluated. It is done following a geometrical approach by constructing different levels of meshes of smaller number of element. Therefore the main issue is to build coarse meshes.

The use of an element-nested method as in [START_REF] Coupez | Toward large scale F.E. computation of hot forging process using iterative solvers, parallel computation and multigrid algorithms[END_REF] is prohibited because the remeshing issue requires a fully automatic procedure to generate the different levels of meshes. The use of a non-nested meshes as in [START_REF] Feng | A non-nested Galerkin multigrid method for solving linear and nonlinear solid mechanics problems[END_REF] is also leave behind because the R and P operators are more complex to build and more expensive to use than in a node-nested method. Therefore this node-nested method is preferred here.

For a 2D geometry, Guillard [START_REF] Guillard | Node-nested multigrid method with Delaunay coarsening[END_REF] and independently Chan and Smith [START_REF] Chan | Domain decomposition and multigrid algorithms for elliptic problems unstructured meshes[END_REF] proposed a Maximal Independent Subset (MIS) algorithm to compute coarse meshes from fine ones. Adams [START_REF] Adams | Parallel multigrid solvers for 3D-unstructured large deformation elasticity and plasticity finite element problems[END_REF] proposed a 3D extension of these algorithms but another technique [START_REF] Carte | Coarsening technique in multigrid applications on unstructured meshes[END_REF] is prefered here to compute coarse meshes. It is based on an iterative node removal algorithm combined to a local remeshing algorithm, with almost constant coarsening ratio. This automatic coarsening technique is part of the MTC® mesh generator toolbox [START_REF] Coupez | 3D anisotropic mesh generation and adaptation with applications[END_REF]. h M is defined by its set of nodes h I and its set of tetrahedra h T . For any node

h I i ∈ , the average length () i l h
of the edges connected to i is computed. Ideally, in the coarse mesh extracted from h M , for any

H I i ∈ the length ( ) i l H
of the edges will be ()

C i l h
, where C is the targeted average length ratio which is usually chosen equal to 2 for an isotropic mesh. The aim of the algorithm is to build a coarse mesh where the average length () i l H of any node i has the targeted value. Let CR be the coarsening ratio that is equal to the number of nodes of the initial mesh divided by the number of nodes of the coarse one.

Let () i T be the set of tetrahedra that contains the node i , ( ) i ε be the set of nodes of ( ) i T

(not including i ) and () i F be the set of exterior faces of ( ) i T (that do not contain i ). The removal of nodes is carried out by the algorithm described in figure 2. --------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------fig. 2: coarsening algorithm

( ) i current ε , let () i T new be the set of tetrahedra of this new triangulation • if () i T new is acceptable then o {} i I I current current - = o () ( ) i T i T T T new current current + - = o 1 + ← i i
The remeshing of the cavity in step 4 of figure 2 produced by removing node i , described in figure 3, is carried out by successively exploring all possible triangulations made by connecting a node of () i ε to all faces of ( ) i F . Among all these triangulations, only the ones which volume almost equals the volume of ( ) i F are kept and the one with the best element quality is selected.

fig. 3: cavity remeshing

This way the coarse meshes are node-nested with the initial one. Figure 4 show examples of coarse meshes obtained with this algorithm in the cases of a connecting rod and a ring at a certain stage of the forming process. 40,897 nodes 4,857 nodes In fact, the construction of coarse mesh is not as straight forward as it may seem from these examples, because the coarse mesh has to satisfy the geometrical constraint of properly representing the shape of the domain. Therefore, in order to obtain a satisfactory fully connex mesh with no significant volume or shape alteration, it is necessary to reduce the expectation in terms of coarsening ratio. For more complex geometries, as encountered at the end of forging, it is then difficult to obtain a coarse mesh with sufficiently small number of nodes, which can penalize the multigrid efficiency because of the computational cost of the direct system resolution on the coarsest level.

Operators computation

The restriction operator, is obtained by a linear interpolation of fine nodes values onto coarse nodes. Figure 5 describes the linear dependencies in a 2D case. For any node i , the coarse element ( )

3 2 1 , , j j j
that contains i is found and its barycentric coordinates are calculated. If the fine node is also a coarse one, there is only one non-zero value its row. Otherwise, there are only four non-zero values its row. With the used of Lagrangian formulation, the mesh endures large deformations during simulation. This deformation is not linear, so the barycentric coordinates of fine nodes with respect to the coarse mesh change. However, for simplicity these changes are not taken into account and operators are assumed to be constant. In [START_REF] Coupez | Toward large scale F.E. computation of hot forging process using iterative solvers, parallel computation and multigrid algorithms[END_REF], it has been observed that this approximation does not reduce the efficiency of the multigrid method with the material deformation.

Other multigrid components

The multigrid method is used as a preconditioner for the conjugate residual solver. This new preconditioned solver is developed within the PETSc library, which still need to be provided by coarse meshes, restriction operators and adequate numerical parameters.

In this frame, several smoothers available in PETSc have been tested. Block Jacobi has been selected because it offers the best compromise between high frequency smoothing and low computational cost.

The coarsest grid correction is calculated by direct resolution. Experience shows that the direct solver is regarded as the quickest when the number of dof is less than 6,000, i.e. when the coarsest mesh contains less than 1,500 nodes. Habitual meshes used for standard forging simulations contain between 10,000 and 80,000 nodes so three different levels with a coarsening ratio of 8 approximately each, are necessary and sufficient to obtain an efficient method where the coarsest mesh resolution can be direct. This is why all our numerical results have been obtained with a three levels method.

The number of pre and post smoothing steps are fit by a parametric study on a typical forging problem. Best results are obtained with one smoothing and one post smoothing iteration. It both minimizes the total number of smoothing steps and the computational time of the resolution.

Numerical results

This section investigates the performances of the developed three levels multigrid method (MG3). First of all the linear rate of convergence of the solver is verified. Then the performances are compared to the reference preconditioned conjugate residual solver (CR/ILU), the same convergence criterion being used for both methods.

Linear rate of convergence

The evolution of the computational time resolution is observed for various systems provided by the connecting rod forming at the end of the simulation (non linear material behaviour and friction coefficient equals 0.14). The different used meshes have 11,245 nodes, 31,200 nodes, 47,083 nodes, 61,987 nodes and 91,127 nodes. The left curve of figure 6 shows the evolution of the resolution time with CR/ILU (plain line) and MG3 (dotted line) solvers.

Comparing the resolution time of the two methods for the same system, the larger the number of dof is, the more the speed-up provided by MG3 is relevant. The gain is effective when the mesh contains more than 10,000 nodes.

The multigrid solver is close to its optimal linear rate of convergence () n O , where n is the number of unknowns: see figure 6 right, where the black line corresponds to a perfect linear behavior and the dotted line to the MG3 rate of convergence.

Application to the simulation of the ring rolling process

efficiency of the multigrid solver is now evaluated for a complex forging simulation, the forming of an external ring of a ball bearing system. Figure 7-1 shows the process with the outer die, the workpiece and the inner die, which has an excentric circular movement. Figure 7-2 shows the workpiece at the beginning and figure 7-3 at the end of the process. The viscoplasticity coefficient is 0.1529, the material has a consistency of 2707 MPa and a temperature of 700°C. Contact between workpiece and tools is supposed to be perfectly sliding. Table 2 summarizes the obtained calculation times. The speed up observed with the multigrid solver directly depends on the mesh size. For 44,980 dof the calculation time is reduced by 2.1. On a more substantial example, with 247,958 dof, the recorded speed up is 7.6 and the time saving is about 12h. These performances are really interesting and provide a considerable minimisation of calculation time for numerical simulations.

Conclusion

A node-nested Galerkin multigrid method has been developed for solving linear system derived from 3D forging simulations with a mixed finite element formulation. This method is plugged to Forge3® software and remains see-through for the user. Moreover, it is compatible with frequent remeshings and with the mixed velocity-pressure formulation. A standard algebraic method might encounter difficulties with the contact penalty terms of the system matrix, while the used mixed algebraic / geometric approach avoids it. An automatic coarsening technique allows computing the various coarse meshes required by the construction of the transfer operators. The node-nested characteristic provides sparse and simple interpolation and restriction passing matrices. The construction of a coarse system by the Galerkin relation is not time consuming. Operators are considered to be constant between two remeshing steps, and even if these steps do not occur frequently, no decrease of performances is observed because of this approximation. When the domain is too complex it is not always possible to build coarse meshes with a coarsening ratio close to 8, but the method has showed to be quite robust when applied to an entire ring rolling simulation with 44 remeshing steps.

This multigrid solver shows a promising global reduction of the computational time. It is also more robust that the previously used solver. Some more forging simulations have been carried out to extend these conclusions to other problems. The speed-up and the robustness of the multigrid solver were not denied. The linear asymptotic convergence rate has been noticed in the particular studied frame of a mixed velocity -pressure formulation.

These first results are really encouraging and confirm that the use of the multigrid is advised for large scale simulations. First improvements are currently "taking place". An iterative resolution method is currently tested to solve the coarsest system when the coarsest mesh contains too many nodes. Therefore, if the coarsening technique does not allow obtaining a sufficiently coarse mesh, an iterative solver will provide a less expensive solution than a direct solving. Future works will deal with the study of larger problems with several bodies, where the number of dof will grow up, so the use of the multigrid method will be even more justified. The solver will also have to be adapted to parallel computation using a SPMP strategy with dynamic mesh partionning.

Fig. 4 :

 4 Fig. 4: fine and coarse meshes of a connecting rod (coarsening ratio: 8.4)

fig. 5 :

 5 fig. 5: interpolation of fine nodes values onto the nodes of a coarse element

Fig. 6 :

 6 Fig. 6: computational cost of the multigrid and CR/ILU solver as a function of dof

Fig. 7 - 1 :

 71 Fig. 7-1: ring rolling forging process Fig. 7-2: ring at the beginning Fig. 7-3: ring at the end

Figure 8 Fig. 8 :

 88 Figure8shows the different meshes used for the 31,200 nodes simulation. The computational mesh is on the left hand side. In the centre, the intermediate mesh is showed, and on the right hand side the coarse mesh.

Table 1

 1 summarizes the data of the computational and different coarse meshes used by the multigrid method.

	Number of dof of the simulation	44,980	124,800	247,948
	Number of nodes of the computation mesh	11,245	31,200	61,987
	Number of nodes of the intermediate mesh	1,452	3,233	6,061
	Number of nodes of the coarse mesh	189	509	790
	Tab. 1: number of nodes of the different meshes

Table 2 :

 2 CPU time and speed up of the first three time increment of simulation

	Number of nodes	11,245	31,200	61,987
	Number of dof	44,980	124,800	247,948
	IC -CR	6min 53sec	3h 39min 51sec	14h 23min 2sec
	MG3	3min 15sec	43min 49sec	1h 53min 20sec
	Speed up	2.1	5.0	7.6