Ide-Flore Kenmogne
email: ide-flore.kenmogne-fokam@etu.utc.fr

Danilo Alves De Lima
email: danilo.alves-de-lima@hds.utc.fr

Alessandro Corrêa

A Global Navigation Management Architecture applied to Autonomous Robots in Urban Environments

Keywords: Sensor-based navigation, Global navigation management, Digital maps

This paper presents a global behavioral architecture used as a coordinator for the global navigation of an autonomous vehicle in an urban context including traffic laws and other features. As an extension to our previous work, the approach presented here focuses on how this manager uses perceived information (from low cost cameras and laser scanners) combined with digital road-map data to take decisions. This decision consists in retrieving the car's state regarding the global navigation goal, selecting which local navigation task should be executed (either lane following or intersection maneuvers), providing some round constraints and further defining the reference trajectory to be met by the selected local task. This system was experimented in a real autonomous car and provided satisfactory results.

I. INTRODUCTION

In these last two decades, we have witnessed great advances in autonomous vehicles conception which resulted in a large number of work. The navigation of these vehicles consists in starting from an initial point (q init ) to a final destination (q goal ) autonomously. Lacroix et al. [START_REF] Lacroix | Autonomous navigation in outdoor environment: adaptive approach and experiment[END_REF] separated this mission in three mandatory steps: the perception and localization, the path planning and the decision processes; which is the triptych of the overall navigation task.

The path planning consists in determining which path to take all-through the environments in order to reach the goal. This step is fundamental in view of the static map structure, the robot dimensions and movement constraints, defined in the car's configuration space [START_REF] Choset | Principles of Robot Motion: Theory, Algorithms, and Implementations[END_REF]. Considering the robot's environment as geometric (composed of interconnected segments and shapes [START_REF] Sotelo | Vision-based navigation system for autonomous urban transport vehicles in outdoor environments[END_REF]), several solutions have recently been proposed in the literature. Most of these are waypoints-based approaches, which have shown to be less computational and to provide an entire path (mostly but not always optimal) that must be performed by the robot at each step (see [START_REF] Ibrahim | Way point based deliberative path planner for navigation[END_REF], [START_REF] Guo | Graphical model-based recursive motion prediction planning algorithm in stochastic dynamic environment[END_REF], [START_REF] Fernandes | CaRINA Intelligent Robotic Car: Architectural design and applications[END_REF] for more details). The implementation of these techniques was stimulated earlier in 2005 by the DARPA Grand Challenge participants, who applied them on off-road autonomous vehicles [START_REF] Thurn | Stanley : The Robot that Won the DARPA Grand Challenge[END_REF]. Along with the planning, a precise localization must be guaranteed while following the planned path. This constraint implies the use of closed-to-market sensors which are unfeasible for commercial vehicles.

In this work, we aim to warrant the global autonomous navigation of a vehicle in an urban environment. Contrary to the wide number of solutions seen, while using low-cost sensors, we will propose a behavioral architecture that will render the execution of local navigation tasks more independent from each other and therefore enable the system to react to external events. In this context, the information on Critical Points (CPs) is very useful to efficiently recover the road geometry (used to find the image reference that is to be followed by the local trajectory tracking component). The digital road map information is also important for the whole navigation task.

In the remainder of this paper, we provide in section II an overview of some proposed solutions and the sensors used to achieve the navigation objectives. Section III presents how the route planning has been done. The Global Manager's architecture is outlined in Figure 3 and detailed in section IV, where some of its advantages and possible applications are mentioned. Before concluding, we summarize our experimental results carried out in section V.

II. RELATED WORKS

Several contests such as the DARPA (Defense Advanced Research Projects Agency) Challenges in the USA (2004, 2005 and 2007); the Korean Autonomous Vehicle Competitions (AVC) in Korea (2010, 2012 and 2013), the Grand Cooperative Driving Challenge in Europe 2011 (GCDC) and many others have been organized worldwide to favor the development of autonomous intelligent vehicles. Considerable issues lie behind the design of these vehicles; starting from the problems related to the localization, the perception of the environment, to the establishment of a reference to be followed by the vehicle. Hereinafter, we underpin approaches implemented by the best participants in these competitions. We would like to emphasize that these solutions ensure a comprehensive autonomous navigation task; however, they do not focus on the more that rest behind the formalization of an overall navigation task management component. This choice was made thanks to our approach stems from a fusion of some of the benefits presented by each of those solutions differentiating mainly on the employed sensors.

Urmson et al. [START_REF] Urmson | Autonomous Driving in Urban Environments: Boss and the Urban Challenge[END_REF] (winner of the DARPA Urban Challenge competition) primarily proposed a behavioral subsystem composed of several sub-components such as a "Goal selection", a "Drive down road" and an "intersection handling". The "Goal selection" is responsible for estimating the vehicle local position according to the Waypoints given by a global route planner (RNDF) and further generating the local goal for execution by the motion planner. The "Drive down road", for its part uses the surrounding conditions of the vehicle (provided by visual sensors and the tight coupled laser scanner and radar) to find which road lane should be followed and what distance interval must be maintained between cars. Finally, the "intersection handler" uses data from the pan-head sensors to compute a save set of steering velocities that can be executed by the motion planner while avoiding collision with surrounding vehicles (output speed, acceleration and lane to be tracked). Although the aim has been reached, this is not marketable because it used accurate albeit expensive sensors such as : an inertial tightly coupled GPS bridges and outages (aided with a sub-meter accuracy), a high-dynamic-range camera (PGF) and many others. Miller et al. [START_REF] Miller | Team Cornells Skynet: Robust Perception and Planning in an Urban Environment[END_REF] also proposed a similar solution but did not really focused on the benefits of this component as Urmson et al. did.

Unlike Urmson et al., Bacha et al. [START_REF] Bacha | Odin: Team VictorTangos Entry in the DARPA Urban Challenge[END_REF] (third place in the race) who proposed a more structured behavioral management, used a hybrid deliberative/reactive architecture to analyze, select the appropriate behavior to adopt, and define a safe path. Their strength relies on the ability of this planner to take the right decision. Actually, the system is currently divided into several drivers associated with different types of situations. So, Taking the right decision consists in choosing the right driver to be applied at a specific time and the behavioral Integrator is responsible for this action. It thereafter outputs the corresponding velocities, accelerations and other parameter limits according to the selected behavior. Meanwhile, the sensors used here are not too different from the ones used before. Despite the difference between these two methods, they are both aimed to locally plan paths between the checkpoints (defining waypoints to guide the robot navigation) and correct the robot's lateral position by means of a high precision localization system and the RNDF information.

Ziegler et al. [START_REF] Ziegler | Making Bertha Drive An Autonomous Journey on a Historic Route[END_REF], participant of the GCDC, lately completed the autonomous guidance in urban environments using a preregistered reference path for the vehicle. A Differential GPS-aided inertial navigation system was used during the tracing process, giving more environment information to the digital map (such as lanes numbers, speed limits, stop signs, etc.) in a post-management step. This post-management was carried-out in a user-generated street map interface, called OpenStreetMap (OSM). Consequently, during the experiment, point feature and lane marking based localization was performed with less costly sensors. However, concerning the decision making process, a Behavior Generation component modeled using state chart notation was used. Therefore, a state will be activated if a series of rules are valid at the same time (leading either to Corridor Driving -while avoiding static obstacles -or Yielding&Merging-for intersection crossing).

Despite the fact that these methods were able to guide the robot globally, they were strongly dependent on the localization system's accuracy (in local tasks) . Victorino et al. [START_REF] Victorino | Safe Navigation for Indoor Mobile Robots[END_REF] Fig. 1. Fully actuated electric car APPACHE with some available resources.

proposed a new approach for indoor environments with a topological representation of the robot path. In this case, nodes were linked by corridors characterizing the relations of accessibility (between doors and crossing corridors). When in a corridor, the robot can move only to the next node, by using a sensor-based control strategy with no localization requirements for example. But, if the robot is near a node, it must localize itself with relation to this node and therefore take the appropriate direction (next in the topology).

If we extend the topology in [START_REF] Victorino | Safe Navigation for Indoor Mobile Robots[END_REF] to urban environments, the corridors are the roads, without any other way to go, and the nodes are the points where some action must be taken (intersections or traffic light stops). In our previous work [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF], the nodes were called by critical points (CPs), where, in comparison to the waypoints, they are sparsely defined and do not represent the shape of the path. Thus, while following the road, robot localization is not required (due to the absence of waypoints), and some sensor-based approach is enough to guide the vehicle to the next CP. At the CP, the vehicle localization can be performed by low cost systems, which can be less dependent on GPS sensors.

This navigation principle mimics the use of a simple personal vehicular GPS which guides the human driver to the desired destination by giving instructions only on some specific points. A topological segmentation was also applied, helping to define a high-level abstract world model which supports the decision-making process of an autonomous car. We will then use this model to define a Global Navigation Manager (decision maker) which is able to achieve the same goal as the previous solutions proposed in the literature but using less redundant information and also less costly sensors. In addition to its main goal in the Global Navigation process, this architecture could eventually be applied with different class of robots and further be suited for convoy missions. Moreover, here we validated this entire approach in a real autonomous car.

III. GLOBAL ROUTE PLANNING

