Software Coupling and Orchestration Tool to the Modeling of Multi-physic Problems L. Mokhtari, B. Delinchant, J-L. Coulomb, T. Le duc, and P. Pham Quang

Abstract-We present in this paper Scot, which is a modular software solution for weakly coupling models, methods and orchestration of the simulation. The goal behind Scot is to make easier the phase of modeling and optimize the phase of simulation by tuning the solvers simulation parameters. Specifications and composition of Scot will be described. Scot has been successfully validated by two different applications, the PEEC-MoM coupling to the modeling of an electromagnetic device and the magnetomechanic coupling to the modeling of deformable nano-switch contact NEMS.

Index Terms-Software component, weak coupling, Tuning solvers, Orchestration, PEEC and MoM method, and magnetomechanic coupling.

I. INTRODUCTION

T HE PROBLEM in multidisciplinary systems modeling is that researchers attempt to model all the different parts of the system that don't appear from their research field, so the models they develop can be of poor quality and may contain several errors of modeling.

In this context, it appears more interesting to use specialized tools to carry out the different models and then link them to rebuild the overall representation of the complex system [START_REF] Mokhtari | Comparing Weak and Strong PEEC-MoM Coupling[END_REF].

Scot (Software Coupling and Orchestration Tool) is a software tool with a modular architecture developed for the purpose of coupling weakly models, methods and specialized tools. It includes some interfaces and makes easy to develop others as required. Models written in any language and embedded in Java, can be easily linked to Scot.

Scot has been successfully tested on two applications. The first concerns the modeling of a real electro-magnetic device by coupling the PEEC (Partial Elements Equivalent Circuit method) and MoM (Method of Moments) method with using a Matlab-Scot coupling. The second test example is a magnetomechanic coupling applied to the modeling of a deformable switch beam NEMS (Nano Electro-Mechanical Systems).

II. SPECIFICATIONS OF SCOT

Coupling strategies are described by two ways in Scot; in the first, the models are implemented directly in Scot and both of the control and the simulation are done on Scot, however, in the second, there are two layers as shown in (Fig. 1). The layer of supervision from where we configure and control the coupling and, the layer of simulation where the L. Mokhtari

Simulation layer

Supervision layer

Fig. 1. Scot Coupling scheme of two models using CT (Coupling Terms).

A. Scot elements 1) Variables: State variables, are tokens consumed or produced by an ScotComponent or ScotTransformer, characterized by its type (integer, real, boolean. . . ), size, and value. Among Scot variables, there are also parameters of configuration and control.

2) Components: ScotSolver is the basic component of Scot, it includes a path to load some type of files, an initialization and solving method, and a name of the recovered variable. ScotSolver is composed by sComp Matlab and sComp Java. The first is an interface that allows to couple external Matlab models (Fig. 3(a)), this type of component is used in the application VI-A. The second one calls models written in java language and can use models in ICAr1 (Fig. 3(b)) components norm as done in the application VI-B.

Socket

Model Model Model Model

Configuration Control

Input Output 3) Composite component: The sComp Composite (Scot composite component) is a class of ScotSolver, it is specified in the goal to build easily complex coupling algorithms. It allows us to couple ScotSolvers sequentially with or without feedback (Fig. 4). With using sComp Composite, some parts of a complex algorithm are encapsulated and seen as a simple ScotSolver from the ScotCoplingTool. 4) Transformers: Since the models are not necessarily designed to be coupled thus state variables produced by one model can not be consumed by another. In this case, it is necessary to use coupling terms "ScotTransformer". These ScotTransformers may also be transition matrices to reduce the order of the state variables or, software that performs interpolation between different scale meshes such as OASIS [3], MpCci [4].

5) Coupler:

Represents the manager of the simulation. At this level, we select the coupled ScotSolvers and Scot-Transformers, and specify the orchestration strategy that these selected components will exchange variables. The parameters and functions of orchestration are given by the ScotCoupling-Tool, these parameters define the coupling strategy.

III. UML DIAGRAMS OF THE SCOT SOFTWARE

To facilitate the understanding of the Scot Composition and architecture, we have given in (Fig. 5) and (Fig. 6) the simplified UML diagrams of Scot.

IV. COUPLING OF SCOT AND MATLAB

A Socket based communication is developed to couple and adapt the Matlab models to the ScotComponent (Fig. 7), it contains 2 java classes and a ".dll" library and implements the following methods:

1) Initialization of the Matlab models (build the models according to the model parameters). 2) resolution of the models and getting the residuals to test the convergence of the coupling. 3) getResults to recover the final results when the convergence criteria is achieved.

These three methods are respectively called by the sComp Matlab methods : init(), solve(), and getOutputs(). The tight connection between Scot and Matlab is performed by Matlab Engine.

V. COUPLING STEPS DESCRIPTION WITH USING SCOT Scot has a suitable graphical user interface that allows us to easily build coupling applications by following the steps below: a) Step one: Build the Scot components and Scot coupling terms as explained briefly in section II-A2. These components appear in the work bar in the left side of the main page of Scot.

b)

Step two: Select a ScotCouplingTool, this one provides as an interface where we can find the ScotSolvers built before, we choose the components to be coupled and give the orchestration parameters. VI. APPLICATION

A. First Application : Micro-inductance

We have chosen as a first application the modeling of a micro-inductance (Fig. 8) by coupling PEEC and MoM method to respectively model the electric inductance and magnetic materials in the device. The numerical model of the micro-inductance is given by the system of equations below [START_REF] Mokhtari | Comparing Weak and Strong PEEC-MoM Coupling[END_REF].

Glue Copper

[M M ] • [M] + [BS m ] • [I] = [0] [Z m ] • [I] + [LF m ] • [M] = [V ] (1)
where, The model of the micro-inductance has been validated with the finite element approach in [START_REF] Le Duc | Delinchant Numerical experimentations on the coupling between PEEC and volume integral method[END_REF]. For the same results of magnetizations and electric current in the inductance, we have compared the time resolution of strong coupling to the tuned weak coupling proposed in [START_REF] Mokhtari | Comparing Weak and Strong PEEC-MoM Coupling[END_REF] regarding the number of unknowns (Fig. 9). The curves show that the simulation of the coupling system has found a good interest in orchestration, it has enhanced and sped up the convergence of the coupling.

B. Second Application : Magneto-Mechanic NEMS

The second test case is a model of Magneto-Mechanic nanoswitch, composed by a permanent magnet, an electromagnet and a deformable beam (Fig. 10). This model has been build and detailed in [START_REF] Pham Quang | Semi-Analytical Magneto-Mechanic Coupling With Contact Analysis for MEMS/NEMS[END_REF]. A weak coupling iterative algorithm is used to couple a numerical mechanic model with an analytical magnetic one. The computed variables are the contact forces produced by the electromagnet on the permanent magnet, the beam deformation vector caused by the displacement of the permanent magnet and the contact length with the fixed substrate. We can see that the progressive accuracy approach is enough faster than the fixed accuracy one, it has converged with a thin relative error and minimized the simulation cpu time by 38 %.

VII. CONCLUSION

We have presented Scot which is a solvers coupling and orchestration tool. The solvers are specified under a norm of Scot Components and have provided the simulation parameters to the manager. These parameters allow us to use the coupling strategies proposed in [START_REF] Mokhtari | Comparing Weak and Strong PEEC-MoM Coupling[END_REF]-[4], develop and test new efficient algorithms to optimize the multi-physic systems simulation.

The concept of layers introduced in Scot minimizes a lot exchanging data with the manager which vacates space memory and requests much less the processor. The architecture based software components used in Scot promotes the reusing of models libraries and so facilitates the modeling of multiphysic systems.

In the applications, we have used a static ScotComponent, our future research concerns the specification of dynamic ScotComponent to the co-simulation of temporal models with variable subcycling step time to accelerate and optimize this type of coupling simulation.

Fig. 2 .

 2 Fig. 2. Representation of Scot Coupling component, its parameters and interactions with Scot manager and other components.

Fig. 3 .

 3 Fig. 3. Graphical user interfaces of Scot Components (a) GUI of Scot Component calling Matlab model, and (b) GUI of Scot Component calling Java class.

Fig. 4 .

 4 Fig. 4. Graphical user interface of Scot composite component.

Fig. 5 .

 5 Fig. 5. Simplified UML diagram of the ScotSolver class and its subclasses.

Fig. 6 .

 6 Fig. 6. Simplified UML diagram of the ScotCouplingTool.

Fig. 7 .

 7 Fig. 7. Steps description of coupling with using Scot.

Fig. 8 .

 8 Fig. 8. Micro-inductance, composed by two magnetic regions (µr = 1000) and an electric inductance powered by a source voltage (V = 1volt, f = 1kHz).

[

  M M ] Magnetostatic matrix, obtained with LOCAPI [6] [BS m ] Impact of the electric on magnetic, coupling term, [Z m ] Impedance matrix, obtained with INCA3D [7] [LF m ] Impact of the magnetic on electric, coupling term. The resolution of the system (1) gives us the vector of material magnetizations [M] and the electric current in the inductance [I], according to the sources vector [V ].

Fig. 9 .

 9 Fig. 9. Curves of the time resolution regarding the number of unknowns.

Fig. 13 .

 13 Fig. 13. Residuals of the mechanic and magnetic model, relative errors, and cpu time of the fixed and progressive accuracy approach.

ICAr : Norm of the software component used in CADES environment [2]

To orchestrate the simulation of the coupled model, the available parameters are the accuracy and the number of iteration of the mechanic model, and the accuracy of the integral in the magnetic one since it is solved numerically (Fig. 11). 1) Coupling strategy: In coupling algorithms, the accuracies of the coupled solvers are unchanged during the simulation. We have proposed a progressive accuracy function based on the coupling iteration number as shown in figure 12, We have presented in this paper Scot which is a solvers coupling and orchestration tool. The solvers are specified under a norm of Scot Components and have provided the simulation parameters to the manager. These parameters allow us to use the coupling strategies proposed in [\ref{Mokhtari}]-[\ref{mpcci}], develop and test new efficient algorithms to optimize the multi-physic systems simulation.\vspace{0.1cm}

