Guillaume Sahut

Giovanni Ghigliotti

Philippe Marty

Guillaume Balarac

Evaluation of Level Set reinitialization algorithms for phase change simulation on unstructured grids

EVALUATION D'ALGORITHMES DE REINITIALISATION DE LA LEVEL SET POUR LA SIMULATION DU CHANGEMENT DE PHASE SUR MAILLAGES NON STRUCTURES

Guillaume Sahut (1) , Giovanni Ghigliotti (1) Philippe Marty (1) , Guillaume Balarac (1) (1) Univ. Grenoble Alpes, CNRS, Grenoble INP, LEGI, 38000 Grenoble, France -e-mail: guillaume.sahut@univ-grenoble-alpes.fr Dans cette étude, on présente différentes méthodes numériques pour la simulation du changement de phase liquide-vapeur (ébullition). On utilise un formalisme Level Set pour capturer l'interface liquide-vapeur. Un tel formalisme nécessite une étape de réinitialisation de la fonction Level Set après advection. Cette étape est critique pour la simulation du changement de phase car elle ne doit ni déplacer l'interface, ni introduire de déformations dans le profil de la fonction Level Set, sous peine de détériorer la précision du calcul de la normale à l'interface et de sa courbure, nécessaires pour définir respectivement la vitesse de l'interface due au changement de phase et le saut de pression à l'interface. On présente d'abord les équations résolues et le couplage des équations de Navier-Stokes avec le taux de transfert de masse modélisant le changement de phase. Puis on détaille différents algorithmes de réinitialisation de la Level Set pour la simulation numérique de l'ébullition, sur maillages structurés et non structurés. Ces méthodes sont ensuite validées par un cas-test de croissance de bulle statique à taux de transfert de masse fixé. En particulier, on observe qu'à l'instant correspondant au doublement du rayon de la bulle, ce dernier converge en maillages pour toutes les méthodes présentées. Le test concluant sur maillages non structurés ouvre la voie à la simulation du changement de phase liquide-vapeur dans des géométries complexes.

MOTS CLEFS : écoulements diphasiques, changement de phase, maillages non structurés, Level Set, ébullition

Evaluation of Level Set reinitialization algorithms for phase change simulation on unstructured grids

In this study, we present different numerical methods for the simulation of liquid-vapor phase change (boiling).

We use a Level Set formalism to capture the liquid-vapor interface. Such a formalism requires a reinitialization (aka redistancing) step of the Level Set function after advection. This step is critical for phase change simulation as it must neither move the interface nor induce perturbations in the Level Set function, otherwise the normal vector to the interface and its curvature, two quantities that are crucial to define respectively the interface velocity due to phase change and the pressure jump at the interface, would be in turn too much perturbed. Here we present a comparison of different reinitialization algorithms of the Level Set function for boiling simulation, on structured and unstructured grids. These methods are then validated against the analytical case of a static growing bubble with a fixed mass transfer rate. In particular, we observe that at the time corresponding to a doubled bubble radius, the error on the bubble radius decreases with the grid cell size for all presented methods. KEYWORDS: two-phase flows, phase change, unstructured grids, Level Set, boiling

I INTRODUCTION

Two-phase flows are encountered in a wide range of industrial applications such as heat exchangers, nucleate boiling or spray cooling. The particularity of two-phase flows is the existence of an interface between the two phases. Numerical simulations of two-phase flows require the localization of the interface. Two-phase flow simulations including liquid-vapor phase change are even more challenging, as the motion of the liquid-vapor interface depends also on the mass transfer rate. There exist several numerical methods to keep track of the interface. In this study we describe different versions of the Level Set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] on both structured and unstructured grids. These methods are then validated and compared on the analytical case of a static growing bubble with an imposed mass transfer rate.

II GOVERNING EQUATIONS

We solve the incompressible Navier-Stokes equations

                      1 p u u t u , (1)
where u  is the velocity, 𝑝 the pressure, 𝜌 the density and Σ the viscous stress tensor of the given phase. The discontinuity in the velocity field at the interface is given by

  , 1 n m u               (2)
where   vap liq q q q    denotes the discontinuity, or jump, of the quantity q across the interface  , m  is the mass transfer rate (in kg m -2 s -1 ) responsible for phase change, and n  is the interface normal vector pointing towards the liquid phase. When solving (1) we also have to take into account the pressure jump at the interface given by

              1 2 m p  , ( 3
)
where  is the surface tension and  the curvature of the interface. Both jumps ( 2) and ( 3) are used in the projection method to solve (1) [START_REF] Tanguy | -Benchmarks and numerical methods for the simulation of boiling flows[END_REF].

In this study, as we focus on the different Level Set approaches available to accurately capture the interface, we assume that the mass transfer rate is fixed. These methods have been implemented in the YALES2 solver [START_REF] Moureau | -Design of a massively parallel CFD code for complex geometries[END_REF] and are presented in the next section.

III LEVEL SET METHODS ON STRUCTURED GRIDS

The most challenging task in two-phase flow simulations is the accurate localization and advection of the interface at every time step. Perturbations in interface localization result in a loss of precision and, in the more severe cases, in physical aberrations. With phase change, this well-known problem is even more restrictive.