The Global Route Planner developed in this work is similar to the one presented by Pereira et al. in [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF]; where the roads are represented by corridors and their intersections by critical points (CPs). Based on digital maps (provided by the OSM), a routing table was generated after a path planning by structuring the global navigation from q init to q goal . The main difference from [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF] is that, we this table have been divided in two parts, where the first is related to the road (way) description after the CP, and the second represents the CP geometry with the angle and direction of every intersecting road. The first part has a constant size, once the possible road parameters are always the same, and the second one has the size depending on the number of roads intersecting the CP. The columns description are as follow:

• Node: with the path's CPs ordered from q init to q goal ; • Lat and Long: the latitude and longitude coordinates of the CP in degrees; • V max : the maximum speed limit of the road (in km/h); • R lanes : the number of available lanes in the road (even if it is a two way road); • R way : the road classification in one way [START_REF] Lacroix | Autonomous navigation in outdoor environment: adaptive approach and experiment[END_REF] or not ( 2

IV. GLOBAL NAVIGATION MANAGEMENT

One cornerstone of the navigation task is the decision making process. Therefore, this paper presents a formalized sensor-based navigation management system aiming to deal with the global navigation problem in an urban environment. Considering the autonomous navigation architecture presented in Figure 3, we propose a navigation management responsible of executing the global route generated by the mission-planner (based on the digital map provided by OSM and low cost sensors). This subsystem would be able to perform lane changes, handle transitions between local tasks and even provide information about the vehicle state (x, y, θ) and the surrounding road to a nearby infrastructure or vehicle if necessary.

The global architecture in Figure 3 puts an emphasis on the location of our behavioral component in the global navigation context. This figure also indicates different links and data exchanges between this component and other navigation subsystems. In this section, we present how the navigation state is defined (possible states and their representation). Then, we highlight the methodology (algorithm) employed to determine this state at a given time. And finally, we present the way transitions are being conducted in other to derive the local task reference trajectory that will be tracked by the control component.

Navigation State Transitions

While structuring the navigation manager, we were led to analyze and consider the Global Navigation problem as successive states in which the vehicle could be all along the process of its mission. Figure 4 depicts these different states and the connections between them illustrate the possibility to go from one state to another. START POINT (state0): represents the vehicle at the initial position q init which is the first CP in the Routing Table . To reach the next CP, the vehicle should pass through a corridor. This will imply the navigation task to switch to the State1 (if, of course q init ̸ = q goal ).

ROAD FOLLOWING (state1): represents the vehicle is computing a road lane following. This means the vehicle is between two critical points. As soon as it approaches one of these critical points, the navigation manager detects this intersection and immediately change the navigation state to State2. It could also stop all the process in case this detected CP is the q goal . ROAD INTERSECTION (state2): for the vehicle in making a road intersection maneuver. This state is maintained throughout the intersection traversal and it either switches back to State1 or it has reached the final destination and therefore changes the state to GOAL POINT(state3).

Determining the navigation state at a defined instant

Algorithm 1 describes how the navigation state is being retrieved and updated at each time by the navigation management. It gives details on the conditions of transition from a state to another. First of all, if the robot is at the initial position (meaning that State = START POINT), we find the next critical point to be reached in the RT and automatically change the state to ROAD FOLLOWING. Then, if the state was already on road lane following, the State selector computes both the distance to the critical point dist (Compute distance2CP) and the angle difference δ θ (compute angle diff) between the actual car direction and the desired CP's direction. If this distance is lower than the previous distance, lower or equal to a defined threshold d max and δ θ lower than an angle tolerance (function of street corner angle recovered using the local occupancy grid); then the intersection crossing is activated (State = ROAD INTERSECTION). This whole conditions have been tested by the function Verifydistance2Intersection and is tested true if bool = 1. Otherwise, the state remains the same or passes to GOAL POINT if the next CP is the goal. The possibility to switch from ROAD INTERSECTION to ROAD FOLLOWING is presented when the distance to the actual CP is increasing while the δ θ decreases until it reaches an intersection-living distance threshold smaller than d max ; this marks the end of the traversal. This is verified by the function Verifydistance2Intersection, which outputs bool = 1 (intersection traversal ended, the state move back to ROAD FOLLOWING). Finally, the State = GOAL POINT is reached whenever the vehicle is very close to the last critical point.

Algorithm 1 Retrieve the Global Navigation Task state. Get Navigation State(x, y, θ, CP, d max )

1: dist = Compute distance2CP(x, y) 2: δ θ = compute angle diff(θ, θ CP ) 3: if (CP == q init andState == ST ART P OIN T ) then 4:
find the next Critical point CP If CP == q g oal 14:

State = GOALP OIN T The proposed algorithm is very generic, simple and can therefore be adapted to more complex situations.

Global manager's procedure

After having determined at what level the vehicle is, the global behavioral component computes, according to the founded state, the motion that must be performed by the local task. We have distinguish two main states :

• State 1 : Road Lane following
In case the vehicle is between two critical points, we used the methodology presented by Lima et al. [START_REF] Lima | An image based dynamic window approach for local navigation of an autonomous vehicle in urban environments[END_REF]; where the navigation reference is represented in the sensor frame. It consist on using the perceived information from the camera to define a visual objective (Visual Servoing Control) that must be reached. The local controller will then be able to regulate the actual lateral error based on this input. The approach is known not reactive. Meanwhile, the authors proposed another methodology called Image-based Dynamic Window Approach (IDWA) to avoid obstacles present in the local occupancy grid; but the final approach Fig. 5. Global navigation experiment for the routing table defined on Table II.

that has been adopted is a hybrid controller named as VS+IDWA, using the benefits of both control design.

• State 2 : Road intersection maneuver (traversal)

The behavioral manager, in case of an intersection, categorizes two main situations. The case the CP is a roundabout or a normal intersection (either left turn or right turn). It constructs a reference based on the principle presented by Jochem in [START_REF] Jochem | Vision-based neural network road and intersection detection and traversal[END_REF]. This consist on reconstructing the intersection's instantaneous rotation center to compute the curvature that the car most follow (arc shape formed by the tangents between the actual CP and the car directions). Once the turning radius is computed, the management express it in the sensor frame (approximate the curve by a polygonal shape) and then activate the image-based hybrid controller VS+IDWA.

V. EXPERIMENTAL RESULTS

To validate the proposed global navigation management, we have implemented and integrated it on the experimental car APPACHE (see Figure 1), using its resources for environment perception and local navigation control. This vehicle is an important test bed for our navigation approach, but, due to the current modifications, it is still not allowed to run in the city roads. However, at the Innovation Center of the Université de Technologie de Compiègne there exist a circuit called SEVILLE; used for autonomous navigation experiments. This circuit possesses some challenging scenarios such as: tight road intersections, roundabouts, sharp turns, and up&downhills (≈ 3.5 • ), which was used during the validation of our methodology.

Table II presents the routing table used for our experiment using only the OSM data. The entire system was running at 10 Hz, while the experimentation was carried out in order to present the global navigation management capabilities during the road following and the road intersection maneuvers in full autonomous mode. As SEVILLE is a closed circuit, the path planning was not applied, it only required of the CPs where the vehicle must pass. The localization tolerance for each CP Fig. 6. Some screenshots of the environment perception system illustrate the management during the road following and intersection maneuvers tasks.

was set to 10 m, regarding the GPS and OSM errors.

Figure 5 shows the global navigation movement performed by the APPACHE for each CP in the Global Routing Table II. In this figure, we could predict the state transitions diagram through which our navigation task have been during each step of the process (states given by the management system). This state transition have been validated by the results obtained in Figure 6. So, not only did our system achieve the goal, but at each level of the task, it could perfectly identify intersections and provide a reference trajectory to be followed (yellow line on the screenshorts of Figure 6). In this Figure, the vehicle always follows the yellow line. The Red line means wrong way road and the green line means right way (or feasible way, but not the one to be followed). Each time a road intersection is approached, the Global Management component outputs an image with two lines (in our scenario) representing the way to take and all other ways with different color than yellow according to the feasibility of the direction. Approaching the GOAL POINT, we can see that the line feature outputed delimit the location of the goal by painting in dark all the point of the line situated after this GOAL POINT.

Additionally, the global navigation management could also perform the speed limitation for each local road. Figure 6 shows the linear velocities applied for the complete movement with the limitation performed for each routing table row. In this figure, the global management is evidenced during the CP transitions, where the road nominal speed (defined as 5, 7, and 13 km/h) is limited to 5 km/h; this to allow safe intersection maneuvers. Also, in order to mimic a human driver behavior and prepare the vehicle for the arrival of an intersection, we anticipate (d max ) the speed limitation to 5 m before the CPs. This distance allows a smooth arrival in the intersection for the current speed limitations. However, for high speed applications, this distance must be related to the vehicle maximum deceleration for a smooth speed reduction.

VI. CONCLUSION AND FUTURE WORK

