César Bernardini

Thomas Silverston

Olivier Festor

Using Social Network Information into ICN

Online Social Networks carry extremely valuable information about their users and their relationships. We argue that this knowledge can help to drastically improve the efficiency of Information Centric Networks.

In this paper, we propose a first step to include social information into ICN architectures. We conjecture a small number of users dominate the activity and receive most attention of others users in the social networks and we argue they produce content that will be more likely to be consumed, and in consequence their content must be replicated. We then propose a caching strategy based on prioritizing their content. We simulate a social network model where the proposed caching strategy is evaluated against common ICN caching strategies. Finally, we show that inclusion of social information into ICN networks may help to improve cache performances.

I. INTRODUCTION

Social Networking has exponentially grown in the last years. Nowadays, millions of users interact with each other through Facebook, Twitter, Digg, Google Plus, etc. People tweets about what they experience on their life with their community: president Obama's re-election picture became the single mostretweeted message in Twitter history with more than 800,000 re-tweets 1 ; In one minute, football player Lionel Messi was mentioned 35,000 times in Twitter after he broke the record of goals in a calendar year 2 ; french people issued 420,574 tweets regarding to the election of Miss France3 .Mobile phones' interfaces include functionalities to instantenously share information through the Social Networks. Most if not all the internet services improve your experience through the addition of social features to rapidly spread interesting content. Companies invest strongly into their Facebook pages to promote new products and profit as much as possible user's feedback [START_REF] Song | Detecting opinion leader dynamically in chinese news comments[END_REF]. News agencies dedicate employees to update their Twitter and Digg pages 4 . 90% of american hospitals use social media to attract new clients and one third has a formal social media plan 5 Users organize social events and keep updated on their friends and families pieces of news through the Internet. We believe these interaction model will continue to expand and to evolve. Internet is becoming social networks oriented.

At the same time and despite its remarkable success, several architectures have been proposed to overcome the limitations of the current Internet such as efficiency, availability, security and mobility. Information Centric Networking (ICN) is a promising new paradigm for the future Internet, where the communication depends on named-data, rather than host names. Indeed, content retrieval, content demand and content identification is lead by its name instead of its physical location. ICN features include in-network caching, multicast support, self-verification, encryption and protection against common current-Internet attacks. ICN architectures include Content Centric Networks (CCN) [START_REF] Jacobson | Networking named content[END_REF], NetInf [START_REF] Dannewitz | Prototyping a network of information[END_REF] and Pursuit [START_REF] Särelä | RTFM: Publish/subscribe internetworking architecture[END_REF] among others, compared in [START_REF] Ahlgren | A survey of information-centric networking[END_REF].

As we said, All ICN nodes support in-networking caching to serve future requests which may help to reduce transport cost for the network and to enhance end-users delivery performances. The availability of different replicas depends on several factors such as cache replacement policies, cache size, content's popularity. Several caching schemes has been evaluated [START_REF] Jacobson | Networking named content[END_REF], [START_REF] Chai | Cache "less for more" in information-centric networks[END_REF], [START_REF] Rossi | Caching performance of content centric networksunder multi-path routing (and more)[END_REF], [START_REF] Mathieu | Information-centric networking: a natural design for social network applications[END_REF], [START_REF] Bernardini | Mpc: Popularity-based caching strategy for content centric networks[END_REF] and there are not concensus about the final caching scheme for ICN. We believe the caching scheme for ICN will be based on social network tendencies.

Online Social Networks carry extremly valuable informations about their users and their relationships. We argue that this knowledge can help to drastically improve the efficiency of Information Centric Networks.

In this paper, we propose a first step to include social information into ICN architectures. We conjecture a small number of users dominate the activity and receive most attention of others users in the social networks. We name them Influential users and we argue they produce content that is more likely to be consumed, and in consequence their content must be replicated. We then propose a caching strategy based on prioritizing content from influential users into the social network. We simulate a social network model where the proposed caching strategy is evaluated against common CCN caching strategies. Finally, we show that inclusion of social information into ICN networks help to improve cache performances.

This article is organized as follows. In section II, we give an insight into social networks and ICN. Then in section III we introduce our social network model used across the paper. Section IV details the simulation environment and section V presents the results. We discuss our model and results in Section VI and finally, our conclusions are exposed in Section VII.

II. RELATED WORK

A. Social Network

Social Network Analysis enables the examination of patterns in the relationships among interacting users. It gained importance with the raising of the Web 2.0. Theory of efficient Hubs or influentials says a small number of users dominate the activity and receive most attention of other users. Enterprises have focused their research into opinion leadership in order to improve quality of their products through interpretation of media messages [START_REF] Bodendorf | Detecting opinion leaders and trends in online social networks[END_REF], [START_REF] Song | Detecting opinion leader dynamically in chinese news comments[END_REF]. While their final goal is different, political analysis follows the same approach [START_REF] Agarwal | Online collective action and the role of social media in mobilizing opinions: A case study on womens right-to-drive campaigns in saudi arabia[END_REF].

We support the idea of using Eigenvector and PageRank as metric to detect influentials nodes in a social network because Google determines most important pages across the Internet using PageRank [START_REF] Franceschet | Pagerank: standing on the shoulders of giants[END_REF]; PageRank is a variant of Eigenvector centrality measure [START_REF] Langville | A Survey of Eigenvector Methods for Web Information Retrieval[END_REF].

To the best of our knowledge, [START_REF] Mathieu | Information-centric networking: a natural design for social network applications[END_REF] is the only one to deal with social networks and ICN at the same time; in order to compare IP, CDN and ICN, they use the same architecture Twitter servers have in the current Internet. In this paper, we differs on the fact we use social networks' information to improve ICN.

B. Caching in ICN

The use of caches to increase content availability and to reduce perceived latency time has been deeply investigated in diverse environments such as Operative Systems, Web-browsers and Proxy-servers, and -we believe-it's welladdressed in [START_REF] Podlipnig | A survey of web cache replacement strategies[END_REF].

In the context of ICN, researched topics include different concepts such as replacement policies [START_REF] Jacobson | Networking named content[END_REF], [START_REF] Rossi | Caching performance of content centric networksunder multi-path routing (and more)[END_REF], [START_REF] Luca Muscariello | Accelerated lru replacement with probabilistic caching: A fast algorithm to approximate lfu[END_REF], alternative policies [START_REF] Chai | Cache "less for more" in information-centric networks[END_REF], [START_REF] Bernardini | Mpc: Popularity-based caching strategy for content centric networks[END_REF], different cache sizes [START_REF] Rossi | On sizing CCN content stores by exploiting topological information[END_REF] over a broad range of topologies such as Binary Trees [START_REF] Jacobson | Networking named content[END_REF], [START_REF] Rossi | Caching performance of content centric networksunder multi-path routing (and more)[END_REF] and common ISP structures [START_REF] Rossi | Caching performance of content centric networksunder multi-path routing (and more)[END_REF], [START_REF] Choi | A Survey on contentoriented networking for efficient content delivery[END_REF], [START_REF] Tyson | A trace-driven analysis of caching in content-centric networks[END_REF], [START_REF] Chai | Cache "less for more" in information-centric networks[END_REF]. Most of the work focus on common IP client-server model [START_REF] Rossi | Caching performance of content centric networksunder multi-path routing (and more)[END_REF], [START_REF] Jacobson | Networking named content[END_REF], [START_REF] Rossi | On sizing CCN content stores by exploiting topological information[END_REF], [START_REF] Choi | A Survey on contentoriented networking for efficient content delivery[END_REF], [START_REF] Chai | Cache "less for more" in information-centric networks[END_REF], [START_REF] Bernardini | Mpc: Popularity-based caching strategy for content centric networks[END_REF] and P2P model [START_REF] Tyson | A trace-driven analysis of caching in content-centric networks[END_REF]. To our knowledge, we are the first to use a social networks model and at the same time to include social networks information to create a social-aware caching strategy.