There are several methods to represent the interface. In this work, the Level Set method is used. A Level Set function is a function seen as a set of iso-levels. The liquid-vapor interface is identified as one specific isolevel [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. The Level Set method has two major advantages. First, the advection of the Level Set function in all the computational domain enables the implicit capture of the interface. Second, geometrical properties such as normal vector and curvature of the interface are embedded in the Level Set field. The main drawback of the Level Set method is the mandatory need of a reinitialization step of the Level Set function.

III.1 Case 1 -Signed Distance Function reinitialized by a Hamilton-Jacobi equation

We start our investigations by the method described by [START_REF] Tanguy | -Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. In this method, the Level Set function is the Signed Distance Function  to the interface  defined for any x  in the computational domain

 by        x x x x     min ) (  . ( 4
)
In this case, the interface is identified as the 0 iso-level of the Level Set function [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF].

The Signed Distance Function is advected by solving the standard advection equation

, vap vap m u t             (5)
where vap u is the vapor phase velocity, vap  is the vapor density and the source term on the right hand side is due to phase change. After advection, the function  is reinitialized as a Signed Distance Function by solving the Hamilton-Jacobi PDE

   , 0 1 0            S (6)
where  is a pseudo-time, S is a smoothed sign function defined by [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF], and 0  is the previously advected Level Set that needs to be reinitialized. The equation ( 6

        n and . n        (7)
Note that high precision is needed in the reinitialization of  to avoid perturbations in the normal vector and curvature of the interface. To this purpose, the term    in ( 6) is computed using a high-order scheme such as the Fifth Order WENO scheme [START_REF] Jiang | Weighted ENO Schemes for Hamilton-Jacobi Equations[END_REF]. High-order schemes require higher computational cost, are difficult to implement on unstructured grids and may reduce performance in a parallel code.

III.2 Case 2 -Signed Distance Function reinitialized by the Fast Marching Method

Another method for the reinitialization of the Signed Distance Function is the Fast Marching Method [START_REF] Sethian | A Fast Marching Level Set Method for Monotonically Advancing Fronts[END_REF]. The Fast Marching Method, based on the solution of an Eikonal equation, can be seen as the stationary version of the Hamilton-Jacobi equation ( 6) and is given by .

1    (8)
The solution of equation ( 8) is based on the propagation of the Signed Distance Function values from the interface along the normal direction to the interface. To limit the computational time needed to perform the Fast Marching Method, we solve it only on a narrow band around the interface, large enough to be able to compute the normal vector and the curvature by equations (7).

III.3 Case 3 -Conservative Level Set on structured grids

In order to improve mass conservation in each phase, [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] proposed the Conservative Level Set method. In this method, the Level Set function  is a smeared-out Heaviside function defined for

  x  by           ) ( 2 ) ( tanh 2 1 2 1 ) ( x x x       , ( 9
)
where  is the Signed Distance Function and  is a scale parameter roughly of the order of the grid cell size.

In this case, the interface is identified as the 0.5 iso-level of the Level Set function. The Conservative Level

Set is advected by solving the conservative form of the standard advection equation with source term

  , vap PC vap u u u t                       (10)
where the source term is composed of the phase change contribution and the divergence of the vapor velocity, which should be null, but it is recommended to account for this correction to decrease the sensitivity to errors in the evaluation of the velocity field. In (10), the interface velocity due to phase change

PC u  is given by , n m u vap PC      and the gradient of the Level Set    is computed as , ) ( 2 ) ( cosh ) ( 4 1 ) ( 2 n x x x x                     (11)
where  is again the Signed Distance Function. This technique to compute    is reused from the reinitialization step suggested by [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] and presented below. The Conservative Level Set embeds the interesting property of volume conservation. The normal vector and the curvature of the interface are given by a similar approach as in (7).

The reinitialization equation for the Conservative Level Set proposed by [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] is

       , 1 n n n                        (12)
where the normal vector is given by

        n . ( 13
)
Nevertheless, the necessity to set  to a small value to improve volume conservation produces sharp gradients in  and thus potential oscillations in n  . To avoid this problem, and with the aim of increasing accuracy, [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] built a new reinitialization equation for the Conservative Level Set given by

                                n n map map     1 2 cosh 4 1 2      , ( 14
)
where the inverse of the Conservative Level Set map  is given by

              1 log map . ( 15
)
The normal vector is computed based on the Signed Distance Function, derived from the Conservative Level Set function by using the Fast Marching Method:

FMM FMM n         , (16)
where FMM  is the Signed Distance Function given by the Fast Marching Method using boundary conditions on the closest nodes to the interface.

IV LEVEL SET METHODS ON UNSTRUCTURED GRIDS

IV.1 Case 4 -Signed Distance Function on unstructured grids

To be able to address complex geometries, algorithms are now extended on unstructured grids. The first unstructured case is similar to Case 1. Only the reinitialization step of the Signed Distance Function is replaced by the method developed on unstructured grids by [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. We make use of the MshDist library implementing this method.

IV.2 Case 5 -Conservative Level Set on unstructured grids

Here we want to extend the methodology of [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] described for structured grids in Section III.3 to unstructured grids. The only difference is the replacement of FMM  by the Signed Distance Function given by the method proposed by [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF].

V NUMERICAL RESULTS

We validate the different Level Set reinitialization methods presented in the previous section on the case of a 2D static growing bubble with a fixed mass transfer rate from [START_REF] Tanguy | -Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. The initial bubble radius  10 -5 kg m -1 s -1 . In

VI CONCLUSION AND PERSPECTIVES

In this study, several numerical methods for the advection and reinitialization of the Level Set function have been presented in the context of liquid-vapor phase change. These methods have been validated and compared on the case of a static growing bubble with a fixed mass transfer rate. All methods present a convergence rate with the grid cell size close to one. Two of the five methods presented are designed for unstructured grids. The ability to accurately compute the interface position allows the quantitative simulation of liquid-vapor phase change on unstructured grids. This opens the path to numerical simulations of liquid-vapor phase change on complex geometries.

The main perspective of our work is the simulation of phase change on unstructured grids with a computed mass transfer rate which depends on the thermal fluxes at the interface and on the latent heat of the fluid.

VII

  ) is solved in pseudo-time until convergence, i.e. until 1    , which is part of the definition of the Signed Distance Function. The normal vector to the interface n  and the curvature of the interface  are then given by

is 0 R

 0 = 10 -3 m and the imposed mass transfer rate is m  = 10 -1 kg m -2 s -1 . The simulations are performed until final time f t = 10 -2 s needed for the bubble radius to double the initial radius. The other physical parameters of interest are liq  = 10 3 kg m -3 , vap  = 1 kg m -3 ,  = 7  10 -2 N m -1 , liq  = 10 -3 kg m -1 s -1 and vap  = 1.78

Fig. 1 .

 1 Fig. 1. The results for the Case 5 on four different unstructured grids are plotted at final time with a characteristic cell size of (a) 4  10 -4 m, (b) 2  10 -4 m, (c) 1  10 -4 m and (d) 5  10 -5 m. The initial interface is plotted in blue, the computed interface in black, and the theoretical interface in red. For clarity, only the coarsest grid is represented in (a). The computed liquid and vapor velocity fields are plotted in (b).

Fig 2 .

 2 Fig 2.The normalized L ∞ norm of the error on the bubble radius  is plotted at final time for the five cases on four different grid cell sizes (a). The analogous error for the normal vector (b) and for the curvature (c).

Table 1 .

 1 Table1, the methods used in all cases are summarized. Summary of all cases presented in the previous sections. SDF stands for Signed Distance Function, CLS for Conservative Level Set, HJ for Hamilton-Jacobi and FMM for Fast Marching Method.

	Grid topology		Structured	Unstructured
	Level Set type	SDF	CLS	SDF	CLS
						[Dapogny and
	Reinitialization method	HJ eq.	FMM	[Chiodi and Desjardins, 2017]	[Dapogny and Frey, 2012]	Frey, 2012] + [Chiodi and
						Desjardins, 2017]
	Cases	Case 1	Case 2	Case 3	Case 4	Case 5

ACKNOWLEDGMENTS This work has been supported by the LabEx Tec21 (Investissements d'Avenir-Grant Agreement No. ANR-11-LABX-0030). Vincent Moureau and Ghislain Lartigue from the CORIA lab and the SUCCESS scientific group are acknowledged for providing the YALES2 code. G. B. is also grateful for the support from Institut Universitaire de France. The authors would like to thank Patrick Begou for support on algorithms implementation and Charles Dapogny for support on the MshDist library.


Cedric Deffo Sikounmo
email: cedric.deffo-sikounmo@univ-smb.fr

Stephane Perrin
email: stephane.perrin@univ-smb.fr

Eric Benoit
email: eric.benoit@univ-smb.fr

Uncertainty management of situations in a housing use context

This paper proposes a new contribution to manage uncertainty in the context of a situation recognition. It gives an application to the measurement of housing use situations. It extends a previous communication that proposes an ontological based model of situations. In this model, a situation is a set of entities linked by relationships themselves determined using sensors. However, such model doesn't allow to manage uncertainty. In order to take into account an uncertainty related to the identification of some entities, we propose to model the set of possible eligible entities, called the universe of discourse, by an anonymous instance in the ontology. This anonymous instance is linked with a "same as" relationships with all possible eligible entities. In addition, we propose to complete this model by the association of a mass with each "same as" relationship. This paper presents a probabilistic semantic to these masses. A simple didactic example on a housing use situations measurement illustrates the presented approach.

I. INTRODUCTION

The Iot field sees the raise of sensible objects that are able to perceive their close environment. By the way, in a community of sensible objects, any object might be able to perceive the other close objects or more generally any physical or biological entity close to it. As a projection of the global state of the entity community on a subset of entities, a situation gives a higher vision of a given environment. In this paper, we focus on the recognition of situations in a housing use context. In this case the biological entities are peoples, children or pets and the physical objects are furniture units and household appliances. The sensible objects are then intelligent connected furniture units [START_REF] Sikounmo | States measurement in a context of intelligent connected furnitures[END_REF] and connected appliances. The goal of the situation recognition in a housing use context is to identify the human uses of the house in order to improve their day-to-day experience. The application fields of the presented study are the assistance to people not only for comfort but also for home support of elderly persons and disabled persons. Indeed, the dynamic identification of human uses allows to detect changes in daily routines as it appends when an underlying medical or health issue appears.

In a previous study, we presented the advantage of representing the global state of the entity community with a factual ontology [START_REF] Sikounmo | Modeling Situations in an Intelligent Connected Furniture Environment[END_REF]. With this approach, a situation is considered to be realized when the ontology that defines this situation is aligned with the factual ontology. In other words, it is realized if one or more instances of this situation definition exists in the factual ontology that describes the global state. This approach supposes that the facts that define the global state are perfectly known. But in a context of a community of objects including sensible objects, some facts are deduced from objects perceptions and are rarely perfectly known. Indeed, as a fact is deduced from measurement results given by low cost sensors. In this case the uncertainty introduced by the measurement is significant. The identification of the situation is provided by a fusion process that uses these measurement results as inputs. Therefore the management of uncertainty is a critical point in the situation identification process.

The final goal of the study presented in this paper is to produce an information on the uncertainty on the realization of a situation. A first chapter will remind the concept of ontologies and its application to the definition of situations then its application to the identification of the realized situations. A second chapter will present the limits of the usual approaches used to consider uncertainty in the factual ontologies. The third chapter is devoted to the presentation of our approach including the method to represent a universe of discourse on the factual ontology and the method to model the uncertainty according several theories.

II. SITUATION IN A HOUSING USE CONTEXT

Let us define an environment made of humans and of household objects including smart connected furniture units. Let us note the various entities describing or influencing the different states of our environment. As an example, a living room can be chosen as a context for the environment. In the same example, the room can be made up of a sofa, chairs, a coffee table, a TV cabinet, a television and living occupants. We remind that the intelligent connected furniture units are sensible objects that hold embedded sensors: the sofa and the chairs may have sensors that allow them to detect if a human sit on them or that allow to perceive the proximity of other objects or of humans. It is also assumed that the sofa dialog with the coffee table itself in communication with the television, which allows an occupant to manage the TV from the sofa. Let admit that the following states can be detected: "occupants watch a movie", "an occupant sits on the couch", "an occupant is watching TV", "occupants are discussing". From the above, we can describe our environment as a community of actors of various natures who interact. We first propose to split the actor category into the human category and the object category. The object category generalizes the smart connected object category that generalizes the smart connected furniture unit category. To find the best way to represent the knowledge of our environment, we need a representation that must be structured, portable, scalable and dynamic while keeping the semantics that corresponds to the environment. To perform this representation we choose the ontologies that already apply in different fields of research namely: knowledge engineering, automatic natural language processing, cooperative information systems and many others.

A. Ontology

The notion of ontology associated with the artificial intelligence appeared early 1990s. The most cited definition is the Gruber's one [START_REF] Gruber | A translation approach to portable ontology specifications[END_REF], he defines an ontology as "an explicit specification of a conceptualization". A conceptualization is a common vision of an object of the real world that can be physical or abstract. In 1997, the definition given by [START_REF] Borst | Construction of engineering ontologies for knowledge sharing and reuse[END_REF] is more refined and specifies that the conceptualization must be shared. Conceptualization should express a shared vision of the described object. Therefore that can easily be interpreted by software. The years after present an ontology as the result of a modeling by emphasizing the characterization of primitives for the representation of knowledge. Finally, in [START_REF] Bouquet | C-OWL: Contextualizing Ontologies[END_REF], Bouquet et al give the following definition: "Ontologies are shared models of a domain that encode a view which is common to a set of different parties". The various and complementary views of these definitions can be noted. Thus, an ontology defines a common and shared vocabulary for a particular domain. This allows communication between agents. The ontologies are designed to facilitate the exchange and reuse of vocabularies.

Hence, the parts of an ontology are:

• classes modeling the entities of the real world. Example: (P aul is_instance_of P erson). In addition, an ontology must satisfy the principles of clarity, objectivity, completeness, consistency, extensibility and modularity as described in [START_REF] Gómez-Pérez | Ontological engineering: A state of the art[END_REF]. A language must be adopted to describe the ontologies in a homogeneous way. The most common one are currently RDFS and OWL (Web Ontology Language) [START_REF] Staab | Handbook on Ontologies[END_REF]. In [START_REF] Horrocks | From SHIQ and RDF to OWL: The making of a Web ontology language[END_REF], the authors briefly describe the evolution of these descriptive languages until the completion of the OWL languages family.

B. Situation recognition

In this study, a situation is defined as a state of a given set of physical or biological entities called "actors of the situation". A situation has an abstract view that we call the reference situation. The situation recognition is then the identification of a realization of the reference situation in the real world.

1) Reference situation: In [START_REF] Sikounmo | Modeling Situations in an Intelligent Connected Furniture Environment[END_REF], a reference situation is then defined by the actor types (called concepts), and by a set of relations between these concepts.

2) Realized situation: A realized situation is the instance of a reference situation. It's a state involving a limited set of actors. In other words, it's a state involving the actors defined in the reference situation. A situation is identified as realized if the factual ontology includes the instances of the concepts and of the relationships that define the reference situation.

Thus, the truth level of the realization of a situation directly depends on the truth level of the relationships defined in the reference situation. As the thruth level of the relationships is given by the measurement performed by the sensible objects, it may hold an uncertainty that must be managed in the ontology.

III. A SHORT RESUME OF UNCERTAINTY MANAGEMENT ON ONTOLOGIES

A. Uncertainty management

Given a smart connected furniture unit embedding sensors of different natures, these sensors send various information on quantities like temperature, pressure, vibration, location. However, the reliability of this information is not always assured. The imperfections of the information can be:

• ambiguity when the returned information doesn't enable to distinguish two different states, • noise that holds a non-pertinent part of the information,

• bias that reflects the unknown of potential information sources, • incompleteness when some information is missing, • inaccuracy, • uncertainty, • inconsistency, • redundancy. In order to consider these imperfections, at least the uncertainty must be managed on the representation introduced section II-A.

B. Uncertainty management on ontologies

Regarding the advance of domains such as data acquisition and information fusion, it was necessary to find ways to associate the logical uncertainty management formalisms with the knowledge representation methods. Indeed, there are several formalisms to manage incomplete information as well as information related uncertainty. Some studies aimed at extending the OWL language and the description logic to the mathematical formalisms supporting uncertainty. This is reflected in the advent of languages such as: PR-OWL, FuzzyOWL, BayesOWL. See for example [START_REF] Lukasiewicz | Managing uncertainty and vagueness in description logics for the semanticweb[END_REF] for a survey on the methods and languages. In the literature, the probability theory is the most used for the uncertainty management. In 2006, Da Costa [START_REF] Da Costa | Pr-owl: A bayesian ontology language for the semantic web[END_REF] proposed a first version for PR-OWL: a probabilistic generalization of the OWL language. In PR-OWL the semantics is based on Multi-Entity Bayesian Networks (MEBN) [START_REF] Laskey | Mebn: A language for first-order bayesian knowledge bases[END_REF]. The MEBN are made of fragments denoted Mfrags that represent the probable information related to a set of random variables. However, this version specifically focuses on the MEBN semantics. Hence the proposal of PR-OWL2 in [START_REF] Carvalho | Pr-owl 2.0 -bridging the gap to owl semantics[END_REF] which places more emphasis on expressiveness and on OWL semantics. The PR-OWL approach is to create a meta-ontology where uncertain concepts are defined with superclasses. In the same context, we can find BayesOWL [START_REF] Ding | Bayesowl: Uncertainty modeling in semantic web ontologies[END_REF]. This kind of ontology, based on Bayesian networks, gives a solution to combine knowledge with observations and a statistic learning theory. BayesOWL gives a set of rules and processes to simplify the migration from OWL to RBs while maintaining the semantics expressed on the initial ontology. Finally, FuzzyOWL presented in [START_REF] Stoilos | Fuzzy OWL: Uncertainty and the semantic web[END_REF] are ontologies based on the fuzzy subset theory with an extension to the description logic named fuzzyDL. In this case, each concept of the ontology is associated to a "fuzzy concept". Properties and axioms are also respectively associated to their respective fuzzy versions.

C. Why a new proposal?

The advantages of these ontology languages lies in their ability to reason. Indeed, the theories of mathematical logics on which they based have proved their advance in the field of modeling, representation and uncertainty management. These studies try to make classes and concepts uncertain but in our approach the concept ontology is stable and only its instances, the factual ontologies, are subject to the uncertainty. Indeed, in our approach the facts are permanently moving and the factual ontology reflects our knowledge on these facts including the uncertainty on the facts truth.

IV. A NEW PROPOSAL

In the field of description logics, knowledge modelling of a domain is performed on two levels. The first level corresponds to the TBox statement that describes the domain general knowledge. The second level corresponds to the ABox statement, also called factual level, that represents a fact, a kind of "instance" related of one configuration. A TBox component includes definition of concepts and roles, and an ABox component describes each individuals by naming them and specifying assertions on individuals that are named using concepts and roles. Several ABox can be associated with a same TBox; each ABox represents a configuration based on a set of individuals and uses concepts and roles included in the TBox. In this paper, we propose to introduce uncertainty. In our case, the uncertainty comes from the measurements. Measures are instanced and therefore are modeled inside ABox. This section describes our proposal to model uncertainty into ABox.

A. Problem definition

As shown in section 1, a furniture environment is described by a factual ontology. This ontology describes both static knowledge and dynamic one. Our goal is to retrieve realized situations of our environment and to associate a truth level to each realized situation. According to this goal, the uncertainty management has to be applied on this ontology based representation. Before defining any meaning for this truth level we will simply call it mass. Depending on the theory chosen to model the uncertainty, the mass will represent for example a probability, a possibility or a belief mass.

A situation corresponds to a graph where relations corresponds to services given by nodes. Nodes represent a piece of furniture or other objects that constitute our environment. For example, suppose that the presence of two humans and a chair is detected in a room. Sensors equip chair allowing it to deliver its state: human is on itself, or is not. Fig. 1 illustrates the ontology of such environment. Suppose that a chair detects that a human is sit on it. On one hand, we want to know which human is sitting on it. More precisely, we want to know the mass associated to the situation "a human is sitting on" for each human candidate. On the other hand, the situation "a human is sitting on" can be uncertainty too. This uncertainty is also represented by another mass. In our precedent work [START_REF] Sikounmo | Modeling Situations in an Intelligent Connected Furniture Environment[END_REF], we proposed a model for the representation of situations. But this model does not take into account uncertainty.

B. Representing a universe of discourse in a factual ontology

The preliminary action to take uncertainty into account is to define the so called universe of discourse that specifies the set of eligible candidates. Fig. 1 illustrates that Human_1 is sitting on chair_1 and represents a certain knowledge. In the real case, the only certain knowledge is that Human_1 and Human_2 are possibly sitting on Chair_1. We need to express in the one hand that a human is sitting on Chair_1. In the other hand that Human_1 and Human_2 are candidates. To express the first need, we propose to use the sitting_on relation to link an anonymous instance with Chair 1 as illustrated in fig. 2. At this step, no individual fits to the anonymous instance. We just modeled on the factual ontology the fact "something is sitting on Chair_1". To express the set of candidate individuals, we propose to use the owl:sameAs relation to link the anonymous instance to all candidate individuals as shown in fig. 3. The universe of discourse we propose is then made of all owl : sameAs relationships related to the same anonymous instance. As such knowledge is produced by sensor measurements, the universe of discourse is expressed on the factual ontology (i.e. the ABox).

This universe of discourse can be formalized using a RDF format triplets:

< Inst i rdf : type Concept i , Relation k , Inst j rdf : type Concept j >
where rdf : type is the RDF relation that links an instance to its concept.

This triplet expresses that:

• Inst i is an instance of Concept i ,
• Inst i is linked to Inst j throw the relationship Relation k , • Inst j is an instance of Concept j . An anonymous instance without type, as shown in fig. 3 is represented by (_ :?).

An anonymous instance with a type, as shown in fig. 4 is represented by (_ :? rdf : type Concept i ).

Let's considering example introduced in section I and the corresponding ontology is presented in fig. 1 to illustrate the of universe of discourse usage. When Chair_1 indicates that sitting_on is an active relationships, the ontology illustrated in fig. 1 becomes the ontology illustrated in fig. 2. Note that an anonymous instance appeared and is linked to the instance Chair_1 using the sitting_on relationship. From this step, the anonymous instance has no type.

The universe of discourse corresponding to the expression of this knowledge is limited to a single RDF triplet and its negation that expresses "no Human is sitting on Chair_1":

Ω 1 = {
< _ :?, sitting_on, Chair_1rdf : typeChair >, < _ :?, sitting_on, Chair_1rdf : type Chair >} (1) However, we consider that any instance of ontology is candidate to be the real instance corresponding to the anonymous instance. In our example, Human_1, Human_2, Object_1 and Chair_1 are candidates. As the sitting_on relation is not reflexive, Chair_1 is dropped from the candidates list. The Fig. 3 illustrates this.

The corresponding universe of discourse is:

Ω 2 = {
To illustrate another example of universe of discourse usage, we consider the same previous example except that the type of the anonymous instance is known: anonymous instance is a Human type. In this case, any instance of ontology restricted to the anonymous instance type is candidate to be the real instance corresponding to the anonymous instance. In this example, Human_1 and Human_2 are candidate.

C. The uncertainty model

We propose to use mass to express uncertainty on a knowledge. The wanted result is a set of candidate entities associated with mass for each of them. Each candidate entity can be the real instance corresponding to the anonymous instance. This link is associated with a given mass. In order to model uncertainty into ontology representation, we propose to associate a mass to the owl : sameAs relationship. According to the sources of information, the masses are distributed or rearranged on the related universe of discourse. . Known relationships on the anonymous instance act as constraints that must be respected when the alignment is performed. In our example, the ABox expresses the fact that the actor sitting on Chair_1 is a Human and Object_1 is not of Human type so it is no more a possible candidate for the identification of the anonymous instance.

Fig. 5. The owl : sameAs relationships related to a same anonymous node hold the mass that will represent the uncertainty on the indetification of the human that is sitting on the chair.

D. Managing uncertainty on a factual ontology

Modeling uncertainty using mass does not include semantic. So, in the first subsection below, we present a probabilistic approach to manipulate the masses; this is traditional approach Then, in the second subsection, we present a possibilistic approach. To facilitate understanding for the reader, we use the same previous example illustrated in fig 5. From this examples, the following events can be deduced:

• E1 = "a Human is sitting on Chair_1" ; • E2 = "Human_1 is sitting on Chair_1" ; • E3 = "Human_2 is sitting on Chair_1". At this step of the paper, we already presented the most inovative part, i.e. the knowledge split and the inclusion of the corresponding universes of discourse on the factual ontology. The refinement of the uncertainty management is now quite simple and is detailed below.

1) Probabilistic representation of uncertainty: In the probabilistic modeling, the mass represents the probability of an event. If all possible exclusive events are represented, the sum of each event probability is equal to 1. In the other case, this sum is less than 1.

According that P (E) is the probability that an event E occurred, and with respect of the RDF format, P (E1), P (E2) and P (E3) are defined as follow. P (E1) = P (< _ :?, sittingOn, Chair_1 >) P (E2) = P (< Human_1, sittingOn, Chair_1 >) P (E3) = P (< Human_2, sittingOn, Chair_1 >)

According that P (E1) = 1 -P (E1) is the probability of the event "no Human is sitting on Chair_1".

The most important point to presented here concerns the decomposition of the event "Human_X is sitting on Chair_1" into the 2 events "a human is sitting on Chair_1" and "the human sitting on Chair_1 is Human_X".

Without external knowledge about Human_1 or Human_2, for example the probability of their presence in the room, the probabilities on the owl : sameAs relationships P (< Human_2, owl : sameAs, _ :? >) and P (< Human_2, owl : sameAs, _ :? >) are equals.

In the general case, the ignorance is modeled by a uniform ditribution of the unitary probability over the universe of discourse.

2) Possibilistic representation of uncertainty: In the possibility modeling, masses are representing the possibility that an event occurred. The supremium of the possibilities of exclusive events is equal to 1. Note that a 0 value for possibility means the impossibility for an event to occur and is already modeled in our approach by no presence of the corresponding relationship. So, for any event E modeled in the factual ontology, the possibility Π(E) respects 0 < Π(E) ≤ 1.