This paper proposed a new formalized sensor based navigation management system aiming to deal with the global navigation problem of an autonomous vehicle in an urban environment. We focused on the structure of the behavioral architecture, mainly on the way decision making is being handled to ensure a more robust transition and reference processing between the lane following and the intersection maneuvering. Our final system was integrated in the AP-PACHE, an experimental car equipped with low cost sensors. Although the information available was very imprecise , the system successfully accomplished its mission by providing instantaneous feedback of its state and also retrieving the appropriate reference trajectory corresponding to each current local task between two CPs. The results showed the good performances of the proposed solution which was also able to provide free-to-go trajectory in case of an obstacle detection. Despite the system's efficiency, the problem of discontinuity of the image reference at intersection traversal still needed a deeper consideration. In this context, the methodology for intersection traversal shortly introduced in Section IV will be developed in future works. In addition, more road constraints must be considered, like road lanes, stop signs, the road priority, etc., to allow road following and intersection maneuvers closer to real situations.

  ); • R a : the CP definition as roundabout point or not; • N roads : the number of roads at the CP; • Θ 1..N roads and D 1..N roads : the sets of angle and direction for each road at the CP, ordered counterclockwise. Θ 1..N roads are the formed angles between the current road and the next one in the intersection. D 1..N roads indicate if the road is the path to follow(4), a wrong way (arrival) road (3), a right way (departure) road(2), an arrival/departure road(1), or if it is the goal point(0).

Fig. 3 .

 3 Fig. 3. Autonomous Navigation Global architecture including a Navigation Management component (BehavioralComponent).

Fig. 4 .

 4 Fig. 4. Global navigation state transitions.

5 :

 5 compute car distance to the critical point dist 6: change State = ROADF OLLOW IN G 7: end if 8: if (State == ROADF OLLOW IN G) then 9: bool = Verifydistance2Intersection(dist,δ θ )

TABLE I ROUTING

 I

lanes Rway Ra N roads Θ1 D1 Θ2 D2

  TABLE FOR THE PATH REPRESENTED IN THE FIGURE 2. Example of a road intersection illustrating the routing table definitions, where the yellow arrows indicates the traffic way.

				CP description					CP geometry
	Node Vmax R qinit Lat Lon 49.392 2.798 30	1	1	0	1	180	4	-	-
	1	49.400 2.800	50	2	2	0	2	110	4	290	1
	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .

Fig. 2.

  Table I presents an example of routing table for the condition shown in the Figure I. A complete routing table based on the OSM data will be presented in the experimental results of this work (Section V). The global navigation, based on the routing table, is performed by following the directions indicated as path to follow (4) at each CP. The management of this table during a sensor-based navigation is presented next.

TABLE II ROUTING

 II

  TABLE FOR THE GLOBAL NAVIGATION MANAGEMENT EXPERIMENTS.

			CP description						CP geometry
	Node	Lat	Lon	Vmax R lanes Rway Ra N roads Θ1 D1 Θ2 D2
	qinit	49.4021061 2.7948685	13	1	1	0	1	180	4	-	-
	1	49.4024539 2.7953365	7	1	2	0	2	135	4	180	3
	2	49.4024597 2.7955631	5	1	2	1	2	180	4	330	3
	3	49.4026169 2.7955953	7	1	1	1	2	180	4	210	2
	4	49.4024800 2.7953713	13	1	2	0	2	180	4	315	2
	q goal 49.4021061 2.7948685	13	1	0	0	1	180	0	-	-

ACKNOWLEDGEMENT

This work was partly carried out in the framework of the Equipex ROBOTEX (Reference ANR-10-EQPX-44-01) and funded in the framework of the Labex MS2T (Reference ANR-11-IDEX-0004-02).


Danilo Alves De Lima

Alessandro Corrêa

Sensor-based

may come

Sensor-based control with digital maps association for global navigation: a real application for autonomous vehicles

Danilo Alves de Lima and Alessandro Corrêa Victorino

Abstract-This paper presents a sensor-based control strategy applied in the global navigation of autonomous vehicles in urban environments. Typically, sensor-based control performs local navigation tasks regarding some features perceived from the environment. However, when there is more than one possibility to go, like in road intersection, the vehicle control fails to accomplish its global navigation. In order to solve this problem, we propose the vehicle global navigation based on a topological representation of the environment using only digital maps and low cost sensors. The approach was developed for two main tasks: road following and road intersection maneuvers. The final solution was completely implemented in a real autonomous car and tested in a challenge circuit, showing the viability of our solution.

Index Terms-Sensor-based navigation, Global navigation management, Digital maps.

I. INTRODUCTION

The complete navigation of self-driving cars is the ability to go from an initial position (q init ) to a final destination (q goal ) in the world, considering that these points are not connected in the robot's workspace. To achieve that, motion planning is fundamental to define a feasible path connecting these points. In this context, several approaches were proposed in the past years for nonholonomic car-like robots motion planning [START_REF] Svestka | Motion planning for car-like robots using a probabilistic learning approach[END_REF], [START_REF] Lamiraux | Smooth motion planning for car-like vehicles[END_REF], defined in the car's configuration space. However, these strategies normally presents a high computational cost associated with the size of the urban environments, mainly if it is dynamic.

To deal with these problems, it is common to consider the vehicle as a punctual holonomic robot during the path planning and simplify the environment in geometric forms and segments. These elements edges define the points where the robot must pass (waypoints). Waypoints can be seen as a simplification of the entire robot path, recovered by connecting each consecutive waypoint. The navigation based on waypoints was applied by several participants in the DARPA Grand Challenge of 2005 for off-road autonomous vehicles [START_REF] Thrun | Stanley: The robot that won the DARPA Grand Challenge[END_REF], provided in the route definition data format (RDDF). Although the waypoints were able to guide the vehicles, they were previously defined by the organizers of the challenge with high precision sensors. In addition, the robots have required accurate localization systems while performing the path following tasks, impracticable for commercial cars.

In urban environments, one way to define waypoints is using road network maps, also known as digital maps. Digital

The authors are with Heudiasyc UMR CNRS 7253 Université de Technologie de Compiègne. Danilo Alves Lima holds a Ph.D scholarship from Picardie region. Contact authors danilo.alves-de-lima@hds.utc.fr maps provide rich information about streets, tracks, railways, waterways, points of interest, also detailing the number of lanes and speed limits. There are many approaches using digital maps to improve the global navigation problem, like enhancing the vehicle localization [START_REF] Schindler | Vehicle self-localization with high-precision digital maps[END_REF] and the difficulties of autonomous vehicles navigation as a whole [START_REF] Hentschel | Autonomous robot navigation based on openstreetmap geodata[END_REF]. Based on digital maps, the DARPA Urban Challenge [START_REF] Buehler | Special issue on the 2007 DARPA Urban Challenge, part I-III[END_REF] organizers provided the road network definition file (RNDF) to the participants. However, the RNDF itself was inaccurate and could add further errors if the vehicle blindly follow the road using that information. Thus, the participants locally planned paths between the checkpoints, defining waypoints to guide the robot navigation and correct the robot lateral position by means of a high precision localization system and the RNDF information. Recently, Ziegler et al. [START_REF] Ziegler | Making bertha drive-an autonomous journey on a historic route[END_REF] performed the autonomous navigation in urban environments using a prerecorded reference trajectory (waypoints) for the vehicle. They use a DGPS-aided inertial navigation system during the recording process, adding more environment information to the digital map (like lanes number, speed limits, stop signs, etc.) in a post-processing step. This post processing was performed in a user-generated street map interface, named as OpenStreetMap (OSM) [START_REF] Haklay | Openstreetmap: User-generated street maps[END_REF]. Thus, during the real car navigation, a map based localization was performed with a less costly system.

Although all these previous methods were able to guide globally the robot, locally they were highly dependent on the localization system. Victorino et al. [START_REF] Victorino | Safe navigation for indoor mobile robots,partii: Exploration, self-localization and map building[END_REF] proposed a new navigation approach for indoor environments with a topological representation of the robot path. In this case, each corridor was connected by nodes representing the relations of accessibility in the environment (doors and crossing corridors). Similarly, in urban environments, corridors are roads, without any-other way to go, and nodes are points where some action must be taken (like turn to a next road in the intersection or stop in a traffic light). In our previous work [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF], the nodes were called by critical points (CPs), where, in comparison to waypoints, they are sparsely defined and do not represent the path's form if connected. Thus, the robot localization is not required during the road following navigation (due to the absence of waypoints), and some sensor-based approach is enough to guide the vehicle to the next CP. At the CP, the vehicle localization can be performed by low cost systems, which can be less dependent on GPS sensors [START_REF] Hassan | Vehicle localization system based on ir-uwb for v2i applications[END_REF].

In this work, the goal is to enable the global navigation of a car-like robot using local navigation techniques and low cost sensors, mainly when there is more than one possibility to go. For that, it is fundamental to provide information associated to the CPs, as their connections and approximated geometric shape. Thus, digital maps (like the OSM) are important tools for this purpose. Differently from our previous work [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF], the robot navigation at intersections will be generalized and validated in a real car-like robot. In addition, the final robot control will be performed by a sensor-based strategy, using the benefits of the local navigation controller VS+IDWA 1 [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF].

This work was structured as follows: the Section II shortly describes the real car-like robot used in the experiments and some available resources; the routing table with the topological representation of the robot path is described in the Section III; Section IV presents the sensor-based approach for the global navigation, with the management of the routing table information; some navigation experiments in our real car are presented in the Section V; and the final considerations and future works are listed in the Section VI.

II. SYSTEM DESIGN

