Mouna Tka Mnad

Christophe Deleuze

Ioannis Parissis

Synchronous Programs Testing Language (SPTL)

SPTL is a language designed to test applications developed for synchronous controllers. It makes possible to provide a specification of the software external environment. This specification can then be processed to generate test input sequences guided by directives such as profiles of use and scenarios. We introduce a definition and an overview of the language through a simple example of a reactive system that we present in this paper.

Introduction

A "reactive system" [START_REF] Harel | Logics and models of concurrent systems[END_REF] is a system that continuously interacts with its environment. Synchronous systems [START_REF] Benveniste | The synchronous approach to reactive and real-time systems[END_REF] is a class of reactive systems that are based on the synchrony hypothesis: they are supposed to be fast enough to take into account any external event. They are widely used in industry for the automation of various processes.

Several programming languages have been defined in the past for developing synchronous systems such as Esterel, Signal or Lustre [START_REF] Halbwachs | The synchronous data flow programming language Lustre[END_REF]. In this research work, we are especially interested in a specific class of synchronous controllers, produced by a leading manufacturer of automation control components. These controllers are used in applications such as heating, air conditioning, access control, water and air treatment, waste treatment or pump management. Unlike usual synchronous systems, these devices are intended to be programmed by a large scope of users, including non-experts. A user-friendly graphical programming environment is used to build and integrate software applications as data-flow networks made of already developed components.

Even if the applications defined by such non-expert users are not highly critical, it is however necessary to adequately test them. The objective of this work is to design an efficient and easy to use environment for testing programs developed with this environment. Users must be able, as they write programs, to express various test scenarios that will be used for testing purposes. These scenarios must be compiled into test data generators that will be used to simulate and test the devices.

In this paper, we present the first component of such a testing environment, an adequate and easy to use test description language.

This document is structured in 6 sections. In the following section, we present related work and the motivations for our research while in the third section we define the main concepts of the testing language. A overview of the language through examples is provided in section 4. In section 5, we describe the semantics of the language.

Related Work

Automatically testing synchronous systems has been addressed in the past in several research works. We focus on investigations on testing Lustre programs [START_REF] Halbwachs | The synchronous data flow programming language Lustre[END_REF], because of the similarity between this synchronous data-flow language and the programming environment of our synchronous controllers.

-Lutess [START_REF] Seljimi | Using clp to automatically generate test sequences for synchronous programs with numeric inputs and outputs[END_REF] has been designed as a black-box testing environment that automatically generates test input sequences from a formal specification of the software external behavior provided in a language similar to Lustre. Lutess requires three components: the software environment description (1), the executable code of the system under test (2), and a test oracle (3) describing the system requirements (see Figure 1). The system under test and the oracle are both synchronous executable programs. The environment description is composed of a set of constraints, provided in a Lustre-like language, that the environment of the software application is assumed to satisfy. Lutess translates the above constraints into a random test generator operating on a single action-reaction cycle. The generator randomly selects and sends a valid input vector to the program under test which reacts with an output vector and feeds back the generator with it. The generator proceeds by producing a new input vector and the cycle is repeated. It must be noted that Lutess constraint resolution mechanism is implemented using constraint logic programming (CLP). A case study using Lutess to test a steam-boiler controller is described in [START_REF] Papailiopoulou | Revisiting the Steam-Boiler Case Study with LUTESS : Modeling for Automatic Test Generation[END_REF].

-Lurette [START_REF] Jahier | Case studies with lurette v2[END_REF] is an automatic test generator for Lustre programs the main concepts of which are similar to Lutess: a specification of the environment behavior is provided and translated into a test generator. The specification language is ad hoc and the constraint resolution mechanism is based on convex polyhedra. -GATeL [START_REF] Marre | Test selection strategies for lustre descriptions in gatel[END_REF] TTCN-3 as a modeling language can support test generation tools by means of some annotations (see for example [START_REF] Wu-Hen-Chang | A new approach in model-based testing: Designing test models in ttcn-3[END_REF]). It applies to a variety of application domains and types of testing. TTCN-3 is a popular test definition language mostly used by experts for communications protocols testing.

The above tools and approaches are either research prototypes or industry tools intended to be used by domain and test experts. Lutess and Lurette use rather complex test modeling languages making the specification of test scenarios difficult while GATeL requires the user to determine a test objective and to assist the tool during the constraint resolution. Similarly, TTCN-3 seems too complex for non-expert users and is not designed for synchronous applications. However, we use the main concepts of these tools to define an easy-to-learn test modeling language, described in the next sections.

Basic Concepts

To make the new testing language easy to use, we define four simple concepts, defined in the next subsections: profiles, categories and groups, and scenarios.

Profile

A system performance and functioning are significantly dependent on the environment in which it operates. A "profile" represents a context in which the system is used. ZigBee [START_REF] Zucatto | Zigbee for building control wireless sensor networks[END_REF], for instance, uses profiles defined as message formats and processing actions to allow devices to exchange data in a given application domain. Profiles are developed by companies to help providing specific needs. For example: commercial building automation, telecom applications, hospital care ... More generally, the notion of "operational profile" has been introduced by Musa [START_REF] Musa | Operational profiles in software-reliability engineering[END_REF] and defined as a quantitative characterization of how the system will be used. Let us see an example of program to illustrate the notion of profile and its utility in our case: an air conditioner. It will not be used in Tunisia in a small house in winter the same way it will be used in Canada in a big building during the summer, because temperatures and time slot of use will be very different. So, to every profile correspond constraints of use that will necessarily affect testing.

Categories and Groups