Magnetic

The concept of layers introduced in Scot minimizes a lot exchanging data with the supervisor which vacates space memory and requests much less the processor. The architecture based software components used in Scot promotes the reusing of models libraries and so facilitate the modeling of multi-physic systems.\vspace{0.1cm}

In the applications, we have used a static ScotComponent, our future research concerns the specification of dynamic ScotComponents to the co-simulation of temporal models with variable subcycling step time to accelerate and optimize this type of coupling simulation. Posteriori parameters analysis can be peformed to determine the optimal parameters of the simulation and then automatically tune the coupled model.

A comparison has been made between the fixed and progressive accuracy approach. We have plotted in (Fig. 13) the force and deformation residuals, and given the relative errors and the simulation cpu time of the two approaches.


Mohamed Lamine Lamali
email: mohamedlamine.lamali@u-bordeaux.frjohanne.cohen@lri.fr

Simon Lassourreuille

Stephan Kunne

Johanne Cohen

A stack-vector routing protocol for automatic tunneling

Keywords: Tunneling, encapsulation, path computation, routing protocol, distributed algorithms

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Routing is one of the most important tasks in any network, and particularly in the Internet. Routing protocols are generally (distributed) versions of path computation algorithms in graphs. For instance, RIP relies on the distributed Bellman-Ford algorithm, while OSPF uses Dijkstra's algorithm. These routing protocols were developed in the early years of the Internet. Thus, they work only in networks using the same protocol 1 (generally IP). However, Internet encompasses now several protocols. For instance, IPv4 and IPv6 coexist, and a lot of other situations involve interoperability of different protocols within the Internet. Thus, a path between a source and a destination may contain several portions using different protocols (for example, if the path crosses several domains or Autonomous Systems). The mapping from a protocol to another one along a path is generally done in two ways: i) Conversion: a packet of some protocol is converted into a packet of another one (e.g., NAT-TP [START_REF] Tsirtsis | RFC 2766 -Network Address Translation -Protocol Translation (NAT-PT)[END_REF]); ii) Encapsulation: a packet of some protocol is encapsulated (or nested) in a packet of another one, thus being transparent to the crossed nodes until its decapsulation (the reverse operation). All these operations are called adaptation functions further in the paper.

A tunnel is a subpath beginning at an encapsulation and ending at the corresponding decapsulation. Tunneling is a ubiquitous feature in the Internet nowadays. It is involved in virtual networks, security (an encryption can be modeled as an encapsulation), interoperability, etc. Tunnels can be nested to achieve different tasks. In this context, topology connectivity is not enough to ensure communication. Even if the network is (strongly) connected, there may not exist a suitable adaptation function sequence to reach some destination from some node. If a path has such a suitable sequence, it is said to be feasible. However, current routing protocols are unable to handle protocol heterogeneity. Nowadays, the endpoints of tunnels are manually configured or determined from a precomputed list. Several works try to automate a part of the process (e.g., TSP [START_REF] Blanchet | RFC 5572: IPv6 Tunnel Broker with the Tunnel Setup Protocol (TSP)[END_REF], ISATAP [START_REF] Blanchet | RFC 5572: IPv6 Tunnel Broker with the Tunnel Setup Protocol (TSP)[END_REF], etc.). They are referenced as automatic tunneling. However, even if the negotiation of some parameters of the tunnel establishment are automated, such as tunnel type and DNS registration, the determination of the tunnel endpoints is not fully automated. Moreover, these approaches cannot handle an arbitrary number of nested tunnels. The main reason is that the underlying algorithms of routing protocols cannot handle adaptation functions. The authors of [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF] propose a centralized polynomial algorithm to solve path computation problem in a multiple protocol context. However, their approach is based on language theory, and is unlikely to be distributed. To the best of our knowledge, no distributed algorithm for fully automated tunneling, i.e., with automatic determination of the (possibly nested) tunnel endpoints was proposed before this work.

Our goal is to design a distributed algorithm for path computation with automatic tunneling. This corresponds to the classical All-Pairs Shortest Path (APSP) problem, but in a network involving adaptation functions. To achieve this goal, we generalize the distributed Bellman-Ford algorithm, in order to take into account possible encapsulations, and to automatically establish (nested) tunnels. The distance-vector is replaced by a stack-vector. The route to follow for a packet is then determined by its destination and its protocol stack.

Our contributions:

1) We design the first distributed algorithm for routing with automatic tunneling; 2) We prove nontrivial bounds on the maximum protocol stack height and on the algorithm convergence: we show that the maximum protocol stack height (number of encapsulated protocols at the same time) of the shortest path (involving tunnels) between two nodes is at most λn2 , where n is the number of nodes and λ is the number of protocols in the network. This implies that the maximum message size is also polynomial, despite the fact that the shortest path may be of superpolynomial length in this context; 3) We design a proto-protocol for routing with automatic tunneling and we evaluate its efficiency through simulations.

The paper is organized as follows: Section II describes the problem and discusses the related work. Section III details the model used through this paper and formalizes the notion of path feasibility. Section IV describes the proposed algorithm and study its convergence and its message size. Section V describes a basic implementation of a routing protocol based on our algorithm, while Section VI presents the simulation results. Finally, Section VII concludes the paper.

II. THE PROBLEM

A. Problem definition and illustration

We aim to illustrate the path computation problem through an example. Figure 1 depicts a network with 8 nodes 2 . We want to compute a path from node S to node D. Both of them use protocol a. However, the nodes connecting S and D use different protocols. And some nodes are able to map a protocol to another one. A packet of any protocol comprises a header field and a data field. The structure of the header field is specific to a protocol. Converting a packet of protocol a to protocol b consists in rewriting the header and converting all the information (for example, the node identifiers, the routing options, etc.). Encapsulating a packet of a in b consists in considering the whole packet (header and data) of a as the data field of b, and then adding a header of protocol b at the top of the packet. For example, node U 1 is only able to encapsulate (a packet of) protocol a in protocol b, while node U 2 is only able to receive and send protocol b without any change of the possible encapsulated protocols (which are transparent to it). The capabilities (adaptation functions) of each node are listed above it. Each step of the path is associated with a protocol stack, i.e., the current sequence of encapsulated protocols, the current protocol being at the top of the stack. In Figure 1a, node S emits (a packet of) protocol a, so at this stage the protocol stack contains only a. When the packet crosses node U 1 it is encapsulated in protocol b. Thus, the current protocol stack is ab (from bottom to top). At U 2 the protocol stack remains the same, since node U 2 passively retransmits packets with current protocol b. Node U 4 then decapsulates protocol a from b, thus the current stack is again a. Node U 5 receives protocol a and it is only able to decapsulate protocol b from a; however, there is no b encapsulated in a at this stage. This path is said to be not feasible. The same problem would appear if a node receives a protocol while it is only able to handle another one. In contrast, a path is feasible if at each stage, the current node can handle the current protocol stack, and the destination node receives a packet with no encapsulated protocols. Figure 1b illustrates a feasible path that involves a loop. This loop is necessary to collect the required encapsulations to cross the subpath U 4 . . . U 6 . This example illustrates some specificities of feasible paths. There is not always a feasible path between two nodes, even if the graph is connected. Moreover, the shortest feasible path (if any) can involve loops. And optimal feasible paths do not exhibit optimal substructure (a subpath of a feasible path is not necessarily feasible, since the protocol stacks at the start and the end of the subpath may be different). In practice, the protocol stack is the sequence of nested headers of a packet that underwent multiple encapsulations.

The problem we address is to compute the shortest feasible path (if any) between each pair of nodes in a distributed way, i.e, the All Pairs Shortest paths (APSP) problem in a network involving adaptation functions. Our algorithm builds routing tables allowing a packet with some protocol stack to reach the destination through the shortest feasible path. The main motivation of this work is to propose a generic algorithm and a basic protocol for routing. Given a network with different protocols and some adaptation functions on some nodes, the goal is to compute the shortest path between any pair of nodes. Tunnels should be automatically established during the routing table construction. This can have important applications on several fields: IPv4/IPv6 interoperability, VPNs, secured tunneling, emulation of lower layer protocols over IP [START_REF] Bryant | RFC 3985: Pseudo Wire Emulation Edge-to-Edge[END_REF], etc.

Note that our problem is different from path computation in multi-layer networks where the number of layers is fixed [START_REF] Jiang | On a class of stochastic multilayer networks[END_REF], while our problem is to resolve the APSP in a network where the number of needed layers is unknown, the layers being induced by the node encapsulation and decapsulation capabilities. Moreover, in our setting, even two layers with the same stack height may not be comparable, since two stacks such as abb or aba cannot be routed in the same way, for example.

B. Related Work

1) Centralized approach: In [START_REF] Kuipers | Path selection in multi-layer networks[END_REF], The authors prove that the path feasibility problem under bandwidth constraint is NP-hard. They propose an exponential algorithm using a breadth-first search approach to solve the problem. The authors of [START_REF] Lamali | Path computation in multi-layer multi-domain networks: A language theoretic approach[END_REF] prove that the problem is polynomial3 if the bandwidth constraint is relaxed. They model the network as a Push-Down Automaton (PDA), and use tools from language theory to find the shortest feasible path between two nodes. In [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF], they widely generalize the previous algorithms, and prove that the feasible path problem under bandwidth constraint remains NPhard even with only two protocols and on symmetric directed graphs. The authors of [START_REF] Iqbal | Technology-aware multidomain multi-layer routing[END_REF] propose a matrix based model and an algorithm to compute feasible paths. However, the algorithm takes exponential time, and can only compute loopless feasible paths. Unfortunately, the language theoretic approach seems unlikely to be distributed. Some works have been devoted to distributed PDAs (see [START_REF] Csuhaj-Varjú | Parallel communicating pushdown automata systems[END_REF] for example). However, in these works, the nodes are themselves PDAs. In contrast, our model is equivalent to a PDA where the states collaborate in a distributed way to find an accepting transition sequence.