III. SOCIAL NETWORK MODEL

In this section, we introduce our model. First, we explain the social network model and how the content is produced and consumed. We then map the social users into an ICN network and present a small example to analyze model behavior. Finally, we introduce a social-aware caching strategy based on prioritizing content from influential users in the social network model.

A. Social Network Model

Social Networks allow users to publish new content at their own will. Their friends are always available to watch their friends updates through a news feed system. Every time, users find interesting content they may share it with their friends spreading and expanding the visible scope of information. This is what we call social network model: users publishing, consuming and sharing.

We define a social network model where users have personal preferences such as topic of interest and groups of friends. In the model, each user has two actions: publish and retrieve.

• Publish: production of new content, all the content issued is self-describing. I.e.: /userC/music/song1 means user C has published a song entitled song1. After retrieving messages, users may share their friends content if they find it interesting enough through a new publish message (retweet). • Retrieve: gets last content issued by their friends. For example, in Figure 1 the user A has a friendship relationship (red dashed line) with users B, E and F. Everytime, A issues a Retrieve message A gets last B, E and F content. We define a sequence of actions as a finite combination of publish and retrieve users' messages.

B. ICN Model: Content Centric Networks

We select Content Centric Networks (CCN) as the ICN Model due to his wide acceptance in the community. CCN communication architectures relies on two named primitives: Interest and Data. The name is self-describing and userreadable. A consumer requests content by broadcasting its interest all over the network; any node hearing the request and having the data can issue a response with a data message. As ICN architectures, nodes cache all the Data messages that have passed by them. Another fact worthy of mention is that caches have finite space which means replacement and acceptance policies must manage them.

As we said, every user has two actions: Publish and Retrieve. Publish consists of generating new personal content, it does not issue any message in the CCN network, it just update information of CCN nodes to future retrievals of the content. The Retrieve operation deals with getting last updates from his friends. The retrieve message generates an Interest for every of his friends. The Interest are then sent across the ICN sorthest path.

C. Example

We illustrate the model by an example where 6 users are distributed over a 9-CCN-nodes network depicting their relationships (red dashed lines) and physical connections (black lines) in Figure 1.

In Figure 1, we can see User E only cares about content from A because it's his only friend and B is concerned by content from Nodes A, C and D. Once B issues a retrieve message, last content from A, C & D users will be requested; when User E issues an retrieve message, only user A is requested by his last updates. In a more concrete example, we define the sequence of actions as [A.publish(1), B.publish(2), C.publish(3), D.publish(4), E.publish(5), F.publish(6), A.retrieve(), E.retrieve()]. After every user publish new content sequentially (A:1, B:2, C:3, D:4, E:5, F:6), internal CCN nodes information is updated and after next-retrieve message, the CCN Data message is issued.

Then, user A and E issue a retrieve message respectivelly. First, A request for last updates of his friends, in Figure 1, users E, F and B. In the CCN architecture, it means an Interest message is issued from node A to node E across the network path [A, D, g, F, E] which stores Content 5; then, A issue another Interest across [A, E, g, F] to receive Content 6, which is stored across the path; finally last Interest for B, stores 2 across the path [A, C, h, B]. Then, E issue a retrieve message which means that path [E, F, g, D, A] caches Content 1. Finally, the post-execution content of the caches is shown in Table I. We can see in the table that non-influential users content (5 and 6) is stored 5 times in the caches while influential content (1, 2) is stored 4 and 5 times. We aim at increasing the availability of influentials contents and we define our social-aware strategy in the incoming section.

D. Eigenvector Social-aware Caching Strategy

In the model, we mentioned ICN nodes' caches are managed with a caching strategy. To compare caching strategies with social-aware caching strategies, we propose a social-aware 1).

caching strategy which prioritizes influential-nodes content. We call our strategy Eigenvector.

The aim of Eigenvector strategy is to replicate content published by influential users before it's requested in order to improve availability and to reduce number of interest to get the content. We believe users with more friends are more influential than those with less number of friends and they produce content that will be more likely to be consumed.

We detect influential nodes using Eigenvector centrality metric. Eigenvector centrality metric is a measure of the influence of a node in a network; the nodes receives a score according to their importance into the graph; For instance, we first calculate the Eigenvector centrality metric into the social graph (red dashed lines in Figure1) and then we define a node/user as Influential if its eigenvector is greater than the average.

In other words, publish messages from non-influential users do not generate replication of content while publish messages from Influential users copy its brand-new content to his friends across the CCN caches in the shortest path of the topology.

Eigenvector works in combination with a common caching strategy (e.g. LRU, FIFO, RAND). For example, LRU+Eigenvector refers to Eigenvector with Last Recently Used (LRU) replacement policy. The common caching strategy manages the cached elements: everytime it receives a Data messages, it's stored and when the cache is full, the strategy decides which item to discard to make room for the new ones. In addition, every time influential nodes publish new content it is automatically replicated into the shortest ICN path to their Social Neighbours.

In Table II, we calculate the Eigenvector centrality metric over the graph of social relationships (red dashed lines) in Figure 1. Users A and B are influentials because its Eigenvector value is greater than Eigenvector averaged of all users (0.58 > 0.38) while users C, D, E and F are non-influentials (0.29 < 0.38). It means when User A publish new content it's copied to his friends across the topology path to E, F and B ( [A, D, g, F, E], [A, D, g, F] &[A, C, h, B] respectivelly). In the other hand, E publish message does not generate copies due his non-influential nature.

IV. SIMULATION ENVIRONMENT

We introduced a model where users publish and request their friends content. Users are distributed across a CCN network. In this section, we show the environment where the simulations were carried out. We developed a simulation tool in Python to represent our model.

To simulate the model, we generate a social network graph and topology structure for CCN. Then, we randomly map social network nodes to CCN nodes. We configure different caching strategies into the CCN nodes. Finally, we analyze the caching performance using a sequence of publish/retrieve messages.

Social networks are well-modeled by small world graphs. That's the reason we represent social networks structures with randomly generated graphs using small-world Newman-Watts-Strogatz algorithm [START_REF] Newman | Random graph models of social networks[END_REF] with 50 nodes.

The topological structure of the ICN Internet has not been decided yet. In this context, we choose randomly generated connected graphs to model the topology. The number of nodes is fixed at 100 and the number of edges is randomly selected. All the graph manipulation is handled with Networkx library [START_REF] Hagberg | Exploring network structure, dynamics, and function using NetworkX[END_REF].

We consider three well-known replacement strategies: Last Recently Used (LRU), Random replacement (RAND) and FIFO. We then compared the results of caching performance with social-aware Eigenvector caching strategies: LRU+Eigenvector, Rand+Eigenvector, FIFO+Eigenvector.

Every user has personal preferences and a random location in the CCN network, and according to his personal preferences: he requests, publishes and shares his thoughts in a twitter-like manner. The sequence of actions is generated according to a uniform distribution of 40 events per node (in average) to simulate caching performance of one day. In 2010, Twitter revealed that it has 105 millions registered users and 3,000 millions requests a day; which in simple calculations stand for almost 29 requests per day per user. We argue 2 years has passed from there, and we estimate 40 requests per day per user. 6 . After analyzing their friends' content he might decide to publish again his friend's content. This decision follows an uniform distributed variable.

V. RESULTS

In this section, we assess the performance of social-aware caching strategies against common-used caching strategies. As we pointed out in previous section, we defined a sequence of publish/retrieve messages and then we analyze the performance of CCN caches using several strategies and cache sizes.