According that Π(E) is the possibility that an event E occurred, and with respect of the RDF format, Π(E1), Π(E2) and Π(E3) are defined from the same example than before. = min(Π(< _ :?, sittingOn, Chair_1 >), Π(< Human_2, owl : sameAs, _ :? >)) [START_REF] Staab | Handbook on Ontologies[END_REF] Without other knowledge about Human_1 or Human_2 i.e. possibility of their presence in the room, the possibility on owl : sameAs relationship Π(< Human_2, owl : sameAs, _ :? >) and Π(< Human_2, owl : sameAs, _ : ? >) are equal to 1.

V. CONCLUSION

The main contribution of this paper is a new proposal to model the universe of discourse on factual ontologies in order to express the uncertainty coming from the measurements performed by sensible objects. The introduction of an anonymous instance to help modeling the universes of discourse is an elegant solution that preserves the original ontology and prevents from the augmentation of its complexity. The proposed approach can be refined to use the usual uncertainty management approaches. The refinements with the probability theory and with the possibility theory are presented in the paper. The application of the evidence theory will be presented later.

The proposed approach helps to model complex interactions between smart connected objects and gives a good support for the identification of situations involving many actors. The significant level of uncertainty produced by the sensible furniture units is now managed in our ontology based model. Indeed the fusion process is able to propagate the given uncertainty in order to identify the situations. Actually the proposed ontological model lays the groundwork of an ambient intelligence that allows us to supervise the evolution of daily routines in order to prevent the consequences of underlying medical or health issues.