The experimental car is an electric Renault Zoe, hereby appointed APACHE and presented in the Figure 1. It was modified for computational actuation of throttle, brake, and steering, receiving commands over the CAN bus. Although these are important modifications, the vehicle still preserving the original human driver conduction capabilities, allowing e.g. a cooperative working mode between the human driver and the computer system. Besides the low level control (steering, throttle, and brake), several proprioceptive information can be retrieved from the CAN bus. Some of them are illustrated in the Figure 1, like the ABS and odometry from the wheels, the linear and steering velocities, the steering angle, etc.

This vehicle was equipped with two PointGrey cameras (Flea3 and BumblebeeXB3) for environment perception and a GPS for global localization (Haicom HI-204E), chosen as a low cost alternative for exteroceptive sensors. The monocular camera (Flea3) has a focal length of 1.8mm and large field of view (FOV) of ≃ 140 • to detect the road lane center features. The stereovision camera (BumblebeeXB3) has a focal length of 6.0mm and FOV ≃ 43 • , which allow the vehicle to detect road and obstacles up to 80m away from it. Both cameras were installed on a rigid structure on the car roof only for 1 VS+IDWA is an acronym for a hybrid controller composed by Visual Servoing (VS) and Imaged-based Dynamic Window Approach (IDWA). validation purposes. The GPS works at 1Hz and produce position errors in a radius of ≈ 5m.

To allow the implementation and execution of this navigation strategy, some embedded computers for data processing and vehicle control were installed in the APACHE (Figure 1). These computers implement:

• The vehicle localization based on GPS and CAN data fusion, by an extended Kalman Filter (EKF), with a final error of ±5 m at 10 Hz;

• The low level velocity control, providing setpoints for the throttle, brake, and steering actuators;

• The environment perception, with the road/obstacles detection and their representation in a local occupancy grid;

• The local navigation control, implementing the VS+IDWA approach [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF], used as a sensor-based control of our global navigation management; and

• The global navigation management, presented in the Section IV, to address the vehicle state in the routing table and provide the right parameters to the local navigation control.

III. GLOBAL ROUTE PLANNING

To validate the present navigation approach, roads were considered as corridors and their intersections as critical points (CPs), similarly to [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF]. Based on digital maps, provided by the OSM, a routing table was created by structuring the global navigation from q to q goal after a path planning. Differently from [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF], we divided this table in two parts, where the first one is related to the road (way) description after the CP, and the second part represents the CP geometry with the angle and direction of every intersecting road. The first part has a constant size, once the possible road parameters are always the same, and the second one has the size depending on the number of roads intersecting the CP. The columns are:

• Node: with the path's CPs ordered from q init to q goal ; • Lat and Long: latitude and longitude coordinates of the CPs in degrees;

• V max : maximum speed limit of the road (in km/h);

• R lanes : number of available lanes in the road (even if it is a two way road); navigation, based on the routing table, can be performed by following the directions indicated as path to follow (4) at each CP. The management of this table during a sensor-based navigation is presented next.

IV. SENSOR-BASED GLOBAL NAVIGATION

For the global navigation problem of guiding the autonomous car from q init to q goal , in urban environments, two local tasks were considered: road following (corridors) and road intersection maneuvers (at CPs). In this context, we proposed the association of a sensor-based controller for local navigation in order to accomplish both tasks, using the environment perception data (from exteroceptive sensors) and the path information from the routing table (previously defined in the Section III). This section presents the navigation management to switch between these two tasks, the local navigation control principle, and its application in the road following and road intersection maneuvers tasks.

A. Navigation management

The role of the navigation management is to determine which local task must be performed (road following or road intersection maneuver) during the car-like robot movement from q init to q goal , with relation to the routing table for its planned path. These tasks can be distinguished by defining how near/far from a CP is the robot. For that, we performed the APACHE movement in a track with CPs defined in all road intersections, presented in the Figure 6. The resulting trajectory, provided by its localization system, always intersects the black circles centered in the CPs. The circles consider the maximum robot localization and OSM errors, resulting in 10 meters of radius. Thus, the vehicle is near to a CP when the distance between them is lower than 10 meters and faraway otherwise. Approaching a CP, the maximum vehicle speed can be reduced to improve the navigation behavior and also apply some localization technique based on the local information [START_REF] Hassan | Vehicle localization system based on ir-uwb for v2i applications[END_REF]. After selecting the local task, the navigation management must define the image setpoint and the control parameters based on the robot perception system for the local navigation control.

B. Local navigation control

The sensor-based approach used for the local navigation control is a hybrid solution for road lane following named as VS+IDWA [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF]. Using cameras as the main sensors, it combines a Visual Servoing (VS) technique as a deliberative control, for road lane following, and an Image-based Dynamic Window Approach (IDWA), for reactive obstacle avoidance. The focus of this controller is to guide the image features set s = [X Y Θ] T to the desired configuration X = X * , Y = Y * , and Θ = Θ * , while avoiding obstacles. At the Figure 3, the features' set was extracted from a line formed by the tangent Γ of the path P at the point D = (X, Y ). Θ represents the angular offset from Γ to the axis -Y (positive counterclockwise). The path P is defined by the road lane center between the boundaries δ 1 and δ 2 .

In this work, the camera was installed on the car's sagittal plane and the desired features configuration represents the car in the center of the road lane, which means the line vertically in the center of the image frame. Thus, for the global navigation, the right image features must be provided to the VS+IDWA controller when performing the road following or the road intersection maneuver. Moreover, the environment information must be modified to avoid invalid movements and the robot speed limit adjusted during the navigation.

C. Local navigation: road lane following

Once the robot is on the road following state, the VS+IDWA controller [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF] can be applied as it was conceived to accomplish the current local navigation task. In this case, the environment perception capabilities of our experimental car APACHE is responsible to detect and provide the correct image feature, as defined in the Figure 3. With the previous knowledge about road speed limit (V max ), number of lanes (R lanes ), and traffic way (R way ), provided by the routing table, the speed setpoint can be adjusted and the lane detection algorithm optimized for each road. In addition, with the traffic way, the obstacle overtaking can be better controlled, avoiding, e.g., movements in discordance with the traffic laws.

D. Local navigation: road intersection maneuvers

Approaches for road following do not traditionally deal with road intersection maneuvers, requiring some intervention at these situations. With this problem in mind, we propose a complement for the VS+IDWA controller [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF], generalizing the idea presented in [START_REF] Pereira | Autonomous vehicle global navigation approach associating sensor based control and digital maps[END_REF]. The first part of this approach is described in the Algorithm 1, with the detection of the road intersection corners and the main steps illustrated in the Figure 4. It was inspired by the laser detection of tunnel intersections [START_REF] Larsson | Laser based intersection detection for reactive navigation in an underground mine[END_REF], where a similar laser reading in the local occupancy grid defines the obstacles position in polar coordinates (r obst , θ obst ) in a maximum perceived distance d max (see the pink dots in the Figure 4b). It also uses the car's width C width , to eliminate false positives, and the routing table information of N roads , with the number of roads at the intersection and their geometry (Figure 4c). These parameters define the path possibilities for the current robot state (Figure 4a), represented by the arrows (red and yellow colors) in the navigable area (black pixels) of the occupancy grid.

Algorithm 1 Find the road intersections. find intersection corners(r obst , θ obst , d max , C width , N roads ) 1: i l , i r = find main road limit(r obst , θ obst , d max ) 2: I polar = Plot(r obst , θ obst ) 3: for all I polar row between the columns i l and i r do 4: find the n pairs (θ obst,l [1..n] , θ obst,r [1..n]), representing the grey pixels limited by the white pixels 5:

if n == N roads then 6: if E dist (Cart(row, θ obst,l [1..n]), Cart(row, θ obst,r [1..n])) > C width then 7: return row, θ obst,l [1..n] , θ obst,r [1..n] 8: end if 9:
end if 10: end for 11: return false Following this algorithm, initially the main road limit is calculated by the Algorithm 2, which return i l and i r . They represent the left and right θ obst indexes with the smallest distance to the vehicle. These points reduce the search algorithm to the region defined between them. The final search is performed in the polar coordinates graphic I polar (Figure 4d), with the radius measurements linked by the white line dividing the area among the vehicle (grey Steps to find the road intersection corners for the car pose (a). Its respective local occupancy grid (b) and routing table state (c) illustrate the road structure at the intersection and the end of the navigable area given by pink dots with polar coordinates (r obst , θ obst ). These polar coordinates were used to detect the main road limits and to generate the graph (d), which defines the road intersections corners (in blue). They are better visualized in the grid point of view in (e), with the main road in yellow, the intersections in blue, and the destination road center in green.

color) and faraway (black color). Thus, the algorithm search for the row discontinuities, which represent the white-grey and grey-white variations between i l and i r . These points are the candidates for the road intersection corners and must be checked by the following two filters:

• N roads : verify if the number of n candidates is equal to the amount of expected roads at the intersection. Each candidate is equal to the pair formed by the whitegrey (row, θ obst,l [k]) and grey-white (row, θ obst,r [k]) variation, for k = [1.

.n];

• E dist : test if the Euclidean distance of each candidate, in Cartesian coordinates (Cart), is higher than the car's width (C width ).

Algorithm 2 Find the main road limit. find main road limit(r obst , θ obst , d max )