Profiles are composed of categories of constraints provided by the tester. As shown in Figure 2, each category contains groups of constraints. The constraints in each group define limits related to a particular environment or users.

Fig. 2. Categories

A profile is obtained by combining groups from different categories. For example in Figure 3 profile 1 is composed of Group1,1, Group3,2 and Group2,3: Profile1 = Group1,1 ∪ Group3,2 ∪ Group 2,3.

Scenarios

When testing reactive systems, we need to express how they are constrained by the context and the user activity. A direct approach is to explicitly provide typical and significant user activities in the testing process. Such descriptions, often called "scenarios", support reasoning about situations of use. A scenario can be seen as a story [START_REF] Carroll | Five reasons for scenario-based design[END_REF]. For example, at the beginning all the buttons of a lighting system are not pressed then the user presses rapidly button1, then presses button2, then keeps pressing button1 until the permanent lighting mode is on, etc.

Overview of SPTL

In this section we present the syntax of SPTL using an example of a reactive system.

Example: Air Conditioner

Let's consider the following simple reactive system, an air-conditioner controller, shown in Figure 4. -IsOn (Boolean): indicates the state of the air-conditioner; -Tout (integer): indicates the temperature of the air that the air-conditioner blows.

Variables

The first part of the test model is a set of variable declarations: the system inputs and outputs are declared as input and output variables. A variable declaration associates a type with an identifier and optionally an initial value. Each identifier is unique. The basic types that are considered are boolean, integer, and time. The word "random" means that the value of the variable will be randomly generated, if there is no other specification about it in the program.

var

) } }

A constraint involves one or more variables, and restricts values that these variables can take simultaneously. These relationships are given by the environment of the system and the user profile. Expressions are defined by combining sub-expressions using arithmetic, logical and temporal operators as in Tuser = pre(Tamb) -3. pre is a temporal operator used to get the value of the variable in the previous cycle. For this to be meaningful Tamb must have a value before the first cycle. Such an initial value is to be specified in the variable declaration. Every output variable must appear within the pre operator since in the calculation of the test values we consider the system outputs generated in the previous cycle. The prob operator is used to define probabilities, useful to guide test data selection. The expression prob(v,p) has the value v with probability p, otherwise it takes a random value different from v.

Scenarios Specification

A scenario is composed of point-wise steps and interval steps. A point-wise step, enclosed between "{" and "}", is a set of constraints to obey in one test generation time. It can be used either in the first step for initializing environment variables or in any step of the execution. An interval step, enclosed between "[" and "]", is a set of constraints on the test values to generate that remains valid until a condition holds. The condition is a boolean expression enclosed between "(" and ")". Below we describe two example scenarios. This first scenario means that Tuser is initialized to 8. Then this temperature will be increased until reaching 10 • C. After this Tuser is set to 12 • C . Finally, for the rest of the test data that will be generated Tuser will be increased by 2 at each cycle until it reaches 22 • C.

Table 1 shows a trace that could be generated from scenario 1. This trace can be examined to check for expected properties from the air-conditioner such as:

-IsOn t ∧ (Tamb t < Tuser t ) ⇒ Tout t > Tamb t : when it's cold, the air conditioner heats and conversely. -OnOff t ⇒ IsOn t-1 = IsOn t : the OnOff button turns on / off the air conditioner.

Table 1. An extract of execution traces of scenario 1 t0 t1 t2 t3 t4 t5 t6 t7 t8