First, we introduce the metrics to evalute the performance of the caches:

• Cache Hit Ratio: the probability to obtain a cache hit all along the path from a requester to a cache node; • Stretch: the number of hops that the data chunk has travelled in the network with respect to the server storing original copy; • Diversity: express the ratio of distinct content stored across all the caches. • Expired-chunks: express the ratio of non-longer-valable content stored across the caches. 6 http://www.businessinsider.com/twitter-stats-2010-4?op=1

We then analyze the results. From Figure 2a, we observe that all Eigenvector social-aware strategies perform 20% -in average-better than common-used caching strategies in terms of Cache Hit. It is noteworthy that bigger cache sizes do not improve performances due to the highly large amount of new publish data.

Stretch is further analyzed in Figure 2b, social caching schemes improves common replacement strategies. Commonused strategies reduces distance to reach content by 33% (stretch of 66%), but our social-aware caching strategy reduce it by 66%. In other words, our social-aware strategy requires to traverse in average 33% of the path to the content while common-used caching strategies need to pass through 66% of the path. Even more, the distance is shortened 50% in comparison with common-used caching strategies. For example using common-used strategies, some user requires to traverse 7-hops in average in order to access content located 10-hops away. while using our social-aware strategy the user requires just to pass through 3-hops.

In the Figure 2c, we depict that social-aware strategies reduce by 6% the ratio of expired elements stored. At first sight, it may seem an slight improvement but every outdate elements represent wasted space into the caches which means we are taking advantage of 6% more of the cache sizes using our approach.

In terms of diversity, common-used strategies keep diverse contents in their caches in comparison with social-aware caching strategies. Even though in Figure 2d common-used strategies are sligthly better in diversity terms (2% in average), our approach drastically increases the performance of 20% and 30% regarding cache hit and stretch respectivelly. For sake of clearness, we only illustrate LRU strategy against LRU+Eigenvector in the Figure 2d.

To sum up, social-aware caching strategies helps to improve caching performance.

VI. DISCUSSION

We consider Social Networking has arrived to completly change the way we use the Internet.

The aim of the work, it's to show that inclusion of social networks' information is beneficial for ICN. In this context, we introduce a caching strategy based on social properties. We give an first step into social-aware caching strategies and start discussing about new ways of including social patterns and preferences into the ICN architecture.

We detail some limitations of the study. Shortest path algorithm was chosen as routing algorithm because there are not agreement about the final routing protocol for ICN, and in particular por CCN. We insist on the fact that we aim to keep the work as simple as possible, some statistical values represent simplifications of the complex reality of the interaction across the social networks.

VII. CONCLUSION & FUTURE WORK

In this paper we present a social networks-like model: users publish and retrieve content in an ICN environment. We then propose to include social networks information into caching strategies in order to improve in-network caching performances to the best of our knowledge, this is the first attempt to merge social networks and ICN.

The social-aware caching strategies show a significant improvement in terms of caching performance. In particular, Eigenvector caching strategy improves results in terms of cache hit (20%), stretch (33%) and at the same time they reduce the number of expired elements into the caches (6%). We consider the Eigenvector strategy as the first step in the research of a social-aware caching strategy for ICN.

Finally, we point out that the inclusion of social information should not be limited to ICN caching related problems. Routing alternatives may as well consider social features to select the best path. Analysis of communities might be a key-point to point ICN routing strategies, we intend to deploy routing and caching mechanisms community-oriented.

Fig. 1 :

 1 Fig. 1: Social Graph and Topological ICN Structure. The black lines represent connections of ICN nodes while red dashed lines symboline ICN nodes users' relationships.

Fig. 2 :

 2 Fig.2: Comparison between social-aware strategies and common-used ICN strategies according to cache hit, stretch, expiredelements and diversity using different cache sizes.

TABLE I :

 I Content of the Caches after publish and retrieve messages' execution.