1: l dmin = r dmin = d max 2: for i = 1 to Size(θ obst )/3 do 3: if l dmin > r obst [i] then 4: l dmin = r obst [i] and i l = i 5: end if 6: if r dmin > r obst [Size(θ obst ) -i + 1] then 7: r dmin = r obst [Size(θ obst ) -i + 1] and i r = i 8:
end if 9: end for 10: return i l , i r

The valid candidates are illustrated by blue dots in the Figure 4d. The main road limits and the intersections found are represented in the grid point-of-view in the Figure 4e. By the intersection geometry, the wrong vehicle directions are obstructed with virtual obstacles (also known as geofencing) at the distance equal to row, which results in the Figure 5a. Once the VS+IDWA requires an image features set for the control calculation, the line connecting the vehicle to the next valid road center (green point at the Figure 5a) was projected in the image plane and defined as the new image features, as seen in the Figure 5b. For more details about the projection principle, see [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF]. This approach is always executed when the vehicle is on the road intersection. Thus, the center of the destination road (green point) can be used as a condition to check if the robot has reached the destination road. If this point is closer enough to the front of the car, the road intersection maneuver was successfully performed and we can switch to the next local navigation method.

V. EXPERIMENTAL RESULTS

To validate the proposed global navigation management, it was completely implemented in the experimental car APACHE (see Section II), using its resources for environment perception and local navigation control. This vehicle is an important test bed for our navigation approach, but, due to the current modifications, it is still not allowed to run in the city roads. Thus, at the Innovation Center of the Université de Technologie de Compiègne was constructed a circuit called SEVILLE for autonomous navigation experiments. It presents some challenging scenarios, like tight road intersections, roundabouts, sharp turns, and up and downhills (≈ 3.5 • ), which was used during the validation of our navigation approach.

With the system running at 10 Hz, a complete experiment was carried out to present the global navigation management capabilities during the road following and the road intersection maneuvers in full autonomous mode. The Table II describes the routing table in this situation using only the OSM data. Once this is a closed circuit, the path planning was not applied, requiring only the order of the CPs where the vehicle must pass. The localization tolerance for each CP was set to 10 m, regarding the GPS and OSM errors.

The Figure 6 illustrates the global navigation movement performed by the APACHE for the CPs of Table II. It Fig. 6. Global navigation experiment for the routing table defined on Table II. The initial position is defined by the blue arrow and the destination by the green point in the map. The robot positions are illustrated in red and the localization tolerance around each CP are represented by the black circles. Some screenshots of the environment perception system illustrate the management during the road following and intersection maneuvers tasks. also presents some screenshots from the environment perception system to illustrate the car navigation during the road following and intersection maneuvers tasks. In this figure, the results for the local navigation control based on the VS+IDWA [START_REF] Lima | A visual servoing approach for road lane following with obstacle avoidance[END_REF] were confirmed, with the robot moving smoothly during the path, trying to bring the image features error to the zero condition (vertically in the center of the image), and avoiding the obstacles detected in the occupancy grid (pink dots for real obstacles and blue dots for virtual ones). Note that the localization system variations make path following applications using only GPS points unfeasible, but they are well fitted in the CPs tolerance radius.

On this previous movement, the global navigation management performs the speed limitation for each local road. The Figure 7 shows the linear velocity profile for the complete movement with the limitation performed for each routing table row. In this figure, the global management is evidenced during the CP transitions, where the road nominal speed (defined as 5, 7, and 13 km/h) is limited to 5 km/h to allow safe intersection maneuvers. In order to mimic a human driver behavior and prepare the vehicle for the arrival of an intersection, we anticipate the speed limitation to 5 m before the CPs. This distance allows a smooth arrival in the intersection for the current speed limitations. However, for high speed applications, this distance must be related to the vehicle maximum deceleration for a smooth speed reduction.

In addition, the Figure 8 presents in detail one of these maneuvers, with some sequential frames. In these frames it is possible to see the robustness of this approach during the movement, detecting the intersecting road with low variations as well as the next road center (green dot). Note that the roads geometry at the intersection do not matter on fact, but can be considered to enhance the detection capabilities. The most important element here is a good representation of the road area by the perception system. If the next road surface is not well detected for example, or it is not on the sensor's field of view, the vehicle will not be able to exit from the intersection situation.

VI. CONCLUSIONS AND FUTURE WORKS

This work presented a new sensor-based approach to deal with the global navigation problem of car-like robots in urban environments, using only low cost sensors and digital map data. We focused on two main tasks, road following and road intersection maneuvers, performed with our local navigation controller VS+IDWA. The final system integrated all the available resources of the experimental vehicle APACHE and successfully accomplished the two main tasks, following the image reference and avoiding obstacles, real or not. The results showed the viability of the solution, even in the presence of a reduced number of sensors for environment perception. Nevertheless, increasing the field of view and, consequently, the number of sensors used is fundamental to assure safe navigation in more situations. Future works will enhance the actual perception configuration for better intersection detection. Moreover, better localization techniques should be tested to avoid GPS problems in urban canyons when approaching an intersection. In addition, more road constraints must be considered, like road lanes, stop signs, the road priority, etc., to allow road following and intersection maneuvers closer to real situations.

Fig. 1 .

 1 Fig. 1. Fully actuated electric car APACHE with some available resources.

Fig. 3 .

 3 Fig. 3. Image frame {I} representation for the camera installed on the sagittal plane of the car's roof. The road lane center projection P (in red) is related to the boundaries δ 1 and δ 2 (in yellow). Its tangent Γ (in blue) at the point D with angle offset Θ from Γ to the axis -Y , defines the image features X, Y , and Θ.

Fig. 4 .

 4 Fig. 4.Steps to find the road intersection corners for the car pose (a). Its respective local occupancy grid (b) and routing table state (c) illustrate the road structure at the intersection and the end of the navigable area given by pink dots with polar coordinates (r obst , θ obst ). These polar coordinates were used to detect the main road limits and to generate the graph (d), which defines the road intersections corners (in blue). They are better visualized in the grid point of view in (e), with the main road in yellow, the intersections in blue, and the destination road center in green.

Fig. 5 .

 5 Fig. 5. Final r obst and θ obst with the virtual obstacles defined by the method described in the Algorithm 1 (a) and the new image features (b), oriented to the center of the next road (green point).

Fig. 7 .

 7 Fig. 7. Velocity profile for the global navigation experiment I of Figure 6.

Fig. 8 .

 8 Fig.8. Vehicle movement during a road intersection maneuver, emphasizing the virtual obstacles (blue dots) and the destination road center (green dot) used to define the new image feature in these conditions (yellow line).

  Θ 1..N roads and D 1..N roads : sets of angle and direction for each road at the CP, ordered counterclockwise. Θ 1..N roads are formed angles between the current road and the next one in the intersection. D 1..N To illustrate these parameters, the Table I presents the routing table for the condition shown in the Figure2. A complete routing table based on the OSM data will be presented in the experimental results of this work (Section V). A global

				CP description					CP geometry
	Node	Lat	Lon	Vmax R lanes Rway Ra N roads Θ1 D1 Θ2 D2
	qinit	49.392 2.798	30	1	1	0	1	180	4	-	-
	1	49.400 2.800	50	2	2	0	2	110	4	290	1
	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .
												roads indicate
							if the road is the path to follow(4), a wrong way
							(arrival) road (3), a right way (departure) road(2), an
							arrival/departure road(1), or if it is the goal point(0).

• R way : road classification in one way

[START_REF] Svestka | Motion planning for car-like robots using a probabilistic learning approach[END_REF]

or not (2); • R a : CP definition in roundabout point or not;

• N roads : number of roads at each CP; • Fig.

2

.

Example of a road intersection to illustrate the routing table definitions, where the yellow arrows indicates the traffic way.

TABLE II ROUTING

 II TABLE FOR THE GLOBAL NAVIGATION MANAGEMENT EXPERIMENTS.

			CP description						CP geometry
	Node	Lat	Lon	Vmax R lanes Rway Ra N roads Θ1 D1 Θ2 D2
	qinit	49.4021061 2.7948685	13	1	1	0	1	180	4	-	-
	1	49.4024539 2.7953365	7	1	2	0	2	135	4	180	3
	2	49.4024597 2.7955631	5	1	2	1	2	180	4	330	3
	3	49.4026169 2.7955953	7	1	1	1	2	180	4	210	2
	4	49.4024800 2.7953713	13	1	2	0	2	180	4	315	2
	q goal 49.4021061 2.7948685	13	1	0	0	1	180	0	-	-

ACKNOWLEDGEMENT

This work was carried out and funded in the framework of the Equipex ROBOTEX (Reference ANR-10-EQPX-44-01). It was equally supported by the French Picardie project VERVE, French Government, through the program "Investments for the future" managed by the National Agency for Research, and the European Fund of Regional Development FEDER. The authors wish to thank the helpful assistance of Gérald Dherbomez, Giovani B. Vitor, Pierre Hudelaine, and Thierry Monglon during the experiments.


Edwin Chacon-Golcher

Sever A Hirstoaga

Mathieu Lutz

Optimization of Particle-In-Cell simulations for Vlasov-Poisson system with strong magnetic field