OnOff 0 1 0 0 0 0 0 0 0 Tamb 13 14 13 12 11 [START_REF] Wu-Hen-Chang | A new approach in model-based testing: Designing test models in ttcn-3[END_REF] This second scenario shows the use of a variable (t) of type time. Time values are positive numbers followed by a unit identifier (s for seconds, ms for milliseconds, c for cycles). (A cycle in the tested program is the necessary time for a full run from starting entries, calculating outputs' values, to finally generating them). t.start is a pseudo constraint having the side effect to start counting time in t. In the example t is used in the condition of the interval step. Table 2 represents a possible trace for scenario 2.

Semantics of the Language

In this section, we give a description of the language semantics. In order to precisely describe how test sequences are generated, we start by defining the semantics of three basic SPTL constructs: expressions, constraints and scenarios.

Expressions

The value of an expression depends on the values of the variables it contains. We call environment and note σ a mapping from variables to their value. Thus σ(x) is the value of variable x in environment σ.

We note e the semantic function of expression e. The (denotational) semantics of most expressions is rather straightforward as shown in Figure 5: the semantics of a constant value is this very value, whatever the environment; the semantics of a variable expression is the variable value in the considered environment; arithmetic and boolean operators map to the obvious mathematical operators. The probabilistic operator introduces non determinism: the expression prob(v,p) has the value v with probability p, or a value w randomly chosen among the other possible values in the type.

cst (σ) = cst id (σ) = σ(id) e1 + e2 (σ) = e1 (σ) + e2 (σ) e1 and e2 (σ) = e1 (σ) ∧ e2 (σ) ... prob(v, p) (σ) = v with
The temporal operator pre allows to refer to the value of a variable (or more generally an expression) at the previous cycle. To define it we extend the environment to include memory variables: such a variable v m always has the value variable v had at the previous cycle. As noted in section 4, initial values must be provided for variables appearing in a pre so that memory variables can be initialized. We show in section 5.4 how these variables are updated.

The set of memory variables is finite and can be built from a static analysis of the model. Let's start with an empty set φ m . At each occurrence of the pre operator in the model, we consider its expression argument e and add v m to φ m for each variable v occurring in e. Nested pre expressions can also easily be dealt with by adding a "second order" memory variable (e.g. v mm for variable v m .)

We define prev as the function mapping an environment σ 1 to an environment σ 2 where each variable v has the value of the corresponding memory variable v m in σ 1 if it exists. We can now state that the semantics of pre(e) in environment σ is the semantics of e in environment prev(σ).

Constraints

When considering constraints we need to distinguish the several kinds of variables that exist in an SPTL model:

output variables are set by the system under test, -input variables are to be computed by the SPTL model, -memory variables remember the value of variables in the previous cycle.

Since output variables can only appear in a pre operator, a constraint only contains occurrences of input and memory variables. 2 It has the form of a boolean expression but its semantics is very different: a constraint defines the possible values of its input variables, according to the values of its memory variables.

Let's split the environment σ into the environment of memory variables σ m and the environment of input variables σ i . We note σ = σ m ⊕ σ i . It is clear that a constraint defines the possible values of σ i from a given value of σ m . Noting c the semantic function of constraint c, Figure 6 shows these possible σ i are those that, when associated with σ m , form a complete environment σ in which the boolean expression has the value true. Also shown on the figure is the rule for constraint composition: the set of possible values for the input variables is the intersection of the sets of possibles values for each constraint.

Scenarios

Recall that a scenario is a sequence of point-wise steps: they specify constraints that apply at one cycle; -interval steps: they specify constraints that apply until some exit condition is verified.

We can describe a scenario by an automaton where each scenario step is a state. A point-wise step has a single true-labeled transition to the next state, while an interval step has a transition to the next state, labeled by its exit condition, and a transition to itself, labeled by the negation of its exit condition.

The scenario 1 in section 4.4 would thus be described by the automaton of Figure 7 where states 1 to 4 hold the sets of constraints of the four scenario steps.

Programs

Recall that a set of profiles have been generated as combinations of groups from the categories. To generate a test sequence we need to consider:

variable declarations, -a profile, -a scenario.

Let's call this triple an SPTL program (as opposed to the whole model ) since that's what will be used for generating a test sequence. Figure 8 shows an informal operational description of the semantics of such a program. The basic idea is that at each cycle a current set of constraints cset is in force (step 1): it is the combination of the profile constraints and the set of constraints of the currently active step in the scenario (scen step).

Solving cset in environment σ m , possible values for the input variables are computed. One of these values σ i is randomly selected (2) and fed to the inputs of the system under test [START_REF] Halbwachs | The synchronous data flow programming language Lustre[END_REF]. System outputs are read (4), an oracle can optionally be given current inputs and outputs (4b), scenario active step is changed (5) if its exit condition is true (evaluated in the environment σ = σ m ⊕ σ i ). 3Memory variables are updated [START_REF] Marre | Test selection strategies for lustre descriptions in gatel[END_REF]. The test terminates if scenario ends, else next cycle begins [START_REF] Musa | Operational profiles in software-reliability engineering[END_REF] Initially, scen step is set to 1, the first step of the scenario. Initial values for the memory variables are provided with variable declarations (0).

For clarity, we haven't included the handling of time variables so far. A time variable counts a number of executed cycles, thus it holds an integer value initialized at 0. The program state includes an initially false internal boolean variable t active for each such time variable. Solving the pseudo constraint t.start has the effect of setting t active to true. At end of cycle all time variables whose associated t active boolean variable is true are incremented. Time variables can not appear in a pre.

Conclusion and Future Work

We have proposed a language to write models to generate test sequences for synchronous reactive systems. It is based on describing both the constraints on the system environment and the profile of the users of the program under test. One of the main design objectives is to allow non testing experts to make realistic tests. This has led us to design a language with few simple concepts.

Moreover, the notion of categories, groups and profiles will reduce the work needed to prepare realistic tests. Indeed, some usual profiles may be predefined, by the manufacturer's engineers or by a third party, to further assist the user in its testing design. Test scenarios make possible to take test objectives into account during test generation.

For the next step in our work, we will focus on the implementation of a prototype of test data generator based on models written in the SPTL language and on its evaluation on case studies. We intend to extend the language to include test oracles, for comparing the actual system outputs with the expected ones, and for observing functionalities and properties of the program under test.

Fig. 1 .

 1 Fig. 1. The Lutess environment

Fig. 3 .

 3 Fig. 3. Profiles

Fig. 4 .

 4 Fig. 4. Air conditioner

Scenario 1 6

 1 ≤ Tamb ; Tamb ≤ 14 ; OnOff = prob(true, 0.1 ) ; begin { Tuser = 8 }| [ Tuser = pre(Tuser) +1 ( Tuser = 10 ) ]| { Tuser = 12 }| [ Tuser = pre(Tuser) + 2 ( Tuser = 22 ) ] end

  Tuser=2; t.start } | [ Tuser=pre(Tuser)+3; Tamb=10 (t = 3c) ] | [ Tuser=pre(Tout)-1; Tamb=pre(Tamb)-1 (Tamb=6) ] | { Tuser=10; Tamb=6 } | end

Table 2 .OnOff

 2 An extract of execution traces of scenario 2 t0 t1 t2 t3 t4 t5 t6 t7 t8

probability p w with probability 1 -Fig. 5 .

 15 Fig. 5. Semantics of expressions

cFig. 6 .

 6 Fig. 6. Semantics of constraints

Fig. 7 .

 7 Fig. 7. Automaton for scenario 1

