Sylvain Kubler
email: s.kubler@univ-lorraine.fr

William Derigent

Alexandre Voisin

André Thomas

Éric Rondeau

Method for embedding context-sensitive information on "communicating textiles" via Fuzzy AHP

Keywords: Product lifecycle management, Data dissemination, Fuzzy AHP, Decision theory, Intelligent product

The amount of data output in our environment is increasing each day, and the development of new technologies constantly redefines how we interact with information. In the context of PLM (Product Lifecycle Management), it is not uncommon to use intelligent products to ensure an information continuum all along the product life cycle (e.g. for traceability purposes). However, it is not that easy to identify what information should be stored on the product. To answer this question, this paper proposes a data dissemination process to select contextsensitive information from the database, that must be stored/replicated on the product (especially if the product is made of "communicating material"). Our approach uses the fuzzy AHP theory for aggregating points of view from different actors. The data dissemination process is then applied on a case study to embed context-sensitive information on a "communicating textile".

Introduction

Intelligent products

Challenges and opportunities arise with concepts such as the Internet of Things (IoT) and ubiquitous/pervasive computing [START_REF] Weiser | Hot topics-ubiquitous computing[END_REF]. Through these concepts, real-world objects are linked with the virtual world. Thus, connections are not just people to people or people to computers, but people to things and most strikingly, things to things. During its product life cycle (PLC), a product moves through numerous companies and a variety of core business sectors. Technical, semantic and organizational interoperability between these companies is not always ensured, thereby contributing to information loss. Considering the product as an information vector (to which information can be linked) should contribute to improved interoperability throughout its PLC [START_REF] López | Taxonomy, technology and applications of smart objects[END_REF]. Research has been conducted on intelligence in objects [START_REF] López | Taxonomy, technology and applications of smart objects[END_REF][START_REF] Meyer | Intelligent products: A survey[END_REF], i.e., products that carry their own information and intelligence. These are also referred to as intelligent or communicating products. However, information is often deported through a network and these solutions do not help to identify the information that must be stored on the product. Indeed, some products are only given an identifier (e.g., via an RFID tag) that provides a network pointer to a linked database and a decision-making software agent [START_REF] Meyer | Intelligent products: A survey[END_REF]. It should be noted that these products have limitations in some respects [START_REF] Mcfarlane | Intelligent Products in the Supply Chain -10 Years On[END_REF]:

-Discrete reading: it is necessary to read a specific zone of the product, 0000-0000/09/$00.00 c 2009 -IOS Press and the authors. All rights reserved -Risk of tag damage: if the unique RFID tag is damaged, all information related to the product will therefore be lost, -Problem of information transfer : when a product is transformed (e.g., cut), some resulting parts will not contain information, -Small memory capacity: RFID tags currently available are memory constrained.

After many years of considering a communicating product as a physical product associated with an informational product (realized via autoidentifying technologies such as RFID), we have proposed a new paradigm, which changes drastically the way we see the material. The concept was first proposed by Thomas [START_REF] Thomas | RFID et nouvelles technologies de communication; enjeux économiques incontournables et problèmes d'éthique[END_REF] and aims to give the ability to the material to be intrinsically and wholly communicating, as discussed in [START_REF] Kubler | Problem definition methodology for the "Communicating Material" paradigm[END_REF][START_REF] Kubler | Prototyping of a communicating textile[END_REF]. The powerful and futuristic idea behind this concept is to imagine that the material is communicating at the cellular/molecular level. Obviously, this vision is far from being possible today, especially due to the technological limitations, but some current research seems to be great on promises. For instance, in [START_REF] Kubler | Prototyping of a communicating textile[END_REF], a prototype of "communicating textile" was designed in which a huge quantity of RFID µtags is spread. Although current technologies do not allow building a "communicating material" up to the molecular level, the designing of new solutions, methods providing new services to users may be addressed. For instance, a product made of "communicating material" could have special abilities like data storage, copy/redundancy/backup information solutions.

Nonetheless, there are still some open questions. These include a specification for the information to be gathered, stored and distributed over the PLC, how information should be managed during the PLC, the linkage of new hardware and software systems with current systems [START_REF] Mcfarlane | Intelligent Products in the Supply Chain -10 Years On[END_REF]. To address some of these matters, this paper develops a data dissemination process to allocate information between databases and "communicating products". This process is described in the next section.

Data dissemination process

The data dissemination process proposed in this paper consists of three steps:

-Process step 1 implements the database system architecture of each actor/company tak-ing place on the PLC. Many studies exists that can help designers to choose the most suitable systems according to their expectations and the application constraints [START_REF] Ozsu | Principles of Distributed Database Systems[END_REF], -Process step 2 selects data appropriate to the expected situation from the database system (i.e., the data that should be stored on the product) thanks to an indicator of relevance which gives the degree to which information might be useful for the subsequent actors (e.g. manufacturers, recyclers) and activities, -Process Step 3 deals with the storage of the data selected in the previous step on the product (made of "communicating material") and with their retrieval. In our study, we propose to use the "communicating textile" developed in [START_REF] Kubler | Prototyping of a communicating textile[END_REF]. When the product is read again (after having stored data on it), data retrieved from the product can be updated in the database thanks to classical protocols (e.g. synchronous/asynchronous protocols) [START_REF] Ozsu | Principles of Distributed Database Systems[END_REF].

This paper mainly focuses on Process step 2 and 3. During its life cycle, a product moves through numerous companies to various core business sectors. The information relevance is therefore dependent upon a variety of factors at each stage of the PLC such as the user concerns, the product environment, etc. An approach relying on a Logical Data Model (LDM) is developed by Chan and Roddick [START_REF] Chan | Context-sensitive mobile database summarisation[END_REF] to select context-sensitive information by formalizing the relationship between context and data. This approach assesses the relevance of all data items, from all tables. Fig. 1(a) shows a part of a LDM, where one entity corresponds to a relational table as depicted with the entity/table MaterialDefinition in Fig. 1(b). The attributes listed in each entity correspond to the table columns and each row is referred to as a tuple. For instance, MaterialDefinition has 3 attributes and 4 tuples. One data item corresponds to one table cell (cf. Fig. 1(b)).

The Chan and Roddick's model uses the notion of priorities which are numerical values (either supplied or generated via observation and experimentation) assigned using a multicriteria evaluation. The higher the relevance value, the higher the necessity that this data item should be stored on the product. In our context, only product-related information is assessed 1 . For instance, in Fig. 1(b), only tuple 3 is assessed (see dashed background) and the relevance value of the data item located at row 3, column 1, noted T MD {3, 1}, is equal to 0.2. Finally, data items are classified in order of relevance (e.g. T MD {3, 3} is the most relevant with a value of 0.36). The approach proposed by Chan and Roddick has been applied in our previous work [START_REF] Kubler | Information dissemination process for context-aware products[END_REF] but their model does not offer possibilities to take into account several opinions and may turn out to be inappropriate to our context. Indeed, at a given moment of the PLC, several experts are concerned by selecting information that must be stored on the product, therefore generating different points of view. All points of view are legitimate and must be taken into account. Accordingly, this paper proposes an approach which combines the use of the analytic hierarchy process (AHP) with the fuzzy logic. AHP enables to handle the multicriteria decision making (MCDM) problem while fuzzy logic enables to aggregate multiple points of view.

Section 2 introduces different approaches used to handle MCDM problems and justify the use of fuzzy AHP. Section 3 details the stages that compose our fuzzy AHP method. Finally, the data dissemination process is applied on a case study in section 5 in order to embed context-sensitive information on a "communicating textile".

Fuzzy AHP/ANP and related work

There are a number of publications existing in the literature which use MCDM methods such as 1 An algorithm to identify all product-related tuples from the database is developed in [START_REF] Kubler | Information dissemination process for context-aware products[END_REF].

AHP (analytic hierarchy process), ANP (analytic network process), TOPSIS (technique for order preference by similarity to ideal situation), ELEC-TRE to solve the selection problem [START_REF] Greco | Multiple criteria decision analysis: state of the art surveys[END_REF]. There are no better or worse techniques, but some techniques better suit to particular decision problems than others do [START_REF] Zheng | Application of a trapezoidal fuzzy AHP method for work safety evaluation and early warning rating of hot and humid environments[END_REF].

AHP/ANP methods have been widely used in the literature to handle the MCDM problems [START_REF] Vaidya | Analytic hierarchy process: An overview of applications[END_REF]. These methods, originally proposed by Saaty [START_REF] Saaty | Decision making with dependence and feedback: The analytic network process[END_REF], have the advantage to organize the critical aspects of the problem in a hierarchical structure, similar to a tree, thus facilitating the decision making process. In these methods, decision makers use linguistic variables rather than expressing their judgments in the form of exact numeric values, which make them usually feel more confident and facilitates the valuation process. Moreover, the use of AHP/ANP does not involve cumbersome mathematics, thus it is easy to understand and it can effectively handle both qualitative and quantitative data [START_REF] Zheng | Application of a trapezoidal fuzzy AHP method for work safety evaluation and early warning rating of hot and humid environments[END_REF]. The ANP extends the AHP to problems with dependencies and feedback among the criteria by using a "super matrix" approach [START_REF] Saaty | Decision making with dependence and feedback: The analytic network process[END_REF]. Like AHP, pairwise comparisons in ANP are performed in a matrix framework, and a local priority vector can be determined as an estimate of the relative importance of the elements being compared 2 , but in the case of interrelated criteria, a network replaces the hierarchy. In this study, we make the assumption that no dependencies among criteria occur, and we therefore propose to implement the AHP method. Nonetheless, we highlight in conclusion that some interdependencies among criteria might be identified and, accordingly, the applica-tion of ANP should be examined in future work so as to compare the results with those of this paper.

Although AHP facilitates the decision making process, some shortcomings can be noted [START_REF] Yang | Key quality performance evaluation using fuzzy AHP[END_REF] such as (i) it does not take into account the uncertainty associated with the mapping of human judgment to a number, (ii) it does not integrate an aggregation mechanism to allow the expression of multiple points of view, (iii) the ranking of the AHP method is rather imprecise. To overcome these problems, several researchers integrate fuzzy theory with AHP [START_REF] Van Laarhoven | A fuzzy extension of Saaty's priority theory[END_REF].

Fuzzy set theory, introduced by Zadeh [START_REF] Zadeh | Fuzzy sets[END_REF], mimics human reasoning in its use of approximate information to generate decisions. It was designed to mathematically represent uncertainty and vagueness and provide formalized tools for dealing with the imprecision intrinsic to many problems. This formalization is gradual rather than crisp. In a classical set A, an element belongs entirely or not to A. In a fuzzy set Ã, an element has a degree of membership in [0, 1].

The earliest work in fuzzy-AHP appeared in [START_REF] Van Laarhoven | A fuzzy extension of Saaty's priority theory[END_REF] in which triangular fuzzy numbers were compared according to their membership functions. Later, [START_REF] Chang | Applications of the extent analysis method on fuzzy AHP[END_REF] introduced a new approach for handling fuzzy AHP, with the use of triangular fuzzy numbers for pairwise comparison scale and the use of the extent analysis method (FEAHP). Many researchers have attempted to use and to develop fuzzy AHP for selecting and solving problems. One of the significant literature surveys is presented in [START_REF] Vaidya | Analytic hierarchy process: An overview of applications[END_REF]. Compared with the other MCDM methods, the use of fuzzy AHP methodology offers a number of benefits. Firstly, the other MCDM method experiences difficulty in capturing uncertain and imprecise judgment of experts [START_REF] Das | Selection of business strategies for quality improvement using fuzzy analytical hierarchy process[END_REF]. Fuzzy AHP can overcome such inability by handling linguistic variables. Secondly, some factors are non-physical and qualitative, measuring qualitative factors by using fuzzy numbers instead of using crisp numbers helps both making decisions easier and obtaining more realistic results, therefore the fuzzy methodology is an excellent tool to handle qualitative assessments [START_REF] Zheng | Application of a trapezoidal fuzzy AHP method for work safety evaluation and early warning rating of hot and humid environments[END_REF]. As the AHP method, fuzzy AHP has been employed in numerous sectors such as political [START_REF] Kuo | A decision support system for locating convenience store through fuzzy AHP[END_REF], hospital [START_REF] Büyüközkan | Strategic analysis of healthcare service quality using fuzzy AHP methodology[END_REF] or still in the transportation sector [START_REF] Vildan Özkir | A fuzzy assessment framework to select among transportation investment projects in Turkey[END_REF].