2) Distributed algorithm related work: In distributed literature, adaptation functions are generally not considered. The distributed Bellman-Ford algorithm (on which relies RIP) solves the APSP problem in O(n) rounds in a synchronized setting if the message size is unbounded (e.g., the LOCAL model [START_REF] Peleg | Distributed computing: A locality-sensitive approach[END_REF]). However, the APSP problem mostly received attention in the CONGEST model, where the message size is restricted to be logarithmic [START_REF] Abboud | Near-linear lower bounds for distributed distance computations, even in sparse networks[END_REF], [START_REF] Elkin | Distributed exact shortest paths in sublinear time[END_REF]. Unfortunately, these works cannot apply in our context since they do not take into account the adaptation functions and the interoperability issues.

3) Networking related work: Feasible path computation is a major challenge in networking. It underlies several technical problems: establishing tunnels, unifying control planes, etc. IPsec is a classical example. The endpoints of IPsec tunnels are often manually determined. The Tunnel Endpoints Discovery mechanism [START_REF] Bang-Li | A Study on IPsec Tunnel Endpoint Discovery Technology[END_REF] allows to automatically determine IPsec tunnel endpoints. However, this can be only done through a single path. This path corresponds to the shortest path without encryption. IPv4/IPv6 interoperability is mainly achieved through conversion of protocols or by tunneling. In the first case, an IPv4 (resp. IPv6) header is converted into an IPv6 (resp. IPv4) header. This operation is costly, and it requires to know the nodes able to perform the conversion. NAT-TP [START_REF] Tsirtsis | RFC 2766 -Network Address Translation -Protocol Translation (NAT-PT)[END_REF] uses this approach. Some other mechanisms (6over4 [START_REF] Carpenter | RFC 2529: Transmission of IPv6 over IPv4 Domains without Explicit Tunnels[END_REF], 6to4 [START_REF] Carpenter | RFC 3056: Connection of IPv6 Domains via IPv4 Clouds[END_REF], etc.) use tunneling. 6to4 is referenced as automatic tunneling. Any 6to4 encapsulation capable router has an anycast address [START_REF] Huitema | RFC 3068: An Anycast Prefix for 6to4 Relay Routers[END_REF]. Discovering the tunnel endpoints is done by sending a packet to the anycast address. However, the end-hosts should know that they need a tunnel. Moreover, this mechanism is specific to IPv4/IPv6 interoperability, and cannot handle nested tunnels with arbitrary protocols. While our goal is to design a generic protocol without these limitations. Some other works attempt to automate a part of the process. For example, the TSP protocol [START_REF] Blanchet | RFC 5572: IPv6 Tunnel Broker with the Tunnel Setup Protocol (TSP)[END_REF] proposes the (automatic) negotiation of the tunnel parameters (e.g., keep alive duration) between the endpoints. However, no algorithm or process is provided to select the optimal (and may be nested) tunnel endpoints. The ISATAP protocol [START_REF] Templin | RFC 5579: Transmission of IPv4 Packets over Intra-Site Automatic Tunnel Addressing Protocol (ISATAP) Interfaces[END_REF] proposes to select the endpoints from a predetermined list via a tunnel Broker. But, again, no algorithm is provided to automatically compute the tunnel endpoints. Moreover, these mechanisms are unable to handle an arbitrary number of nested tunnels. Fully automatic tunneling is still an open problem. The RFC 7059 [START_REF] Steffann | RFC 7059: A Comparison of IPv6-over-IPv4 Tunnel Mechanisms[END_REF] proposes a comparison between the different IPv4/IPv6 tunneling mechanisms.

C. Our approach

We propose to generalize the distributed Bellman-Ford algorithm. Since the path to follow by a packet depends upon the destination, but also upon the packet protocol stack, the distance-vector should be replaced by a stack-vector. Each node advises its neighbors that it can reach some destination, with some protocol stack, and at some cost. The rows of the routing table of each node should contain a next hop per destination and per protocol stack. The termination of such an algorithm is an important issue. Since feasible paths may involve loops, the same destination can be reached with different protocol stack heights. The question is to know at what stack height to stop. We resolve the problem by showing a bound on the stack height of any shortest feasible path. The convergence speed of the algorithm depends upon the length of the shortest feasible paths. Our algorithm is generic: it can handle any number of protocols, and it does not limit the number of nested tunnels to find a feasible path.

III. MODEL AND DEFINITIONS

This section describes the mathematical model used in this paper and formalizes the notion of path feasibility. We use the same notations and definitions as in [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF].

A. Network model

A network is modeled as a 4-tuple N = (G, A, F, w). The network topology is represented by a symmetric directed graph G = (V, E). Each node in V corresponds to a router, and each pair (U, V ) ∈ E is a unidirectional link. The number of nodes is denoted by |V| = n, and the number of links is denoted by |E| = m. The set of protocols available in the network is denoted by A = {a, b, . . . }. Its cardinality is denoted by |A| = λ. Each node U ∈ V has a set F(U ) of adaptation functions that it can perform. These functions are:

• Conversion: A protocol a is converted into a protocol
b without any change of the possible encapsulated protocols. It is denoted by (a → b) (example: IPv4/IPv6 conversion through NAT-PT). Observe that a classical retransmission without protocol change is a special case of conversion. It is denoted by (a → a).

• Encapsulation: A protocol a is encapsulated in b. It is denoted by (a → ab) (e.g., IPv4/IPv6 encapsulation).

• Decapsulation: A protocol a is decapsulated from b.

It is denoted by (a → ab).

The set of all adaptation functions available in the network is denoted by F. We denote the set of protocols that a node can receive by In(U ), and the set of protocols that a node can send as Out(U ).

Finally, performing an adaptation function on a node U has a cost defined by the function w : V ×F ×V → + . The value w(U, f, V ) (where U, V ∈ V and f ∈ F(U )) corresponds to the cost of using link (U, V ) with adaptation function f on node U . Hence, function w represents any additive metric associated with both links and adaptation functions. So the cost of a path is the sum of the costs of each triple (U, f, V ) involved in it in the network. For example, if one wants to minimize the number of encapsulations in the path, the cost function should be w(U, f, V ) = 1 if f is an encapsulation and 0 otherwise.

B. The protocol stack

A sequence of adaptation functions induces a protocol stack. For example, the sequence (a → a)(a → ab)(b → b) induces the stack ab (from bottom to top). For each position i in a path, H i denotes the protocol stack at this position, i.e., the protocol stack induced by f 0 . . . f i , and h i denotes the protocol stack height. The protocol at the top of a stack H is denoted by T op(H) and the protocol just below T op(H) in the stack is denoted by T op -1 (H). The height of a stack H is denoted by h(H). The "forbidden" stack is denoted by ∅ (note that it should not be confused with the empty word ).

More formally, let f be an adaptation function, and let H be a stack and H = H .T op(H) 4 , where H is eventually empty. We will also denote by f the function taking as argument a stack and performing the adaptation function on this stack:

• if f = (x → y) and T op(H) = x, then f (H) = H .y • if f = (x → xy) and T op(H) = x, then f (H) = H.y • if f = (x → xy) and T op(H) = y and T op -1 (H) = x, then f (H) = H • f (H) = ∅ otherwise.
We also denote by f the reverse function of f , i.e., if H = f (H) and H = ∅ then f (H ) = H. Note that if f is an encapsulation, then f is the corresponding decapsulation. f (H) = ∅ means that the adaptation function f cannot handle the stack H (e.g., f encapsulates a in b while T op(H) = a). Note that f (∅) = ∅ for any function in our context.

Thus, the protocol stack H i induced by a sequence of adaptation functions f 0 . . . f i , is recursively defined as following:

• H 0 = x if f 0 = (x → x) and H i = f i (H i-1 )

C. Path feasibility

In our context, a path should take into account the adaptation function capabilities of the network nodes. A path should contain the list of adaptation functions involved in it. Thus, a path from node S to node D in a network N is a sequence of nodes and adaptation functions Sf 0 U 1 f 1 U 2 f 2 . . . U f D where each U i , i = 1, . . . , , is a node, and each f i is an adaptation function 5 . When needed, a directed path in the graph without taking into account the adaptation functions (i.e., just a sequence of nodes) is referred to as a classical path. Definition 1:

A path P = Sf 0 U 1 f 1 U 2 f 2 . . . U f D is feasible if and only if: 1) The sequence SU 1 U 2 . . . U D is a classical path in G = (V, E) and each f i ∈ F(U i ); 2) H = x and x ∈ In(D).
Actually, the protocol sequences of feasible paths can be characterized as a context-free language [START_REF] Lamali | Path computation in multi-layer multi-domain networks: A language theoretic approach[END_REF]. In Figure 1a, the depicted path is:

S(a → a)U 1 (a → ab)U 2 (b → b)U 4 (a → ab)U 5 (b → ba)U 6 f D.
It cannot be feasible for any adaptation function f ∈ F(U 6 ) since (a → ab) appears before any encapsulation of protocol b. In contrast, Figure 1b depicts the feasible path:

S(a → a)U 1 (a → ab)U 2 (b → b)U 3 (b → ba)U 1 (a → ab) U 2 (b → b)U 4 (a → ab)U 5 (b → ba)U 6 (a → ab)D.
The corresponding protocol stacks are below the links in Figure 1.