Fig. 8 .

 8 Fig. 8. Execution cycle of an SPTL program

  A subprogram, like ComputeTemperature in the example, connects an input with a set of outputs together and admits intermediate variables. The declaration of a subprogram contains a header, an optional block for variables declaration and the body of the subprogram that contains the equations.

	input bool OnOff;
	input random int Tamb = 10;
	input int Tuser = 7;
	output bool IsOn;
	output int Tout;
	4.3 Categories Specification
	categ CountrySeason
	{ group FranceSummer
	{ 20<Tamb ; Tamb<44 ;
	Tuser = ComputeTemperature(Tamb) }
	group TunisiaWinter
	{ 6 <= Tamb ; Tamb <= 14 }
	sp ComputeTemperature(int Tamb) returns(int Tuser)
	{ Tuser = pre(Tamb) -3 }
	}
	categ TypePlace
	{
	group Company { OnOff=prob(true, 0.001) }
	group House	{ OnOff=prob(true, 0.1

Each category contains a set of groups of constraints to be respected and possibly sub-programs used for the calculation of values of certain variables. Thereafter, a set of profiles can be generated combining groups of different categories. The tester can then choose a profile consisting of multiple constraints (environment, use ...), modify it or add test cases.

  . cset ← profile ∪ scenario[scen step].constraints 2. compute σi by solving cset (using σm) 3. apply σi to system inputs 4. read system outputs to σo (4b) feed oracle with σi and σo 5. evaluate exit condition from scenario step, increment scen step if true 6. update σm (from σi and σo) 7. if end of scenario then stop else goto[START_REF] Benveniste | The synchronous approach to reactive and real-time systems[END_REF]

	(0) scen step ← 1
	initialize σm
	1.

http://www.ttcn-3.org/

Of course, memory variables do not appear explicitly in the model. We call memory variable any occurrence of a variable in a pre expression.

A point-wise step can be defined as an interval step with a true-labeled exit condition.


Hervé Guillard
email: herve.guillard@inria.fr

Mapped Fourier Methods for stiff problems in toroidal geometry

Keywords: Spectral methods, Fourier approximation, mapped Fourier method, pellet injection, Tokamaks Méthodes Spectrales, Approximation de Fourier, méthode de Fourier transformée, injection de glaçon, Tokamaks

Fourier spectral or pseudo-spectral methods are usually extremely efficient for periodic problems. However this efficiency is lost if the solutions have zones of rapid variations or internal layers. For these cases, a large number of Fourier modes are required and this makes the Fourier method unpractical in many cases. This work investigates the use of mapped Fourier method as a way to circumvent this problem. Mapped Fourier method uses instead of the usual Fourier interpolant the composition of the Fourier interpolant with a mapping in such a way that in the computational space, the functions to represent are not stiff. This work gives some examples of the usefulness of this method and apply it to a simple model of pellet injection in tokamaks as an example of its potential interest for complex multi-dimensional problem.

Introduction

Many problems in science and engineering involve one or several periodicity directions. In these directions, the numerical approximation of the mathematical models can be performed in an efficient way using Fourier decomposition of the solutions in term of trigonometric polynomials. For infinitely smooth solution, the convergence of Fourier methods is more rapid than any finite power of 1/N where N is the number of Fourier modes : this property is called spectral accuracy and is one of the main advantage of Fourier methods versus finite difference (FD) or finite element (FE) local approximation methods since it allows to reach a very good accuracy using a small number of modes. Moreover, the existence of fast algorithms (Fast Fourier Transform) commonly available on vector and parallel computers allows an efficient implementation of these method in term of CPU cost. Therefore in different branches of science and engineering, a large number of simulation codes rely on the use of Fourier techniques. This is for instance the case in astrophysics or in plasma physics and tokamak research where mature simulation codes designed, validated and tuned against experimental and observational results are used in a routine or production way in many different teams worldwide. However, all these advantages of Fourier methods break down if the solution possesses localized internal layers and region of rapid changes which are small in comparison of the computational domain. For these cases, even if the solution is infinitely smooth 1 , spectral accuracy is only recovered for a large number of Fourier modes, making the use of these methods unpractical. To simulate these cases, one possibility is to use local FD/FE methods. However, this can imply a large re-coding effort that is not always possible due to time and money constraints. Another possibility is to alter the physics of the problem by enlarging artificially the region of large variations or modifying the governing parameters of the problem to obtain a less stiff problem. For complex problem governed by non linear partial differential equation, the correctness of the results obtained in this way is not guaranteed. In this work, we experiment another possibility whose advantage is to not require large modification of existing codes. Only the parts of the code performing the computations of the spatial derivatives in the periodicity directions need to be changed and moreover these modifications can be light (depending of the structure of the existing code however). The basic idea is to change the global trigonometric interpolant by functions much more adapted to the representation of functions exhibiting rapid internal variations. This is done by composing the Fourier modes with an appropriate change of variables. Another way to present the same idea is to write the initial problem in a computational domain where the solution is not stiff. This idea to use a change of variable to map a physical domain to a periodic one for computational purpose is not new. The term mapped Fourier methods to designate these methods have been used in quantum computations studies [START_REF] Fattal | Phase space approach for optimizing grid representations: The mapped fourier method[END_REF] almost 20 years ago and it is for instance discussed at length in the well-known Boyd's text book [START_REF] John | Chebyshev and Fourier spectral methods[END_REF]. In the cases where the localization of the internal zones of rapid variations is not know but is part of the problem, algorithms to construct an adapted mapping from the knowledge of the current solution or of the residual have been designed see [START_REF] Guillard | On the use of spectral methods for stiff problems[END_REF] and [START_REF] Bayliss | An adaptive pseudo-spectral method for reaction diffusion problems[END_REF]. For some theoretical development on the convergence of these type of methods in a non-periodic context, one can consult [START_REF] Shen | Error analysis for mapped legendre spectral and pseudospectral methods[END_REF][START_REF] Shen | Error analysis for mapped jacobi spectral methods[END_REF]. The structure of this paper is as follows. In section 2, we give some examples of the behavior of Fourier methods when the solutions exhibit localized internal layers to emphasize the loss of fast convergence of Fourier methods in these problems. Section 3 recall the principle of mapped Fourier methods and revisits the examples given in section 2 to show the potentialities of the method. Section 4 apply this technique for the simulation of a crude model of pellet injection, a technique used to refuel the core plasma in tokamaks. Finally we conclude by some remarks.

Spectral accuracy and spatial resolution 2.1 Accuracy of the numerical derivatives

We begin with the very simple problem of computing the successive derivatives of a function possessing a region of rapid variation. The test function considered is the Maxwellian :

u(x) = e -(x-π) 2 ∆ 2
Computing the derivatives of this function by a Fourier method is done by sampling u on a set of discrete points x j = 2π/N, j = 0, . . . , N -1 and replacing u by its Fourier interpolant I N u that is a member of S N the set of trigonometric polynomials of degree N . The Fourier interpolant is defined by

I N u = N/2-1 k=-N/2 ûk e ikx
where the coefficients û are computed by the Discrete Fourier Transform (DFT):

ûk = 1 N N -1 j=0 u(x j )e -ikxj
Then computing the pseudo-spectral derivative of u consists in computing the exact derivative of the Fourier interpolant I N u. Fig. 1 compares the exact first derivative of u with the first derivative of I N u while Fig. 2 plots the results for the second derivative. It is seen that if the spatial resolution is not sufficient and not able to capture the zone of rapid variation of the function, then the results are extremely inaccurate and that the Fourier interpolant (the green line) is an extremely oscillatory function. This phenomena (the Gibbs phenomena) is well known and studied in details for instance in [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF]. Fig. 3 that displays the discrete L2 norm of the error vs ∆ is another way to consider the same results. It shows that either the resolution is sufficient and spectral methods are extremely accurate, either the spatial resolution is not sufficient and the accuracy of the methods methods is extremely bad. Spectral methods have thus a "binary" character : They are then either very good or very bad. for ∆ = 0.1

Laplace equation

Here we consider in some sense the inverse problem of the one considered in section 2.1 and solve by a Fourier Spectral method the problem :

-u xx = f (1)
The right hand side f is chosen such that the exact solution of this problem is u = e -(x-π) 2 ∆ 2 . In the following numerical experiments, ∆ = 0.1 and the Laplace problem can be considered as moderately stiff. The figure 4 displays the trigonometric interpolating polynomial of the right hand side of equation ( 1) sampled respectively on 32, 64 and 128 collocation nodes. The collocation nodes are displayed by blue dots on this figure and it is seen that if the spatial resolution is not sufficient, the trigonometric interpolating polynomial passing through these points is highly oscillatory. For this quite large value of the width parameter ∆, 128 nodes at least are necessary to avoid large oscillations.

Solving 1 by a Fourier pseudo-spectral method on a grid of N nodes is equivalent to solve the problem : Find u N ∈ S N such that where I N stands for the Fourier interpolant. Due to the highly oscillatory character of the trigonometric interpolating polynomial when the spatial resolution is not sufficient, one can expect the error to be extremely large in this case. This is what is shown in figure 5 that displays the solution of 2 for 32, 64 and 128 collocation nodes respectively. The red dots displays the solutions on the collocation points where the solution is effectively computed while the red solid lines display u N the interpolating polynomial passing through these points. For N = 32, the error of the numerical solution is extremely large and the Fourier method is clearly useless in this case. Even for N = 64 where there are 5 collocation points in the width ∆ of the Maxwellian, Inria small oscillations can be noticed in the solution and at least 128 are needed to solve this problem with a sufficient accuracy.

-u N = I N f (2)
3 Mapped Fourier methods

Principle of the method

Consider a partial differential equation on [0, 2π] that we write symbolically as

L(u, ∂u ∂x , ∂ 2 u ∂x 2 , . . . ) = f (3)
The principle of Fourier methods is to look for the solution u of this PDE in the space S N of the trigonometric polynomials of order N . Although there are many variations of spectral methods (Galerkin, Tau, collocation), a common characteristics of Fourier methods is to compute the approximate solution by sampling u (or the RHS f ) on a set of discrete points x j = 2π/N, j = 0, . . . , N -1 and to replace it by its Fourier interpolant

I N u = N/2-1 k=-N/2 ûk T k (x)
where the basis functions T k = e ikx are the Fourier modes and where the coefficients û are computed by the Discrete Fourier Transform (DFT):

ûk = 1 N N -1 j=0 u(x j )e -ikxj
As shown in section 2 a drawback of this type of method is that the Fourier modes are unable to represent a function displaying small localized features when the spatial resolution is not sufficient : either the resolution is sufficient and Fourier methods are extremely efficient and accurate, either the spatial resolution is not sufficient and the accuracy of Fourier methods is extremely low. Since the Fourier modes are not adapted to the representation of stiff functions with localized features, a simple remedy is simply to change the Fourier modes for functions that are better adapted to the characteristics of the functions that are to be computed. This can be simply done by introducing a smooth change of variables h : ξ → x and to consider instead of the Fourier modes T k (x), mapped Fourier modes defined by T k o h -1 (x). The function u(x) is then developed on the basis of the mapped Fourier modes :

I N (h)u = N/2-1 k=-N/2 ûk T k o h -1 (x)
and standard (Galerkin or collocation) methods can be used as with usual methods.

Another way to describe this technique is to consider that to each function u(x) defined on Ω x = [0, 2π] is associated another function U(ξ) defined on a computational space Ω ξ by :

U(ξ) = u(x(ξ)) = u o h(ξ)
The mapping h is chosen to provide an adequate resolution in the computational space Ω ξ or in other words, such that the function U(ξ) is not stiff. The solution u(x) of (3) can then be obtained from the solution of the transformed pde

L ξ (u, ∂U ∂ξ , ∂ 2 u ∂ξ 2 , . . . ) = f ξ (4)
obtained from repeated application of the chain rule to the partial derivatives present in (3). Equation ( 4) is then discretized using a standard Fourier method. This procedure can look difficult to implement as the pde (4) can be considerably more complex than (3) especially if (3) contains high order derivatives. However if a collocation method is used, to solve (3), only the expression of the derivatives at the collocation points x j are needed. Then, let D be the matrix representing the Fourier derivative at the collocation points, i.e

(I N u) (x i ) = D ij u(x j )
and denote u (resp U ) the N-dimensional vector of the values of u at the collocation points x i i.e u i = u(x i ) (resp. the N-dimensional vector of the values of U at the collocation points ξ i i.e U i = U(ξ i ) ) then the application of the chain rule give :

u (x j ) = U (ξ j ) h (ξ j ) = 1 h (ξ j ) (DU ) j = 1 h (ξ j ) (Du) j (5)
since u(x j ) = U(ξ j ) and the vectors U and u are identical. The mapped Fourier technique thus reduces to multiply by a diagonal matrix the result of the standard Fourier derivative. This last operation is usually done without assembling the matrix D and using instead Fast Fourier Transform.

Let us now consider the evaluation of higher order derivatives, one possibility is to evaluate them using the Faà di Bruno formula. For instance, we have

u (x j ) = 1 (h (ξ j )) 3 [h (ξ j )U (ξ j ) -h (ξ j )U (ξ j )] (6)
and U (ξ j ) can be evaluated as U (ξ j ) by FFT. However, an alternate way to evaluate high order derivatives is simply to apply in a repeated way the expression [START_REF] Fattal | Phase space approach for optimizing grid representations: The mapped fourier method[END_REF]. This is the strategy that has been followed in this paper.

In the following two sections, we repeat the numerical experiments of sections 2.1 and 2.2 to demonstrate the benefits of using mapped Fourier methods. The mapping considered is

x = h(ξ) = atan(Ltanξ) (7)
Use of this mapping has been proposed and advocated in [START_REF] Boyd | The arctan/tan and kepler-burgers mappings for periodic solutions with a shock, front, or internal boundary layer[END_REF] that gives good arguments to prefer it against alternate choices in the periodic case.

Accuracy of the numerical derivatives

We repeat here the numerical experiments of the section 2.1 but use for the evaluation of the derivatives the mapped Fourier method. The value of the parameter L in [START_REF] Guillard | On the use of spectral methods for stiff problems[END_REF] has been chosen as L = ∆. Fig. 6 and 7 have to be compared with the figures 1 and 2 and illustrates the improvement that can be expected from the use of the mapped Fourier method. On these figures, note that not only the values of the derivatives on the collocation points are extremely accurate but that also the interpolating function (that is not a trigonometric interpolant but a member of the space S N (h) of the function composing a Fourier interpolant with the mapping h) is not oscillatory and superpose exactly with the exact derivatives. for : left plot ∆ = 0.1, right plot ∆ = 0.01, note that the right plot is a zoom on the central region

Laplace equation

In this section, we again consider the solution of equation [START_REF] Bayliss | An adaptive pseudo-spectral method for reaction diffusion problems[END_REF]. The exact solution is the Maxwellian

u = e -(x-π) 2 ∆ 2
and as ∆ → 0 this solution becomes harder and harder to catch with conventional Fourier methods. However, the use of transformed Fourier series as explained in section 3.1 can cure this problem in an efficient way. We again use the change of variables [START_REF] Guillard | On the use of spectral methods for stiff problems[END_REF] with L = ∆. Figure 8 shows the results obtained with 32 nodes for values of the width parameter ∆ decreasing from 0.1 to 0.01. As in figure 5 the red dots displays the solutions on the collocation points where the solution is effectively computed while the red solid lines display u N the interpolating function passing through these points. Note that in contrast with conventional Fourier method, the interpolating function is not a trigonometric polynomial in the physical space : it is a member of the space S N (h) and results from the composition of a trigonometric polynomial with a change of variables.

As shown in these figures, the accuracy obtained with the MF method is very good. Not only the solution computed on the collocation points coincides with the exact solution but the interpolating function matches the exact solution almost exactly. Figure 9 presents a zoom of From a computational point of view, the computational cost of solving (1) with a transformed method is identical to the cost of solving a problem of the type :

- d dx [a(x) du dx ] = f (8)
This requires essentially 4 FFT while the use of conventional Fourier methods require only 2 FFT for the Laplace problem. However, in practice this additional cost is largely balanced by the reduction in the number of necessary nodes. Moreover, in many practical application, the diffusion coefficients are non-constant and the real problems to solve are of the type (8) instead of (1). The cost of the transformed method is thus almost identical to the cost of standard methods.

Pellet Injection

In magnetic confinement experiments for fusion research, the injection of cryogenic DT pellet in the plasma is considered as one of the most interesting way to fuel the core plasma. Pellets which are composed of neutral gas are rapidly heated by the plasma electrons and ionize to refuel the plasma. Moreover present tokamaks based on the so-called H-mode regime are known to be subject to some kind of MHD instabilities known as ELMs for Edge Localized Modes that can cause an excessive erosion of the plasma facing component. The triggering of ELMS by pellet injection has been demonstrated experimentally to be a possible technique to reduce the amplitude of naturally occurring large ELMs. However, a lot of uncertainties remains for the practical application of this technique in the future ITER tokamak and reliable numerical simulations can be extremely useful to reduce these uncertainties. The simulation of this problem is complex and involves the coupling of a non-linear MHD model with a pellet ablation model. This is usually done by introducing in the MHD model a density source term to model the ablation of the pellet. One of the difficulties of this type of simulation comes from the fact that the pellet size is much smaller than the tokamak size. A typical ratio of these dimensions is of order 10 -3 . Consequently the number of numerical simulations for this type of problem is relatively small : Among the few works devoted to this question one can notice the work of Samtaney et al [START_REF] Samtaney | 3d adaptive mesh refinement simulations of pellet injection in tokamaks[END_REF] that use AMR technique to deal with the size problem or Futatani et al [START_REF] Futatani | Non-linear MHD modelling of ELM triggering by pellet injection in DIII-D and implications for ITER[END_REF] that artificially increase the pellet size.

In this preliminary work, we do not intend to use a complex MHD model to deal with this question but rather we want to study how the use of mapped Fourier methods can help to deal with this problem.

One dimensional problem

In this preliminary section, we consider a very simple diffusion problem with a source term localized in space. More specifically the mathematical model writes :

∂u ∂t -u xx = S(x, t) (9)
where the source term is S(x, t) = 5ue

-( t-0.3 0.15 √ 2 ) 2 × e -( (x-π)) ∆ ) 2 (10)
As ∆ → 0, the solution of this equation becomes more and more difficult to compute accurately by Fourier methods. We have used here ∆ = 0.03 that corresponds to a difficult problem for Fourier methods. Equation ( 9) is discretized in space using the mapped Fourier method described in the previous section. Again the mapping is given by [START_REF] Guillard | On the use of spectral methods for stiff problems[END_REF]. The results of the mapped Fourier technique will be compared with those obtained by the standard Fourier pseudo-spectral method (that corresponds here to L = 1 in the formula ( 7)). A simple explicit first-order time advancing scheme have been used. More sophisticated and accurate schemes can be used but this was unnecessary for these preliminary experiments. We have found that a time step constraint of the form :

∆t ∼ 1 max j (h (ξ j )) 2 n 2 (11)
where n is the number of collocation points, is sufficient to ensure the stability of the scheme.

The results for standard Fourier methods are displayed in Fig. 10 for different times between t = 0 where the solution is uniform and equal to 2 and t = 0.5 for different numbers of Fourier modes. Due to the specific form of the source term, [START_REF] Shen | Error analysis for mapped legendre spectral and pseudospectral methods[END_REF] this one is maximum for t = 0.3 and almost zero for t > 0.5 and thus after t > 0.5 the solution simply relaxes to a constant. Fig. 10 shows that the computed solutions are characterized by the presence of huge oscillations whose amplitude exceed the maximum value of the true solution (Fig. 10 d.). As the number of Fourier modes increases, the amplitude of the oscillations are reduced. However this reduction is slow and 512 modes are necessary to obtain an oscillation free solution. In contrast, mapped Fourier methods can capture the localized source term quite accurately with a very moderate number of Fourier modes. Fig. 11 compares the results obtained with 64 modes and L = 0.1 with the result obtained by the standard Fourier technique with 512 modes. It is seen that the results obtained with the mapped Fourier technique is not only free from oscillations but also as accurate as the results obtained with standard techniques and 512 modes.

Two dimensional problem

In this section, a crude model of pellet injection in tokamaks is considered. This simulation do not intend to represent accurately the physical problem but tries to solve some of the numerical difficulties present in a complete model of pellet injection. For more complete modeling and numerical simulation, we refer to [START_REF] Samtaney | 3d adaptive mesh refinement simulations of pellet injection in tokamaks[END_REF], [START_REF] Futatani | Non-linear MHD modelling of ELM triggering by pellet injection in DIII-D and implications for ITER[END_REF]. The computational domain is a disk (r, φ) ∈ [1, 2] × [0, 2π] with r the radial direction and φ the toroidal one that corresponds to the midplane of a tokamak. The temperature2 evolution is governed by the heat equation and the ablation rate of the injected pellet during its passage through the background plasma is expressed in terms of a source term localized in time and space. More specifically the mathematical model writes :

∂u ∂t -∇.q = S(x, t) (12)
where q is the heat flux. Due to the anisotropic character of magnetized plasma, the heat flux in the parallel direction to the magnetic field far exceeds the one in the perpendicular direction.

Considering in our problem that the magnetic field is purely parallel to the toroidal direction, we have used :

q = κ // r ∂u ∂φ e φ + κ ⊥ ∂u ∂r e r (13)
with κ // = 10(u/u 0 ) 5./2 and κ ⊥ = 0.01. The initial temperature varies linearly between u = 2 at r = 1 and u = 1 at r = 2. Due to the highly non-linear behavior of the parallel heat diffusion coefficient, this one varies between 1.76 and 10 and is larger in the center of the disk. where the temperature is large. The ratio between parallel and perpendicular diffusion is between 176 and ∼ 1000. The source term modeling pellet injection is given by : S(x, t) = -250ue

-( t-0.3 0.15 √ 2 ) 2 × e -(φ-π) 2 +(r-1.5) 2 ) ∆ 2 (14)
More sophisticated models do exist and are available, a commonly used model for these studies being the one by Parks and Turnbull [START_REF] Parks | Effect of transonic flow in the ablation cloud on the lifetime of a solid hydrogen pellet in a plasma[END_REF]. However the simple model (14) contains the essential numerical difficulties that we want to address and is sufficient for our purpose. Note that the source term is negative as the injection of cryogenic pellets in a very hot plasma gives rise to a local drop of the plasma temperature. The parameter ∆ is related to the pellet radius that is by several orders of magnitude smaller than the tokamak size. In this work, we have taken ∆ = 0.03. A mapped Fourier method has been used in the toroidal direction and a simple finite difference one in the radial direction. In Fig. 13 are compared results obtained from standard Fourier methods using from 64 to 256 modes with the results obtained with a mapped Fourier method with 32 points. In all these computations, the number of nodes in the radial direction is 32. The standard computations done with 64 and 128 Fourier are clearly unusable due to the presence of large oscillations on the opposite side of the disk. This is particularly true for tokamak researches since pellet injection is considered as a possible way to mitigate the occurrence of large scale ELMs (Edge Localized modes) that can potentially affect the integrity of the vacuum chamber of tokamaks by triggering small amplitude instabilities. The presence of large scale oscillations on the opposite side of the tokamak with respect to the location of the pellet source can then lead potentially to erroneous conclusions on the triggering of instabilities by pellet injections. For a number of mode equal to 256, these oscillations almost disappear (actually some small amplitude fluctuations are still present in these computations). However this large number of modes is far beyond what can be done for these computations considering the large number of nodes that must be used in the poloidal section. In contrast, mapped Fourier methods are particularly efficient in this context and only 32 nodes have been required to obtain oscillation free computations. Fig. 14 presents an enlargement of the source region showing that the results obtained with mapped Fourier methods compare very favorably to the results obtained by standard methods with a much larger number of points.

Conclusion and remarks

This preliminary work intends to demonstrate the possibilities of mapped Fourier methods for improving the computation in periodic domains of the solution of PDE when these solutions possesses internal layers of rapid variation. We have concentrated this study on the simple case where the localization of this zone is a priori known and do not evolve in time. This is a simplifying assumption but it already covers many practical situations. In particular, for Fusion Plasma, we believe that this technique can be applied to pellet injection simulations in tokamak with great benefit. Of course, this remains to be proved by using more sophisticated models than the one used here. This will be the purpose of a next series of investigations.

We also note that problem of adapting the mapping to the characteristics of the solution is also of great interest for the development of new numerical methods. The design of adaptive meshes linked with a priori or a posteriori error estimates is a mature field for local methods (FD/FE ones). This is not true for global (spectral) methods. In this area some works do exist, Figure 14: Left : Standard Fourier with 256 modes, Right : Mapped Fourier with 32 modes see [START_REF] Guillard | On the use of spectral methods for stiff problems[END_REF][START_REF] Bayliss | An adaptive pseudo-spectral method for reaction diffusion problems[END_REF] for instance. Boyd [START_REF] John | Chebyshev and Fourier spectral methods[END_REF] table 16.2 gives a list of works performed in this direction. However, the majority of these works are in a simple 1D setting and real large scale applications using these techniques are still missing. We hope that this work can help to show that this area of investigation may have some interest.

Figure 1 : 2 ∆ 2 for ∆ = 0. 1 Figure 2 :

 12212 Figure 1: Pseudo Spectral representation of the first derivative of e -(x-π) 2 ∆ 2 for ∆ = 0.1

Figure 3 :Figure 4 : 1 Figure 5 :

 3415 Figure 3: Evolution of the error vs ∆

Figure 6 :

 6 Figure 6: Mapped Fourier representation of the first derivative of e -(x-π) 2 ∆ 2

Figure 7 :

 7 Figure 7: Mapped Fourier representation of the second derivative of e -(x-π) 2 ∆ 2 for : left plot ∆ = 0.1, right plot ∆ = 0.01, note that the right plot is a zoom on the central region

Figure 8 :

 8 Figure 8: Numerical solution of equation (1) with 32 nodes : left ∆ = 0.1, middle ∆ = 0.05, right ∆ = 0.01

Figure 9 :

 9 Figure 9: Zoom of the Numerical solution of equation (1) with 32 nodes : ∆ = 0.01

Figure 10 :

 10 Figure 10: 1D diffusion problem. Standard Fourier : a) left top n = 64, b) right top n = 128, c) left bottom n = 256, d) right bottom n = 512

Figure 11 :

 11 Figure 11: 1D diffusion problem Mapped Fourier (MF) compared to standard method : a) left Standard Fourier with 512 modes, b) right MF using 64 modes and L = 0.1

  Fig 12 presents the distribution of collocations points in the computational domain. The right plot shows that the points are highly concentrated at the location of the source.

Figure 12 :

 12 Figure 12: Left : Distribution of collocation points in the physical space, Right : Zoom on the source region

Figure 13 :

 13 Figure 13: Top Left : Standard Fourier with 64 modes, Top Right Standard Fourier with 128 modes, Bottom Left : Standard Fourier with 256 modes, Bottom Right : Mapped Fourier with 32 modes

The case where the solution possesses shocks or discontinuities is different : for shock, spectral accuracy is never recovered, even if N is very large. The present work does not address this problem and concentrates on problems with smooth solutions.RR n°8566

RR n°8566

alternatively this equation can be considered as a density evolution equation but in this case the sign in front of the source term must be changed RR n°8566

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399