We study the dynamics of charged particles under the influence of a strong magnetic field by numerically solving the Vlasov-Poisson and guiding center models. By using appropriate data structures, we implement an efficient (from the memory access point of view) particlein-cell method which enables simulations with a large number of particles. We present numerical results for classical one-dimensional Landau damping and two-dimensional Kelvin-Helmholtz test cases. The implementation also relies on a standard hybrid MPI/OpenMP parallelization. Code performance is assessed by the observed speedup and attained memory bandwidth.

Introduction

This work focuses on the numerical solution of a multiscale, four-dimensional Vlasov-Poisson equation which statistically describes the evolution in time of charged particles under the action of an external magnetic field and a self-generated electrostatic field. This system models a strongly magnetized plasma constituted by electrons and an uniform and motionless background of ions. When tracking the electrons on a long enough time scale, a drift phenomenon occurs due to the self-consistent electric field in the plane orthogonal to the magnetic field. In order to study this multiscale behaviour we consider the particular configuration where the strong magnetic field is aligned with the x 3 direction: B = (0, 0, 1/ε) T , where ε is a small parameter and where we assume that the electric field created by the particles evolves only in the plane orthogonal to the magnetic field. It is well-known (see [START_REF] Bostan | The Vlasov-Maxwell System with Strong Initial Magnetic Field: Guiding-Center Approximation[END_REF]9]) that when ε goes to zero, this Vlasov model may be approximated by a reduced model: the guiding center equation. However, when solving the Vlasov equation, the small parameter induces fast time oscillations in the solution, thus making classical numerical schemes very expensive.

In this context, the aim of this work is twofold. Firstly, we want to compute the numerical solution of a four-dimensional Vlasov-Poisson system and to illustrate a convergence result. Secondly, we aim at implementing classical numerical methods in a sufficiently efficient way that enables us to produce simulations executing in reasonable time. The first aim deals with verifying the results of the computations, i.e. ensuring that the code solves the mathematical model in question. To accomplish this we use known analytical approximations or known conserved quantities. The second aim is related with code implementation features designed to reduce the computation time. In this regard we deal with issues like optimized data structures, streamlined memory access and parallelization. In addition, we notice that, with the further aim to improve the recent time scheme in [8] and to extend it to more general problems than the aforementioned multiscale model, we need to have at our disposal a fast code yielding the reference solution. In the following, we understand by reference solution of a model its numerical solution calculated with classical numerical schemes and small enough numerical parameters.

Since a Vlasov model has the structure of a transport equation in phase space, most of the numerical methods are based on the study of its characteristics. In this paper we use a standard Particle-In-Cell (PIC) method to solve the model (see [START_REF] Birdsall | Plasma physics via computer simulation[END_REF][START_REF] Sonnendrücker | Approximation numérique des équations de Vlasov-Maxwell, notes de cours de M2[END_REF]). We are thus led to advance in time a number of macro-particles following these characteristics. Although particle methods are known to be noisy (see [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]), methods exist to improve their precision, as the linearly-transformed particle scheme recently studied in [START_REF] Campos-Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF] for Vlasov-Poisson simulations. An alternative approach is to use semi-Lagrangian methods (see [START_REF] Sonnendrücker | Approximation numérique des équations de Vlasov-Maxwell, notes de cours de M2[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF]) which show to be very accurate when solving Vlasov-like equations. Lastly, another powerful method to address noise in PIC simulations is to increase the number of particles in the simulation, but this is only practical when the implementation is efficient. We pursue this method here.

In this work we follow a standard PIC framework for solving a Vlasov-Poisson system. In order that phenomena occurring at small scales can be captured, one needs a very large number of macro-particles and fine enough grids in time and space. This yields a very computationally demanding problem. To make the resulting simulations execute in an acceptable time, we choose standard methods of parallelization (using the well-known multiprocess and multithreading libraries MPI and OpenMP) and less widely employed methods of code design using specialized data structures to optimize memory access patterns. See the approach developed in [START_REF] Bowers | Accelerating a Particle-in-Cell Simulation Using a Hybrid Counting Sort[END_REF] and [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF]. We have thus performed simulations able to process, at their best, 265 million particles per second in average, on one node with 16 cores. The implementation was done in the library SeLaLib1 . For verification of the code we have performed simulations of the wellknown Landau damping problem. Additionally, we simulated the Kelvin-Helmholtz instability test case. Finally we have compared the numerical solution of the multiscale four-dimensional Vlasov-Poisson system against the one for the guiding center model: we show the convergence result and assess it in quantitative and qualitative ways (see Section 2.3).

The rest of the paper is organized as follows. In Section 1 we write the equations to be solved and we briefly describe the PIC method for the numerical solution. In Section 2 we detail the PIC implementation, discuss about its efficiency, and show the numerical results concerning the Vlasov-Poisson and the guiding center models. We conclude in Section 3 and draw some lines of future work.

1 The numerical approach for the Vlasov-Poisson model

The Vlasov-Poisson and the guiding center models

The Vlasov-Poisson system in which we are interested in this paper is

                 ∂ t f ε + 1 ε v • ∇ x f ε -e E ε (x, t) + 1 ε v ⊥ • ∇ v f ε = 0, E ε (x, t) = -∇ x φ ε (x, t), -∆ x φ ε = ρ i + ρ ε e , ρ i = 1 and ρ ε e (x, t) = -e R 2 f ε (x, v, t)dv, f ε (x, v, t = 0) = f 0 (x, v) , (1)
where the unknown f ε = f ε (x, v, t) represents the distribution of electrons in phase space (position and velocity coordinates) and in time in the plane orthogonal to the magnetic field. We denote by -e = -1 the electron charge,

x = (x 1 , x 2 ) the position variable, v = (v 1 , v 2 ) the velocity, v ⊥ = (v 2 , -v 1 )
, and by E ε the electric field induced by the particles. The electric field is obtained by solving a Poisson equation coupled with the particle's equation through the source terms: ρ i denotes the constant ion charge density and ρ ε e stands for the electron charge density. The term (1/ε)v ⊥ corresponds to the projection of the magnetic force in the plane orthogonal to the magnetic field. The other singularity in the equation above is due to a time rescaling that allows us to study the long-term behaviour of the electron distribution function and thus, to capture the drift phenomenon. The regime in which the Vlasov-Poisson system (1) is written corresponds to a time rescaling of the classical drift-kinetic model. The scaling procedure leading to this regime is given in [START_REF] Frénod | Two-scale expansion of a singularly perturbed convection equation[END_REF] and the time rescaling is done in [9] (see also Section 5.2.1 in [8]).

The limit model of (1) is derived in [START_REF] Bostan | The Vlasov-Maxwell System with Strong Initial Magnetic Field: Guiding-Center Approximation[END_REF]9], under some hypotheses for f 0 . The authors have established that the electron charge density in (1) converges in some ways, when ε goes to zero, towards the unique solution to the guiding center model:

         ∂ t f GC + E ⊥ • ∇ x f GC = 0, E (x, t) = -∇ x φ(x, t), -∆ x φ(x, t) = ρ i -f GC , f GC (x, t = 0) = R 2 f 0 (x, v) dv. (2)
As already mentioned, the numerical solution of ( 1) is based on the characteristics of the Vlasov equation:

     dX(t) dt = 1 ε V(t), dV(t) dt = - 1 ε E ε (X(t), t) - 1 ε 2 V ⊥ (t). (3)
The main difficulty when solving these ODEs is in the stiff term 1/ε 2 which requires that the explicit numerical schemes use a very small time step. This drives the interest to the asymptotic model in (2) whose solution does not contain high oscillations in time. However, the aim of this work is to produce a reference solution of the Vlasov-Poisson equation ( 1) for several values of ε. We notice that this asks for long-term simulations since the final time is of order 1 while the numerical schemes need to resolve the 1/ε 2 term.

The Particle-In-Cell method

We solve the Vlasov-Poisson system (1) on a doubly periodic physical domain by a standard PIC method (see [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]). This consists in discretizing the distribution function at every time step t n = n∆t by a collection of N p macro-particles with coordinates (x n k , v n k ) and then applying a regularizing step with a convolution kernel S

f ε S (x, v, t n ) = 1 N p Np k=1 S (x -x n k ) S (v -v n k ) . (4)
As usual on cartesian meshes, we used B-splines as the convolution kernel (see [START_REF] Sonnendrücker | Approximation numérique des équations de Vlasov-Maxwell, notes de cours de M2[END_REF]). More precisely, tests were carried out with B-splines of degree 1 and 3, but all the simulations in this paper were mainly performed with a degree 1 B-spline, which in one dimension of the physical space writes

S(x) =      1 ∆x 1 - |x| ∆x if |x| < ∆x 0 otherwise.
After a randomized initialization for positions and velocities, the macro-particles are advanced in time with a leap-frog discretization of the characteristics in ( 3)

             x n+1 k = x n k + ∆t ε v n+ 1 2 k v n+ 1 2 k = v n-1 2 k - ∆t ε   E n (x n k ) + 1 ε   v n+ 1 2 k + v n-1 2 k 2   ⊥    . (5)
The electric field is computed by solving the Poisson equation on a uniform cartesian grid by a Fourier method. To this aim, at every time step the particle charge density is computed on the grid by