D. Problem formalization

The cost of a path

P = Sf 0 U 1 f 1 U 2 f 2 . . . U f D from node S to node D is defined as w(P) def = i=0 w(U i , f i , U i+1
) where S = U 0 and D = U +1 . An optimal (i.e., shortest) feasible path P between two nodes is a feasible path that minimizes w(P). When needed, we refer to the number of hops in a path P by |P|. Note that |P| may be different from w(P).

The problem we aim to solve is the following: compute the routing tables at each node, such that for each pair of nodes U and V , a packet from U to V following these routing tables follows an optimal path from U to V . This corresponds to an APSP in a network involving adaptation functions. Note that there is not always a feasible path between two nodes even if the underlying graph G is strongly connected.

IV. STACK-VECTOR ROUTING ALGORITHM

A. Different approaches to distributed routing

There are two main approaches in distributed routing: link-state routing and distance/path vector routing. The first one consists in spreading the whole topology then computing locally the path in a centralized way. However, this approach has a drawback: computing the whole shortest path at a node is not enough to route the packets. A node should know the next neighbor at each stage of the route. This is not possible unless the header of the packet contains the whole path to follow. The second approach is to use a Bellman-Ford algorithm (as RIP) where each node shares its current routing table with its neighbors. We cannot use this approach since the path leading to the destination depends upon the protocol stack of the packet to route. Storing the cost and the destination in the routing table is not enough.

Path vector routing consists in sharing not only the destination and the cost, but also the whole path to reach the destination. The main goal to do so is to avoid loops (if a node receives a path where it already appears, then it discards it). BGP uses this approach. This would be a possible solution for our problem, since a complete feasible path (the node sequence and the adaptation function sequence) characterizes the route to take for a packet. However, the shortest feasible path can be superpolynomial [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF]. A lower bound of Ω((n/λ) λ ) is known for an arbitrary number of protocols λ. Even for two protocols, there are shortest paths of length at least Ω(

√ n2 Ω( √ n)
). The best available upper bound is 2 O(λ 2 n 2 ) -1 (see [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF] for more details). Using a path vector protocol would lead to a superpolynomial message length. We opted for a stack vector algorithm (by analogy with distance vector protocols), where the destination, the cost and the protocol stack are stored in the routing table.

B. The algorithm

The initialization algorithm (Algorithm 1) allows each node U to share with all its neighbors N ei(U ) the set of protocols In(U ) that it can receive. Sending a message (U, H, c) means that the sender can reach destination U at cost c if it receives a packet with protocol stack H. Thus, the initialization phase consists of informing its neighbors that it can receive any packet with protocol x ∈ In(U ) (without nested encapsulations) at cost 0. for all V ∈ N ei(U ) do Send(U, H, 0) to V Each row of the routing table T of each node is a 5-tuple (D, H, c, V, f ) where D is the destination to reach, H is the protocol stack needed to reach D, c is the cost of the remaining path to reach the destination, node V is the next neighbor (next hop) to reach the destination, and f is the adaptation function to perform at the current node. The row is indexed by the pair (D, H). Thus, T (D, H) returns the row corresponding to destination D, and stack H. When a node tries to add a new row (D, H, c, V, f ) to its routing table T , it checks if the tuple (D, H) is already in its routing table. If it is not, then it inserts the new row to the table. Otherwise, it compares the cost of the new route with the old one. The table is updated by replacing the old route by the new one if the new cost is lower. This step is done according to Algorithm 2.

Algorithm 2 Add a row to a routing table T of node U

Input: A row (D, H, c, V, f ) 1: if (D, H) / ∈ T then 2:
Add (D, H, V, c, f ) to T 3: else if T (D, H).cost > c then

T (D, H).f unction ← f

Algorithm 3 is the main algorithm. When receiving a new message (D, H, c) from a neighbor V , node U determines which of its adaptation functions it can apply to stack H (lines 3-5). Then it computes the new cost of the remaining path by adding the cost of (U, f, V ) to the old cost (line 6). It tries to add the new row corresponding to the message and the chosen adaptation function to its routing table T , according to Algorithm 2 (line 7). If the row corresponds to a new route in the table, then it sends the new message (D, H, c) to its neighbors, indicating that it can reach the destination with the new stack H. Observe that if f is a decapsulation (a → ab), and U can reach D with some protocol stack H such that T op(H) = a, then U can handle a packet with stack H.b to reach D. Similarly, If f is a conversion (a → b), and if U can reach D with some protocol stack H such that T op(H) = b (i.e., H = H .b for some stack H ), then U can handle a packet with stack H a to reach D. Thus, if node U can reach a destination D with some stack H, for each f ∈ F(U ), it should apply the reverse adaptation function f (line 4) to H before sending a message to the other nodes. Note that only stacks with heights less than or equal to λn 2 are kept and shared (line 5). In the next section, we prove that this height is sufficient to reach any destination with a shortest feasible path if there exists one.

Algorithm 3 Routing table construction algorithm of node U

1: loop 2:

Receive (D, H, c) from V

3:

for all f ∈ F(U ) do 4:

H ← f (H) 5: if H = ∅ and h(H) ≤ λn 2 then 6: c ← c + w(U, f, V ) 7:
Add row (D, H, c, V, f ) to T

8:

# according to Algorithm 2

9:

if T has been modified then 10:

for all W ∈ N ei(U ) do 11:

Send (D, H, c) to W

We assume that the packets to route are of the form (destination, protocol stack, payload); for example, (D, H, data). Once all routing tables are computed, if node U receives a packet with protocol stack H and destination D, it first searches for tuple (D, H) in its routing table. If there is no corresponding row, then node U has no route for the destination with the received protocol stack. Otherwise, if the corresponding row is (D, H, c, V, f ), then it sends the packet (D, f (H), data) to V . The routing procedure is illustrated by Algorithm 4.

Algorithm 4 Routing a packet. Send the packet (D, f (H), data) to V .

C. Correctness and complexity

1) Stack height upper bound and message size: First, we prove that for any shortest feasible path, the maximum stack height reached along the path is polynomially bounded.

Lemma 1: The shortest feasible path (if any) between two nodes reaches a maximum stack height of at most λn 2 protocols.

Proof: See Appendix A.

Corollary 1:

The maximum message size is in O(λ log λ n 2 ).

Proof: A message contains node identifiers, that are in O(log n) if the number of nodes is n. A protocol stack, its height is at most λn 2 according to Lemma 1. The identifier of a protocol is in O(log λ) for a number λ of protocols. We assume that the cost is bounded.

2) Convergence and correctness: We will prove that in a synchronized setting (for example the LOCAL model in distributed computing [START_REF] Peleg | Distributed computing: A locality-sensitive approach[END_REF]) where all nodes receive messages at the same time t, then process them and send them at the same time t+1, Algorithm 3 correctly converges in polynomial time 6 according to the network size and its diameter.

We define the diameter of a network N , denoted by diam N , as the length (in number of hops) of the shortest feasible path that maximizes the number of hops:

diam N def = max P shortest feasible |P| .
Proposition 1: Algorithm 3 computes the correct routing tables after O(λn 2 diam N ) rounds.

Proof: See Appendix B.

V. A PROTO-PROTOCOL

A. Practical limitations and algorithm adaptation

The main issue for the implementation of the proposed algorithm is that the stack height and the message length bounds are too large, even if they are polynomial. For example, in a network of 100 nodes with only 2 protocols, the stack height may reach 2 × 10 4 protocols. This is not sustainable for real applications. Moreover, for a stack height h, there are λ h possible stacks, i.e., in the worst case, there may be n 1-λ λn 2 +1 an exponential number of exchanged messages. Finally, as demonstrated by the authors of [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF], diam N can be superpolynomial in the size of the network N . These bounds are tight, i.e., it is possible to exhibit a network where the shortest feasible path reaches a stack height of λn 2 . Thus, these limitations are not due to our algorithm but are inherent to the problem.

However, the simulations performed in [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF] show that such instances are extremely unlikely to appear. Thus, we propose to set the maximum stack height as a parameter of the protocol. A small value of this parameter is enough to compute the shortest path in most cases. We should bound the maximum stack height by a constant h max for any real implementation. Thus, the condition h(H) ≤ λn 2 in Algorithm 3 (line 5) should be replaced by h(H) ≤ h max . In such a case, the maximum number of different stacks of height h max would be λ hmax , and the maximum number of rows in a routing table would be n 1-λ hmax+1 1-λ -1 , which is linear in n and polynomial in λ.

B. Routing proto-protocol specification

The main requirements are that the nodes and the protocols must have a unique identifier. For example, 1 byte for the protocol identifiers and 16 bytes for the node identifiers. The routers involved in the routing protocol should have a specific multicast address. The routing protocol messages that advise a route must contain: the destination of the route, the protocol stack needed for a packet in order to reach the destination, and the cost of the route. In addition, some classical parameters can be exchanged: keep alive duration, emitted message timestamps, etc. The maximum stack height should be set at the start of the process. The routing table must contain entries indexed by the destination and the protocol stack of a received packet. Each entry must contain the cost, the next hop, and the adaptation function to perform before sending the packet to the next hop.

C. Routing a packet

The packet to route contains a stack of headers of different protocols. It can be seen as a meta-header. It must contain:

• The unique identifier of the destination node;

• The source (even if it does not impact the routing process, since it is a per destination/stack routing);

• The stack height: the current stack height of the packet

• The protocol stack: the stack of identifiers of the protocols corresponding to the nested headers. This will speed up the routing process, since it avoids to access to the whole header stack;

• The header stack: each encapsulated header preceded by its length;

• The payload of the inner packet.

Figure 2 illustrates such a meta-header.

VI. SIMULATIONS

In order to evaluate the efficiency of our algorithm, we performed simulations with different parameters.