Future studies will consist in further expressing various uncertainties, especially when information is provided from several instanced entities, creating a multi-source context. These sources of information can deliver various uncertainty information, that could be contradictory. Evidence theory is an appropriate candidate to manage this kind of information in order to facilitate decision process.

Fig. 1 .

 1 Fig. 1. This figure presents a concept ontology (named Tbox), a factual ontology (named Abox or instance ontology) and the rdf : type relation that links them. In this factual ontology, we can recognize a situation where a human is sitting on a chair.

<

  _ :?, owl : sameAs, Human_1 rdf : type Human >, < _ :?, owl : sameAs, Human_2 rdf : type Human >, < _ :?, owl : sameAs, Object_1 rdf : type Object >}

3 )Fig. 2 .

 32 Fig. 2. An anonymous instance is used in the factual ontology to represent the fact that something is sitting on the chair chair_1.

Fig. 3 .

 3 Fig. 3. We propose to use the sameAs relationship to define the set of all instances that possibly be the anonymous instance. This example modelizes the fact that either Human_1, Human_2 or Object_1 are possibly sitting on Chair_1

  Fig 5 illustrates this example.

Fig. 4

 4 Fig.4. Known relationships on the anonymous instance act as constraints that must be respected when the alignment is performed. In our example, the ABox expresses the fact that the actor sitting on Chair_1 is a Human and Object_1 is not of Human type so it is no more a possible candidate for the identification of the anonymous instance.

  Whe have: P (E2) = P (< Human_1, sittingOn, Chair_1 >) = P (< _ :?, sittingOn, Chair_1 > ∩ < Human_1, owl : sameAs, _ :? >) = P (< _ :?, sittingOn, Chair_1 >) .P (< Human_1, owl : sameAs, _ :? >) (4) and P (E3) = P (< Human_2, sittingOn, Chair_1 >) = P (< _ :?, sittingOn, Chair_1 > ∩ < Human_2, owl : sameAs, _ :? >) = P (< _ :?, sittingOn, Chair_1 >) .P (< Human_2, owl : sameAs, _ :? >)

  Π(E1) = Π(< _ :?, sittingOn, Chair_1 >) Π(E2) = Π(< Human_1, sittingOn, Chair_1 >) Π(E3) = Π(< Human_2, sittingOn, Chair_1 >) Whe have: Π(E2) = Π(< Human_1, sittingOn, Chair_1 >) = Π(< _ :?, sittingOn, Chair_1 > ∩ < Human_1, owl : sameAs, _ :? >) = min(Π(< _ :?, sittingOn, Chair_1 >), Π(< Human_1, owl : sameAs, _ :? >)) ) = P (< Human_2, sittingOn, Chair_1 >) = Π(< _ :?, sittingOn, Chair_1 > ∩ < Human_2, owl : sameAs, _ :? >)