ρ ε S (x, t n ) = -e 1 N p Np k=1 S (x -x n k ) . (6)
In order to get conservation of the total momentum, we use the same convolution kernel S in order to compute the field at the particle positions:

E n (x, t) = ∆x 1 ∆x 2 j E n j S (x -x j ) , (7)
where x j j stand for the grid nodes and E n j j denote the grid values of the electric field at time t n .

PIC implementation and simulations

The multiscale features of the model above require a large number of particles and small time steps in order to resolve the fine structures in phase space, which makes this an expensive computation. In order to make simulations run in reasonable times, code optimization and parallelization are needed. We test the performance of the PIC implementation only on a standard Landau damping test case, briefly recalled in section 2.1. Speedup due to parallelization and access to memory curves of the code are given in section 2.2. Finally, in section 2.3, we perform four-dimensional simulations and show results that illustrate the convergence result of the Vlasov-Poisson equation [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] to the guiding center model [START_REF] Bostan | The Vlasov-Maxwell System with Strong Initial Magnetic Field: Guiding-Center Approximation[END_REF].

All the simulations in this paper were done at the computing center of the Université de Strasbourg (Mésocentre), on Intel(R) Xeon(R) E5-2670 @ 2.54 GHz processors. The typical node used at the Mésocentre contains 2 sockets, each with 8 cores, 32 GB of local memory and an L3 cache size of 20 MB.

The PIC code uses standard parallelization libraries for both, the distributed memory paradigm (using multiple software processes with MPI libraries), and the shared memory paradigm (multiple threads with OpenMP). While using these libraries, care was taken to associate both processes and threads to the underlying hardware to reduce the requests to remote memory. For this purpose, software processes were placed in single processor sockets (the hardware unit that groups a set of cores) and software threads in single cores. In this way, we attempted to minimize the negative impact of non-uniform memory access (NUMA) which can happen when a software process uses cores placed in different sockets, for example. Thus we employ parallelization with multiple processes running on individual sockets and within these processes we parallelize with multiple threads. Details of what is parallelized in the code with MPI and with OpenMP and other optimization techniques are discussed below.

Landau damping simulations