A. Simulation methodology

All the networks used in the simulations are generated according to the following steps:

1) We generate a random undirected graph of a given size according to a preferential attachment mechanism (the Barabàsi-Albert model [START_REF] Barabási | Emergence of scaling in random networks[END_REF]), where each added node is attached to 5 existing nodes; 2) The graph is then converted into a symmetric directed graph. Each undirected link is converted into two directed links; 3) For a given number λ of protocols, there are 3λ 2 possible adaptation functions. Each adaptation function is available on a node with a given probability p.

The algorithm is implemented in Python 3.4.5, using the NetworkX package 7 . The implementation is done in a distributed fashion: each node is simulated by a thread, and a directed link (U, V ) is implemented as a queue where U can only write, and V can only read. The simulations were performed on a multi-core server with CPUs 1.59GHz.

The input parameters of the algorithm are: the number of nodes, the probability p of availability of an adaptation function, the number of protocols, and the maximum stack height h max . The main output results are the convergence time, and the percentage of times where the algorithm finds the shortest path between the network extremities. Note that, if the algorithm does not find the shortest feasible path, it may be because there is no feasible path in the given network. The probability of existence of a feasible path according to different parameters can be found in [START_REF] Lamali | Algorithmic and complexity aspects of path computation in multi-layer networks[END_REF]. All the results are averaged over 1000 runs.

B. Convergence time

Figure 3 shows the convergence time of the algorithm according to different parameters. Figure 3a shows the convergence time according to the network size (number of nodes) in three cases: p = 0.05, p = 0.10, and p = 0.20. The maximum stack height h max is set to 3, and the number of protocols is set to 2. It appears that the convergence time hugely depends on the parameter p. For p = 0.05, the convergence time is 2.74s in a network of 50 nodes; while for p = 0.20, the convergence time is 113s. However, note that if p = 0.20, the average number of adaptation functions per node is 2.4, which is unrealistic, since only a few number of nodes should be able to perform conversions/encapsulations. Actually, Figure 3b shows the impact of the probability p (and the average number of adaptation functions per node) on the convergence time (with h max = 3 and 2 protocols). For p = 0.05, the convergence time is around 0.7s (resp. around 3s) in a network of 30 (resp. 50) nodes; while for p = 0.30, the convergence time is around 50s (resp. around 3min) in a network of 30 (resp. 50) nodes. We can see that the processing time hugely depends upon the number of adaptation functions per node. Figure 3c shows the impact of the parameter h max on the convergence time. The number of protocols is set to 2 and the network has 50 nodes. Note the logarithmic scale on the y-axis. The maximum stack height hugely impacts the convergence time. With p = 0.05 (resp. p = 0.20), the convergence time is 0.7s (resp. around 30s) if the maximum stack height is set to 2. However, the convergence time is around 21s (resp. 11min) if the maximum stack height is set to 5.

C. Algorithm efficiency

Figure 4 shows the average number of times where Algorithm 3 found a shortest feasible path between the network extremities, according to different parameters. For each run, we check if a shortest feasible path is computed between the network extremities. Note that, if there is no path, it can be due to two different reasons: i) there is no feasible path between these nodes; ii) the shortest feasible path reaches a stack height larger than h max .

Figure 4a shows the impact of probability p on the number of found paths. For 30 nodes (resp. 50 nodes) and h max set to 3, the number of runs where a shortest feasible path was found is 1% (resp. 3%) if p = 0.05. While it is 80% (resp. 73%) if p = 0.30. Note that the probability of existence of a feasible path is much smaller if p = 0.05 than if p = 0.30. Figure 4b shows the impact of the maximum stack height h max on the number of found paths. The probability p is set to 0.20 and there are 2 protocols. Obviously, the algorithm finds more feasible paths if it is allowed to explore paths with higher stacks. However, the difference between stack heights 3 and 5, for example, is negligible. It seems that almost all the shortest feasible paths are of maximum stack height at most 3 (i.e., two nested tunnels).

These simulations shows that the problem is complex, and the convergence time can be prohibitive, especially if the number of adaptation functions per node is large. However, a convergence time of few minutes is sustainable in a network with infrequent topology changes.

VII. CONCLUSION

Nowadays, the Internet encompasses several protocols. The interoperability between theses protocols is an important issue, and is ensured thanks to adaptation functions. A path between two nodes may involve different protocols at different stages. The current routing protocols are not able to automatically compute such paths due to the lack of distributed algorithms taking into account the adaptation functions. In this paper, we design the first fully distributed algorithm taking into Fig. 4: The % of found feasible paths by Algorithm 3 according to different parameters. account these functions. Our algorithm builds at each node a local routing table that allows to route a packet following the optimal path. Moreover, we prove that our algorithm converges polynomially in the size of the network and its diameter, and it uses messages of polynomial size, despite the fact that the shortest feasible paths can be of superpolynomial length. We propose a basic implementation of our algorithm as a stackvector routing protocol, and we evaluate its efficiency through simulations. We believe that this work can have an important impact on protocols such as TSP, on IPv4/IPv6 interoperability, and on automatic tunneling more generally. As a future work, we plan to study the application of our algorithm to secured communication, more precisely to nested encrypted tunnels. In this paper, we assumed that each node had access to the full protocol stack of a packet that it receives. This implies that the data can be encrypted several times but that the headers must not be encrypted. This is not suitable for security reasons. Thus, the main possible improvement to our algorithm would be to adapt it to the case in which all the encapsulated headers are also encrypted. This means that any node has only access to the top (outer) protocol. This would have applications to distributed onion routing.

guages by means of right-congruences," Theor. Comput. Sci., vol. 70, no. 2, pp. 213-232, 1990.

[21] A. Amarilli and M. Jeanmougin, "A proof of the pumping lemma for context-free languages through pushdown automata," arXiv preprint arXiv:1207.2819, 2012.

APPENDIX A

Proof of Lemma 1

Proof: We will prove that if there is a feasible path from S to D that reaches a maximum stack height h max > λn 2 , then there is a shorter feasible path from S to D that reaches a maximum stack height h max < h max . The main ideas of the proof are results from language theory, and are related to the pumping lemma. They can be found in [START_REF] Sénizergues | A characterisation of deterministic context-free lan[END_REF] and [21] for example.

Let P be a feasible path of length between two nodes. Suppose that P reaches a maximum stack height h max > λn 2 at some position j (i.e., h j = h max ). For each stack height h such that h max -λn 2 ≤ h ≤ h max , let i h (resp. k h ) be the last position before (resp. the first position after) j reaching the stack height h. More formally:

• i h = max{i ≤ j | h i = h} • k h = min{k ≥ j | h k = h}
And let h be the 3-tuple (U i , U k , a i ) where:

• U i (resp. U k ) is the current node at position i h (resp. k h ),
• a i is the current protocol at position i h (i.e., a i = T op(H i )). Note that T op(H i ) = T op(H k ) since f i is followed by an encapsulation, and f k follows a decapsulation.

For any h, there are λn 2 possible values of h . However, since h max -λn 2 ≤ h ≤ h max , h can take λn 2 + 1 different values. Thus, there are two stack heights h = h such that h = h . Let h < h for convenience. By definition of i h and k h , the subpath from position i h to k h does not involve a protocol stack smaller than h. Thus, it depends only on the part of the stack above h, on the top of the stack, and on the current node at position i h (recall that i h is the last position before j reaching stack height h, the following positions until k h involve stacks higher than h). Likewise, the subpath between positions i h and k h does not involve a stack smaller than h . The part of the stack below height h (resp. h ) is transparent to the nodes between i h and k h (resp. i h and k h ). The same sequence of adaptation functions and nodes leading from i h to k h can also lead from i h to k h . It is then possible to shortcut the part of the path between i h and i h together with the part from k h to k h . Thus, the path:

P = Sf 0 U 1 f 1 . . . U i h -1 f i h -1 U i h f i h . . . U j f j . . . U k h f k h U k h +1 f k h +1 . . . D
is feasible, shorter than P, and reaches a maximum stack height h max -(h -h ).

APPENDIX B

Proof of Proposition 1

First we prove the following lemma:

Lemma 2: Let N be a network. If after a number of rounds t, all the shortest feasible paths reaching a maximum stack height at most h max are already computed, then all the shortest feasible paths reaching a maximum stack height at most h max + 1 will be computed after t + O(diam N ) rounds.

Proof: Computing a shortest feasible path of maximum stack height h max +1 from shortest feasible paths of maximum stack height h max can be done in two steps:

1) The new path is a path of maximum stack height h max following an encapsulation and followed by a decapsulation; 2) Concatenation of several new paths of maximum stack height at most h max + 1.

Suppose that after some rounds, all the shortest feasible paths that reach a maximum stack height h ≤ h max are computed. Let P = Sf 0 U 1 f 1 U 2 f 2 . . . U f D be a shortest feasible path from node S to D that reaches a maximum stack height h j at some position j, such that h j > h max . Let i h (resp. k h ) be the last position before (resp. the first position after) j that reaches stack height h = h j -h max , and follows an encapsulation (resp. is followed by a decapsulation). More formally:

• i h = max{i < j | h i = h j -h max and h i-1 = h i -1} • i k = min{k > j | h k = h j -h max and h k+1 = h k -1}
Then the path U i h f i h . . . U j f j . . . f k h -1 U k h is feasible and is already computed, i.e., node U i h knows that it can reach U k h by emitting a packet of some protocol x, and (U k h , x, c, U i h +1 , f ) is in its routing table for some cost c and some protocol x.