ACKNOWLEDGMENT

This research is supported by the Universite Savoie Mont Blanc, the Association des Pays de Savoie and Miliboo Corporation.


D Tromeur-Dervout

LID-DRIVEN CAVITY BENCHMARK WITH PETSc

Keywords: High performance computing, Navier-Stokes, lid-driven cavity problem, Navier-Stokes system

published or not. The documents may come

Introduction

This paper is devoted to the use of the Portable, Extensible Toolkit for Scientific Computation (PETSc) for solving the unsteady three-dimensional Navier-Stokes equations written in a conservative form of the velocity-vorticity formulation applied to the liddriven cavity of spanwise aspect ratio 3:1 at a Reynolds number Re = 3200 on uniform non-staggered grids covering all the span.

From the numerical viewpoint, the three-dimensional flows in a cavity serve as ideal prototype non-linear problems for testing numerical codes. Geometry simplicity and weIl defined flow structures make these flows very attractive as test cases for new numerical techniques and also provide benchmark solutions to evaluate differencing schemes and problem formulation. We use this test case to study the PETSc features in the teaching of high performance computing to the fifth year engineer students in applied mathematics at Polytech Lyon engineering school (graduate students that already having 240 European Credit Transfer and Accumulation System (ECTS) and an M1 level in the bachelor's master's doctorate system ). These classes follow the teaching of the Message Passing Interface library, which facilitates the understanding of the PETSc data distribution and communication through MPI communicators.

We show how different time discretizations of this Navier-Stokes system of equations obtained with writing implicit in the time formulation more and more terms can illustrate on this problem the use of linear and nonlinear solvers features of PETSc software. Further investigations on the French national resources permitted some comparisons (up to 2048 cores) in terms of parallelism speed-up and efficiency of the time discretization strategies .We also investigated the impact of the time discretizing strategies on the flow behavior, showing that having some term explicit in time on the boundary conditions leads to a delay in time in the flow behavior. From our knowledge there no simulation of velocity-vorticity 3D formulation of Navier-Stokes with fully implicit time discretization in the literature. Some works can be found for the fully time implicit 3D lid-driven cubic cavity in primitive variables with Reynolds 1000 [START_REF] Elman | A parallel block multi-level preconditioner for the 3D incompressible Navier-Stokes equations[END_REF] or also [START_REF] Loghin | Schur complement preconditioners for the Navier-Stokes equations[END_REF]. Comparison between the flow behaviors in a cubical cavity and in a cavity of spanwise aspect ratio 3:1 is over the scope of this paper. We are mainly focus to see how the time discretisation strategies (semi implicit or fully implicit) with the same set order for the time and the space discretizations impact the flow.

The plan of this paper is as follows: section 2 describes the 3D Navier-Stokes equations governing the flow written in vorticity-velocity formulation, and the special care needed to define the boundary conditions on the vorticity. Section 3 focuses on the time discretization and its programming counterpart in the PETSc coding framework. We show three time discretizations strategies with taking implicitly more and more terms in the system of equations. Section 4 gives the results in term of parallelism of the three discretizations in time strategies while section 5 exhibits the flow behavior with respect of the time strategies. Section 6 concludes this paper.

1 Lid-driven cavity test case and space discretization

Governing equations

The velocity-vorticity formulation of the unsteady three-dimensional NavierStokes equations is mathematically equivalent to the primitive variables velocity-pressure formulation as demonstrated in [3] [4]. It is often chosen, as opposed to the primitive variable velocity-pressure Vp formulation, as transport is the critical physical phenomenon of unsteady viscous flows, it leads to a more natural decoupling of the governing equation by separating the spin dynamic of a fluid particle (represented by the vorticity transport equation) from its translation kinematics (represented by the elliptic velocity problem) [START_REF] Gatski | Review of incompressible fluid-flow computations using the Vorticiy Velocity formulation[END_REF] [6] and is completely independent of the pressure. Notably, as pointed by [START_REF] Huang | Velocity-Vorticity Simulation of Unsteady 3-D Viscous Flow within a Driven Cavity[END_REF], the vorticity transport equation is quasi-linear in vorticity and independent of pressure, whereas the velocity transport equation (momentum equation) is nonlinear in velocity and coupled to the pressure.

The lid-driven cavity is a classical test case in fluid mechanics. Several papers investigated the flow behavior with respect to the Reynolds numbers [START_REF] Bhaumik | A new velocity-vorticity formulation for direct numerical simulation of 3d transitional and turbulent flows[END_REF][START_REF] Guj | A vorticity-velocity method for the numerical solution of 3d incompressible flows[END_REF]. The flow in a lid-driven cavity of spanwise aspect ratio 3:1 at a Reynolds number Re = 3200 is a challenging problem as fully transient solutions are expected to show up. The difficulties for meaningful calculations come from both space and temporal discretizations which have to be sufficiently accurate to resolve detailed structures like Taylor-Görtlerlike vortices and the appropriate time development. We expect to exhibit in the plan parallel to the flow a primary eddy at the center and three secondary eddies that rotate in the opposite way than the primary eddy. In the plan orthogonal to the flow some Taylor-Görtler (TG) vortices in the flow appear after some time. The number of these TG vortices seems to be 7 at time T = 100 and 9 at time T = 200 [START_REF] Tromeur-Dervout | Parallelization via Domain Decomposition Techniques of Multigrid and ADI Solvers for Navier-Stokes Equations[END_REF]. This flow structures are depicted in the figure 1.

Figure 1: Lid-driven cavity with aspect ratio 3:1

The dimensionless unsteady incompressible Navier-Stokes equations in velocityvorticity ( Vω) conservative form are formulated as folIows, neglecting body forces:

∂ ω ∂t -∇ × ( V × ω) = 1 Re ∆ ω + B.C. + I.C (1)
∆ V = -∇ × ω + B.C. + I.C. (2)
The Reynolds number is defined as Re = U ∞ L ν where ν is the kinematic viscosity, L = 1 is the characteristic length of the cavity and U ∞ = 1 is the velocity modulus of the lid driven. We have a transport equation with a time derivative for the vorticity ω de f = ∇ × V while an elliptic equation gives the velocity. These equations have the advantages to avoid the difficult computing of the pressure that guarantees the incompressibility in primitive variables and also to compute directly the vorticity for rotating flows. The drawback is the definition of the vorticity boundary condition that is not naturally defined and must be deduced from the vorticity definition. The elliptic equation must be accurately solved as it guarantees the incompressibility. The other drawback is that we have to solve six scalar equations: three for the velocity and three for the vorticity although some mathematical manipulations can reduced the number of equation for the velocity.

Space discretization