The purpose of this section is to solve numerically the Landau damping test case in order to verify the code (see [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Sonnendrücker | Approximation numérique des équations de Vlasov-Maxwell, notes de cours de M2[END_REF] and the references therein). The relevant equations are

             ∂ t f + v • ∇ x f -eE(x, t) • ∇ v f = 0, E (x, t) = -∇ x φ(x, t), -∆ x φ = ρ i + ρ e ρ i = 1 and ρ e = -e R 2 f (x, v, t)dv, f (x, v, t = 0) = f 0 (x, v) (8)
where -e = -1 stands for the electron charge and the initial condition is given by

f 0 (x, v) = 1 2πv 2 th 1 + α cos (kx 1 ) exp - v 2 1 + v 2 2 2v 2 th , (9)
where v th = 1, k = 0.5, and

Ω x = 0; 2π k × [0; 1] (10)
is the domain in the physical space. We use classical periodic boundary conditions for the physical domain. Even if the data stand in four-dimensional phase space, this test case is mainly two-dimensional since the perturbation α occurs only in the x 1 direction of the physical space.

Before showing the implementation and performance issues, we discuss the numerical results. It is known that the total energy of the Vlasov-Poisson model is constant in time:

d dt Ωx R 2 f (x, v, t) |v 2 | 2 dvdx + 1 2 d dt Ωx |E (x, t) | 2 dx = 0. ( 11
)
The numerical conservation of this quantity is important for the verification of the simulation.

In Fig. 1 we plot the time history of the total energy for different values of the perturbation α.

These runs used 2 million particles, a grid of 256 x 2 cells and a time step of 0.1. We observe a good behaviour of the time evolution of the total energy which is preserved up to 0.05% when α = 0.5 (this is comparable to the results in [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF]). Clearly, the smaller the perturbation α, the smaller the error in the total energy. Then, we are interested in the time evolution of the logarithm of the L 2 -norm of E. When α is small, by doing a linear approximation to the Vlasov-Poisson equation one may compute an analytical solution approximating the electric field in (8) (see [START_REF] Sonnendrücker | Approximation numérique des équations de Vlasov-Maxwell, notes de cours de M2[END_REF]). More precisely, after a short time, we have a very good approximation of E 1 : E 1 (x 1 , t) ≈ 4 α 0.3677 e -0.1533 t sin(kx 1 ) cos(1.4156 t -0.536245).

Thus, we can compare the numerical results with the analytical approximation by computing the logarithm of the electric field's L 2 -norm. We have done tests with α = 0.5, α = 0.1, α = 0.01 and α = 0.005 (see Fig. 2). For the smallest values of α (0.01 and 0.005, the so-called weak Landau damping test case) we tested that increasing the number of particles per cell makes the computed electric field energy match on longer runs the analytical linear approximation. In the case of the strong Landau damping, when the linear approximation is no longer valid, we obtained when α = 0.5 (see Fig. 2) the two benchmarked rates for the initial damping and for the growth between times t = 20 and t = 40 (see [START_REF] Campos-Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF] and the references therein). As for the case of α = 0.1 we can remark that the analytical approximation is valid only until time t ≈ 15. A drift in the oscillation frequency with respect to the linear approximation appears after this time. Finally, we need to remark that the small values of α pose a challenge to the PIC simulation, as they need a large number of particles in order to achieve a low-noise electric field energy in the long run.

Code optimization and performance

For the sequential implementation of this code we have focused on the same general type of data structures and algorithms as introduced in the work by Bowers (see [START_REF] Bowers | Accelerating a Particle-in-Cell Simulation Using a Hybrid Counting Sort[END_REF] and [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF]) and also implemented in other codes (see [START_REF] Tskhakaya | Optimization of PIC codes by improved memory management[END_REF]). This approach is tailored to the modern commodity CPU architecture as it attempts to minimize and order the memory access as its guiding principle. This acknowledges that memory access is relatively expensive compared with computations. Thus, a minimization of the number of memory accesses and furthermore, an increase of the amount of sequential memory access results in very large performance gains. These principles are accomplished through several key design points:

-particle lists organized as arrays of structures: the array of structures obeys the principle of "data that are used together should be stored together"; all the particle information needed to represent a particle can be thus accessed in a single memory stream. Contrast this with an approach based on a structure of arrays in which the information of position, velocity, etc. are stored independently in their own arrays. The explicit representation of the particle's cell enables us to implement an efficient particle loop, but this can only be seen when considering all of the following items together.

-reduced memory footprint of the particle structure: the choice of representing the particle coordinates as a combination of the cell index plus the local (normalized) offset within the cell, as opposed to using the particle's global coordinates permits the reduction of the memory needed for each particle. Furthermore this representation reduces the amount of arithmetical operations needed during the interpolation and charge deposition procedures. There are other advantages to this approach related to the flexibility of the meshes that can be modeled, but these advantages were not exploited in this work.

-redundant, cell-based data structures to store field and deposited charge values: the field quantities that need to be accessed during the particle push are stored on a cell-basis (all the field values at the corners of the cell are stored in a single data structure which is then laid out as a single array). Together with the data structure used to represent the particle (see below), this maximizes the amount of memory that is accessed sequentially. What makes these data structures redundant is the fact that the field values for a single point must be stored in multiple cell-based structures, but the cost of accessing more memory is offset by ensuring that this access is done sequentially. Furthermore, the execution time is dominated by the size of the particle arrays, not the field arrays.

-use of particle sorting techniques to maximize sequential memory access: With the previously described items, all the information needed to push the particles is stored in three different arrays in memory (one for particles, two for fields), the field information is indexed by cell number and the particle's cell is explicitly represented. If only the particles are sorted by cell number, this enables us to stream through these arrays in sequential form. Hence the use of a particle sorting algorithm. The push therefore looks like this: take a particle, read its cell, summon the necessary field data and complete the push. Read the next particle, read its cell (most likely the same as the previous particle), summon the field data (most likely already in the cache due to the previous particle, etc. As the simulation evolves, the particles get unsorted, so this algorithm needs to be applied periodically, but far from every time step.

-integrated inner loop: the main loop is responsible for advancing the particles, checking if the particle is out of bounds or needs further post-processing, and depositing the charge on the grid. All of these operations are condensed in a single pass over the particle list.

One necessary step, the conversion of the charge data structue into a two-dimensional array usable by the Poisson solver, is done outside of the inner loop. At least one more memory stream is needed in the more general version of this approach, to store references to particles which may need postprocessing outside of the main loop.

The above general features are enough to implement an optimal PIC loop from the perspective of memory access patterns. As mentioned above, during the natural progression of the simulation, the particles fall out of their cell-based ordering, thus yielding some performance degradation over time. For this reason some periodic re-sorting of the particle lists needs to be applied. In the tests reported here, the sorting frequency has been chosen to maximize the overall particle throughput. Now, we describe the parallelization issues. First, by means of a distributed memory paradigm we have followed a simple strategy that permits scaling up to a few nodes: parallelization of the particle list processing while using interprocess communications to gather the charge density information and then redundantly computing the electric fields over the whole domain. No domain decomposition is involved. In other words, each process takes responsibility of the whole life-cycle of a particle list with particles defined over the whole region: from creation, initialization and time-stepping. Only at the time where the electric fields are needed the processes share their individual contributions to the total charge distribution. Second, within each process, we use the multithreading for assigning different segments of the particle list to the different threads. More precisely, we parallelize the loops over particles for the pusher in time, the charge accumulation, and the interpolation of the electric field. Then, we take into account the thread/core affinity control features of OpenMP to maximize the chances that a given core will access the memory that is closest to it. Thus, management of thread/core affinity through Intel's KMP AFFINITY environment variable (= compact) yielded clear benefits. In addition, this setting yielded uniform run times on multiple runs on 16 cores.

We run 2 MPI processes, each on a different socket and with a number of threads from 1 to 8. This leads to simulations using from 1 to 16 threads in all, but distributed over 2 processes as described above. In order to assess the efficiency of our code we performed simulations with three data sizes: 2 million, 20 million, and 200 million particles. Thus, the data (particles and fields) size is bigger than the L3 cache size (thus the memory access speed is determined by the access speed to the main memory). The reference grid was set at 512 x 16 cells, the time step to ∆t = 0.01 and the number of iterations to 1000, when α = 0.1.

The results of using thread affinity are illustrated in Fig. 3 and Table 1 where the best run times among several runs are given. The best speedup of 9.43 was reached on 16 cores for the biggest data size. As expected, Fig. 3 (at right) shows the linear growth of the run time with respect to the increase of the data size. This occurs when using one or several cores. It is interesting to get an idea of the performance of the memory subsystem in these simulations. We have computed at each time step the obtained memory bandwidth during the particle push. This number is N • f /T , where N is the number of particle loads and stores, f the memory-size of one particle (in our case, 32 bytes), and T the time to process the particle list. In Fig. 4 we illustrate the time evolution of this measure in GB per second when running on a single CPU. We can see that the highest memory bandwidth speeds are obtained when particles are sorted, as expected. In Fig. 5 we observe that, for two data sizes, similar values were obtained, when using one or several cores.

Maybe the most interesting measure of the performance of the whole simulation (not only remark that the value of the relative error for ε = 0.5 is much larger than those for smaller ε. This illustrates that the guiding center model is not a good approximation anymore for the Vlasov-Poisson system when ε = 0.5. A more precise view of these facts can be observed in Fig. 8, where cuts of the densities are presented at 3 values of the x 1 direction and at two different final times, t = 12 and t = 22. Then, going further, in Fig. 9 we show contours of the densities for ε = 0.1 and ε = 0.5. Other simulations for ε = 0.05 like those in Fig. 9 (but not shown here) display a qualitatively very good match between the Vlasov-Poisson reference solution and the guiding center solution. Thus, our simulations show that smaller values of ε make the Vlasov-Poisson solution closer to the guiding center solution. We finally remark that the guiding center numerical solution is a good approximation of the Vlasov-Poisson numerical solution when ε = 0.1.

Conclusion and outlook

The aim of this paper was to efficiently implement a Particle-In-Cell method for solving some nonlinear Vlasov-like models in a four-dimensional phase space. The organization of the code follows an existing work (see [START_REF] Bowers | Accelerating a Particle-in-Cell Simulation Using a Hybrid Counting Sort[END_REF]) and mainly focuses on optimization of the memory access. We demonstrated the performance of our code within some numerical experiments on the Landau damping test case. We also verified the code on a simulation of the Kelvin-Helmholtz instability. In addition, this test case allowed us to numerically illustrate some convergence result from the Vlasov-Poisson model with a strong magnetic field towards the guiding center model.

In future work, we intend to implement the time stepping scheme introduced in [8]. This method allows to solve multiscale Vlasov-Poisson models by treating the stiff terms as an exponential integrator. In the framework of an optimized PIC implementation, this scheme potentially provides a solution to perform efficient simulations of highly oscillating six-dimensional models without solving the smallest scale.

Acknowledgment

The first author's work was supported by the project ELI: Extreme Light Infrastructure (CZ.02.1.01/0.0/0.0/15 008/0000162) from European Regional Development.

Figure 1 :

 1 Figure 1: Time evolution of the relative error on the total energy for three values of the perturbations. At the left α = 0.5, in the middle α = 0.1, and at the right α = 0.005.

Figure 2 :

 2 Figure 2: Time evolution of the log of the electric field energy for four values of α. The time step is 0.1 and the grid has 128 x 2 cells.When α = 0.5, the slopes for the initial damping and for the growth between t = 20 and t = 40 are -0.2920 and 0.0815 respectively (see[START_REF] Campos-Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF]).

Figure 3 :

 3 Figure 3: Speedup and best times by core number for simulations with 3 data sizes

Figure 4 :

 4 Figure 4: History of the instantaneous memory bandwidth in GB per second (as measured during the particle push loop only) for the 2 million particles run on 1 core. At right a zoom at the beginning of the simulation. In green the steps where sorting is done.

Figure 6 :

 6 Figure 6: The time evolutions of the relative errors in the electric field energy (left) and in the enstrophy (right) for the Kelvin-Helmholtz instability.

Figure 7 :

 7 Figure 7: The global relative error at time t = 5 of the density of the Vlasov-Poisson solution with respect to the guiding center solution.

Figure 8 :

 8 Figure 8: The Kelvin-Helmholtz instability with 10 million particles: Cuts in different x 1 directions of the densities of the guiding center model (in black), the Vlasov-Poisson model with ε = 0.05 (in blue), with ε = 0.1 (in red), and with ε = 0.5 (in green). The final times are 12 (at left) and 22 (at right).

Figure 9 :

 9 Figure 9: The Kelvin-Helmholtz instability with 10 million particles: contour plots of the densities given by the guiding center model (top row), by the Vlasov-Poisson model with ε = 0.1 (middle row) and ε = 0.5 (bottom row), at times 12 (left) and 22 (right).

Table 1 :

 1 Best speedup and best times by thread number for the simulation with 20 million particles, 1000 iterations and a grid of 512 x 16 cells.

	Cores Time (seconds) Speedup
	1	703.2	1.0
	2	348.2	2.02
	4	191.5	3.67
	6	139.2	5.05
	8	112.9	6.23
	12	87.7	8.02
	16	75.3	9.34

Table 2

 2 for the three considered data sizes. Thus, we achieved a simulation with 265 million particles advanced per second on a 2-socket node with 16 cores, which is a very satisfactory result.

	Simulation \ Cores	1	2	4	8	16
	2 million particles	28.5 56.8 101.7 165.2 224.3
	20 million particles	28.4 57.4 104.4 177.1 265.6
	200 million particles 28.0 56.9 103.5 174.4 263.8

Table 2 :

 2 Millions of particles processed per second for 3 runs on a single or several cores of 1000 iterations with a grid of 512 x 16 cells.

http://selalib.gforge.inria.fr

Figure 5: History of the instantaneous memory bandwidth in GB per second (as measured during the particle push loop only) for the 2 million (left) and 20 million (right) particles runs, when using several cores.

Vlasov-Poisson vs. guiding center simulations

The aim of this section is to compare a reference solution of the Vlasov-Poisson system in (1) when ε decreases to a reference solution of the guiding center model in [START_REF] Bostan | The Vlasov-Maxwell System with Strong Initial Magnetic Field: Guiding-Center Approximation[END_REF]. For both systems we implement standard PIC methods as described previously. Next we discuss the numerical parameters and the results.

We consider the initial condition for a Kelvin-Helmholtz instability (see [START_REF] Sonnendrücker | Approximation numérique des équations de Vlasov-Maxwell, notes de cours de M2[END_REF])

defined in Ω x × R 2 , where Ω x = [0; 2π/k] × [0; 2π] and k = 0.5. For the guiding center model, the initial condition consists obviously in the x part in [START_REF] Tskhakaya | Optimization of PIC codes by improved memory management[END_REF], defined in Ω x . First, for the numerical solution of the guiding center system, we push the particles in time with an explicit second order Runge-Kutta scheme. The time step is ∆t = 0.01. We put 10 million particles and a cartesian grid with 256 x 128 cells. The energy (L 2 -norm of E) and the enstrophy (L 2 -norm of f GC ) are known to be conserved quantities in time. Therefore, in order to verify the code, we compute these values. In Fig. 6 we plot in time

and the same quantity for the enstrophy. We observe a good conservation of these values (this is comparable to the results in [START_REF] Christlieb | A High Order Time Splitting Method Based on Integral Deferred Correction for Semi-Lagrangian Vlasov Simulations[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF]).

Second, we illustrate the convergence result mentioned in Section 1.1. With this aim, we recall that we need to investigate long-term simulations (of order 1) since the fastest scale of the system (1) is of order ε 2 . Thus, we compute a reference solution for the Vlasov-Poisson model by using the leap-frog scheme in [START_REF] Campos-Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF] with the time step dt = 2πε 2 /80, in order to resolve the time period of the solution. The same numerical parameters as for the guiding center simulation are used. We recall that periodic boundary conditions are used in physical space. Now, for a quantitative verification of the convergence result we compute some norm in space of the densities. Thus, we first plot in Fig. 7, as a function of ε, the maximum in time until t = 5 of the quantity ||f GC -ρ ε || L 2 /||f GC || L 2 where || • || L 2 stands for the L 2 norm on Ω x . We