TABLE II :

 II Eigenvector values calculated over the graph of social relationships (red dashed lines into Figure

http://mashable.com/2012/11/06/obama-wins-twitter/

http://www.fcbarcelona.com/club/barca-2-0/detail/article/leo-messi-srecord-reverberates-throughout-social-networks

[START_REF] Dannewitz | Prototyping a network of information[END_REF] http://www.peopleinside.fr/twitter-miss-france-signe-unrecord/news/

52782 4 http://www.clarin.com/politica/cuenta-Twitter-Camioneros-diapicante 0

814118751.html 5 http://www.healthleadersmedia.com/content/TEC-245418/Few-Hospitals-Use-Social-Media-Effectively-Says-Study


Eric Dalissier

Cindy Guichard

Pascal Have

Roland Masson

P Havé

C Yang

ComPASS : a tool for distributed parallel finite volume discretizations on general unstructured polyhedral meshes

come

Introduction

ComPASS, standing for Computing Parallel Architecture to Speed up Simulations, is a project focusing on the solution of time-evolution problems on general polyhedral meshes using Finite Volume (denoted FV in the following) discretizations which has been the subject of active researches over the last decades. The long-term ambition of this project is to develop a parallel prototype code in order to implement and compare on relevant problems a large class of promising FV schemes with degrees of freedom (d.o.f.) typically located at cells, faces or vertices as exhibited in the recent 3D benchmark for diffusion problems [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF]. The targeted applications are those of the INRIA project team COFFEE (COmplex Flows For Energy and Environment), such as the modelization of multiphase flows in porous media applied to oil and gas recovery in petroleum reservoirs, oil exploration in sedimentary basins, the geological storage of CO 2 in saline aquifers or the disposal of high level radioactive waste. These applications imply the simulation of multiphase Darcy flows in heterogeneous anistropic porous media, strongly coupling the mass conservation laws with the local hydrodynamical and thermodynamical closure laws and involving a large range of time and space scales. The time integration schemes are usually fully implicit in • general meshes including polyhedral cells with an arbitrary number of faces and possibly non planar faces,

• various and open request connectivities for cells, faces, and vertices, and one full layer of ghost cells (including all its cell, face and vertex d.o.f.) in order to easily compute and assemble locally on each processor a large class of compact FV discretizations,

• in order to simplify the first implementation, limit in a first step the size of the problem up to say 10 7 cells on a maximum number of one thousand processes.

One possibility would have been to fully embed the code on an existing open source software such as the library PETSc [START_REF]PETSc[END_REF] which offers a tool-kit for scientific computing with emphasis on linear systems solutions, or the Trilinos Project [START_REF]The Trilinos Project[END_REF] which is an object-oriented software framework for the solution of large-scale complex multi-physics problems, or also the object oriented distributed framework DUNE [START_REF]DUNE -Distributed and Unified Numerics Environment[END_REF], which is a modular C++ library for the solution of partial differential equations using grid-based methods. However, it was not clear that these libraries were easy to adapt in order to deal with general polyhedral meshes, and our team did not include advanced users of these libraries. Alternatively, in order to have an easy full control on the core of the code (mesh format, connectivity and distribution, FV schemes, physics, assembly of the non-linear and linear systems) we have decided to develop it from scratch and to rely on existing up to date scientific computing libraries only for specific tasks such as the partitioner, the linear solver and the visualization tools. Among well known programming languages for High Performance Computing (HPC), Fortran 90 combined with MPI has been chosen for our implementation rather than an object oriented language like C++ mainly because our team is already experienced with Fortran 90 and not with C++. The research session of the CEMRACS 2012 summer school, devoted to HPC, and joining at the CIRM around 100 researchers during 5 weeks, has been clearly a perfect framework to start the development of the ComPASS project.

Models and discretization

As stated in the introduction, the final goal of the ComPASS project is to simulate complex flows in geosciences such as compositional multiphase flows in porous media used in many applications. For example, in oil reservoir modelling, the compositional triphase Darcy flow simulator is a key tool to predict and optimize the production of a reservoir. In sedimentary basin modelling, such models are used to simulate the migration of the oil and gas phases in a basin saturated with water at geological space and time scales. The objectives are to predict the location of the potential reservoirs as well as the quality and quantity of oil trapped therein. In CO 2 geological storage, compositional multiphase Darcy flow models are used to optimize the injection of CO 2 and to assess the long term integrity of the storage. The numerical simulation of such models is a challenging task, which has been the subject of a lot of works over a long period of time, see the reference books [START_REF] Aziz | Petroleum reservoir simulation[END_REF] and [START_REF] Peaceman | Fundamentals of numerical reservoir simulation[END_REF]. Let us briefly sketch the framework of these models.

We recall here the basis of the Coats' formulation [START_REF] Coats | Implicit compositional simulation of single-porosity and dual-porosity reservoirs[END_REF] for compositional multiphase Darcy flow models. Let P denote the set of possible phases α ∈ P, and C the set of components i ∈ C. Each phase α ∈ P is described by its non empty subset of components C α ⊂ C and is characterized by its thermodynamical and hydrodynamical properties depending on its pressure P α , its molar composition C α = (C α i ) i∈C α , and on the volume fractions or saturations denoted by S = (S α ) α∈P . The phases can appear or disappear due to thermodynamical equilibrium, which leads to the introduction of the new variable Q ⊂ P representing the number of present phases at each point of the space time domain. The set of unknowns of the system in the Coat's formulation is defined by U = (P α , S α , α ∈ P, C α , α ∈ Q) and the flow is governed by the following conservation equations,

φ(x) ∂ t n i (U ) + α∈Q div (-m α i (x, U ) Λ(x) ∇P α ) = 0, for all i ∈ C, (1)
where n i denote the number of moles of the component i ∈ C per unit volume, m α i is roughly speaking the mobility of the component i in the phase α, φ(x) and Λ(x) are respectively the porosity and the permeability tensor and both are functions of the position x in space. The system is closed by the hydrodynamical and thermodynamical laws as well as the volume balance equation α∈P S α = 1. We refer to [START_REF] Eymard | Vertex-centred discretization of multiphase compositional darcy flows on general meshes[END_REF] and the references therein for a more detailed discussion about this formulation. There are two main difficulties to solve such models coupling parabolic/elliptic and hyperbolic/degenerate parabolic unknowns. First, degeneracies appear in the set of equations (1) typically at the end points of the hydrodynamical laws, and at phase transitions. Second, the heterogeneous and anisotropic nature of the media leads to possibly strong jumps of the permeability tensor Λ(x), of the porosity φ(x), and also of the mobilities m α i (x, U ) and of the capillary pressures P c,α,β (x, S) = P α -P β . FV discretizations are particularly adapted to the discretization of such models due to their conservative properties and their ability to capture discontinuous solutions. They are usually combined with a forward Euler integration in time coupling implicitly the closure laws and the conservation equations in order to obtain the stability of the scheme even for large time steps. The extension of FV discretizations to general polyhedral meshes and to anisotropic heterogeneous media has been the object of intensive research in the last two decades and has lead to numerous schemes using various d.o.f. as exhibited in the recent 3D benchmark for diffusion problem [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF]. Nevertheless, finding the 'best' FV scheme is still a problem dependent issue and the main objective of the ComPASS project is to build an appropriate tool to test promising schemes on problems with realistic physics, mesh size, and time scales, using parallel distributed computations.

A 3D mesh format for unstructured general grids

In geological applications, in order to capture the geometry and the main heterogeneities of the media, the grid generation is constrained to match the main structural features of the reservoir or basin geometry defined by the geological sedimentary layers possibly cut by a network of faults. For that purpose, the most commonly used mesh format in reservoir simulation is still the so called Corner Point Geometry (CPG) which is a globally structured hexahedral grid, locally unstructured due to erosion leading to degenerate hexahedra (collapse of edges, or of entire faces), and locally non-matching due to faults as illustrated in Figure 1. In addition, non-conforming Local Grid Refinement (LGR) can be used in nearwell regions. Although CPG grids are widely used for full field reservoir simulations, they are not adapted to modelize complex fault networks such as Y-shaped faults at basin scales, nor to discretize complex geometries at finer scales such as in the neighbourhood of complex wells or fracture networks. This justifies the growing interest to use unstructured meshes for such applications. For instance, let us consider a 3D nearwell mesh exhibited in Figure 2. The first step of the discretization is to create a radial mesh that is exponentially refined down to the well boundary to take into account the singular behaviour of the pressure close to the well. Then, this radial grid is matched to the reservoir grid using tetrahedra as well as pyramids at the interface, as can be seen in Figure 2(b). This transition zone enables to couple the nearwell simulation to a global reservoir simulation. To conclude, geological applications require to be able to use general polyhedral meshes with possibly non planar faces. Our 3D mesh format implements the mathematical definition given below. It is a conforming mesh format with a description of the cells by their set of faces and of the faces by their set of vertices given in cyclic order. It is assumed that the possible non-conformities have been solved by computing the co-refinement of the surfaces on both sides of the non-conformities generating a conforming mesh with general polyhedral cells. Our format is very close to the one used for the 3D benchmark [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF] without the description of the set of edges which is expensive in terms of memory requirement and usually not used in practice for most FV schemes. For the sake of simplicity, it will be assumed in the following that the edge data structure is not used. It could always be computed and added if required by the discretization scheme. Let us also specify that the rock properties such as the permeability tensor and the porosity will be defined as usual as cellwise constant.

Definition 1.1 (Admissible space discretization) Let Ω be an open bounded subset of R 3 and ∂Ω = Ω \ Ω its boundary. An admissible space discretization of Ω is defined by assuming that a primary polyhedral mesh M is given where each element ("commonly named cell") K ∈ M is strictly star-shaped with respect to some point x K , called the centre of K. We then denote by F K the set of interfaces of non zero two dimensional measure among the inner interfaces K ∩ L, L ∈ M, and the boundary interface K ∩ ∂Ω possibly split in several boundary faces. Each σ ∈ F K ("face of K") is assumed to be the union of triangles τ ∈ S σ . We denote by V σ the set of all the vertices of σ, located at the boundary of σ ("vertices of σ"), and by V 0 σ the set of all the internal vertices of σ. Therefore, the 3 vertices of any τ ∈ S σ are elements of V 0 σ ∪ V σ . Let us denote by

F = K∈M F K
the set of faces, and by

V = σ∈F V σ ,
the set of vertices of the mesh M.

In practice, it will be assumed that the set V 0 σ reduces to the isobarycentre of the set of vertices V σ . Hence, the mesh format reduces to a numbering of the cells M, of the faces F and of the vertices V, to the coordinates of the vertices, to the definition of the sets F K for all cells K ∈ M, and of the sets V σ in cyclic order for each face σ ∈ F. In addition, the set of neighbouring cells M σ is given for each face σ ∈ F. It is well known that, for a conforming mesh, the set M σ contains either two cells for an inner face and one cell for a boundary face.

Example of the VAG scheme on a parabolic problem

In the short term of the CEMRACS, it has been decided to implement the following simple linear parabolic equation,

∂ t u(x, t) -△u(x, t) = f (x, t), (2)
in the space-time domain Ω × (0, T ), with the initial condition u(x, 0) = u ini (x), for a.e. x ∈ Ω, and an homogeneous Dirichlet boundary condition, under the following assumptions:

• Ω is an open bounded connected polyhedral subset of R 3 , and T > 0,

• f ∈ L 2 (Ω × (0, T )) and u ini ∈ L 2 (Ω).
The implementation of this simple model already includes most of the features required for more complex models such as the compositional Darcy flow model (1): a time loop for time integration, the discretization of a second order elliptic operator that will be implemented using FV fluxes, and the solution of a linear system at each time step. The ComPASS framework aims to implement a large class of FV schemes for the discretization of such fluxes on general meshes. It includes for example the cell centred Multi-Points Flux Approximation (MPFA) schemes with in particular the well-known MPFA O-scheme introduced in [START_REF] Aavatsmark | Discretization on non-orthogonal, curvilinear grids for multi-phase flow[END_REF][START_REF] Edwards | A flux continuous scheme for the full tensor pressure equation[END_REF]. These cell centred schemes are consistent on general meshes but they are non-symmetric, and their coercivity and convergence are mesh dependent. In order to recover an unconditional symmetry and coercivity while keeping a compact flux stencil, one needs to use additional d.o.f. For instance, Mimetic schemes [START_REF] Brezzi | A family of mimetic finite difference methods on polygonal and polyhedral meshes[END_REF] or similarly the Sushi scheme [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes, sushi: a scheme using stabilisation and hybrid interfaces[END_REF] use face unknowns in addition to the cell unknowns. The DDFV (Discrete Duality Finite Volume) scheme [START_REF] Andreianov | Benchmark 3d: a version of the ddfv scheme with cell/vertex unknowns on general meshes[END_REF] uses both cell and vertex unknowns. For our test case, we have chosen to implement the VAG (Vertex Approximate Gradient) scheme, also using both cell and vertex d.o.f. It has been introduced in [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF] in the framework of the gradient schemes. We will briefly describe the main characteristics of this scheme and we refer to [START_REF] Eymard | Vertex-centred discretization of multiphase compositional darcy flows on general meshes[END_REF] for its application to multiphase Darcy flows.

Gradient schemes for a parabolic problem

The gradient schemes were recently introduced in order to obtain a uniform framework for the definition and the convergence analysis of a large family of schemes including conforming discretizations (such as finite element methods) and non conforming discretizations (such as mixed finite element and finite volume methods). The gradient schemes framework has already been applied to several problems, we mention especially here the case of the Stefan problem [START_REF] Eymard | Gradient schemes for the stefan problem[END_REF] since it includes our simple model (2) as a particular case. Gradient schemes are based on the weak primal formulation of the model which in the case of the parabolic equation ( 2) leads to the following variational formulation: find u ∈ L 2 (0, T ;

H 1 0 (Ω)) such that T 0 Ω ( -u(x, t)∂ t ϕ(x, t) + ∇u(x, t) • ∇ϕ(x, t)) dxdt -Ω u ini (x)ϕ(x, 0)dx = T 0 Ω f (x, t)ϕ(x, t)dxdt, (3)
for all ϕ ∈ C ∞ c (Ω × [0, T [).
The spirit of the gradient schemes is to replace the continuous operators by their discrete versions leading to the following definition of a space gradient discretization. 3 is the discrete gradient operator. It must be chosen such that

∇ D : X D,0 → L 2 (Ω)
• D := ∇ D • L 2 (Ω) 3 is a norm on X D,0 .
The time interval (0, T ) is discretized using, for the sake of simplicity, a constant time step denoted by δt and defined by δt = t (n+1) -t (n) with n = 0, . . . , N such that t (0) = 0 < t (1) . . . < t (N ) = T . The problem (3) is then approximated by a space-time gradient discretization using the following Euler implicit time integration scheme. Given u (0) ∈ X D,0 with Π D u (0) an approximation of u ini , the discrete solutions u

(n) ∈ X D,0 satisfy        Ω Π D u (n+1) -u (n) δt (x)Π D w(x) + ∇ D u (n+1) (x) • ∇ D w(x) dx = 1 δt t (n+1) t (n) Ω f (x, t)Π D w(x)dxdt,
∀w ∈ X D,0 , ∀n = 0, . . . , N -1.

(4)

Definition the VAG scheme

Following the Definition 1.1 of the space discretization of Ω, we also define V K as the set of all elements of V which are vertices of K. For any K ∈ M, σ ∈ F K , τ ∈ S σ , we denote by S K,τ the tetrahedron with vertex x K and basis τ as illustrated in the Figure 3.

Let us then define the gradient scheme space and operators of Definition 1.2 for the VAG scheme.

• We define X D,0 as the space of all vectors u = ( (u

K ) K∈M , (u v ) v∈V ), u K ∈ R, u v ∈ R such that u v = 0 for all v ∈ V ∩ ∂Ω.
• Then, the mapping Π D is simply defined, for any u ∈ X D,0 , by Π D u(x) = u K , for a.e. x ∈ K.

• The mapping ∇ D is defined, for any u ∈ X D,0 , by ∇ D u = ∇ Π D u, where Π D u is the continuous reconstruction which is affine in all S K,τ , for all K ∈ M, σ ∈ F K and τ ∈ S σ , with the values u K at x K , u v at any vertex v of τ which belongs to V σ , and x∈Vσ α x v u x at any vertex v of τ which belongs to V 0 σ , where the coefficients α x v , x ∈ V σ are the barycentric coordinates of v.

v ∈ V 0 σ x K Figure 3: Tetrahedral submesh.
Using the above definition of the VAG scheme, and precisely thanks to the result that for all x ∈ K, for K ∈ M, the discrete gradient ∇ D u(x) is expressed as a linear combination of (u K -u v ) v∈V K , the Problem (4) can be rewritten as follows,

K∈M K u (n+1) K -u (n) K δt w K + K∈M v∈V K F K,v (u n+1 )(w K -w v ) = K∈M 1 δt t (n+1) t (n) K f (x, t)w K dxdt, ∀w ∈ X D,0 , ∀n = 0, . . . , N -1, (5)
where F K,v (u n+1 ) is interpreted as the flux from a cell K ∈ M to one of its vertices v ∈ V K and is defined by,

F K,v (u n+1 ) = v ′ ∈V K T v,v ′ K (u (n+1) K -u (n+1) v ′ ). ( 6
)
Remark 1 Thanks to the definition of the above fluxes, the VAG scheme can be seen as a FV method. For this, we need to define some disjoint arbitrary control volumes V K,v ⊂ v∈τ S K,τ for all K ∈ M and v ∈ V K .

It is important to notice that it is not in general necessary to provide a more precise geometric description of V K,v than its volume. The mapping Π D is then defined for any u ∈ X D,0 by Π D u(x) = u K , for a.e.

x ∈ K \ v∈V K V K,v , and by Π D u(x) = u v for a.e. x ∈ V K,v . Finally, related to the Definition (6) of the fluxes, the VAG scheme can be seen as a MPFA FV method for multiphase Darcy flow models in the set of control volumes defined as the union of the cells and of the vertices. We refer to [START_REF] Eymard | Vertex-centred discretization of multiphase compositional darcy flows on general meshes[END_REF] for a more detailed discussion.

Implementation and discussion

Finally, setting in ( 5), on the one hand w v = 1 with all other values equal to zero, and, on the other hand w K = 1 with all other values equal to zero, one obtains the following discretization of Problem (2): given

u (0) ∈ X D,0 , find u (n+1) ∈ X D,0 at each time step n ≥ 0 such that K∈Mv v ′ ∈V K -T v,v ′ K (u (n+1) K -u (n+1) v ′ ) = 0 ∀v ∈ V\∂Ω, . (7a)
|K| δt (u

(n+1) K -u (n) K ) + v∈V K v ′ ∈V K T v,v ′ K (u (n+1) K -u (n+1) v ′ ) = |K|f (x K ) ∀K ∈ M, ( 7b
)
where |K| is the volume of the cell K and M v is the set of cells sharing the vertex v, defined by

M v = {K ∈ M | v ∈ V K }.
The advantage of this scheme is that it allows to eliminate without any fill-in all the cell unknowns (u K ) K∈M by using the equation (7b). The system reduces to the Schur complement on the vertex unknowns (u v ) v∈V leading to a compact vertex centred scheme. Its stencil involves only the neighbouring vertices of a given vertex, with a typical 27 points stencil for 3D topologically Cartesian meshes. Let us also point out that the VAG scheme applies to general meshes with possibly non planar faces. In addition, the VAG scheme is consistent and unconditionally coercive and has exhibited a good compromise between accuracy, robustness and CPU time in the recent FVCA6 3D benchmark [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF].

In particular, the VAG fluxes can be used for the discretization of multiphase Darcy flows leading to a very efficient scheme on tetrahedral meshes or on hybrid meshes combining typically tetrahedra, hexahedra and pyramids. For such problems on tetrahedral meshes, the VAG scheme is say 15 times more efficient in comparison with cell centred discretizations such as the MPFA O-scheme. This result can be explained by the ratio of around 5 between the numbers of cells and the number of vertices, and by the very large stencil of the MPFA O-scheme on such meshes (about 3 times larger than the stencil of the VAG scheme). In addition, the usual difficulties of vertex centred approaches to deal with multiphase Darcy flows in heterogeneous media can be solved thanks to the fact that the cell unknown is kept in the discretization and thanks to the flexibility provided by the VAG scheme in the choices of the volumes V K,v (see [START_REF] Eymard | Vertex-centred discretization of multiphase compositional darcy flows on general meshes[END_REF] for more details).

ComPASS parallel distributed processing

In this section, assuming that a mesh satisfying Definition 1.1 is given, the main steps of the algorithm for the parallel distributed solution of equation ( 7) is detailed (see also Figure 4). It will be assumed that the topology of the mesh does not change during the time evolution. Consequently, the partitioning and the distribution of the mesh on the set of processes P will be done in a preprocessing step once and for all. Also in order to define the overlaps between processes (ghost cells, faces or vertices), the FV schemes are assumed to be compact in the sense that their stencil involves only d.o.f. located at the neighbouring cells of a given item. This remark leads to our definition of the overlap as one layer of cells.

ComPASS fits into the SPMD (Single Program, Multiple Data) paradigm. To simplify the implementation, it is assumed that a master process reads and stores the mesh, computes the partitioning, and distributes the local meshes to each process. After distribution of the mesh by the master process, all the computation are done in parallel, and synchronisations of the d.o.f. at the ghost items can be performed in order to match the overlapping d.o.f. between the processes. Following our example of equation ( 7), we will describe the parallel computation of the FV scheme as well as the assembling and the solution of the linear system at each time step using the PETSc [START_REF]PETSc[END_REF] library. Thanks to our general definition of the overlap, the assembling step will be done locally on each process without communication.

• mesh reader

• mesh distribution

• messages ordering

• memory preallocation

• VAG scheme computation

• matrix and RHS assembling

•

Distributed mesh

The aim of this subsection is to distribute the mesh M on each process p ∈ P. Since our objective is to deal with a large class of FV schemes with d.o.f. possibly located at various mesh items such as cells, faces or vertices, we have chosen to partition the set of cells

M = p∈P M p . ( 8
)
The overlapping decomposition of the set of cells

M p , p ∈ P,
is chosen in such a way that all compact finite volume schemes can be assembled locally on each processor. Hence M p is defined as the set of all cells sharing a vertex with M p .

The overlapping decompositions of the sets of vertices and faces follow from this definition

V p = K∈M p V K , p ∈ P,
for the set of vertices V, and

F p = K∈M p F K , p ∈ P,
for the set of faces F. These overlapping decompositions are sketched in Figure 5. All these decompositions are computed on the master process from the initial mesh. The partitioning of the set of cells M is computed by the subroutine METIS PartGraphRecursive of the Metis library [START_REF] Metis | [END_REF] which uses a multi-constraint multilevel recursive bisection algorithm to minimize the number of edge cuts. It is provided with the number of processes card(P), and with the sets M K , for K ∈ M of cells connected to K by a vertex defined by

M K = v∈V K M v .
The Metis subroutine output is the partitioning of the set of cells [START_REF] Aziz | Petroleum reservoir simulation[END_REF]. The next step is to choose the partitioning of the remaining sets of items, i.e. the set of vertices V and the set of faces F. Let us start first with the vertices, the methodology is similar for the faces. The partitioning of the set of vertices is obtained by a rule of belonging of a given vertex to a process p ∈ P, chosen to ensure a well balanced partitioning. To be more specific, let us denote by # M K the unique global index of a cell K ∈ M.

We have the following definition.

Definition 2.1 (Partitioning of the set of vertices) A vertex v ∈ V belongs to the same process p ∈ P as the process of the cell with the smallest global index # M K among the set M v of cells sharing v. Following this rule, we obtain the partitioning V = p∈P V p with

V p = {v ∈ V | K(v) ∈ M p }, p ∈ P,
where, for all v ∈ V, K(v) = argmin L∈Mv {# M L}.

Similarly the partitioning of the set of faces F is given by the following definition.

F p = {σ ∈ F | K(σ) ∈ M p }, p ∈ P,
where, for all σ ∈ F, K(σ) = argmin L∈Mσ {# M L}.

Let us emphasize that one clearly has from the above definitions that

V p ⊂ K∈M p V K and F p ⊂ K∈M p F K .
It results, thanks to our choice of the overlaps, that the computation of the FV schemes will always be done locally on each process without communication, for all compact schemes requiring only one layer of neighbouring cells.

Communication and synchronization between processes

A basic ingredient for the parallelization of the solution algorithms is the synchronization of the d.o.f defined at the ghost items I p \ I p , for all p ∈ P, and for I = M, V or F. This synchronization is achieved using MPI (Message Passing Interface, see [START_REF] Gropp | A high-performance, portable implementation of the mpi message passing interface standard[END_REF][START_REF] Gropp | Using MPI: portable parallel programming with the message passing interface[END_REF] for basic references) Send and Receive communications between processes. To that end, the following subsets of P are computed in a preprocessing step once and for all: first the subsets of processes from which p must receive ghost items

P p,← I = {q ∈ P\{p} | I p ∩ I q = ∅}, I = M, F, V, p ∈ P, (9)
and second, the subsets of processes to which p needs to send ghost items

P p,→ I = {q ∈ P\{p} | I p ∩ I q = ∅}, I = M, F, V, p ∈ P. (10)
For each mesh item I = M, V, F, we also store once and for all, for each process p ∈ P, the partitioning

I p \ I p = q∈P p,← I I q ∩ I p , (11)
of the ghost items of p (see Figure 6 (b)), and the partitioning

I p ∩   q∈P,q =p I q   = q∈P p,→ I I p ∩ I q , (12)
of the own items of p belonging to the ghost items of q = p, q ∈ P (see Figure 6 (a)).

I p

Proc p

I p ∩ I q0 Proc q 0 ⋯ I p ∩ I qi Proc q i ⋯ I p ∩ I q (n → p,I -1)
Proc

q (n → p,I -1)
(a)

m i I p I p Proc p I q0 ∩ I p Proc q 0 ⋯ I qj ∩ I p Proc q j ⋯ I q (n ← p,I -1) ∩ I p Proc q (n ← p,I -1)
(b)

m j
Figure 6: Communication of processor p ∈ P with other processors. (a) the processor p sends message m i = (Send, p, q i ) to the processor q i ∈ P p,→ I , i = 0, . . . , n → p,I -1, with n → p,I = card(P p,→ I ); (b) the processor p receives message m j = (Receive, p, q j ) from the processor q j ∈ P p,← I , j = 0, . . . , n ← p,I -1, with n ← p,I = card(P p,←

I

).

These partitionings define the set of Send and Receive MPI messages that need to be sent and received to synchronize a given vector w p , p ∈ P with w p ∈ V I p , where V is a given vector space of d.o.f (typically R n for n a given number of unknowns per item). The synchronization provides the values of the vector at the ghost items on each process p ∈ P defining the overlapping vector w p , p ∈ P with w p ∈ V I p . For each mesh item I = M, V, F, the set of all Send and Receive messages for a given process p ∈ P is denoted by M p I = {(Receive, p, q), q ∈ P p,← I } ∪ {(Send, p, q), q ∈ P p,→ I } ,

where, for the message m = (Send, p, q), the process p sends to q ∈ P p,→ I the values of w p at the set of items I p ∩ I q , and for the message m = (Receive, p, q), the process p receives in w p from q ∈ P p,← I the values of w q at the set of items I q ∩ I p .

Definition 2.3 (Complementary message) For a given mesh item I = M, V, F, and a given message m = (mtype, p, q), we denote by m -1 its complementary message which is defined by m -1 = (mtype -1 , q, p) where Send -1 = Receive, and symmetrically Receive -1 = Send.

We have chosen to write the synchronization algorithm in a way that the execution will be deterministic and reproducible. This choice requires to use blocking point-to-point communications and an ordering of the messages to avoid possible deadlocks. This strategy has already proven its efficiency in the Arcane Framework [START_REF] Grospellier | The arcane development framework[END_REF] achieving a good scalability of the synchronization algorithm on more than 10000 processes. Following [20, Subsection 2.1], there are two situations of deadlock that can occur when using blocking pointto-point MPI communications:

(d1) for a given message, its complementary message does not exist, (d2) a Send-Receive cycle occurs, due to a wrong coordination of the Send and Receive messages.

The first case (d1) is avoided thanks our definition of the sets ( 9)-( 10)- [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes, sushi: a scheme using stabilisation and hybrid interfaces[END_REF]. Indeed, considering a mesh item I and a process p ∈ P, then, for a given message m ∈ M p I , there clearly exists a unique process q ∈ P p,← I P p,→ I such that m -1 ∈ M q I . A Send-Receive cycle (d2) arises when two (or more) processes are all waiting for their complementary messages. Let us illustrate this case in the example exhibited on Figure 7 taken from [START_REF] Luecke | Deadlock detection in mpi programs[END_REF]Subsection 2.1].

In this example P = {p i , i = 0, 3} and the set of messages [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes, sushi: a scheme using stabilisation and hybrid interfaces[END_REF] of each process is defined by M pi I = {(Send, p i , p i+1 ) , (Receive, p i , p i+1 )} for i ∈ 0, 3 , p 4 = p 0 , with I a fixed mesh item. On the subfigure 7(a), deadlocks occur since all processes are sending their send message first. The reordering of the sets (M pi I ) i=0,3 proposed on the subfigure 7(b) allows to avoid such occurrence. We have to define a proper ordering of the set of messages M p I in order to avoid deadlocks of type (d2). Roughly speaking our ordering rule will prescribe that:

(a) each process p performs its sendings and receptions in ascending order of the process ranks q ∈ P p,← I P p,→ I , (b) when two processes p and q communicate, the one with lowest rank must sends its message first.

To be more specific, let us denote by #p the rank of a process p ∈ P. Given two messages m i = (mtype i , p, q i ), i = 1, 2, m 1 = m 2 belonging to M p I , they are ordered according to the rule in the following Table 1.

1: if q 1 = q 2 = q then 2:

if mtype 1 = Send then else if mtype 1 = Receive then 9:

if #q < #p then end if 15: else if #q 1 < #q 2 then 16:

m 1 ≺ m 2 17: else 18: m 2 ≺ m 1 19: end if Table 1: Ordering of two messages m 1 = m 2 ∈ M p
I with m i = (mtype i , p, q i ), i = 1, 2 and where m i ≺ m j means that the message m i has to be executed before m j for i = j. Theorem 2.4 (Synchronisation without a deadlock situation (d2)) For a given mesh item I = M, V, F, the execution of the messages in the ordering defined in the Table 1 applied to each set M p I , for all p ∈ P, leads to a synchronisation without a deadlock occurrence in the sense of definition (d2).

Proof. Let P = {p i } i=1,n the set of processes, with card(P) = n > 2, and # p i = i for all i = 1, • • • , n. We assume, for the sake of simplicity, that all processes communicate, possibly with empty messages. This assumption simplifies the writing of the proof and corresponds in practice to a worst case assumption. Indeed, an empty message implies that its complementary message is also empty, which means in practice that no communication arises between both processes. The proof will be shown to hold even if a communication is considered in such a case. The proof is done by induction over the number of processes n, and for a given mesh item I. Thus, let us first prove the result for n = 2. The ordering of the two sets of messages, M p1 I and M p2 I is done according to the point (b) of our ordering rule (corresponding to the lines (1 -14) of the algorithm described in the Table 1). It results that the ordered sets of messages are

M p1 I = [ (Send, p 1 , p 2 ) ≺ (Receive, p 1 , p 2 ) ] and M p2 I = [ (Receive, p 2 , p 1 ) ≺ (Send, p 2 , p 1 ) ].
It is clear that at the first execution step, the process p 1 executes the message (Send, p 1 , p 2 ) and the process p 2 executes the complementary message (Receive, p 2 , p 1 ). Then, at the second step the process p 2 executes the message (Send, p 2 , p 1 ) and the process p 1 executes the complementary message (Receive, p 1 , p 2 ) which shows that the communications are achieved without deadlock. Let us now assume that the result is satisfied for n processes, and check that it is still true for n + 1 processes. We denote by M pi I , for i = 1, • • • , n, the ordered sets of messages between the first n processes which are, by assumption, well ordered to avoid a deadlock. The sets M pi I , for a given i ∈ 1, n , are the union of M pi I and of the messages (Send and Receive) between the processes p i , i ∈ 1, n and p n+1 . Since the process p n+1 has the highest rank, applying point (a) of our ordering rule (corresponding to the lines (15 -19) of the algorithm presented in the Table 1), each set of messages M pi I , for i = 1, • • • , n, is ordered starting with the ordered set M pi I . Then, using the point (b) of our ordering rule, we obtain the ordering,

M pi I = [ M pi I ≺ (Send, p i , p n+1 ) ≺ (Receive, p i , p n+1
) ], ∀i ∈ 1, n . Similar arguments lead to the following ordering of the set of messages of the process p n+1 ,

M pn+1 I = [ (Receive, p n+1 , p 1 ) ≺ (Send, p n+1 , p 1 ) ≺ • • • ≺ (Receive, p n+1 , p n ) ≺ (Send, p n+1 , p n ) ].
It results from our induction assumption that the messages ( M pi I ) i=1,n will be executed first without deadlock. Next, from the above ordering of the messages, for i = 1, • • • , n in this order:

• the process p i executes the message (Send, p i , p n+1 ),

• the process p n+1 executes the complementary message (Receive, p n+1 , p i ),

• the process p i executes the message (Receive, p i , p n+1 ),

• the process p n+1 executes the complementary message (Send, p n+1 , p i ).

We conclude that all the communications are achieved without deadlock.

Computation of the Finite Volume scheme

From our choice of the overlaps, the computation of the FV scheme can always be done locally on each process without communications, for all schemes requiring only one layer of neighbouring cells. Let us consider our example of the computation of the VAG fluxes [START_REF] Aavatsmark | Discretization on non-orthogonal, curvilinear grids for multi-phase flow[END_REF]. For all process p ∈ P, we compute the transmissivities T v,v ′ K , v, v ′ ∈ V K for all K ∈ M p . In our case, the diffusion coefficient is fixed to 1 and the transmissivities can be computed only once in the simulation, before the time-loop. Similarly, the cell volumes |K| and the cell centres x K are also computed locally for all K ∈ M p and for each process p ∈ P. speaking a minimum of 50000 cells for each process. This explains the decrease of the parallel efficiency in our case over a few tens of processes. All together, the results of the two test cases exhibit a rather good strong scalability of the code ComPASS which validates our implementation.

Conclusion and perspectives

The CEMRACS has offered an ideal framework to start in a collaborative project the development of the new code ComPASS adapted to the distributed parallelization of Finite Volume schemes on general unstructured meshes using various d.o.f. including typically cell, face or vertex unknowns. The parallelization is based on a classical domain decomposition approach with one layer of cells as overlap allowing to assemble the systems locally on each process. The code is adapted to fully implicit time integrations of evolution problems thanks to the connection with the PETSc library in order to solve the linear systems.

To simplify the implementation, at least in a first step, it has been decided to read and partition the mesh on a master process which distributes it to all processors in a preprocessing step. All the remaining computations are done in parallel. The code has been tested on a parabolic toy problem to evaluate its parallel efficiency on up to 500 processes. The results exhibit a rather good strong scalability of the CPU time with a typical behaviour of the parallel efficiency for implicit discretizations on a fixed mesh size. This is clearly only the first step of development of a more advanced version of the code. In the near future, we intent first to improve the parallelization of the post-processing step of the code, and of course to include more physics starting with multiphase porous media flow applications. Last but not least, this CEMRACS project has also offered to the youngest authors, the opportunity to learn MPI parallel programming and to work in a collaborative project using appropriate software engineering tools.

Figure 1 :

 1 Figure 1: Examples of general cells due to geological erosion (left part) or faults (the last at the right).

  (a) exponentially refined radial mesh (b) hybrid mesh with hexahedra, tetrahedra and pyramids

Figure 2 :

 2 Figure 2: Space discretization of a deviated well.

Definition 1 . 2 (

 12 Space gradient discretization) A space gradient discretization of the problem (3) is defined by D = (X D,0 , Π D , ∇ D ) where, 1. the set of d.o.f. X D,0 is a finite dimensional vector space on R suited to the Dirichlet homogeneous boundary condition, 2. the linear mapping Π D : X D,0 → L 2 (Ω) is the reconstruction of the approximate function, 3. the linear mapping

Figure 4 :⋯Proc p n- 1 I pn- 1 I pn- 1 I

 4111 Figure 4: Main steps of the algorithm for the parallel distributed solution of equation (7).

Figure 5 :

 5 Figure 5: The overlapping decompositions of the sets of cells, vertices and faces. I = M, V, or F; n = card(P).

Definition 2 . 2 (

 22 Partitioning of the set of faces) A face σ ∈ F belongs to the same process p ∈ P as the process of the cell with the smallest global index # M K among the set M σ of cells sharing σ. Following this rule, we obtain the partitioning F = p∈P F p with

Figure 7 :

 7 Figure 7: Example of communications between four processes.

  (a) Cartesian grid (b) Randomly perturbed grid

Figure 9 :

 9 Figure 9: Speed-up study using two types of topologically Cartesian grids with 128 3 cells and up to 500 processes. Speed Up : Sp (n proc) = Time (1 proc) /Time (n proc) ; linear Speed up Sp (n proc) = n.

Numerical validation on a parabolic modelIn this section, the strong scalability of ComPASS is studied on the discrete problem described in subsection 1.2 using the parallelization of the algorithm described in section 2. The mesh is a topologically Cartesian grid of size fixed to 128 3 cells, and the code is run on a number of processes ranging from 1 to 500. To be more specific, two types of grid have been tested. The first one is a Cartesian mesh and the second one is obtained

Acknowledgement

The authors would like to thank the organisers of CEMRACS 2012 and especially Pierre Navaro for its precious help in using the package PETSc and running our code on the cluster of Marseille.

this author is partially supported by European Research Council ERC starting Grant 2009, project 239983-NuSiKiMo.

1

Assembly of the linear system and linear solver

Thanks to our choice of the overlap, and to the above computation of the tranmissivities, cell volumes and cell centres for all cells K ∈ M p , the following subset of the set of equations ( 7) is assembled on each process p ∈ P without any communications between processes:

-

|K| δt (u

Note that it includes the equations of all own vertices v ∈ V p , and the equations of all own and ghost cells K ∈ M p . This set of equations can be rewritten as the following linear system (the time superscript (n) is removed for the sake of clarity),

denote the vectors of discrete unknowns at vertices and cells, b p ∈ R M p the right hand side, and where the matrices have the following sizes:

Moreover, since the block D p is a non singular diagonal matrix, the Schur complement can be easily computed without fill-in to reduce the linear system to the vertices unknowns only as mentioned in the presentation of the VAG scheme in the subsection 1.2. We obtain

and the Schur complement system

Finally on each process p ∈ P, we construct the matrix

This computation requires no MPI communications thanks to our choice of the overlap. It is transferred line by line to PETSc. Next, the parallel linear system obtained from [START_REF] Gropp | A high-performance, portable implementation of the mpi message passing interface standard[END_REF] for all p ∈ P, is solved using the GMRES PETSc Krylov linear solver with a block Jacobi ILU2 preconditioner. PETSc provides the solution vector U V p ∈ R V p at own vertices for each process p ∈ P. The solution at own and ghost vertices U V p , p ∈ P is obtained by a synchronization of the vertices, and the solution at own and ghost cells U M p is computed from ( 16) locally on each process p ∈ P.

from the previous one by a random perturbation of its nodes leading to hexahedral cells with non planar faces.

For the sake of clarity, a coarser randomly perturbed grid of size 32 3 is exhibited on the Figure 8(a), and its partitioning is presented on the Figure 8(a). The time interval is given by [0, 0.05] and is discretized using a constant time step of 0.001.

This numerical example has been run using the OpenMPI library version 1.4.5 on the meso-centre of Aix-Marseille University which provides a cluster architecture of 96 nodes each equipped with 12 cores and connected by an InfiniBand QDR communications link (see [START_REF]Meso-centre of Aix-Marseille university[END_REF] for more details). The Figure 9 exhibits the speed-up obtained with both types of grids with the following splitting of the run times:

(i) computation of the VAG transmissivities as discussed in the subsection 2.3, (ii) creation of the PETSc objects, assembly of the matrix and of the right-hand side (denoted RHS), and computation of the Schur complement as described in the subsection 2.4, (iii) computation of the discrete solution, including both the call to the linear solver of PETSc, and the calculation of the cell unknowns using [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF].

We observe on the Figure 9 that the computation of the VAG transmissivities (i) satisfies as expected the optimal linear speed-up since all the computations are done without communication between the processes. For the two other steps (ii) and (iii), the parallel efficiency decreases clearly with the number of processes. This can be explained by the cost of the communications needed for each step. Indeed, for step (ii), one synchronization of the global index of the local vertices is required for the PETSc objects creation, and the synchronization of the discrete solution at ghost vertices is required at the end of each time step. For step (iii), the linear solver involves communications for the matrix vector product and for the scalar products at each solver iteration. It results that for that type of implicit algorithm, the speed up can only be optimal for roughly