On the other hand, node U k h already knows that it can reach node U k h +1 with some protocol x, since it received the message (U k h +1 , x, 0) from U k h +1 at the first round, thanks to Algorithm 1. Thus, at round t, node U i h +1 knows that it can reach node U k h +1 with some protocol stack. After two rounds, through node U i h , node U i h -1 will know that it can reach node U k h +1 with some protocol x. Now it remains to prove that if P 1 , . . . , P d are shortest feasible paths already computed, and the last node of P i (call it D i ) is the first node of P i+1 , then the shortest feasible path P = P 1 . . . P d (if any) is computed after O(diam N ) rounds. It means that the first node of P 1 knows that it can reach the last node of P d by emitting a suitable protocol. By the same argument, the first node of P 1 will know that it can reach the last node of Proof: We prove by induction that after O(λn 2 diam N ) rounds, all the shortest feasible paths that reach a maximum stack height h ≤ h max are computed, i.e., the routing tables are able to route packets following these paths.

Basis: Suppose that there is a feasible path P = Sf 0 U 1 f 1 U 2 f 2 . . . U f D from S = U 0 to D = U +1 that keeps stack height of 1 (i.e., there is neither encapsulations nor decapsulations along the path). Since, by definition, |P| ≤ diam N , after at most diam N rounds, the routing table of S contains the row (D, x, c, U 1 , w(P)) for 0 ≤ i ≤ and some protocol x ∈ In(S), since P is feasible.

Induction: Suppose that after t rounds, all the shortest feasible paths that reach a maximum stack height at most h max are computed. By Lemma 2, all the shortest feasible paths that reach a maximum stack height of h max + 1 are computed after t + O(diam N ) rounds. By Lemma 1, the maximum stack height of a feasible shortest path (if any) is λn 2 . This concludes the proof.

  (a) Unfeasible path without loop. (b) Shortest feasible path with loop.

Fig. 1 :

 1 Fig. 1: A network requiring a shortest path with loop.

Algorithm 1

 1 Initialization algorithm of node U 1: for all x ∈ In(U ) do 2:

4 :T

 4 (D, H).cost ← c 5:T (D, H).next hop ← V 6:

1 :

 1 Receive a packet (D, H, data) 2: if (D,

Fig. 2 :

 2 Fig. 2: A header of a packet to route.

  Convergence time according to hmax. Note the logarithmic scale on the y-axis.

Fig. 3 :

 3 Fig. 3: The convergence time of Algorithm 3 according to different parameters.

Proposition 1 :

 1 P d after at most |P 1 , . . . , P d | rounds. Since, by definition, |P 1 , . . . , P d | ≤ diam N , this concludes the proof. Recall now Proposition 1: Algorithm 3 computes the correct routing tables after O(λn 2 diam N ) rounds.

This paper will appear in the proceedings of IEEE INFOCOM 2019.[START_REF] Tsirtsis | RFC 2766 -Network Address Translation -Protocol Translation (NAT-PT)[END_REF] In this paper, communication protocols and routing protocols should not be confused. For simplicity, we refer to the first ones as protocols, and to the second ones explicitly as routing protocols.

Note that, in this example, the underlying topology is not a symmetric directed graph, while in our model (see Section III) the directed graph must be symmetric. In this context, a symmetric directed graph is defined as a graph where a link (U, V ) exists if and only if the reverse link (V, U ) exists. However, the example on Figure1is simplified for readability purpose. It can be easily converted into a symmetric directed graph with the same properties.

More precisely, the decision version, i.e., deciding if a feasible path exists between two nodes, is polynomial. However, as stated in Section IV-A, explicitly computing the path can be superpolynomial.

The notation "." stands for a simple concatenation. For example, if H = abab then H.b = ababb.

The adaptation function f 0 is dummy. If node S emits packets of protocol a, then f 0 is denoted by (a → a) by convention.

https://networkx.github.io/

Acknowledgement. The authors would like to thank Géraud Sénizergues for his invaluable help and explanations about the proof of Lemma 1. The first author was partially supported by the H ÉRA project, funded by The French National Research Agency. Grant no.: ANR-18-CE25-0002.


Christophe Rubeck

Bertrand Bannwarth

Olivier Chadebec

Benoît Delinchant

Jean-Paul Yonnet

Jean-Louis Coulomb

JCuda vectorized and parallelized computation strategy for solving integral equations in electromagnetism on a standard personal computer

The paper presents a computation strategy for solving integral equations in electromagnetism. Nowadays, powerful programmable Graphic Processing Units (GPU) can be found in any standard computer. The paper investigates the benefits of the use of GPUs in addition to the CPU one in order to improve computation speed by using integral methods. Java language and the JCuda library, not often used in speed calculation by the computing community, has been used here. A 100 time speed-up is reported in matrix assembly between an optimized traditional CPU computation and a CPU+GPU one. FMM… Index Terms-Fast Multipole Method on GPU, JCuda computing, Pure Java ...

I. INTRODUCTION

NTEGRAL equation methods (IEM) are currently widely used in electromagnetic modeling. Unlike the finite element method (FEM), they do not require the meshing of non-active materials like air. However, they are based on the computation of electromagnetic interactions between all elements (i.e. full interaction). Therefore, they lead to fully dense systems of equations. They are well known to be easily parallelizable because of the independence of interactions. Moreover, the interest of IEM has considerably increased since the emergence of acceleration methods such as the Fast Multipole Method (FMM). In this kind of algorithm, near and far field interactions are separated [START_REF] Greengard | A fast algorithm for particle simulations[END_REF]:

far near V V V   (1)
While far fields computations are highly accelerated by FMM, the near field interaction is treated classically. In particular, full near field matrices have to be computed with a high performance strategy in order to keep the advantage of using FMM.

In this work, a parallelized and vectorized full matrix interaction computation strategy is implemented, first using a classical CPU on a standard personal computer then a CUDA capable GPU [START_REF] Nvidia | NVIDIA CUDA programming guide[END_REF]. The main software is developed thanks to Java language so the use of the JCuda library enables GPU interfacing directly from Java [START_REF] Hutter | JCuda (Java bindings for CUDA[END_REF]. The choice of using Java can seem to be surprising for speed computations, but the use of this language enables easy portability, robust and fast software developments, and performances are respectable in comparison with most commonly used language like C++ [START_REF] Reinauer | Object-Oriented Development and Runtime Investigation of 3-D electrostatic FEM problems in Pure Java[END_REF].

II. CHARGE DENSITY COMPUTATION

Let consider a perfect conductor be in free space associated to a known potential V0. To compute the charge density in electrostatics, the following integral equation has to be solved:

dS r V S     0 0 4 1 ( 2
)
Where S is the surface of the conductor, σ is the charge density, ε0 is the vacuum permittivity and r is the distance between the point where the potential is expressed (on the conductor) and the integration point.

III. NEAR FIELD COMPUTATION

The near potential (1) is computed as presented in [START_REF] Nvidia | NVIDIA CUDA programming guide[END_REF]. Let mesh the surface into a set of triangle patches. The system of linear equations is generated using a point matching approach with 0-order shape function. This method is very simple but has already shown its accuracy.

A. High performance matrix assembly

In (2), if S is meshed into N cells, N integrals on N cells have to be computed. This is why the computation time increases in N 2 . Moreover, the numerical evaluation of the integrals of ( 2) is sometime numerically singular (in particular for the computation of the interaction of an element on itself so when r is close to zero). The use of analytical formulae [START_REF] Lezar | GPU acceleration of method of moments matrix assembly using Rao-Wilton-Glisson basis functions[END_REF] to evaluate the kernel of ( 2) is then preferred but these computations can be time-consuming. The chosen approach is a mix of numerical and analytical integral computation in order to get the best ratio between accuracy and speed.

The integrals are first computed thanks to a numerical Gauss integration technique. In this approach, there are three overlapped loops in the algorithm (fig. 1).

I

Manuscript received January 1, 2008 (date on which paper was submitted for review). Corresponding author: F. A. Author (e-mail: f.author@nist.gov).

Digital Object Identifier inserted by IEEE // Loop 1: on all the N elements For i = 0,1,…N // can be parallelized if multi-CPU // Loop 2: interaction of element i with all elements For j = 0,1,…N // Loop 3: Gauss integration For k = 0,…number of Gauss points Integral(i,j) += 1/r(i,j,k) *weight(k) * jacobian(k) End End End Working on larger and continuous sets of data improves memory access speed and therefore the computation speed. That is why the two last loops are switched in order to have a high number of indexes in the last loop. It is allowed because all interactions are independent. Furthermore vectorized operators can be used as shown in the fig. 2. The benefit of the vectorized algorithm toward the classical one is a speed-up of almost 30 to 50 times for the matrix assembly.

B. Fixing singularities

After the matrix assembly, the artificial singularities which have been introduced with the numerical integration are fixed. The diagonal coefficients of the matrix are corrected by the analytical solution of the corresponding integrals.

C. Pure Java computation

An optimized vectorized Java matrix package has been developed in the G2Elab [START_REF] Coulomb | Numerical Design Of Experiments and Optimization[END_REF]. It is based on contiguous memory storage of the matrix, adapted indexes and macro matrix manipulation commands.

Before assembling the matrix, a pre-processing is needed. A table containing all the coordinates of elements Gauss points in the main referential is generated. The process is repeated for the Gauss weights, the Gauss jacobians and the matching points. Then the matrix is assembled line per line thanks to vectorized operators and the diagonal is fixed.

D. JCuda implementation

Thanks to the JCuda library, it is possible to call CUDA kernels directly from Java. A Java GPU matrix library has been developed. It enables the management of the GPU memory allocation, the data transfer between the GPU device and the host, matrix manipulations and the call of kernels.

1) GPU computing accuracy It is well known that GPU are much faster in single precision computing than in double precision. Furthermore, only latest GPU are double precision capable. Therefore the good computation accuracy of (2) in single precision must be checked before developing a CUDA method. Computations of (2) on a triangle in single and double precision are compared for analytical and numerical method. Errors are noticed in the order of magnitude of 1e-7. So the quality of the evaluation of (2) is only a few influenced by the single precision computation. For the same calculations, results between CPU and GPU are compared. A difference on the 1-3 last digits is noticed. The explanation is that in CUDA most functions are implemented in non-standard-compliant way [START_REF] Nvidia | NVIDIA CUDA programming guide[END_REF].

2) JCuda matrix assembly Performances in CUDA programming are better if the algorithm is massively parallelized, therefore the matrix has to be computed with a high number of threads. The chosen approach allocates one thread to each interaction. So, N 2 threads are defined. Each thread contains only the Gauss integration loop. A simplified version of the matrix assembly CUDA kernel is presented in fig. 3.