We discretize the computational domain using the same regular step size in the three dimensions. Then we approximate the differential operators using second order finite differences. This leads to the standard discretization stencil with 7 points for the Laplacian and we use centered second order discretization for the first order partial derivatives of the convection term (for the forthcoming time discretization strategies 1 and 3 and first order decentered finite differences discretization for time discretization strategy 2). This space discretization gives regular data dependencies that links each points to its six neighbor points if they exist (see figure 2. Our approach has the advantage to avoid the corner singularities in the computation but misses the effect on the flow of this corner singularities [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF]. This effect should be taken in account in further developments. Here, we are mainly focused on the effect of the time discretization on the parallelism efficiency of its PETSc implementation and on the flow behavior.

Boundary conditions

For t ≥ 0 we impose no-slip boundary condition for the velocity:

V = (V x = 0,V y = 0,V z = 0) for (x = ± 3 2 , y = ± 1 2 , z = -1 2
), V = (0, 1, 0) for (z = + 1 2 ). For the vorticity boundary condition we use the vorticity definition, ∇ × V associated to the values of the velocity derivatives at the walll. For example, at the wall z 0 = ± 1 2 , we have :

V x = V y = V z = 0, and ∂V x ∂x = ∂V x ∂y = ∂V y ∂x = ∂V y ∂y = ∂V z ∂x = ∂V z ∂y = 0 (3)
Using the definition of the vorticity we obtain:

z 0 = ± 1 2 : ω x (x, y, z 0 ) = - ∂V y ∂z ; ω y (x, y, z 0 ) = ∂V x ∂z ; ω z (x, y, z 0 ) = 0
The same considerations lead to the other vorticity boundary conditions:

x 0 = ± 3 2 : ω x (x 0 , y, z) = 0; ω y (x 0 , y, z) = - ∂V z ∂x ; ω z (x 0 , y, z) = ∂V y ∂x y 0 = ± 1 2 : ω x (x, y 0 , z) = ∂V z ∂y ; ω y (x, y 0 , z) = 0; ω z (x, y 0 , z) = - ∂V x ∂y
In order to have a second order approximation for the vorticity boundary conditions, we use a limited development of the velocity components at the neighborhood of the wall of the cavity. Considering, the Poisson equations for the velocity on the wall z = -1 2 , we have :

∂ 2 V x ∂z 2 = ∂ω y ∂z , ∂ 2 V y ∂z 2 = - ∂ω x ∂z
Using a Taylor series of V x and V y at the viscinity of the wall z 0 = -1 2 , we have:

V x (x, y, - 1 2 + ∆z) -V x (x, y, - 1 2 ) = ∆z ∂V x ∂z (x, y, - 1 2 ) + ∆z 2 2 ∂ 2 V x ∂z 2 (x, y, - 1
2 ) + O(∆z 3 ) V y (x, y, - 1 2 + ∆z) -V y (x, y, - 1 2 ) = ∆z ∂V y ∂z (x, y, - 1 2 ) + ∆z 2 2 ∂ 2 V 2 y ∂z 2 (x, y, - 1
2 ) + O(∆z 3 )
With replacing the partial derivatives of V x and V y with respect of those of ω x and ω y :

V x (x, y, - 1 2 + ∆z) -V x (x, y, - 1 2 ) = ∆zω y (x, y, - 1 2 ) + ∆z 2 2 ∂ω y ∂z (x, y, - 1 2 ) + O(∆z 3 ) V y (x, y, - 1 2 + ∆z) -V y (x, y, - 1
2 ) = -∆zω x (x, y, - 1 2 ) - ∆z 2 2 ∂ω x ∂z (x, y, - 1
2 ) + O(∆z 3 )
With discretizing ∂ω y ∂z and ∂ω x ∂z at the first order, we deduce the boundary condi-

tions at z = -1 2 : ω x (x, y, - 1
2 ) + ω x (x, y, - 1 2 + ∆z) = - 2 ∆z (V y (x, y, - 1 2 + ∆z) -V y (x, y, - 1 2 )
)

ω y (x, y, - 1 2
) + ω y (x, y, -

1 2 + ∆z) = 2 ∆z (V x (x, y, - 1 2 + ∆z) -V x (x, y, - 1 2 )
)

ω z (x, y, - 1 2 ) = 0
We see that the boundary condition on the vorticity connects the point on the boundary and the neighbor point within the computational domain. This is an important consideration as all the flow dynamics come from the lid-driven. The other point to highlight is the introduction of a certain delay in time between the vorticity and the velocity if these boundary conditions are not implicitly treated (considering the vorticity and the velocity at the same time step).

Time discretization strategies and their coding in PETSc

This lid-driven cavity problem has been implemented in a research parallel code, requiring three years man of development, and using ADI for time marching on vorticity (second order in time) and multigrid accelerated by schur dual domain decomposition with generalized conjugate residual Krylov method for the velocity [START_REF] Tromeur-Dervout | Résolution des Equations de Navier-Stokes en Frmulation Vitesse-Tourbillon sur Systèmes Multiprocesseurs à Mémoire Distribuée[END_REF]. This code implementation was efficient numerically and in elapsed time in the benchmark of codes dedicated to solve this test case in the workshop [START_REF] Tromeur-Dervout | Parallelization via Domain Decomposition Techniques of Multigrid and ADI Solvers for Navier-Stokes Equations[END_REF].

A direct implementation of this code in the PETSc framework would not be possible without devoting a lot of effort as there is no ADI solvers implemented in PETSc. The two difficulties would have been to define Krylov solvers associated to the tridiagonal operator in each direction of space and to split the vector of unknowns defiined on the 3D distributed mesh in a set of vectors of unknowns associated to a 2D distributed mesh for each space direction. This can be an opportunity of development of PETSc solvers.

Instead, we investigate different features of the PETSc [START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF] parallel software to implement three time discretization strategies for (1)-( 2) that take more and more terms implicitly in the equations and boundary conditions. This benchmark constitutes the core content to teach PETSc software to engineer students in applied mathematics in few hours. We must notice that the computational domain is simple with structured data. If we need to deal with more complex geometries and still use this structured data we should use some domain decomposition strategies to split the computational domain in block structured data and use some transformation mapping of the block structured meshes to fit the complex geometry. The other solution should be to work with unstructured data, and use some graph partitioning such as par-metis available in PETSc to distribute the vector components between processors.

As the problem geometry is a cavity and the space discretization is with second order finite differences then the velocity and vorticity fields can be represented by structured data like the PETSc distributed array (DA). This is done by the DMDACreate3D function of Table 1. The User can let PETSc choose how to distribute the data between the MPI processes. He has to give the dimensions in each direction, the type and the length l of the computing stencils (stencil box for data dependencies in (i±l, j ±l, k ±l) or Stencil star for data dependencies in (i ± l, j, k), (i, j ± l, k), (i, j, k ± l)) ) and the degree of freedom per point DOF. This DA defines the data dependencies and consequently the maximum number of coefficient entries per row in the discretization matrices.

Table 1: PETSc coding to distribute the fields between processors

Strategy 1 : time explicit convection term

The first strategy solves (1) with taking the convection term explicitly and then solves (2) with taking the curl of the computed vorticity as right hand side:

(I - ∆t Re ∆) ω n+1 = ω n + ∆t ∇ × ( V n × ω n ) (4)
∆ V n+1 = -∇ × ω n+1 , after solving (4) (5)
We use a distributed array with 3 DOF and we define two Krylov Space Projection Table 2: PETSc coding of the KSP set up for velocity solvers KSP one for the velocity (KSPV) and one for the vorticity (KSPW) (see Table 2). The KSP and its preconditioner PC have to be defined with the KSPSetType and PCSetType functions. The KSP is associated to the DM through the KSPSetDM function. This defines the potential matrix structure of the operator associated to the KSP. The user provides two functions, one for building the matrix (here MatrixVelocity) and for building the right hand side (here RHSVelocity) that are associated to the KSP through the KSPSetComputeOperators and KSPSetComputeRHS functions respectively (see Table 2).

Table 3: Part of the PETSc coding for the RHS vorticity function

The main programming effort is to define the functions giving the matrix and the RHS (see Table 3). The DM is retrieved from the KSP through KSPGetDM, then DMDAGetLocalInfo gets the sizes (xm,ym,zm) of the subdomain managed by the process and the size (mx,my,mz) of the global computational domain aand also the starting point (xs,ys,zs) at the left down corner of the subdomain. The data are extracted from the vector by DMDAVecGetArrayDOF. Extracting the data field with its ghost points in the subdomain needs to create a local vector using DMGetLocalVector, then to extract a localArray with again DMDAVecGetArrayDOF and finally to update the ghost point with neighbor subdomains with DMGlobalToLocalBegin and DMGlobalToLo-calEnd. The data of the RHS are then computed using the 4-index array

[k][j][i][dof]
The data associated to the RHS have to be restore in the vector using DMDAVecRe-storeArrayDOF and must be assembly with VecAssemblyBegin and VecAssemblyEnd. Local vectors have to been restored before leaving with DMRestoreLocalVector.

Then for strategy 1 one time step consists to solve the linear system with KSPSolve and to extract the solution from KSP with KSPGetSolution. We tested several KSP Table 4: PETSc coding of the time loop for strategy 1 solvers associated to different preconditioners. It appears that the best pair KSP/PC for the elliptic equation is the GMRES [START_REF] Saad | GMRES -A Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF][START_REF] Kuznetsov | Numerical methods in subspaces[END_REF] associated to the HYPRE Algebraic multigrid [START_REF] Falgout | hypre: A library of high performance preconditioners[END_REF] preconditioner that converges in 5 (respectively 3) iterations per time step for the velocity (respectively the vorticity) for the 192 × 64 × 64 mesh.

As both operators for the vorticity solution and the velocity solution do not depend on the time iteration, we also studied the choice of KSP composed by KSPPRE-ONLY where only the PC is applied with the choice of Gauss factorizaton with the SuperLU dist package [START_REF] Li | SuperLU DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems[END_REF]. This is done simply by passing the "-ksp type preonly" option in the execution order.