Existing approaches for fuzzy AHP weight derivation can be classified into two categories:

approaches which derive a set of fuzzy weights from a fuzzy pairwise comparison matrix, -approaches which derive a set of crisp weights from a fuzzy pairwise comparison matrix.

Our study relies on the second category which includes the FEAHP method [START_REF] Chang | Applications of the extent analysis method on fuzzy AHP[END_REF]. Since fuzzy weights are not as easy to compute as crisp ones, the literature survey made by [START_REF] Wang | Fuzzy analytic hierarchy process: A logarithmic fuzzy preference programming methodology[END_REF] shows that the vast majority of the fuzzy AHP applications uses FEAHP for simplicity of calculation [START_REF] Vildan Özkir | A fuzzy assessment framework to select among transportation investment projects in Turkey[END_REF]. Based on this review, our attention has turn to FEAHP which is, in part, implemented and described in the next section.

Process step 2: Computation of data item relevances using fuzzy AHP

In our study, fuzzy AHP consists of 5 stages as depicted in Fig. 2: 1. breakdown the MCDM problem into a hierarchical structure: criteria/alternatives, 2. collection & aggregation of information to determine: i) fuzzy sets of alternatives with respect criteria, ii) the criteria importance, 3. creation of the fuzzy judgment matrix Ã, 4. computation of the fuzzy performance matrix H by synthesizing à with criteria importance, 5. alternative ranking: necessity of aggregating the alternative multi-criteria performance in a fuzzy vector R and then, to rank alternatives. These five stages are respectively detailed in section 3.1 to 3.5. To make the understanding easier, a scenario is considered in the rest of the paper whose parts are preceded by the symbol "➫". This scenario relies on 19 entities which come from crisp performance matrix which gives the crisp performance score of all alternatives with respect to all criteria under a degree of confidence α and a risk level β R α β (l) the final crisp score which indicates the relevance of the alternative l. The higher the relevance value, the higher the necessity that the data item l should be stored on the product the LDM of the B2MML standard 3 . Fig. 1(a) illustrates 4 of the 19 entities.

Stage 1: AHP structure

Our MCDM problem is broken down into the hierarchical structure depicted in Fig. 3. The alternatives are the data items (cf. level 3) which must be assessed and ranked in term of relevancy (cf. level 1). Four criteria are defined at level 2: Enumeration, Contextual, Data Size and Model-Based which are respectively abbreviated C e , C c , C s , C m and are detailed in stage 2 (next section).

Level 1

Data item classification (in order of relevance)

Level 2 Enumeration Contextual Data Size Model-Based Level 3 T MD{3,1} T MD{3,2} • • • data item n

Stage 2: Information collection and expert's point of view aggregation

This section details how a decision maker evaluates each alternative with respect to each criterion (C e , C c , C s , C m ) and how he evaluates the criteria importance. In our study, there are p max experts who give their evaluation in a criterion (p max may vary according to the criterion). The main notations used in this paper are synthesized in Table 1.

Enumeration -C e

Through this criterion, the user may enumerate information he judges important to store on the product. To do so, each decision maker p enumerates attributes from tables. Let t a table from the database (t ∈ T ) and v an attribute of t. If the attribute v is enumerated by an expert p, the crisp enumeration score s p (v, t) = 1, otherwise 0, as noted in equation 1. Then, all decision makers' opinion s p (v, t) are aggregated in equation 2 through a triangular fuzzy number s(v, t). Since an attribute corresponds to a column of a table, which in turn, consists of one or many data items, therefore, if a data item l belongs to v, its fuzzy score with respect to the criterion C e , noted φe (l), will therefore be equal to s(v, t).

s p (v, t) = 1 enumerated 0 not enumerated , ∀ p = {1..pmax} (1)
s(v, t)= L M U min(s p (v)) pmax p=1 s p (v) pmax max(s p (v)) (2)
➫ Two decision makers D 1 , D 2 enumerate respectively the attributes listed in Table 2. By applying equations 1 and 2, the triangular fuzzy number of each attribute of each table t is com-puted. Let MaterialDefinition be the table t and Value be the attribute v. Table 2 shows that D 1 enumerates v (i.e. s 1 (v, t) = 1) and D 2 does not enumerate v (i.e. s 2 (v, t) = 0). s(v, t) can therefore be computed as in equation 3. Let us consider the data item T MD {3, 3} (cf. Fig. 1) as the alternative l. Since T MD {3, 3} belongs to the attribute Value, its fuzzy score with respect to C e , noted φe (l), will therefore be equal to s(v, t) = [0, 0.5, 1].

s(v, t) = L M U min(0, 1) 0+1 2 max(0, 1) (3)
= 0 0.5 1 v = Value ∈ MaterialDefinition

Contextual -C c

The previous criterion allows users to specify information that must be stored on the product. However, they could omit important information. Indeed, they might not be aware of all the data needed by the downstream actors (along the PLC). As a result, a new criterion referred to as contextual -C c is integrated in order to moderate and to balance C e . First, let us note that a multitude of information systems exist over the PLC (e.g. ERP, PDM, MES) which are not concerned by the same data (i.e. the same entities from the LDM). The idea is to identify specific "entity groups" through the LDM according to, for instance, the information systems and, therefore, to evaluate their importance over the PLC. In our study, each decision maker p performs pairwise comparisons between entity groups as in equation 4, with z the number of groups defined trough the LDM. The importance of entity group i over entity group j evaluated by the decision maker p is noted s p ij . This evaluation is based on the 1 to 9-point scale from [START_REF] Saaty | The Analytic Hierarchy Process[END_REF]: {1, 3, 5, 7, 9}. s p ij =1 means that groups i and j are equal in importance and s p ij =9 means that group i is strongly favored over j. Then, all crisp scores s p ij are aggregated through a triangular fuzzy number, noted sij , in equation 5. Equation 6 is finally applied in order to acquire a unique fuzzy score (triangular fuzzy number) φc (G i ) for each entity group i. φc (G i ) actually indicates the relative importance of information involves in G i . If a data item l is contained in a table included in G i , its fuzzy score with respect to C c , noted φc (l), is therefore equal to φc (G i ).

    G 1 • • • Gz G 1 1 • • • s p 1z . . . . . . . . . . . . Gz s p z1 • • • 1     , s p ji = s p ij i = j (s p ij ) -1 i = j (4) sij = L M U min(s p ij ) pmax p=1 s p ij pmax max(s p ij ) (5) φc(Gi) = z j=1 sij z k=1 z j=1 skj (6)
➫ In our case study, the experts cluster the 19 entities in 4 groups as detailed in Table 3. The Equipment and Personal data groups report information about equipments and persons which/who are somehow related to the product (e.g. equipments used for manufacturing it). The Product and Production data groups relate respectively information about the product composition (e.g. raw material) and operations (e.g. production rule). Two decision makers D 1 , D 2 perform the pairwise comparison matrices as given in equation 7. For instance, D 1 strongly favors information from G 1 over G 2 (s 1 12 = 9), while D 2 favors slightly G 1 on G 2 (s 2 12 = 3). Therefore, the vector s12 can be determined as in equation 8. All vectors sij ∀ i, j = {1, 2, 3, 4} are synthesized in equation 9. The importance value φc (G i ) of each group i is then computed as detailed with φc (G 1 ) in equation 10 4 . Since a triangular fuzzy set φc (G i ) represents the preference of all experts about storing on the product information belonging to G i then, the higher the triangular fuzzy numbers, the higher the importance/necessity. Accordingly, it can be observed that information from the entity groups G 1 and G 3 are more important compared to information from G 2 and G 4 .

-Ã ⊕ B = [a + d, b + e, c + f ] -Ã ⊗ B = [a × d, b × e, c × f ] -Ã B = [ a f , b e , c d ]
Let us consider T MD {3, 3} as the alternative l. Since T MD {3, 3} belongs to MaterialDefinition, which is a table included in G 1 , its fuzzy score with respect to C c , noted φc (l), will therefore be equal to φc (G 1 ) = [0.29, 0.52, 0.97].

D 1 =      G 1 G 2 G 3 G 4 G 1 1 9 5 5 G 2 1 9 1 1 5 1 G 3 1 5 5 1 3 G 4 1 5 1 1 3 1      , D 2 =     1 3 5 3 1 3 1 1 3 1 3 1 5 3 1 1 5 1 3 3 5 1     (7) s12 = L M U min(9, 3) 9+3 2 max(9, 3) (8) = 3 6 9     [1, 1, 1] [3, 6, 9] [5, 5, 5] [3, 4, 5] 1 9 , 2 9 , 1 3 [1, 1, 1] 1 5 , 4 15 , 1 3 1 3 , 2 3 , 1 1 5 , 1 5 , 1 5 [3, 4, 5] [1, 1, 1] 1 5 , 8 5 , 3 1 5 , 4 15 , 1 3 [1, 2, 3] 1 3 , 8 3 , 5 [1, 1, 1]     (9) φc(G1) = [1, 1, 1] ⊕ [3, 6, 9] ⊕ . . . ⊕ [3, 4, 5] [1, 1, 1] ⊕ . . . ⊕ 1 3 , 8 3 , 5 ⊕ [1, 1, 1] (10)

Data Size -C s

This criterion favors the storage of information on the product according to the data size. Since products are often memory-constrained, the data relevance should decrease when the data size in-creases. Such a behavior can be obtained via equation 11, with b the size in bytes of a data item l and k a constant adjusted by an expert p. Fig. 4 shows two functions according to the data size b and for two different k. It can be observed that the smaller k is, the bigger the data authorized to be stored on the product is (e.g., the expert who fixes k = 1.08 does not wish to store on the product data > 60 bytes).

k -b k ∈ R + ; b ∈ N + (11)
The coefficient k must be adjusted by each expert p as in equation 12. Then, all coefficients s p k are aggregated through a triangular fuzzy number, noted sk in equation 13. The fuzzy score of an alternative l with respect to C s , noted φs (l), is therefore computed in equation 14.