Only the mesh is sent to the GPU. Gauss points, Gauss jacobians, Gauss weights and matching points tables are generated on the GPU thanks to CUDA kernels. The tables stay in the GPU global memory. The access to the data from the assembling CUDA kernel is done continuously to preserve high speed computation (i and j indexes are inverted). Finally the diagonal of the matrix is fixed with analytical solutions.

3) JCuda linear system solving Once the set of linear equations is obtained and stored in the global memory of GPU, the full interaction problem is solved by a GMRES iterative method [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. The iterative method developed here is not a fully GPU algorithm because it is hardly parallelizable, only the matrix product for the computation of the residual is done using the JCublas library.

Single precision

Double precision

Fig. 4. Relative error of the charge distribution depending of integration technique and variables precision (lower is better).

Therefore at each iteration a vector is transmitted between the host computer and the CUDA device. A simple Jacobi preconditioning method is used.

4) Hardware adaptability

The main limiting factor of GPU computing is the graphic memory amount. That is why two strategies have been developed. First if the GPU can contain the matrix then the problem is solved with the iterative method presented before. If there is no memory enough, then the matrix is assembled by blocks and transferred to the host computer. The system is solved by a classical CPU method.

E. Results

1) Benchmark problem In our example, a spherical iso-potential conductor (radius of 10 cm) is modeled. The potential is set to 1 mV. This example is simple but the theoretical charge density is known, so the accuracy of computing can be precisely evaluated.

2

) Hardware specifications

The CPU results are obtained with an Intel Xeon 2.67 GHz in single core. The CUDA devices are a Geforce 320M (48 cores, 0.95 GHz, 250 MB of shared memory) and a Tesla C1060 (240 cores, 1.30 GHz, 4 GB of RAM).

3

) Integration method accuracy

The sphere is meshed into 9068 triangular patches. The GMRES convergence break is set to 1e-9 and the Krylov subspace is extended to 30. Therefore the errors are mainly due to the meshing and the integration technique. The fig. 4 shows the relative error of the charge distribution between the theoretical values and the computations, depending of integration technique and variable precision. As expected, a double precision computation is better than a single one, and a pure analytical computation [START_REF] Rao | A simple numerical solution procedure for statics problems involving arbitrary-shaped surfaces[END_REF] is better than a numerical one. However, the differences are really poor (0.004%) and it is reasonable to use the numerical method in single precision.

4) Matrix assembly speed

The mesh is refined to 20602 patches and the matrix is computed in single precision. As shown on fig. 5, numerical technique is faster than the pure analytical one (3-4 times). There is one order of magnitude in computation times between the CPU computation and the low cost Geforce 230M one, for which the matrix is assembled in 9 parts and transferred to the host computer. Finally, two orders of magnitude are noticed between the CPU and the Tesla C1060 computations. The speed up in using GPU is not the number of GPU cores because a GPU core is less powerful than an Intel CPU one.

5) Iterative solver speed

A gain of an order of magnitude between CPU and GPU in the solving of the linear system is noticed in spite of the data transfers between the CUDA device and the host computer. The computation is done in single precision and the meshing and solving parameters are the same as previously.

6) JCuda limitations Some commonly optimizations, like the use of pinned memory in order to speed up data transfers, are not attractive. Indeed the data transfer is faster but the access to these data from Java is slower. In fact JCuda is based on C library and performances decrease in using C memory scheme from Java.

F. Conclusion

A very efficient integration method has been developed, mixing analytical and numerical approach. Matrix assembly is faster without real losing of precision. Furthermore, the use of GPU thanks to JCuda library speeds up the computation of two orders of magnitude. The linear solving is also speeded up of one order of magnitude.

IV. FAR FIELD COMPUTATION

A. Fast Multipole Method B. Single precision computation C. JCuda implementation D. Results

E. Conclusion

V. CONCLUSION AND PERSPECTIVES

We have reported in this paper a strategy for computing electromagnetic fields on Java platform with the JCuda library on a standard computer. ….

Wavelets compression [START_REF] Scheiblich | Compression of the Radiative Heat Transfer BEM Matrix of an Inductive Heating System Using a Block-Oriented Wavelet Transform[END_REF], [START_REF] Scheiblich | Parallel Hierarchical Block Wavelet Compression for an Optimal Compression for 3-D BEM Problems[END_REF] which seems more parallelizable than the FMM will be investigated.

Fig. 1 .

 1 Fig. 1. Classical matrix assembly algorithm.

Fig. 3 .

 3 Fig.3. Simplified CUDA kernel for matrix assembly using Gauss integration technique.

Fig. 5 .

 5 Fig. 5. Comparison of matrix assembly times depending of architecture and integration technique (lower is better).

JCuda vectorized and parallelized computation strategy for solving integral equations in electromagnetism on a standard personal computer

  . Rubeck, B. Bannwarth, O. Chadebec, B. Delinchant, J-P. Yonnet and J-L. Coulomb Grenoble Electrical Engineering Laboratory, INP/UJF/CNRS UMR 5269, 38402 Saint Martin d'Hères Cedex, France

C


Keywords: 2D magnetostatic, NRM, MEC, BIM, TCM

Hybridizing the network reluctance method (NRM) and the boundary integral method (BIM) aims to take advantage of both methods. First for ferromagnetic materials in order to take non linearity into account, and second for surrounding air in order to include fringing and leakages accurately. The automation of such modeling has been done with dedicated software and is compared to the tooth contour method (TCM) and finite element (FEM) simulations, which is applied to an E-core actuator.

I. INTRODUCTION

Several methods are available to simulate electromagnetic actuators. The choice is made usually for purposes of accuracy or speed. This diversity is explained by the fact that there is no universal method which treats all problems and meets all needs. In fact, all these programs are complementary. For example, the finite element method aims to model objects with complex shapes (2D or 3D) including numerous physical behaviours (saturation, skin effect ...), but it requires significant computing resources. The method of equivalent magnetic circuit (MEC), also called reluctance network method (NRM) is well suited for quick calculations, however the model development is time consuming and requires significant expertise from the user.

The objective of our work is to develop a numerical simulation method, which, like the reluctance network method, is well suited for quick calculations, but whose implementation requires much less time and expertise.

Hybridizing techniques already exists like NRM-MoM [START_REF] Janet | Magnetic Moment and Reluctance Network Mixed Method Applied to Transformer's Modeling[END_REF] but the idea is now to combine NRM and the boundary integral method (BIM). On the one hand NRM is relatively easy to implement in the absence of magnetic leakage in the air, but becomes tedious if leakages are significant. On the other hand, BIM is well known for its effectiveness in fields modelling in the air, but taking into account the method of non-linear materials is more difficult. Hybridizing is then to couple these two methods, ferromagnetic materials will be modelled by NRM and the surrounding air by BIM.

A first method using such a hybridizing is the tooth contour method (TCM) [START_REF] Kuznetsov | Numerical modelling of electromagnetic process in electromechanical systems[END_REF] in which leakage flux of a NRM model is identified using a BEM numerical software package. A second method is FVM-BIM, which automates both NRM creation in iron (using finite volumes method) and air (using BIM).

II. MODELING

A. Application description

Fig. 1 shows the E-core actuator that is chosen because of its large leakage in order to compare each method. The mobile part is parameterized by its horizontal and vertical translations respectively x 0 and y 0 in order to compute horizontal and vertical forces (F x and F y ) in several configurations. Fig. 2 shows the flux lines obtained with a FEM simulation, in which y 0 has been defined to 10mm. A 2 nd order mesh is used in Flux2D™ with about 3000 nodes for x 0 = 0 mm and 4000 nodes for x 0 = -45 mm. A simulation takes about 1 second. FEM results will be our reference values.

B. FEM simulation

C. Network Reluctance Method (NRM)

Reluctool software [START_REF] Du Peloux | Automatic generation of sizing static models based on reluctance networks for the optimization of electromagnetic devices[END_REF] is used to define and simulate a reluctance network. In order to model translation in both x and y directions, the network include parameterized reluctances using analytical formulas. The computation time with linear material is about 10 milliseconds for each position, which is significantly faster than the FEM model.

D. Tooth Contour Method (TCM)

The TCM is a numerical method to obtain accurate airgap reluctances including leakage and fringing effects.

Hybridizing network reluctance and boundary integral methods: comparisons on an E-core actuator

Tran Van T. 1 , Martins Araujo D. To obtain the airgap reluctances, the analogy between electric fields and magnetic fields is used. With boundary element (BEM) software (ELECTRO2D) the reluctances in the airgap are obtained. When the reluctances are obtained in-and outside the magnetic material, the calculation time is comparable to original NRM models. Fig. 3 shows the evaluated reluctance network, in which only 12 of the 161 airgap reluctances are shown.

E. Finite Volumes Method -Boundary Integral Method (FVM-BIM)