Strategy 2 : semi-implicit convection term

Strategy 2 solves (1)-( 2) with taking the convection term semi-implicitly:

(I -∆t Re ∆). + ∆t∇ × ( V n × .) 0 ∇ × . ∆. ω n+1 V n+1 = ω n 0 ( 6
)
Let us notice that for strategy 2, we still use the boundary condition for the vorticity using the V n as we split the solution of (6) in two stages, the solution of the velocity and then the solution for the vorticity.

The PETSc implementation of strategy 2 involves two distributed array (DA) with three degree of freedom (DOF) with two Krylov solvers (KSP) with updating a part of the matrix vorticity at each time step. The vorticity boundary condition would have been taken implicitly if we used a DA with six DOF and on KSP taking in account the six equations. One time step for strategy 2 differs from those of strategy 1 with calling Table 5: PETSc coding of the time loop for strategy 2 KSPSetComputeOperators function to update the matrix vorticity before KSPSolve associated to the vorticity equation (see Table 5).

Again the best pair KSP/PC for strategy 2 is the GMRES/HYPRE pair , that converges in 5 (respectively 3) iterates per time step for the velocity (respectively the vorticity) for the 192 × 64 × 64 mesh and in 6 (respectively 4) iterates per time step for the velocity (respectively the vorticity) for the 384 × 128 × 128 .

Strategy 3 : totaly implicit formulation

The third strategy solves (1)-(2) totally implicitly:

(I -∆t Re ∆)( ω n+1 -∆t ∇ × ( V n+1 × ω n+1 ) -ω n ∆ V n+1 + ∇ × ω n+1 = 0 (7)
In this strategy the boundary condition for the vorticity is totally implicit as well as the convection term as the six equations are solved together.

Table 6: PETSc coding of the SNES setup for strategy 3 and the time loop associated For strategy 3 the DM has six DOF and we create a SNES (system of Nonlinear Equations Solver) object. This DM is associated to the SNES through the command SNESSetDM. Then the user has to provide the function to minimize, implementing the six discretized equations and boundary conditions with SNESSetFunction. As the SNES may needs the Jacobian of the function to minimize we can provide the jacobian to the SNES with the SNESSetJacobian function or we can let the Jacobian be evaluated with differentiating with providing a Index Set coloring associated to the data dependencies in the function to minimize in order to optimize the number of call to the function. This is done by the MatFDColoringCreate creating a MatFDColoring object that has been associated to the SNES through the SNESSetJacobian command.We also can change the KSP from the SNES through SNESGetKSP. One time step iteration consists in applying SNESSolve and then get de solution with SNESGetSolution.

We also tested several SNES solver the two candidates have been newtonls (for Newton based nonlinear solver that uses a line search [START_REF] Deuflhard | Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms[END_REF]), and ngmres (for Nonlinear Generalized Minimum Residual method [START_REF] Oosterlee | Krylov subspace acceleration of nonlinear multigrid with application to recirculating flows[END_REF] [21] that combines m previous solutions into a minimum-residual solution by solving a small linearized optimization problem at each iteration, very similar to the Anderson mixing [START_REF] Walker | Anderson acceleration for fixed-point iterations[END_REF]). While newtonls was performing better to reach a relative tolerance of 1e-6 ,ngmres with up to 300 previous solutions was prefered to reach a relative tolerance of 1e-8 as we will discuss in the flow behavior result. Figure 3 gives the comparison of the number of ngmres iterations with respect to the simulaton Time (s) for the 192 × 64 × 64 mesh and for two time steps dt = 0.005 and dt = 0.02. We see that for the two simulations the number of ngmres oscillates with quite the same behavior until T = 40. The mean of the number of iteration is 245.19 (respectively 259.17) for dt = 0.005 (respectively dt = 0.02). We see that the computational effort for ngmres is slightly greater with the time step dt = 0.02 than for dt = 0.005 but not in the ratio of 4 that the no limitation on the CFL condition allows.

3 Parallelism results

Computational resources and time measurement

Table 7 presents the characteristics of computational resources used to obtain the parallelism results. We dispose of a local computational resource that allows us to define the methodology and to obtain the flow behavior and some results on parallel efficiency in a shared memory environment. The national allocated computing resource occigen is finite 50000 hours and is used to obtain parallel speed-up in an dedicated environment with up to date software framework by an expert dedicated team. The results on occigen will use as much as possible the load balancing between memory and CPU with using, when is possible, not all CPUs of each node. PETSc provides a smart mechanism to monitor the performance with a command line option -log view that provides counters on events defined by the user. The user has to declare PetscLogStage and use the PetscLogStagePush(stage) and PetscLogStage-Pop(); before and after the portion of code that must be monitored as illustrated in Table 8. Then when the code terminates, the log of selected events with their time measurement are displayed as shown in table 9.

Table 9: PETSc log of the code for 10+1(first) time steps of strategy 1 on glenan for the mesh 192 × 64 × 64

We made the observation that the first solve always takes more elapse time, this why we separated the first solve from other iterations in the time measurement. Figure 4 gives the speed-up for strategy 1 for two meshes 192 × 64 × 64 and 384 × 128 × 128 on glenan and occigen with the GMRES+Hypre solver. We see that the speed-up on glenan is limited to 8 for 64 threads while the KSPW scales better on occigen with a speed-up of 512 for 1024 processes for the 192 × 64 × 64 mesh and 741,86 for 2048 processes for the 384 × 128 × 128 mesh. We also note that the KSPW scales better than the KSPV on occigen due to the nature of the operator nearby the identity compared to the Laplacian operator, that should give less complex algebraic multigrid projection and restriction operators and so more efficient in term of parallelism preconditioner. Nevertheless the speed-up on glenan for the KSPW is strongly degraded on the 384 × 64 × 64 with nearby 4 for 32 threads. This degradation is less sensible on occigen. Figure 5 gives the speed-up for strategy 2 for two meshes 192 × 64 × 64 and 384 × 128 × 128 on glenan and on occigen with the GMRES+Hypre solver. Again the speedup is better for the KSPW than for KSPV, even if the stagnation of the speed-up comes more early with the introduction of convective term in the vorticity operator that leads to an operator with the diagonal dominance reinforced . A speed-up of 426.67 for 1024 processes is reach on the 192 × 64 × 64 on occigen for KSPW. This speed-up stagnates on 2048 processes. The speed-up is again around 786.36 for 2048 processes on occigen for the 384 × 128 × 128 mesh.

Speed-up for strategy 1

Speed-up for strategy 2

Speed-up for strategy 3

Figure 6 gives the speed-up for the strategy 3 on occigen for a mesh of 192 × 64 × 64 with Nonlinear GMRES with up to 300 previous solutions. A speed-up of is obtained for 1024 CPUs. We see that the SNES's initialization (Snes init), that corresponds to the Figure 6: Speed-up for strategy 3 with ngmres and 192 × 64 × 64 buliding of the data structures for the snes, scales until 64 and the scaling performances are degraded after. As this initialization time is small compared to the time associated to the number of time iterates, this result is not significant. We still spearate the first solve to the other iterations as it takes more time. The SNES solving scales well until 64 processes on occigen, and a speed-up of 298 ( 357 respectively) is obtained on 1024 processes (respectively 2048 processes).

Elapsed time comparison for the three strategies

Table 10 gives the elapse time per iteration for all the strategies for the 192 × 64 × 64 mesh with dt = 0.005 on occigen computer.

• The strategy 1 with superLU didn't run onto occigen due to lack of memory with 110Gb. Strategies 1 and 2 take quite the same time on one processor while strategy 3 took 74.76 times more time than strategy 2 on one process. This ratio drops to 20 on 256 processes. Strategy 3 still reduces the elapsed time after 256 processes.

• Strategy 2 is quite always better than strategy 1 even if the matrix for the vorticity change. This is quite surprising. This may due to the decentered discretization of the convection term for the vorticity equation that leads to more diagonal dominant matrix than for the Laplacian and consequently has a better Krylov convergence. The velocity equations are solved identically in both strategies.

• For 256 processors the strategy 3 becomes reachable as there is no CFL condition on the time step. To be complete, we must notice that the elapsed time per point and time step is around 1.60e -6 to 3e -7 for strategy 2 that is not impressive compared to the elapsed time of the optimized research code with ADI and multigrid and domain decomposition done on the CRAY YMP 25 years ago with 2.510 -5 s per point and time step for a 81 × 41 × 41 mesh [START_REF] Tromeur-Dervout | Parallelization via Domain Decomposition Techniques of Multigrid and ADI Solvers for Navier-Stokes Equations[END_REF].

Nevertheless the programming effort was low (once understanding the PETSc philosophy of implementation); we also let PETSc deciding of the data distribution and the equations have been solved with a better accuracy at each time step.

Flow behavior

For the flow behavior comparison between the three strategies, we restrict ourselves to compare the flow behavior in the xz plane at y = 1/6 where the Taylor-Gortler vortices should appear. We extract the distributed component ω y from the distributed mesh through MPI implementation although we could use some vecGather function from the PETSc library. The figures have been generated with matlab contourf function for isovalues equal to [-10, -5, -3, -1, -0.5, 0, 0.5, 1, 3, 5, 10] 4.1 Underresolved problem for strategy 3 Figure 7 exhibits a problem that occurred for the strategy 3 when the solver a newtonls with a rtol = 1e -6 was used. It presents the flow behavior at T = 38 that keeps its symmetry while it loses its symmetry for T = 60. This is due to some underresolved solution at some time after T = 38. We advocate that as the equations are symmetric as the boundary conditions also, we should keep the symmetry of the flow. The problem was solved by changing the solver with the ngmres with up to 300 previous solutions that allows to solve the problem with a rtol = 1e -8.

Figure 7 gives the flow behavior for T = 38 and T = 60 for strategy 3 with newtonls solver with a tolerance of 1e -6. We see that the flow loses its symmetry due to an underresolving during some time step. This why we easely changed the solver with ngmres with 300 fixed-point iterations in order to reach a tolerance of 1e -8 . Figures 8 ,9, and 10 compare the flow behavior for the three strategies at time T = {20.28, 45.00, 73.48} for the mesh 192 × 64 × 64 and a time step value of dt = 0.005 (which is imposed by the CFL condition of strategy 1). The solvers used are su-perlu+dist, gmres+hypre and ngmres 300 for strategies 1, 2 and 3 respectively.

Figure 8 exhibits some quite similar results for strategy 1 and strategy 2 with slightly more strong vortex in position i=45 and k=50. The difference coming from taking the vorticity term implicit in time in the convection term but the treatment of the vorticity boundary conditions are unchanged between the two strategies. For the strategy 3, it shows that the vortices attached to the boundaries are more strong, and the flow is slightly in advance compared to the two other strategies as the vortex in position i=20, k=55 is already attached to the vortex on the boundary. We conclude that taking explicit in time the convection term and with more impact taking explicit in time the right hand side of the vorticity boundary conditions introduce a delay on the flow.

The flow behavior at T=45 on Figure 9 still exhibits the quite same flow structure for the three strategy, with more strong vortices for the strategy 3. The flow keeps its symmetry. This is also the case for Figure 10 where the strategy 3 exhibits more defined vortices structures.

Conclusions

This paper shows the use of different functionalities of PETSc to solve the lid-driven cavity problem with aspect ratio 3:1 for Reynolds 3200 from a semi-implicit to a fully implicit in time formulations. We notably exhibits that writing all the terms of the vorticity boundary condition at time step t n+1 leads to stronger vortices much more than writing all the terms involved in the convective term at same time step t n+1 . This conclusion is somewhat foreseeable as all the dynamics of the flow comes from the boundaries. The PETSc library permits us to tests different solvers and to choose the solver and its parameters the best adapted for the computation. It was somewhat surprising in a first approach that taking the convective term implicit gives better results on the parallelism for the solver associated to the vorticity. This is mainly due to the HYPRE preconditioner behavior, that creates algebraically the restriction and projection operators, which performs well for the operator where we reinforce the diagonal dominance.

The fully implicit formulation has the advantage to not have CFL condition and allows us to violate at least 4 times this CFL of the strategy 1 with keeping the right flow behavior. One output of this work, is to demonstrates that the fully implicit in time strategy can be competitive (reachable) with a sufficient number of processes. Nevertheless, we must mention that we could provide to the PETSc library with the KSPRegister instruction our owner implementation of gmres where we can stored the directions of descent from different consecutive time steps in a cumulative Krylov space since the operator for strategy 1 does not change. Then we can initialize the new time step solution with projecting the new right hand side with respect to this Krylov space [START_REF] Tromeur-Dervout | Résolution des Equations de Navier-Stokes en Frmulation Vitesse-Tourbillon sur Systèmes Multiprocesseurs à Mémoire Distribuée[END_REF][START_REF] Erhel | An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems[END_REF][START_REF] Tromeur-Dervout | Choice of initial guess in iterative solution of series of systems arising in fluid flow simulations[END_REF]. We now have a non-trivial PETSc base test case that can allow us to develop non-linear acceleration methods by further entering in the PETSc coding to define our own SNES solvers, using for example singular value decomposition of previous time step solutions to accelerate the Newton [START_REF] Tromeur-Dervout | Choice of initial guess in iterative solution of series of systems arising in fluid flow simulations[END_REF]. The discretization could also be improve by several way, with at least better space discretization with compact schemes [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] for example and better discretization in time. The perspective of this work could be also to go further in controling the error on the solution with using for example BDF schemes with adaptive time step of the SUNDIALS [START_REF] Hindmarsh | Sundials: Suite of nonlinear and differential/algebraic equation solvers[END_REF] library.

Figure 2 :

 2 Figure 2: Stencil of the data dependancies.

Figure 3 :

 3 Figure 3: Comparison of the number of ngmres iterations with respect to the simulaton Time (s) for the 192 × 64 × 64 mesh and for two time steps dt = 0.005 (full line) and dt = 0.02 (dash line)

Figure 4 :

 4 Figure 4: Speed-up (strong scaling) of strategy 1 on glenan (left) and occigen (right) for domain sizes 192 × 64 × 64 (top) and 384 × 128 × 128

Figure 5 :

 5 Figure 5: Speed-up (strong scaling) of strategy 2 on glenan (left) and occigen (right) for domain sizes 192 × 64 × 64 (top) and 384 × 128 × 128

Figure 7 :

 7 Figure 7: Effect on the undersolve time step solution on the flow behavior: up the Newton LS algorithm for strategy 3 with a tol = 1e -6 at T = 38 and bottom the lost of the flow symmetry at T = 60 .

Figure 8 :

 8 Figure 8: Comparison of the flow behavior for the three strategies (strategy 3 up, strategy 2 middle, strategy 1 bottom) at time t=20.

Figure 9 :

 9 Figure 9: Comparison of the flow behavior for the three strategies (strategy 3 up, strategy 2 middle, strategy 1 bottom) at time t=45.

Figure 10 :

 10 Figure 10: Comparison of the flow behavior for the three strategies (strategy 3 up, strategy 2 middle, strategy 1 bottom) at time t=74.

4. 3 3 Figure 11

 3311 Figure11gives the flow behavior at time T = 76 for the strategy 3 with the ngmres solver for two time steps dt = 2.10 -2 and dt = 5.10 -3 . We see a good agreement between the two computations and see the advantage of strategy 3 that has no CFL condition. The limitation for the time step is only due to physics consideration, in order to catch the right dynamics and we see that we can take a 4 times greater time step for this flow behavior computation.

Figure 11 :

 11 Figure 11: Comparison for strategy 3 (with no CFL time step limitation) for two time steps dt = 2.10 -2 (top) and dt = 5.10 -3 (bottom).

Table 7 :

 7 Computing resources characteristics

	Computer	glenan	occigen
	Model	Dell PE920	bulx DLC
	Processor	Xeon	Xeon
		E7-4880v2	E5-2690v3
		4x15C 2.5Ghz	1x12C 2.6Ghz
	Memory/Node	512Gb	110Gb
	Network	shared memory Infinityband FDR
	Resource	laboratory	national (cines)
		∞ hours	50000 hours
	Use	speed-up tests	speed-up tests
		flow behavior

Table 8 :

 8 Coding of the event log for parallel performances

Table 10 :

 10 Time (s) per iteration: occigen dt=0.005 192x64x64

	P	S1 S1 LU	S2	S3	comments
	1	12.61	-12.59 941.25 too much even with no CFL
	2	9.97	32.89	9.09 479.35
	4	6.69	20.48	6.78 316.63	S1 vs S2 similar (<>tol)
	8	4.56	10.77	3.82 311.45
	16	2.91	8.50	2.52	99.18
	32	1.83	4.54	1.58	39.20
	64	0.99	4.16	0.83	17.84
	128	0.68	2.37	0.56	10.03
	256	0.40	2.42	0.34	6.04	S3 reachable with no CFL
	512	0.30	1.83	0.28	4.11
	1024	0.23	5.05	0.28	3.15
	2048	-	-	0.26	2.63

Acknowledgements

This work was granted access to the HPC resources of CINES under the allocation1 2018-AP01061040 made by GENCI (Grand Equipement National de Calcul Intensif). Author thanks also the Center for the Development of Parallel Scientific Computing (CDCSP) of University of Lyon 1 for providing local computing resources to design the methodology.