s p k = k (12) sk = L M U min(s p k ) pmax p=1 s p k pmax max(s p k ) (13) φs(l)= L M U sk (U ) -d sk (M ) -d sk (L) -d (14

Model-based -C m

This criterion is based on the relationships implied through the LDM. Chan and Roddick [START_REF] Chan | Context-sensitive mobile database summarisation[END_REF] explain that the shorter the distance between tables, the higher the data correlation. In our context, the network pointer carried by the product refers to a given table of the database (i.e. from the LDM) which is called "product table". This table is at the centre of our concerns and thus, the data relevance would decrease as the modeled distance between the product table and a distant table increases. The distance corresponds to the shortest path (i.e. the number of relations which separate them). Let, for instance, MaterialLot be the product table in Fig. 1(a) and ManufacturingBill be the distant table. The distance between both tables equals 2 (relations ①-②). Once again, the expert can adjust his preference about storing distant data on the product through the same equation than previously, as detailed in equation 18 with d the distance and k a constant adjusted by the expert p.

k -d k ∈ R + ; d ∈ N + (18)
In this case, the lesser the coefficient k, the more distant information is favored. It is therefore necessary to study the entire LDM to fix k (depending on the maximal distance through the LDM). The fuzzy score of an alternative l with respect to C m , noted φm (l), is also obtained through equations 12 to 14.

➫ Two decision makers D 1 , D 2 , adjust the coefficient k. Since the model is not so large (maximum distance through the LDM equals 18), they respectively fix k at 1.08 and 2 as in equation 19. Let us consider the data item T MD {3, 3} as the alternative l. The fuzzy score of alternative l, noted φm (l) is obtained via equations 20 and 21. The power value in equation 21 is equal to -1 because the table containing T MD {3, 3} is distant of 1 relation from the product table 5 .

s 1 k = 1.08 s 2 k = 2 (19) sk = L M U min(1.08, 2) 1.08+2 2 max(1.08, 2) (20) = 1.08 1.09 2 φm(l) = L M U 2 -1 1.09 -1 1.08 -1 (21)

Criteria importance

Experts specify the criteria importance via pairwise comparisons as in the Contextual criterion. Each decision maker D p performs the pairwise comparison matrix given in equation 22. Then, all decision maker's grades are integrated in the fuzzy vector sij in equation 23. Finally, equation 24 is applied in order to acquire a unique fuzzy score/fuzzy vector φ(C x ) which represent the relative importance of criterion x. ➫ The evaluation of the importance between each criterion x is performed via pairwise comparisons but only the resulting importances (triangular fuzzy vectors) of the four criteria are given in equation 25. The results show that the criterion C e is the most important at this moment of the PLC. In other words, the experts give freedom to the user to decide what information must be stored on the product.

Stage 3: Fuzzy judgement matrix Ã

After getting all alternative fuzzy scores with respect to all criteria, the fuzzy judgment matrix à is formed as in equation 26, where φx (l) denotes the fuzzy judgement score of alternative l with respect to criterion x.

à =      Ce Cc Cs Cm A 1
φe( 1) φc( 1) φs( 1) φm( 1)

A 2 φe(2) φc(2) φs(2) φm(2) . . . . . . . . . . . . . . . An φe(n) φc(n) φs(n) φm(n)      (26)
➫ All along the scenario considered in the previous section, the fuzzy scores of T MD {3, 3} with respect to C e , C c , C s and C m have been computed, which are respectively equal to: φe(TMD{3, 3}) = [0, 0.5, 1] φc(TMD{3, 3}) = [0.29, 0.52, 0.97] φs(TMD{3, 3}) = 1.08 -4 , 1.08 -4 , 1.08 -4 φm(TMD{3, 3}) = 2 -1 , 1.09 -1 , 1.08 -1

Stage 4: Fuzzy performance matrix H

At this stage, only fuzzy scores of alternatives with respect to criteria are taken into account, without considering the relative importance of criteria. As a result, we synthesize the fuzzy matrix à with the criteria importance in a fuzzy performance matrix H. The fuzzy performance score hx (l) consists in multiplying the fuzzy set φx (l) by the criterion importance itself φ(C x ) as in equation 27. Finally, the fuzzy performance matrix H is obtained as depicted in equation 28.

hx(l) = φx(l) ⊗ φ(Cx) (27)
H =           Ce Cc Cs Cm A 1 he(1) hc(1) hs(1) hm(1) A 2 he(2) hc(2) hs(2) hm(2) . . . . . . . . . . . . . . . An he(n) hc(n) hs(n) hm(n)           (28)
➫ The equation 29 details the calculations of he (T MD {3, 3}) and hc (T MD {3, 3}).

Stage 5: Defuzzification/Alternative ranking

This stage aims at ranking alternatives according to their fuzzy sets. To do so, defuzzification is first applied to transform each fuzzy set into a crisp value. Then, alternatives are ranked by considering the crisp values based on the TOPSIS method [START_REF] Hwang | Multiple attribute decision making: methods and applications: a state-of-the-art survey[END_REF]. Defuzzification and TOPSIS methods are respectively presented in section 3.5.1 and 3.5.2.

Defuzzification step

Defuzzification is executed (i) by using the αcut method on the fuzzy performance matrix H and (ii) by taking into account the risk index β to compute the crisp performance matrix H α β . i. The α-cut level of an alternative l, noted h α x (l), is computed via equations 14 to 16. h α xa (l) and h α xb (l) denotes respectively the left and right point of the level of the triangle as shown in Fig. 5 (0 ≤ α ≤ 1). The α-cut value reflects the degree of confidence of the decision makers when they subjectively evaluate alternative and criteria scores for the MCDM problem. If decision makers set up a high degree of confidence, it means that they have gathered enough knowledge to support their decisions. ii. In our approach, the risk index β is applied as a defuzzifier. The crisp performance matrix H α β is calculated in equation 18 based on equation 17, where h α xβ (l) denotes the crisp performance score of alternative l with respect to the criterion C x under a degree of confidence α and a risk level β. In equation 17, when β continuously increases, the crisp performance score progressively approaches the left point of the interval. A risk index β = {0, 0.5, 1} respectively indicates an optimistic, moderate and pessimistic viewpoint of the decision maker choice.

α-cut Lx(l) h α xa (l) h α xb (l) ➘ ➘ Ux(l) Mx(l) x 0 0.5 1 u à (x) h α xa (l) = Lx(l) + α (Mx(l) -Lx(l)) (14)
h α xb (l) = Ux(l) -α (Ux(l) -Mx(l)) (15)
h α x (l) = [h α xa (l) h α xb (l)] (16)
h α xβ (l) = βh α xa (l) + (1 -β) × h α xb (l) (17)
H α β =             Ce Cc Cm Cs a1 h α eβ (1) h α cβ (1) h α mβ (1) h α sβ (1) a2 h α eβ (2) h α cβ (2) h α mβ (2) h α sβ (2) . . . . . . . . . . . . . . . an h α eβ (n) h α cβ (3) h α mβ (n) h α sβ (n)             (18)

Alternative ranking

In TOPSIS, the ideal solution h α+ xβ and the negative ideal solution h α- xβ are respectively defined as the best and the worst crisp performance scores among all alternatives on a criterion, as specified in equation 19. Subsequently, the distance between the ideal and the negative ideal solution for each alternative is respectively calculated in equations 20 and 21, where S α+ β (l) and S α- β (l) represent the distances between the crisp performance scores h α xβ (l) of an alternative l with respect to all criteria C x and all the ideal and negative ideal solutions h α+ xβ , h α- xβ , respectively.

S α- β (l) = x h α xβ (l) -h α- xβ 2 l = 1, 2, .., n (20)
A prior alternative has a longer distance to the negative ideal solution and a shorter distance to the ideal solution. Consequently, the relative closeness to the ideal solution for each alternative can be formulated as in equation 22, where R α β (l) denotes the final performance score. The larger the score R α β (l), the more relevant the alternative l. All alternatives/data items can therefore be ranked from the more relevant to the less relevant.

R α β (l) = S α- β (l) S α+ β (l) + S α- β (l) l = 1, 2, .., n (22)

Synthesis

In resume, when the machine/operator needs to store data on the product, the relevance of each data item from the database is computed thanks to the fuzzy AHP method. Data items are then ranked in order of relevance and information of the highest relevance/importance is therefore stored on the product thanks to the process step 3 which is detailed in the next section 4. Data item storage/retrieval using "communicating materials": Process step 3

As mentioned in the introduction, a prototype of a "communicating textile" was designed in our previous work [START_REF] Kubler | Prototyping of a communicating textile[END_REF], throughout which a large quantity of RFID µtags are spread (up to 1500µtags/m 2 ). This "communicating textile" is shown in Fig. 6. In this section, an appropriate architecture combined with a protocol of communication is developed to disseminate data items all over the material and then, to retrieve them.

RFID tag

Fig. 6. Prototype "communicating textile" designed in [START_REF] Kubler | Prototyping of a communicating textile[END_REF] As the RFID tags are memory-constrained, the idea is to split the set of data items among several tags. Fig. 7 depicts the architecture that must be implemented for storing the list of data items (identified thanks to process step 2) all over the "communicating textile" (split thanks to process step 3). A specific application protocol, named splitting protocol, is developed in [START_REF] Kubler | Storage of information on manufactured products using "communicating materials[END_REF] to split the data items. Indeed, a RFID tag may store more or less information according to the technology and, therefore, one data item may require more memory space than that available in a unique tag. This header enables to know in what order data items have been split. The next section presents the case study in which process step 2 (fuzzy AHP method) and process step 3 are applied so as to disseminate relevant information over the "communicating textile".

5.

Case study: Dissemination of information using a "communicating textile"

All along section 3, we detailed the set of computations needed to determine the fuzzy performance matrix H (see equations 28 and 29). Then, H α β and R α β (l) are calculated based respectively on equations 17-18 and equations 19 to 22 with α = 0.85 (decision makers judge they have gathered enough knowledge to overcome uncertainty in their evaluations) and β = 0.2 (reasonable to match the real situation). Let us remind ourselves that the higher the final score R α β (l), the more relevant the data item l. The results about the relevance R α β (l) of all data items (i.e. alternatives) are commented in section 5.1. Section 5.2 gives insight into the way data items are split on the "communicating textile" and then, rebuilt when reading.

Results about data relevance: Process step 2

Fig. 8(a) provides the resulting list of data items (which are product related data items) ordered from the highest relevant R α β (l) to the smallest. The two more relevant data items come from MaterialLot (R α β (l) = 0.8232 and 0.6524) and the 3 rd is T MD {3, 3}. Due to the large amount of data items included in the list (225 alternatives), the results are presented in the form of diagrams (whisker diagram and pie chart).

First, let us look at the whisker diagram in Fig. 8(b). For each table t ∈ T is given the min, the 1 st -3 rd quartile, the median and the max R α β (l) of the set of data items included in the list and belonging to t. We can see that the most relevant data items come from MaterialLot (which includes the two first data items of the list), but also from MaterialDefinition, ProductSegment and ProductionOrder. This is due to the fact that attributes from these four tables have been enumerated by the decision makers (cf. Table 2). We can see that D 1 and D 2 enumerate the same two attributes in MaterialLot, which explains why this table stands out among the other tables. Moreover, note that these four tables are included in G1 and G3 and that the experts, concerned by C c , have highly recommended to select information from both groups (cf. vectors in equation 10).

In our scenario, it is not possible to embed more than the first 155 data items from the list because the product memory is limited (as highlighted in Fig. 8(a)). The pie chart in Fig. 8(c) shows the percentage of data items among the 155 which belong to each entity group. We can observe that 49% (i.e. ≈ 75 data items out of 155) are included in tables clustered in G3 and 38% in tables clustered in G1. This is largely due to the choices made in the enumeration and contextual criteria as explained previously. In this respect, these results meet the expert specifications. The list, the whisker diagram and the pie chart are displayed to the user and can be used as decision-support tools. The last step consists in storing the first 155 data items on the product thanks to the process step 3 which is the subject of the following section.

Data Storage and retrieval: Process step 3

In our case study, the "communicating textile" designed in [START_REF] Kubler | Prototyping of a communicating textile[END_REF] (see Fig. 6) is used to store the first 155 data items from the list. Only two data items from the table MaterialDefinition must be stored on the product, namely T MD {3, 1} and T MD {3, 3} as shown in Fig. 9, because T MD {3, 2} got the 199 th rank (i.e. > 155). Fig. 9 exemplified the splitting process of these two data items in addition to the remaining 153 others, where the protocol header developed in [START_REF] Kubler | Storage of information on manufactured products using "communicating materials[END_REF] is added to each resulting piece of data item (gray background). These pieces of data items can therefore be stored/disseminated among the RFID tags.

The "communicating textile" continues its life cycle, and let us assume that a downstream user read the information carried by that one. The 155 data items stored in the "communicating textile" are retrieved, rebuilt, and then displayed . hea d.

Communicating textile

RFID tag

Fig. 9. Illustration of the splitting protocol on a mobile device via an application software (JAVA programming). The software actually recomposes all tables based on the set of data items retrieved from the material, enabling users to view each table and product-related tuples. For example, Fig. 10 shows the unique tuple of MaterialDefinition which has been rebuilt. This is the "product-related tuple" and more exactly, the data items T MD {3, 1} and T MD {3, 3} which were included in the first 155 th data items of the list. In this example, the software cannot displayed T MD {3, 2} (see Fig. 10) which may lead to unanswerable queries. Based on the application software, many services can be programmed: e.g. queries may be directly performed via the JAVA software without requiring access to the database.

Conclusion

In the framework of the product life cycle management, it is not uncommon to use intelligent/communicating products for ensuring an information continuum over the product life cycle -PLC (e.g. for traceability purposes). If one considers the product as an information vector (on which information could be stored), it would contribute to improve interoperability all along its life cycle. However, it is not that easy to identify, at a given stage of the PLC, what information should be stored on the product. To answer this question, a data dissemination process is developed in this paper to select context-sensitive information from the database and to store/replicate it on the product.

The approach proposed in this paper uses the fuzzy AHP theory for aggregating points of view from different actors/experts. An assumption is made about independencies among criteria, as required in AHP. Currently, no strong dependencies have been clearly identified among our criteria, even if some may be sensed. Indeed, since the product is memory-constrained, the number of enumerated data might, for instance, impact the maximal data size allowed to be stored on the product. Further experiments should therefore confirm or reveal such interdependencies, which cannot be captured by AHP. Accordingly, in future work, the application of ANP will be examined and results compared with those of this paper.

1 axFig. 2 .

 12 Fig. 2. Fuzzy AHP consisting of 5 stages

Fig. 3 .

 3 Fig. 3. General architecture of the hierarchy

4

  Let à = [a, b, c] and B = [d, e, f ] be two fuzzy numbers:

=

Fig. 4 .

 4 Fig. 4. Adjustment of k performed by an expert p

  j={e,c,m,s} sxj k={e,c,m,s} j={e,c,m,s} skj[START_REF] Yang | Key quality performance evaluation using fuzzy AHP[END_REF]

  φ(Ce) = [0.41, 0.53, 0.71] ; φ(Cc) = [0.21, 0.28, 0.37] (25) φ(Cs) = [0.07, 0.11, 0.18] ; φ(Cm) = [0.05, 0.07, 0.1]

  he(TMD{3, 3}) = φe(TMD{3, 3}) ⊗ φ(Ce) (29) = [0, 0.5, 1] ⊗ [0.41, 0.53, 0.71] = [0, 0.265, 0.71] hc(TMD{3, 3}) = φc(TMD{3, 3}) ⊗ φ(Cc) = [0.29, 0.52, 0.97] ⊗ [0.21, 0.28, 0.37] = [0.061, 0.146, 0.359]

Fig. 5 .

 5 Fig. 5. α-cut on each fuzzy performance score hlx

-h α+ xβ 2 l

 2 = 1, 2, .., n

Fig. 7 .

 7 Fig. 7. Global communication architecture

Fig. 8 .

 8 Fig. 8. Results of the data item's relevance 1 (★) 2 3 ID Material Description Value MDWP0 Wood plank with a nominal"... 4m of...1

Fig. 10 .

 10 Fig. 10. Display of tuples related to MaterialDefinition

  This corresponds to the data item noted T MD{3,1}The relevance value of T MD{3,1} is equal to 0.2(b)Table MaterialDefinition & data item relevance Fig. 1. View of a Logical Data Model (LDM) and a relational table

	ManufacturingBill ID ManBill ★ Description Quantity ID MatDef ✩	2	MaterialDefinition ★ ID MatDef Description Value	1	★	MaterialLot ID MatLot Description Status Quantity	1	1 (★) ID Material MDWP0	2 Description Wood plank with a nominal"...	3 Value 4m of...
	✩	ID ProdDef				✩	ID MatDef	2	MDP-B	Textile with a high developed... 3mm...
			3	★ : Primary Key (PK)	3	MD06	0.2	Textile which is provided...	0.05	15 • C	0.36
	PublishedDate Description ID ProdDef ★ ProductDefinition		x : relation number ✩ : Foreign key (FK)		4	MDH-V1	Vehicle headrests which conf...	2 × ...
		(a) Example of a Logical Data Model (LDM)

Table 1

 1 Main notations used in the paper {1, 2, .., pmax}; pmax varies according to the criterion x φx(l) fuzzy score (i.e. triangular fuzzy number) of the alternative l with respect to criterion x φc(Gi) fuzzy score (i.e. triangular fuzzy number) representing the relative importance of group G i over the other groups φ(Cx) fuzzy score (i.e. triangular fuzzy number) representing the relative importance of criterion Cx over the other criteria à fuzzy judgement matrix which gives all alternative fuzzy scores with respect to all criteria H fuzzy performance matrix synthesizing the decision matrix à and the fuzzy weight vector φ(Cx) ∀ x = {e, c, m, s}

	Notations

T

set of tables in the database. Let t be a table ∈ T , a data item from t is noted Tt{u, v} where u corresponds to the row/instance and v to the column/attribute l alternative l (i.e. a data item l) in the AHP structure, with l = {1, 2, .., n} Cx abbreviation for the criterion x with x = {e, c, m, s} (e: Enumeration, c: Contextual, s: Data Size, m: Model-based) G i abbreviation for a group of entity i defined in the LDM, with i = {1, 2, .., Z} D p decision maker p with p =

Table 2 "

 2 Table attributes" enumerated by D p

	D p Table Name (t)	Enumerated attributes (v)
		MaterialLot	{IDLot, Quantity}
	D 1	MaterialDefinition {IDMatDef, Value} ProductSegment {IDProdSeg, Duration, Unit}
		ProductionOrder	{IDProdOrd, StartTime, EndTime}
	D 2	MaterialLot ProductSegment	{IDLot, Quantity} {IDProdSeg, Description}

Table 3

 3 LDM entities grouped according to the type of information

	Product data group (G 1 ) Personal data group (G 2 )	Production data group (G 3 )	Equipment data group (G 4 )
	MaterialLot	Person	ProductionOrder	Equipment
	MaterialDefinition	PersonClass	ProductSegment	EquipmentClass
	MaterialClass	ActualPersonSegment	ProductDefinition	EquipmentSegmentSpecificat •
	ManufacturingBill	PersonSegmentSpecificat •	SegmentRequirement	ActualEquipmentSegment
	ActualMaterialLot		SegmentResponse
	MaterialSegmentSpecificat •

  List of data items ordered from the highest relevant to the smallest. The relevance of the data item l is given by R α iβ .

		Rang l 1 st 2 th	Table Name MaterialLot MaterialLot	Attribute Name Quantity IDMaterialLot	Primary Key LPB61 LPB61	Cell Value 50 LPB61	R α β (l) 0.8232 0.6524	Set of data items that must be stored on
		3 rd	MaterialDefinition	Value	MD06	15 • C	0.3639	the product
		. . . . . . 224 th 225 th	. . . . . . EquipmentClass PersonClass	. . . . . . Description Product storage limit: 155 th . . . . . . ECBSM01 Description PCPLW01	. . . . . . Bonding machi... Worker Line...	. . . . . . 0.0036 0.0029	Set of data items which cannot be stored on the product
	iβ Data item relevance R α	(a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 P e r s o n P e r s o n C la s s E q u ip m e n t E q u ip m e n t C la s s M a t e r ia lD e fi n it io n M a t e r ia lC M a t e r ia lL la s s P r o d u c t S e g m e n t o t M a t e r ia lS e g m e n t S p e c ifi Table name c a t io n E q u ip m e n t S e g m e n t S p e c ifi c a t io n P e r s o n S e g m e n t S p e c ifi c a t io n P r o d u c t D e fi n it io n M a n u f a c t u r in P r o d u c t io n O r d e r g B il l S e g m e n t R e q u ir e m e n t S e g m e n t R e s p o n s e A c t u a lM a t e r ia lL A c t u a lE q u ip m e n t S e g m e n t o t A c t u a lP e r s o n S e g m e n t	group (G1) Product data Production data group (G3)	49% 38%	9% 4%	(G4) data group Equipment Personal data group (G2)
		(b) Result of the data item's relevance in a whisker diagram

It consists in deriving weights from a pairwise comparison matrix.

B2MML is meant to be a common data format to link business enterprise applications with manufacturing enterprise applications.

S. Kubler et al. / Method for embedding context-sensitive information on "communicating textiles" via Fuzzy AHP

In our study, the network pointer carried by the product refers to MaterialLot which is therefore considered as the product table.


Jean-Philippe Braeunig
email: braeunig@unistra.fr

Nicolas Crouseilles

Virginie Grandgirard

Guillaume Latu

Michel Mehrenberger

Eric Sonnendrücker

Some numerical aspects of the conservative PSM scheme in a 4D drift-kinetic code

Keywords: numerical simulation, conservative scheme, maximum principle, plasma turbulence 2010 MSC: 65M08, 76M12, 76N99

The purpose of this work is simulation of magnetised plasmas in the ITER project framework. In this context, kinetic Vlasov-Poisson like models are used to simulate core turbulence in the tokamak in a toroidal geometry. This leads to heavy simulations because a 6D dimensional problem has to be solved, even if reduced to a 5D in so called gyrokinetic models. Accurate schemes, parallel algorithms need to be designed to bear these simulations. This paper describes the numerical studies to improve robustness of the conservative PSM scheme in the context of its development in the GYSELA code. In this paper, we only consider the 4D drift-kinetic model which is the backbone of the 5D gyrokinetic models and relevant to build a robust and accurate numerical method.

Introduction

The ITER device is a tokamak designed to study controlled thermonuclear fusion. Roughly speaking, it is a toroidal vessel containing a magnetised plasma where fusion reactions occur. The plasma is kept out of the vessel walls by a magnetic field which lines have a specific helicoidal geometry.

However, turbulence develops in the plasma and leads to thermal transport which decreases the confinement efficiency and thus needs a careful study.

Plasma is constituted of ions and electrons, which motion is induced by the magnetic field. The characteristic mean free path is high, even compared with the vessel size, therefore a kinetic description of particles is required, see Dimits [4]. Then the full 6D Vlasov-Poisson model should be used for both ions and electrons to properly describe the plasma evolution. However, the plasma flow in presence of a strong magnetic field has characteristics that allow some physical assumptions to reduce the model. First, the Larmor radius, i.e. the radius of the cyclotronic motion of particles around magnetic field lines, can be considered as small compared with the tokamak size and the gyration frequency very fast compared to the plasma frequency. Thus this motion can be averaged (gyro-average) becoming the so-called guiding center motion. As a consequence, 6D Vlasov-Poisson model is reduced to a 5D gyrokinetic model by averaging equations in such a way the 6D toroidal coordinate system (r, θ, φ, v , v ⊥ , α) becomes a 5D coordinate system (r, θ, φ, v , µ), with v the parallel and v ⊥ the perpendicular to the field lines components of the particles velocity, α the angular velocity around the field lines and µ = m v 2 ⊥ /2B the magnetic momentum depending on the velocity norm |v ⊥ |, on the magnetic field magnitude B and on the particles mass m which is an adiabatic invariant. Moreover, the magnetic field is assumed to be steady and the mass of electrons m e is very small compared to the mass of ions m i . Thus the cyclotron frequency ω i,e = q i,e B/m i,e is much faster for electrons than for ions ω e >> ω i . Therefore the electrons are assumed to be at equilibrium, i.e. the effect of the electrons cyclotronic motion is neglected and their distribution is then supposed to be constant in time. The 5D gyrokinetic model then reduces to a Vlasov like equation for ions guiding center motion:

∂ fµ ∂t + ∇ X • dX dt fµ + ∂ v dv dt fµ = 0 (1)
where fµ (X, v ) is the ion distribution function for a given adiabatic invariant µ with X = (r, θ, φ), velocities dX/dt and dv /dt define the guiding center trajectories. If ∇ (X,v ) •(dX/dt, dv /dt) t = 0, then the model is termed as conservative and is equivalent to a Vlasov equation in its advective form:

∂ fµ ∂t + dX dt • ∇ X fµ + dv dt ∂ v fµ = 0. ( 2
)
This equation for ions is coupled with a quasi-neutrality equation for the electric potential Φ(R) on real particles position, with R = X -ρ L (with ρ L the Larmor radius) :

- 1 Bω i ∇ ⊥ • (n 0 ∇ ⊥ Φ) + e κT e (Φ-< Φ > θ,φ ) = fµ dµdv -n 0 (3)
where n 0 is an equilibrium electronic density, T e the electronic temperature, e the electronic charge, κ the Boltzmann constant for electrons and ω i the cyclotronic frequency for ions.

These equations are of a simple form, but they have to be solved very efficiently because of the 5D space and the large characteristic time scales considered. This work is then a contribution in this direction, following

Grandgirard et al who develops the GYSELA 5D code that solves this 5D gyrokinetic model, see [START_REF] Grandgirard | Global full-f gyrokinetic simulations of plasma turbulence[END_REF] and [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. Looking at the model, one notices that the adiabatic invariant µ acts as a parameter. Therefore for each µ we have to solve a 4D advection equation as accurately as possible but also taking special care on mass and energy conservation, especially in this context of large characteristic time scales. The maximum principle that exists at the continuous level for the Vlasov equation should also be carefully studied at discrete level:

min i (f (x i , t n )) ≤ f (x, t n+1 ) ≤ max i (f (x i , t n ))
with f (x i , t n ) the value at x i in cell i at time t n . There is no physical dissipation process in the gyrokinetic model ( 1) that might dissipate over/undershoots created by the scheme and the loss of this bounding extrema of the solution at t n+1 > t n may even eventually crash a simulation.

Those studies will be achieved in this paper on a relevant reduced model, the 4D drift-kinetic model described in section 4, which has the same structure than equations [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF]. The geometrical assumptions of this model for ion plasma turbulence are a cylindrical geometry with coordinates (r, θ, z, v ) and a constant magnetic field B = B z e z , where e z is the unit vector in z direction. This 4D model is conservative and will be discretized using a conservative semi-Lagrangian scheme, the Parabolic Spline Method scheme (PSM, see Zerroukat et al [START_REF] Zerroukat | The parabolic spline method (PSM) for conservative transport problems[END_REF] and [START_REF] Zerroukat | Application of the parabolic spline method (PSM) to a multi-dimensional conservative semi-Lagrangian transport scheme (SLICE)[END_REF] ). It is a fourth order scheme which is equivalent for linear advections to the Backward Semi-Lagrangian scheme (BSL) currently used in the GYSELA code (see Grandgirard et al [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]) and introduced by Cheng-Knorr [2] and Sonnendrücker et al [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF]). This conservative PSM scheme based on the conservative form of the Vlasov equation will be described in section 4 and properly allows a directional splitting.

In this paper, the BSL and PSM schemes will be detailed with an emphasis on their similarities and differences. We will see that one difference is about the maximum principle. The BSL scheme satisfies it only with a condition on the distribution function reconstruction and the conservative PSM scheme does not satisfy it without an extra condition on the volumes conservation in the phase space. The last condition is equivalent to try to impose that the velocity field is divergence free at the discrete level. A scheme is given to satisfy this constraint in the form of an equivalent Finite Volume scheme. Moreover, we have designed a slope limiting procedure, Slope Limited Splines (SLS), to get closer to a maximum principle for the discrete solution, by at least diminish the spurious oscillations appearing when strong gradients exist in the distribution function profile.

The outline of this paper is the following : in section 2 will be recalled some important properties of Vlasov equations at the continuous level. Then BSL and PSM schemes will be described and compared, according to properties of the discrete solutions. In section 3, a numerical method will be given to improve the respect of the maximum principle by Vlasov discrete solutions when using the PSM scheme and particularly to keep constant the volume in the phase space. In section 4, practical aspects of the PSM scheme use will be described in the context of the 4D drift-kinetic model and at last we will comment on numerical results. with an arbitrary divergence free velocity field:

   ∂ t f + a • ∇ x (f ) = 0, ∇ • a = 0 and f (x, t) ≥ 0 (4)
with position x ∈ R D and a(x, t) ∈ R D the advection velocity field.

The solutions satisfy the maximum principle:

0 ≤ f (x, t) ≤ max x (f (x, t 0 )) (5)
for any initial time t 0 < t.

Since ∇ • a = 0, we can also use an equivalent conservative formulation of the Vlasov equation:

∂ t f + ∇ x • (a f ) = 0, (6)
For more details, see Sonnendrücker Lecture Notes [START_REF] Sonnendrücker | Lecture notes CEA-EDF-INRIA, Modèles numériques pour la fusion contrôlée[END_REF]. One obvious property of this conservation law (Reynolds transport theorem) is to conserve the mass in a Lagrangian volume V ol(t), by integrating the distribution function on each Lagrangian volume element dΩ:

d t m = d t V ol(t) f (x, t)dΩ = 0.
Let us introduce the convective derivative d t (.) = ∂ t (.) + a • ∇ x (.), thus [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF] becomes:

d t f + f ∇ x • a = 0. ( 7
)
Considering a Lagrangian motion of an infinitely small volume V ol(t), we have d t m = d t (f V ol) = 0, thus we obtain:

d t V ol V ol = ∇ x • a. (8)
Obviously, a divergence free flow ∇ x • a = 0 conserves a Lagrangian volume in its motion.

2.2. Maximum principle in the BSL and PSM schemes 2.2.1. Backward semi-Lagrangian (BSL)

Let us consider a Vlasov equation in its non conservative form:

∂ t f + a • ∇ x f = 0, (9)
with f (x, t) a scalar function, position x ∈ R D and a(x, t) ∈ R D the advection field. The BSL scheme, see Sonnendrücker et al [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF], is based on the invariance property of function f along characteristic curves to obtain values f n+1 at time t n+1 from the values f n at t n :

f n+1 X(x n+1 , t n+1 ) = f n X(x n+1 , t n ) , (10)
with x the Eulerian coordinates and the characteristic curves X defined as

dX(x, t) dt = a(x, t) (11)
with the initial position x = X(x, t n ) at t n . Let us locate the discrete function values

f n i = f n X(x n+1 i , t n+1 ) at mesh nodes x n+1 i = X(x n+1 i , t n+1 ).
We solve the following nonlinear system which is a second order approximation of d t X(t) = a(x, t):

X n+1/2 i = X(x n+1 i , t n+1 ) + X(x n+1 i , t n ) /2, X(x n+1 i , t n ) = X(x n+1 i , t n+1 ) -∆t a X n+1/2 i , t n+1/2 , f n+1 X(x n+1 i , t n+1 ) = f n h X(x n+1 i , t n ) , (12)
with ∆t = t n+1 -t n . The function f n h (x) is a reconstruction of the solution f n (x) according known values at nodes x n+1 i using cubic splines basis functions on the domain to obtain the value at x n i = X(x n+1 i , t n ), which is not a mesh node in general.

Properties of the BSL scheme. This scheme is formally fourth order in space.

It is second order in time using for instance a Leap-Frog, Predictor-Corrector or Runge-Kutta time integration. Mass is not conserved by this scheme, because it has no conservative form. However, an approximated maximum principle is satisfied. Let us consider f n h (x) the cubic spline interpolation of the distribution function f (x, t n ) at time t n , we have for any x:

f n+1 X(x, t n+1 ) = f n h (X(x, t n )) .
It then naturally appears a "discrete" maximum principle:

min x (f n h (x, t n )) ≤ f n+1 (x) ≤ max x (f n h (x, t n )). ( 13
)
Comparing with the property (5), we have here min

x (f n h (x, t n )) = 0 and max x (f n h (x, t n )) = max x (f n (x, t n ))
because the cubic spline reconstruction does not satisfy a maximum principle. If we have a manner to enforce this property to this reconstruction, a maximum principle is granted for the BSL scheme. No directional splitting is allowed since the BSL scheme is based on the non conservative form of the Vlasov equation, see [START_REF] Sonnendrücker | Lecture notes CEA-EDF-INRIA, Modèles numériques pour la fusion contrôlée[END_REF].

Semi-Lagrangian Parabolic Spline Method (PSM)

Let us consider a Vlasov equation in its conservative form:

∂ t f + ∇ x • (a f ) = 0, (14)
with f (x, t) a scalar function, position x ∈ R D and a(x, t) ∈ R D the advection field. Notice that with the hypothesis ∇ x • (a) = 0, conservative form (14) and non-conservative form (9) of the Vlasov equation are equivalent.

The PSM scheme, see Zerroukat et al [START_REF] Zerroukat | The parabolic spline method (PSM) for conservative transport problems[END_REF] and [START_REF] Zerroukat | Application of the parabolic spline method (PSM) to a multi-dimensional conservative semi-Lagrangian transport scheme (SLICE)[END_REF], is based on the mass conservation property of function f in a Lagrangian volume V ol to obtain the value f n+1 at time t n+1 :

V ol n+1 f (x, t n+1 )dΩ = V ol n f (x, t n )dΩ, (15)
with the characteristic curves X defined as dX(x, t) dt = a(x, t) and x n = X(x, t n ) with t n the initial time, and the volume V ol n = {X(x n+1 , t n ) such that X(x n+1 , t n+1 ) ∈ V ol n+1 } defined by the Lagrangian motion with the field a(x, t). The important point is that this conservative formalism properly allows a directional splitting without loosing the mass conservation.

Indeed, equation ( 14) may be solved with D successive 1D advections still of conservative form:

∂ t f + ∂ x k (a k f ) = 0, k ∈ [1, D]. (16)
We then approximate a 1D equation for each direction k using the conservation property. Omitting subscript k, the PSM scheme writes in 1D as follows:

x n+1 i+1/2 x n+1 i-1/2 f (x, t n+1 )dx = x n i+1/2 x n i-1/2 f (x, t n )dx, (17)
with x n+1 i+1/2 = X(x n+1 i+1/2 , t n+1 ) settled as the 1D mesh nodes and

x n i+1/2 = X(x n+1 i+1/2 , t n ) the associated foot of the characteristic curve, V ol n i = [x n i-1/2 , x n i+1/2 ] and V ol n+1 i = [x n+1 i-1/2 , x n+1 i+1/2
]. Let us define the unknowns of the scheme as the average of f in cell i

f n+1 i = 1 ∆x x n+1 i+1/2 x n+1 i-1/2 f (x, t n+1 )dx, (18)
and the primitive function

F n (z) = z x 1/2 f (y, t n )dy, (19)
with the uniform space step ∆x = x n+1 i+1/2 -x n+1 i-1/2 and x 1/2 an arbitrary reference point of the domain and for instance the first node of the grid {x i-1/2 } i=1,N +1 . Therefore, one has to solve a nonlinear system, which is similar to the BSL one, to obtain a discrete solution of equation ( 17) that writes:

X n+1/2 i+1/2 = X(x n+1 i+1/2 , t n+1 ) + X(x n+1 i+1/2 , t n ) /2, X(x n+1 i+1/2 , t n ) = X(x n+1 i+1/2 , t n+1 ) -∆t a X n+1/2 i+1/2 , t n+1/2 , f n+1 i ∆x = F n h (X(x n+1 i+1/2 , t n )) -F n h (X(x n+1 i-1/2 , t n )) (20)
with the time step ∆t = t n+1 -t n and the uniform space step ∆x = x n+1 i+1/2 -

x n+1 i-1/2 .
The computation of the reconstructed primitive function F n h (x) is based on values at mesh nodes x n+1 i+1/2 :

F n h (x n+1 i+1/2 ) -F n h (x 1/2 ) = i k=1 f n k ∆x.
Then this set of values is interpolated by cubic splines functions to obtain an approximated value F n h (z) of the primitive function F n (z) at any point z of the domain:

F n h (z) ≈ F n (z) = z x 1/2 f (y, t n )dy. ( 21
)
Properties of the PSM scheme. This scheme is formally fourth order in space and strictly equivalent to the BSL scheme for constant linear advection, see [3]. It is second order in time using for instance a Leap-Frog, Predictor-Corrector or Runge-Kutta time integration scheme. Mass is exactly conserved by this scheme for each 1D step k of the directional splitting. However, no maximum principle does exist for each step k even for the exact solution: in general ∂ x k a k = 0, even if the velocity field is divergence free ∇ • a = 0. Let us consider the scheme in D dimensions of space:

∂ t f + ∇ • (a f ) = 0, (22)
with f (x, t) a scalar function, position x ∈ R D and a(x, t) ∈ R D the advection field. Let us consider a cell i, where the solution is described at time t n+1 by its average in cell i, the PSM scheme then writes:

f n+1 i V ol n+1 i = V ol n+1 i f (x, t n+1 )dΩ = V ol n i f (x, t n )dΩ. ( 23
)
We thus obtain the following relation:

f n+1 i = f n * i V ol n i V ol n+1 i (24)
with the average of the distribution function in the Lagrangian volume at time t n :

f n * i = 1 V ol n i V ol n i f (x, t n )dΩ.
Here clearly appears two conditions, both difficult to satisfy especially in the context of a directional splitting, to have a maximum principle defined as follows:

min j (f n j ) ≤ f n+1 i ≤ max j (f n j ). ( 25
)
1. Maximum principle on the distribution function in V ol n i :

min j (f n j ) ≤ f n * i ≤ max j (f n j ). (26)
2. Conservation of volumes in the phase space at the discrete level:

V ol n i = V ol n+1 i . ( 27
)
The first condition is difficult to ensure in general, because a maximum principle should be satisfied for any average of the distribution function on an arbitrary volume V ol n i . Moreover, in the context of a directional splitting, it is impossible to satisfy a maximum principle for a 1D step k, because it does not exist at the continuous level since in general ∂ k a k = 0. Therefore it is probably impossible to recover a maximum principle of the reconstruction after all steps of the directional splitting.

The second condition is true at the continuous level while ∇ • a = 0, since we have d t V ol = V ol ∇ • a = 0, see equation ( 8) in section 2.1. As well as for the first condition, in the context of a directional splitting, it is difficult to ensure a constant volume evolution V ol n i = V ol n+1 i after all steps of the directional splitting, where compressions or expansions of the Lagrangian volume occur successively.

As a consequence, we will propose a form of the conservative PSM scheme that does not use a directional splitting. However we will not write the Semi-Lagrangian form of the PSM scheme in D dimensions of space because it is costly in computational time, because of the reconstruction step, and it is difficult to handle with arbitrary coordinate systems. The solution we choose is to use an equivalent Finite Volume form of the PSM scheme described in section 3.2, which is locally 1D at each face of the mesh. It is therefore possible to design 1D numerical limiters to try to better satisfy the maximum principle condition (26). Moreover, we will show that this form allows an exact conservation of the volumes in the phase space (27). The maximum principle and therefore the robustness of this scheme will thus be considerably improved.

3. Maximum principle for the PSM scheme

Numerical limiters for the distribution function reconstruction

Enforcing the first condition on the maximum principle of the distribution function reconstruction (26) can be really costly in computational time.

Instead of trying to correct the cubic spline reconstruction, we will reduce the spurious oscillations, generated by high order schemes when strong gradients appear in the distribution function profile, by using a classical Van Leer like slope limiting procedure, see for instance LeVeque [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. We propose here to measure the gradients in the flow and to add diffusion where they are detected. The diffusion is added by mixing the high order PSM flux with a first order upwind flux. The evaluation of the gradient is given by the classical function θ and we estimate the diffusion needed with a function γ(θ) ∈ [0, 1] based on a minmod like limiter function (see Fig. 1). The resulting limiter we propose here is called SLS (Slope Limited Splines), see [START_REF] Guterl | Test of some numerical limiters for the PSM scheme for 4D Drift-Kinetic simulations[END_REF] for details:

φ SLS i+1/2 = γ(θ i+1/2 ) φ P SM i+1/2 + (1 -γ(θ i+1/2 )) φ upwind i+1/2
where

φ upwind i+1/2 = a i+1/2 f n i + f n i+1 2 -sign(a i+1/2 ) f n i+1 -f n i 2 .
We define θ i+1/2 as the classical slope ratio of the distribution which depends on the direction of the displacement:

θ i+1/2 =          f n i -f n i-1 f n i+1 -f n i if a i+1/2 > 0 f n i+2 -f n i+1 f n i+1 -f n i if a i+1/2 < 0
However, the classical limiter minmod, where γ i+1/2 = max(0, min(θ i+1/2 , 1)), set γ to 0 when θ < 0. That means that the scheme turns to order 1 when an extrema exists, i.e. the slope ratio θ < 0. These extrema are thus quickly diffused and that leads to loose the benefits of a high order method. For SLS, the choice is to let the high-order scheme deal with the extrema and We also introduce a constant K in relation to control the maximum slope allowed without adding diffusion, i.e. mixing with the upwind scheme, see figure 1:

γ i+1/2 = max(0, min(K|θ i+1/2 |, 1)). ( 28
)
with the constant K = 5 experimentally settled.

We present in Fig. 2 the results of the linear advection of a step function with the standard PSM scheme with and without the SLS limiter (K=5).

The domain is meshed with 70 cells with periodic boundary conditions and the displacement is set to 0.2 cell per iteration. One can see that as any high order scheme, the PSM scheme produces spurious oscillations at the discontinuity or at a stiff gradient location. The SLS limiter do well with K = 5

to reduce these oscillations without introducing much diffusion. However, a maximum principle is not granted. This limiter has been further studied and compared with other limiters in the report [START_REF] Guterl | Test of some numerical limiters for the PSM scheme for 4D Drift-Kinetic simulations[END_REF].

Finite Volume form of the PSM scheme

Let us consider the 1D conservative advection equation of the form:

∂ t f + ∂ x (a x f ) = 0, (29)
with f (x, t) a scalar function, position x ∈ R and a x (x, t) ∈ R the advection field. We recall that X(x n+1 i+1/2 , t n+1 ) = x i+1/2 is the position of the mesh node i + 1/2. Let us set the notation X(x n+1 i+1/2 , t n ) = x * i+1/2 for the "foot" position at t n on the characteristic curve. Let us rewrite the PSM scheme (20):

x * i+1/2 = x i+1/2 -∆t a x x n+1/2 i+1/2 , t n+1/2 , f n+1 i ∆x = F n h (x * i+1/2 ) -F n h (x * i-1/2 ).
with the primitive function F n h (z) is interpolated by cubic splines at mesh nodes such that

F n h (z) ≈ F n (z) = z x 1/2
f n (y)dy at fourth order in space.

Let us make appear explicitly in 1D the fluxes at cell faces i ± 1/2, by introducing the primitive values at cell faces F n (x i±1/2 ):

f n+1 i ∆x = F n h (x * i+1/2 ) -F n h (x i+1/2 ) -F n h (x * i-1/2 ) -F n h (x i-1/2 ) +f n i ∆x with f n i ∆x = F n h (x i+1/2 ) -F n h (x i-1/2 ). It yields f n+1 i -f n i ∆t + F n h (x i+1/2 ) -F n h (x * i+1/2 ) ∆x∆t - F n h (x i-1/2 ) -F n h (x * i-1/2 ) ∆x∆t = 0. ( 30
)
The PSM fluxes at cell faces i ± 1/2 clearly appear:

f n+1 i -f n i ∆t + Φ P SM i+1/2 -Φ P SM i-1/2 ∆x = 0 (31) with Φ P SM i+1/2 = F n h (x i+1/2 ) -F n h (x * i+1/2 ) ∆t ≈ 1 ∆t x i+1/2 x * i+1/2 f n (y)dy. (32)
A simple Taylor expansion shows that this PSM flux, which consists in a cubic spline approximation of the integral of f n (x) along the characteristic curves at cell faces, is a consistent approximation at node i + 1/2 of the continuous flux Φ = a x f in equation (29), i.e. Φ P SM i+1/2 ≈ (a x f ) i+1/2 . Moreover, this flux is an approximation of the integral of (a x f n on cell faces.

Coming back to (29) and integrating on the cell volume V ol n+1 i = A x ∆x, with A x the bounding faces Γ i±1/2 area transversal to x of cell i:

V ol n+1 i ∂f ∂t + V ol n+1 i ∂ x (a x f )dΩ = 0. ( 33
)
We obtain using Green formula:

∂f ∂t + 1 V ol n+1 i ∂V ol n+1 i =Γ i-1/2 ∪Γ i+1/2 f (a x • n x )dΓ = 0. ( 34
)
with ∂V ol n+1 i = Γ i-1/2 ∪ Γ i+1/2 the surface bounding V ol n+1 i and n x its outgoing normal.

∂f ∂t + 1 V ol n+1 i Γ i+1/2 f a x dΓ - Γ i-1/2 f a x dΓ = 0. ( 35
)
Comparing with formula (31), we see that the PSM flux is an approximation of the flux average at cell faces:

Φ P SM i+1/2 ≈ 1 A x Γ i+1/2 f (a x • n x )dΓ (36) because V ol n+1 i = A x ∆x.
The extension to D dimensions of space is then straightforward, because it only consists in adding the fluxes through the faces of a cell in every direction d. Considering the Vlasov equation in its conservative form in dimension D and for an arbitrary coordinate system:

∂(Jf ) ∂t + ∇ • (Ja f ) = 0
with J the geometric Jacobian of the cell, f (x, t) a scalar function, position

x ∈ R D and a(x, t) ∈ R D the advection field. In a classical way in Finite

Volume methods, we integrate this local equation on the cell C i of volume V ol n+1 i and we use the Green formula:

V ol n+1 ∂f ∂t Jdx + Γ d ∈C i (a • n d )f Jdx = 0, (37)
with Γ d the face of cell i perpendicular to direction d of area A d and of outgoing normal unit vector n d . Using the 1D flux formula (36) in direction n d and a first order time discretisation, we obtain the following Finite Volume scheme:

V ol n+1 i f n+1 i -f n i ∆t + Γ d ∈C i A d Φ P SM d = 0, (38)
with

f n i = 1 V ol n+1 i V ol n+1 i f (y, t n )J(y)dy. ( 39
)
and Φ P SM d the flux that goes through Γ d :

Φ P SM d = F n h,d (x d ) -F n h,d (x * d ) ∆t ≈ 1 A d Γ d (a • n d )f (y, t n )J(y)dy, (40)
with x * d = x d -∆t (a(x d ) • n d
) the foot of the 1D characteristic curve and F n h,d (z) the primitive function at t n reconstructed using cubic splines in the direction of n d :

F n h,d (z) ≈ F n d (z) = z d x d,1/2 f (y, t n )J(y)dy. ( 41
)
We see here a Finite Volume form of the Semi-Lagrangian PSM conservative scheme. This equivalence is however restricted, because this Finite Volume form is submitted to a CFL condition as any scheme of this form:

∆t ≤ min d   ∆x d max x d (a n d (x d ))   .
Moreover, as we will see in section 3.3, the Lagrangian volume evolution is here approximated by the cell faces motion only in their normal direction, instead of the general motion as it is described in the Semi-Lagrangian formalism (15). It is the same volume evolution as the Semi-Lagrangian method with 1D directional splitting. However it is the classical Finite Volume formalism and it is the key point that will permit to enforce a divergence free evolution of the flow.

Notice that the Finite Volume form (38) can be directionally split keeping exactly the same result. Indeed, it only consists in adding the flux in two successive operations instead of in one. As an example, let us consider the 2D (x, y) case of a cartesian mesh:

V ol n+1 i f n+1 i -f n i ∆t +A x i+1/2,j Φ P SM i+1/2,j + A x i-1/2,j Φ P SM i-1/2,j
+A y i,j+1/2 Φ P SM i,j+1/2 + A y i,j-1/2 Φ P SM i,j-1/2 = 0.

(42)

It is stricly equivalent to use the directional directional splitting:

V ol n+1 i f nx i -f n i ∆t +A x i+1/2,j Φ P SM i+1/2,j + A x i-1/2,j Φ P SM i-1/2,j = 0 V ol n+1 i f n+1 i -f nx i ∆t +A y i,j+1/2 Φ P SM i,j+1/2 + A y i,j-1/2 Φ P SM i,j-1/2 = 0. ( 43
)
with the only condition that all fluxes Φ P SM are computed using f n and a n at time t n as in the unsplit scheme (42).

Conservation of volumes in the phase space

The second condition to have a maximum principle for the PSM scheme is to satisfy the multi-dimensional condition (27) of conservation of volumes, i.e. V ol n = V ol n+1 . Equation [START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system[END_REF] showed that at continuous level the volume is constant in its evolution in the phase space if the advection field is divergence free. Therefore we will study the PSM scheme to find out a divergence free condition that should be satisfied at the discrete level ∇ h • a = 0 in such a way V ol n = V ol n+1 , in the same way ∇ • a = 0 at the continuous level. With the idea of making appear the total evolution of volumes between t n and t n+1 , we will use the Finite Volume form of the PSM scheme given in section 3.2 for the 2D polar coordinate system. Let us consider radial r and orthoradial θ directions, with constant space steps ∆r, ∆θ and the volume of the cells V ol i,j = r i ∆r∆θ with r i the mean radius of the cell. We consider that the mesh in polar coordinates has locally no curvature, i.e. each mesh is a trapezium with straight edges. This is important to be noticed to write the Finite Volume scheme and calculate the volume swept by the cell edges, see Fig. 3. Let us set a velocity field (a r (r, θ), r a θ (r, θ)) such that:

∇ r,θ • a = 1 r ∂ r (r a r ) + 1 r ∂ θ (r a θ ) = 0.
Let us write the conservative advection equation in polar coordinates:

∂ t (rf ) + ∂ r (r a r f ) + ∂ θ (r a θ f ) = 0. ( 44
)
Notice that the geometric Jacobian J = r for polar coordinates.

The PSM scheme without directional splitting in the Finite Volume form (38) reads here:

V ol i,j f n+1 i,j -f n i,j ∆t +A r i+1/2,j Φ P SM,r i+1/2,j -A r i-1/2,j Φ P SM,r i-1/2,j +A θ i,j+1/2 Φ P SM,θ i,j+1/2 -A θ i,j+1/2 Φ P SM,θ i,j-1/2 = 0 (45)
with Φ P SM,θ i,j±1/2 and Φ P SM,r i±1/2,j positioned at cell faces center and with cell i, j of volume V ol i,j = r i ∆r∆θ and faces areas A r i±1/2,j = ∆θ and A θ i,j±1/2 = ∆r. The cell averaged values of f used in the scheme are:

f i,j = 1 V ol i,j V ol i,j f (r, θ, t) rdrdθ. ( 46
)
Using the integral form (40) of the fluxes:

r i ∆r∆θ(f n+1 i,j -f n i,j ) +∆θ r i+1/2 r * i+1/2 (r, θ i , t n )rdr -∆θ r i-1/2 r * i-1/2 f (r, θ i , t n )rdr +∆r θ j+1/2 θ * j+1/2 f (r i , θ, t n )r i dθ -∆r θ j-1/2 θ * j-1/2 f (r i , θ, t n )r i dθ = 0. ( 47
)
Let us introduce the volumes swept by each cell face in its normal motion in accordance with the Green formula and the way of computation of feet of characteristic curves normal to cell faces, i.e. without taking into account the tangential motion at the cell faces or the curvature of the mesh, see Fig. 3:

δV ol r i±1/2,j = r i±1/2 ∆θ (r i±1/2 -r * i±1/2 ) δV ol θ i,j±1/2 = r i ∆r (θ j±1/2 -θ * j±1/2 ).
Therefore we obtain:

V ol i,j f n+1 i,j rdrdθ = V ol n i,j f (r, θ, t n ) rdrdθ (48)
with

V ol n i,j = V ol i,j -δV ol r i+1/2 + δV ol r i-1/2 -δV ol θ j+1/2 + δV ol θ j-1/2 .
We here recover a discrete mass conservation formulation. To obtain V ol i,j = V ol n i,j , and thus preserve a constant function, it yields:

δV ol r i+1/2 -δV ol r i-1/2 + δV ol θ j+1/2 -δV ol θ j-1/2 = 0.
Using δV ol k definitions, we thus obtain a discrete divergence formulation in polar coordinates ∇ h • a = 0 to be nullified:

1 r i r i+1/2 a r (r i+1/2 , θ j ) -r i-1/2 a r (r i-1/2 , θ j ) ∆r + 1 r i r i a θ (r i , θ j+1/2 ) -r i a θ (r i , θ j-1/2 ) ∆θ = 0. ( 49
)
with the following first order definition of the characteristic curves feet computation:

r i±1/2 -r * i±1/2 = ∆t a r (r i±1/2 , θ j ) and θ j±1/2 -θ * j±1/2 = ∆t a θ (r i , θ j±1/2 ).
As a conclusion, we have presented a general methodology for any coordinate system to compute the associated discrete divergence free condition, by using the approximation of cell edges by straight lines and by only considering the normal to cell faces motion of the volume as it has to be when invoking the Green formula in this Finite Volume framework. The discrete divergence formulation (49) is independent of the time integration method. It is a discrete consistent relation for the advection field of the form ∇ h • a = 0. It should be satisfied to get the conservation condition on volumes V ol n i = V ol n+1 i in the phase space, which is necessary to obtain a maximum principle for the PSM scheme or actually for any Finite Volume scheme. This condition is also necessary when using the Semi-Lagrangian PSM scheme with directional splitting as described in section 2.2.2 as well as when using the Finite Volume form described in section 3.2. In Fig. 4, we compare the results of a 4D drift-kinetic benchmark (see section 4.3.1

for details) obtained with the Semi-Lagrangian PSM scheme 2.2.2 with an advection field computed: first in such a way the discrete divergence free condition (49) is satisfied and second with an advection field computed by cubic spline interpolation without satisfying this condition. Respecting condition (49) for the advection field not only leads to a better respect of the maximum principle, it is actually necessary to ensure the stability of the scheme. The result in figure 4 diverges from realistic physics.

4. Use of the PSM scheme in a 4D drift-kinetic code

Drift-kinetic model

This work follows those of Grandgirard et al in the GYSELA code, see [START_REF] Grandgirard | Global full-f gyrokinetic simulations of plasma turbulence[END_REF] and [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. The geometrical assumptions of this model for ion plasma turbulence are a cylindrical geometry with 4D coordinates (r, θ, z, v ) and a constant magnetic field B = B z e z , where e z is the unit vector in z direction. The model is the 4D Drift-Kinetic equations described inGrandgirard et al [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]:

dr dt = v GCr ; r dθ dt = v GC θ ; dz dt = v ; dv dt = q i m i E z (50)
with v GC = (E × B)/B 2 and E = -∇Φ with Φ the electric potential.

The 4D Vlasov equation governing this system, where the ion distribution function is f (r, θ, z, v , t), is the following:

∂ t f + v GCr ∂ r f + v GC θ ∂ θ f + v ∂ z f + q i m i E z ∂ v f = 0. ( 51
)
This equation is coupled with a quasi-neutrality equation for the electric potential Φ(r, θ, z) that reads:

-∇ ⊥ Φ • n 0 (r) B Ω 0 ∇Φ + e n 0 (r) T e (r) (Φ-< Φ > θ,z ) = n i -n 0 (52)
with n i = v f (r, θ, z, v )dv and constant in time physical parameters n 0 , Ω 0 , T e and e. Let us notice that the 4D velocity field a = (v GCr , v GC θ , v , q/m i E z ) t is divergence free:

∇ • a = 1 r ∂ r (r v GCr ) + 1 r ∂ θ (v GC θ ) + ∂ z v + ∂ v (q/m i E z ) = 0 (53) because of variable independence ∂ v E z = ∂ v (∂ z Φ(r, θ, z)) = 0 and ∂ z v = 0
and we have v GC = (E × B)/B 2 , with E = -∇Φ and B = B z e z , thus

v GCr = 1 B z - 1 r ∂ θ Φ and v GC θ = 1 B z (∂ r Φ) (54) and ∇ r,θ • a = 1 r ∂ r (r v GCr ) + 1 r ∂ θ (v GC θ ) = 1 r B z (∂ r (r (-1/r)∂ θ Φ) + ∂ θ (∂ r Φ)) = 0. (55)
Therefore, one can write an equivalent conservative equation to the preceding Vlasov equation (51):

∂ t f + ∂ r (v GCr f ) + ∂ θ (v GC θ f ) + ∂ z (v f ) + ∂ v q i m i E z f = 0 (56)

Computation of a divergence free velocity field at the discrete level

We have obtained a discrete form of the velocity field divergence to nullify (49), as a necessary condition to obtain a numerical solution with a maximum principle. We saw in (53) that ∇ • a = 0 is satisfied equivalently if ∇ rθ • a = 0 (55) is satisfied and this is still true at the discrete level (independence of variables). Therefore, the velocity field should nullify the discrete polar divergence (49):

1 r i r i+1/2 a r (r i+1/2 , θ j ) -r i-1/2 a r (r i-1/2 , θ j ) ∆r + 1 r i r i a θ (r i , θ j+1/2 ) -r i a θ (r i , θ j-1/2 ) ∆θ = 0 (57) with a r = dr/dt = v GCr = -1 r B z ∂ θ Φ and a θ = dθ/dt = v GC θ /r = 1 r B z ∂ r Φ,
using definitions given in (54).

Proposition 1.

Let us define the electric potential at the nodes of the mesh Φ i+1/2,j+1/2 , whatever the way it is computed. Let us set the following natural finite difference approximation for the velocity field:

a r (r i+1/2 , θ j ) = -1 r i+1/2 B z Φ i+1/2,j+1/2 -Φ i+1/2,j-1/2 ∆θ a θ (r i , θ j+1/2 ) = 1 r i B z Φ i+1/2,j+1/2 -Φ i-1/2,j+1/2 ∆r . ( 58
)
With this approximated velocity field, the approximation of ∇ rθ • a = 0 (57) is satisfied.

The proof is easy, we just have to put the velocity field (58) in (57) to see that all terms annulate each others.

Remark 2. Notice that the electric potential Φ should be computed at nodes (i ± 1/2, j ± 1/2) of the mesh to obtain velocities at the center of cell faces (i ± 1/2, j) and (i, j ± 1/2). It is well adapted to the PSM schemes, where the displacement should be calculated at cell faces. In this section, we will compare the numerical methods on a 4D driftkinetic benchmark, following the paper Grangirard et al [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. The model is described in section 4.1. We will compute the growth of a 4D unstable turbulent mode. The benchmark consists of exciting the plasma mode (m, n), with m the poloidal mode (θ) and n the toroidal mode (z). The initial distribution function is the sum of an equilibrium and a perturbation distribution function f = f eq + δf . The equilibrium distribution function has the following form:

f eq (r, v ) = n 0 (r) (2πT i (r)/m i ) 1/2 exp - m i v 2 2T i (r) (59)
and the perturbation δf δf (r, θ, z, v ) = f eq (r, v ) g(r) h(v ) δp(θ, z)

with g(r) and h(v ) two exponential functions and δp(θ, z) = cos 2πn L z z + mθ with L z the length of the domain in z direction, m i , T i (r), n 0 (r) physical constant profiles, see [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF] for details. We have set here m = 16 and n = 8.

Algorithm

At the beginning of the time step, the distribution function f (x, v , t n ) is known at time t n , with x = (r, θ, z). The time step ∆t = t n+1 -t n is computed at each step with the CFL like condition: The operator splitting between the quasi-neutral equation and the Vlasov transport equation is made second order using a Predictor-Corrector scheme in time:

1. Time step ∆t computation.

2. Quasi-neutral equation (52) solving at t n using the distribution function f n (actually the density) to obtain the electric potential Φ n (x) at time t n+1/2 . The advection field a(x, t n ) is computed with Φ n (x) according to equation (53) and using formula (58) .

3. 4D Vlasov equation solving at t n with time step ∆t/2 to obtain the distribution function f n+1/2 (x) at time t n+1/2 using the advection field a(x, t n ).

4. Quasi-neutral equation (52) solving at t n+1/2 using the distribution function f n+1/2 (actually the density) to obtain the electric potential Φ n+1/2 (x) at time t n+1/2 . The advection field a(x, t n+1/2 ) is computed with Φ n+1/2 (x) according to equation (53) and using formula (58) .

5. 4D Vlasov equation solving at t n with time step ∆t to obtain the distribution function f n+1 (x) at time t n+1 using the advection field a(x, t n+1/2 ).

In the two following paragraphs, we describe the schemes for the 4D

Vlasov equation (56) solving of the algorithm with ∆t * = ∆t/2 in the prediction step and ∆t * = ∆t in the correction step.

4D Semi-Lagrangian PSM sheme with directional splitting.

• PSM 1D advection of f (x, v , t n ) in direction v with velocity a n v and time step ∆t * /2 to obtain f (x, v , t v /2 ).

• PSM 1D advection of f (x, v , t v /2 ) in direction z with velocity a n z and time step ∆t * /2 to obtain f (x, v , t z/2 ).

• PSM 1D advection of f (x, v , t z/2 ) in direction θ with velocity a n θ and time step ∆t * /2 to obtain f (x, v , t θ/2 ).

• PSM 1D advection of f (x, v , t θ/2 ) in direction r with velocity a n r and time step ∆t * to obtain f (x, v , t r ).

• PSM 1D advection of f (x, v , t r ) in direction θ with velocity a n θ and time step ∆t * /2 to obtain f (x, v , t θ ).

• PSM 1D advection of f (x, v , t θ ) in direction z with velocity a n z and time step ∆t * /2 to obtain f (x, v , t z ).

• PSM 1D advection of f (x, v , t z ) in direction v with velocity a n v and time step ∆t * /2 to obtain f

(x, v , t v ) = f (x, v , t n+1 ).
Each PSM 1D advection is achieved using the standard 1D semi-Lagrangian PSM scheme as described in section 2.2.2. The directional splitting is second order by using a Strang like decomposition. Since we use here a directional splitting and a second order scheme in time for the computation of the characteristic curves, the volumes are not strictly conserved in the phase space, because the scheme does not satisfy the discrete divergence free condition (49).

4D Finite Volume form of the PSM scheme:.

• PSM 1D advection of f (x, v , t n ) in direction v with velocity a n v and time step ∆t * /2 to obtain f (x, v , t v /2 ).

• PSM 1D advection of f (x, v , t v /2 ) in direction z with velocity a n z and time step ∆t * /2 to obtain f (x, v , t z/2 ).

• PSM 2D advection of f (x, v , t z/2 ) in each plane (r, θ) with velocities (a n r , a n θ ) and time step ∆t * to obtain f (x, v , t r,θ ).

• PSM 1D advection of f (x, v , t r,θ ) in direction z with velocity a n z and time step ∆t * /2 to obtain f (x, v , t z ).

• PSM 1D advection of f (x, v , t z ) in direction v with velocity a n v and time step ∆t * /2 to obtain f (x, v , t v ) = f (x, v , t n+1 ).

Each PSM 1D advection is achieved using the standard 1D semi-Lagrangian scheme as described in section 2.2.2. The PSM 2D advection in (r, θ) is achieved with the Finite Volume form as described in section 3.2. Since we use here the scheme 3.2, the volumes are strictly conserved in the phase space, because the scheme does satisfy the discrete divergence free condition (49). Even if we use the semi-Lagrangian PSM 1D advection in z and v directions, the property is kept because the velocity is constant in these directions.

Results

The mesh is 128 × 256 × 32 × 16 cells in r, θ, z, v directions. Boundary conditions are periodic for directions θ and z and Neumann (∂f /∂n = 0) in r and v . We first ran the reference test case with the non-conservative Backward Semi-Lagrangian (BSL) scheme in section 2.2.1 which is currently used in the GYSELA code. Then we ran four test cases to show the influence of each numerical treatment: the standard conservative semi-Lagrangian PSM scheme in section 4.3.2 with 1D directional splitting (PSM Directional Splitting 1D) and the same with the SLS limiter (SLS Directional Splitting 1D), the unsplit Finite Volume form of the PSM scheme in section 4.3.2 (PSM Finite Volume) and the same with the SLS limiter (SLS Finite Volume).

The computed 4D distribution functions are pictured in Fig. 5 at time t = 1800 and in Fig. 6 at time t = 4400. We only present 2D slices (r, θ) of the distribution function at v = 0 and for a given value of z = z 0 . In these figures X stands for r direction and Y for the θ direction. At initial time t = 0, the minimum and maximum values of the distribution function in this slice are (min . = 0.331, max . = 0.4187) and these values should be the same at any time of the computation if the maximum principle would be respected. In Fig. 5, we show pictures of each scheme result at time t = 1800, which corresponds approximately to the beginning of the non-linear turbulent phase saturation. Small structures are appearing and interact with each others. All results are still close qualitatively. However, we already see oscillations in the solution obtained with PSM DS (Directional Splitting), where the minimum and maximum values (min . = 0.3086, max . = 0.4427) are already quite different than the one at initial time. The PSM Finite Volume (PSM FV) form and the SLS DS better keep these extrema, but only SLS Finite Volume keep the extrema unchanged until time t = 1800 with a really similar behaviour of the solution. In Fig. 6, we show pictures of each scheme result at time t = 4400 when turbulence is well developed. We see that the standard PSM scheme creates a lot of unphysical oscillations SLS Directional Splitting SLS Finite Volume (the SLS limiter does not provide a maximum principle) and the solution is smooth. We may say that the added diffusion with the limiter helps the scheme to diffuse subgrid structures without creating oscillations. The divergence free property of the Finite Volume scheme is important to cure to solution from instabilities that can be seen at r values close the average value of r (vertical line at the middle of pictures in Fig. 6) in the SLS Directional Splitting solution compare to the SLS Finite Volume solution. In Fig. 7, we see in the reference BSL solution at time t = 1800 spurious oscillations produced during the reconstruction step of the distribution function, which is the only possibility to break the maximum principle for the BSL scheme:

here extrema are (min . = 0.3124, max . = 0.4430) instead of values at initial time (min . = 0.331, max . = 0.4187). At time t = 4400, we see spurious oscillations as well, but the maximum principle is better satisfied than with the standard PSM DS scheme in Fig. 6, because no conservation of volumes in the phase space has to be satisfied, as it is explained in section 2.2.1.

Conclusion and perspectives

The PSM scheme has been successfully integrated in the GYSELA code and has been tested on 4D Drift-Kinetic test cases. We had first experimentally stated and afterward explained in this paper that the PSM scheme can be unstable without taking care of a velocity field divergence free condition.

The numerical results show that the study of the volume evolution in the SLS Directional Splitting SLS Finite Volume

2 .

 2 Semi-Lagrangian schemes for Vlasov equation 2.1. Basics of the Vlasov equation Let us consider an advection equation of a positive scalar function f (x, t)

Figure 1 :

 1 Figure 1: γ function for the SLS limiter

Figure 2 :

 2 Figure 2: Linear advection of a step function. The exact solution is green and the numerical results are in blue, above the standard PSM scheme and below the PSM scheme with SLS limiter with K = 5.

Figure 3 :

 3 Figure 3: Polar cells in black and the Finite Volume representation of the cells in red as a trapezium with the volumes swept by the cell faces in their normal motion. The direction for all faces motion is drawn outward, but it could be inward as well, function of the velocity field.

Figure 4 :

 4 Figure 4: Result at time t = 60 with (left) the advection field computed in a way (see section 4.3.1) that satisfy the discrete divergence condition (49) and with (right) the advection field computed with cubic splines, which do not satisfy this condition (49).

4. 3 .

 3 Numerical tests 4.3.1. Drift-kinetic 4D model, PSM schemes comparison

  CF L r = CF L θ = 0.5, because the flow is highly nonlinear in (r, θ) planes thus characteristics should not cross each others during one time step, and CF L z = CF L v = 8 because it is linear advection in direction z and v so characteristics can not cross each others and then we allow a maximum displacement of 8 cells. Actually excluding the linear phase, the most restrictive directions for the time step are r and θ, in such a way this last value (8) has a minor importance compare to the leading parameters CF L r = CF L θ = 0.5.

Figure 5 :

 5 Figure 5: Simulation with 128x256x32x16 cells -PSM Directional Splitting 1D (up-left) -PSM Finite Volume (up-right) -SLS Directional Splitting 1D (down-left) -SLS Finite Volume (down-right) -time =1800.

Figure 6 :

 6 Figure 6: Simulation with 128x256x32x16 cells -PSM Directional Splitting 1D (up-left) -PSM Finite Volume (up-right) -SLS Directional Splitting 1D (down-left) -SLS Finite Volume (down-right) -time =4400, with all color tables set to the minimum and maximum value at initial time.

BSL time=1800 BSL time=4400

Figure 7: Reference BSL Simulation with 128x256x32x16 cells -BSL at time=1800 with the real color table values (left) -BSL at time=4400 with the color table set to the minimum and maximum value at initial time (left).