A new hybrid method is developed in Grenoble, and is implemented in software called MAGOT, in which FVM and BIM are coupled automatically. Like NRM and edge element method equivalence [START_REF] Demenko | Reluctance network formed by means of edge element method[END_REF], FVM can be considered as equivalent with NRM. Magot software is then able to define a meshing parameter (number of elements in a flux tube cross section, n = 2 in Fig. 4) corresponding to the reluctances in the iron. Fig. 5 shows convergence and computation time regarding this subdivision parameter. The accuracy reaches less than 4% with 5 subdivisions leading to a computation time of 1 second. But a result with about 7% of accuracy can be reach in 30 milliseconds. This FVM is then coupled with BIM, as FEM and BEM can be coupled together [START_REF] Salon | The hybrid finite element-boundary element method in electromagnetics[END_REF]- [START_REF] Meunier | Hybrid finite element boundary element solutions for three dimensional scalar potential problems[END_REF].

F. Results comparison

Fig. 6 shows that the results for the hybrid methods (TCM and FVM-BIM), which are significantly better than for the NRM. Table I gives the forces for x 0 is 0 and -45mm

The NRM is the fastest method but requires deep knowledge in order to build a correct network topology and to make an accurate parameterization of reluctances. The results for the TCM are better than for the NRM, due to improvement of the leakage reluctances, however a similar network topology must be build in advance. FVM-BIM has been fully automated and reluctances meshing can be configured in order to deal with time/accuracy trade off leading rapidly to good results. The advantage of hybridizing is clear for such an application where leakages and fringing are important. Moreover, the automation of the modeling made with FVM-BIM allows having results very quickly and to do parameterized study with high confidence results since leakages are computed for each new position. In the full paper, details will be given on reluctance networks, modeling procedure of each methods, and comparison results on both vertical and horizontal forces.

Fig. 1 :

 1 Fig. 1 : Magnetic actuator parameters description

Fig. 2 :

 2 Fig. 2 : The flux lines of the E-core actuator with y0 = 10 mm and x0 = [0; -45] mm

Fig. 3 :

 3 Fig. 3 : Reluctance network including the evaluated airgap reluctances for one segment

Fig. 4 :

 4 Fig. 4 : Ferromagnetic part meshing in Magot software

Fig. 5 :

 5 Fig. 5 : Accuracy and computation time regarding the number of subdivisions of the iron

Fy

  1 , Coulomb J.-L. 1 , Delinchant B. 1 , Chadebec O.1 Vrijsen N.H. 2 , Jansen J.W.2 , Lomonova E.A.2

	(1) Grenoble Electrical Engineering	(2) Eindhoven University of
	lab	Technology
	BP 46, 38402 Saint-Martin d'Hères,	Den Dolech 2, 5612 AZ Eindhoven,
	France	The Netherlands
	Jean-Louis.Coulomb@g2elab.grenoble-inp.fr
			x 1 =x 2 =x 3 =30 mm,
			y 2 =y 3 =30 mm,
		x0	y 1 =40 mm,
			z 1 =30 m
			N.I=2000 A.t
			µ

r = 1000 (linear)

TABLE I VERTICAL

 I FORCE SIMULATED WITH EACH METHOD


Abir Rezgui

Laurent Gerbaud

Benoît Delinchant

VHDL-AMS to support DAE-PDE coupling and multilevel modeling

Keywords: AMS, 2D magneto static, PDE-DAE, V shaped Cycle, Multi-domain simulation, multilevel modeling

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The emergence of the modeling language VHDL-AMS, defined by IEEE 1076.1 standard [START_REF]Norme VHDL-AMS VHDL Analog and Mixed-Signal Extensions[END_REF] has proved extremely useful for the defining and modeling of multi physical and multi-technological systems. VHDL-AMS is a powerful hardware description language that allows the modeling of the behavior of mixed-signal (analog and digital), the multiphysics systems such as electrical, magnetic, thermal, and so on, as well as the hierarchical modeling of complex multi-physical systems and phenomena [START_REF]Norme VHDL-AMS VHDL Analog and Mixed-Signal Extensions[END_REF], [START_REF]Norme VHDL-AMS VHDL Analog and Mixed-Signal Extensions[END_REF]. It specifies the possibility to describe the behavior of a complex system by formulating directly the corresponding differential algebraic equations (DAE) of the system, and thanks to its standard language in terms of syntax and semantics, it ensures portable descriptions among different industrial simulators (Portunus™, SMASH™, Simplorer™…). Due to their complexity, partial differential equations (PDE) have been intentionally left out in VHDL-AMS language. This limits the accurate modeling of system blocks that include physics effects. Various approaches to solve PDE models [START_REF] Abd Rahman | Modelling of Mixed Physical-domain System[END_REF] or to extend VHDL-AMS standard are proposed [START_REF] Nikitin | Modeling partial differential equations in VHDL-AMS[END_REF].

The paper deals with the simulation of systems by coupling VHDL-AMS models and others models based on "black box" components. In the approach, a VHDL-AMS model calls external models. This approach should assure the interoperability between models, and should offer, to the designer, the possibility to use VHDL-AMS as a system description language integrating many modeling aspects. Our approach is applied on a classical E-shaped actuator with a coil bounded to its central leg. It shows how to model such a system by invoking the multi physical aspects and the multi-level modeling of VHDL-AMS.

II. COUPLING WITH VHDL-AMS

Since engineering often has to deal with optimization, some models come as functional descriptions which are not coded in a conventional description language (like VHDL-AMS) but rather in pure programming languages such as C/C++. To implement our approach, two schemes for the coupling of models can be defined: • coupling using a language: VHDL-AMS provides a "foreign" attribute to decorate functions [START_REF]Norme VHDL-AMS VHDL Analog and Mixed-Signal Extensions[END_REF].

• coupling using simulation environment functionalities.

A. Foreign VHDL-AMS keyword

A first coupling approach is possible thanks to VHDL-AMS "FOREIGN" attribute [1][5]. Fig. 1 shows the principle of this approach to integrate C functions in VHDL-AMS models. Basically, C/C++ models or dynamic linked library (DLL) are shown as external functions by the VHDL-AMS model, and allow an exchange of data [START_REF] Smash | Foreign C functions in VHDL-AMS and SPICE models[END_REF]. Then, if more complex algorithms are required for calculation, software components can be used thanks to "Interface Adapter" ensuring the communication with attribute "foreign". Fig. 2 shows the principle of the second approach. Many simulators provide a powerful C communication interface which allows dynamic data exchange between VHDL-AMS models and external models or applications on the same simulator. This offers the possibility to create new models or to encapsulate a fine modeling of subsystems using the programming language C/C++. Implementing models using the C-interface is helpful if the model behavior cannot be described by standard models or their combinations within macros models. Recent researches, inspired from the classical V-shaped design approach, tend to propose new design methodologies based on functional virtual prototyping (FVP), as developed by Y.Hervé [START_REF] Hervé | Functional virtual prototyping design flow and VHDL-AMS[END_REF], which increases each step of design by model and simulation studies. The subject of our study is to present the advantage to use VHDL-AMS for different modeling level of the same system, and it was carried out using only the classical V-shaped design.

VHDL-AMS to support DAE-PDE coupling and multilevel modeling

The methodology of modeling by VHDL-AMS contains three modeling description level: "functional", "behavioral and structural" and "physical" [START_REF]Norme VHDL-AMS VHDL Analog and Mixed-Signal Extensions[END_REF]:

-"Functional modeling": denoted by «system-level», consists in translating the characteristics of the system resulting from the functional specifications into functions which can be simulated. -"Structural Behavioral modeling": describes phenomena and interactions like digital-analog one and multidomain interactions. -"Physical modeling": studies the physical behavior of the device. The combination of modeling with VHDL-AMS and the Vshaped design approach ensures consistency of physical performance criteria and the potentials of functional performances.

IV. APPLICATION ON AN ELECTROMAGNETIC DEVICE

A. Description Fig. 3 shows a classical E-shaped actuator with a coil bounded to its central leg. Only the linear vertical movement of the moving part of the actuator is possible, which is subject to different damping sources due to viscous damping and translational linear spring.

B. Modeling

The dynamic model of the device represents a coupling between the electrical circuit feeding the coil and the mechanical part of the moving part. To model this device, magnetic, electrical and mechanical aspects are carried out by using computation tools with either fine modeling or behavioral modeling.

1) Magnetic: fine modeling using FEM FLUX2D software was used to model the magnetic force and the inductance value. On the one hand, FLUX2D solver was controlled and linked to the "Foreign" key word (1 st method), on the other hand, it was linked to the simulation software (2 nd method).

2) Magnetic: behavioural modeling using RNM The modeling of the magnetic behavior of the actuator was made by reluctance network method (RNM) (Fig. 4). For this, a library of all magnetic components, such as linear reluctance, air-gap, magnetic coil, was developed on VHDL-AMS. 3) Electric The state system modeling of the electrical circuit feeding the coil was defined by two differential equations. If state variables are x, their derivatives x'dot are maximal. The inductance value is provided by the magnetic model.

4) Mechanic

The state variable of the mechanical part of this system are the core position z, and its velocity v. So, the state system , was defined by two differential equations..

V. CONCLUSIONS

VHDL-AMS, as a unified modeling language, has to support all the necessary modeling concepts with mathematically founded semantics and to support multidomain aspects of systems. Using "Foreign" key-word or simulation environment functionalities allow the coupling of VHDL-AMS with PDE solvers or other numerical modeling tools. The paper demonstrates the use of this functionality in order to do multi-method modeling in order to perform multiphysical and multi-level simulations. The V-Shaped method will be detailed in the full paper and applied on an E-Core.

Fig. 1 :Fig. 2 :

 12 Fig. 1 : Coupling approach using "foreign" attribute B. Interfacing Solver with Programming Language C/C++ interface

Fig. 3 :

 3 Fig. 3: Magnetic actuator description

Fig. 4 :

 4 Fig. 4: Reluctance network modelling

Fig. 5 :

 5 Fig. 5: VHDL-AMS electrical model