Quan V Nguyen

Francis Colas

Emmanuel Vincent

Franc ¸ois Charpillet

Localizing an Intermittent and Moving Sound Source Using a Mobile Robot

This paper addresses the problem of localizing and tracking one intermittent, moving sound source using a microphone array on a mobile robot. Robot motion provides a solution for estimating the distance to the source and avoiding front-back ambiguity. We propose a mixture Kalman filter (MKF) framework in order to fuse the robot motion information and the measurements taken at different poses of the robot. Experiments and statistical results demonstrate the ability of the proposed method to track one intermittent sound source in a reverberant environment where false measurements of the source angle of arrival (AoA) and the source activity often occur compared to a method that does not consider tracking source activity into account.

I. INTRODUCTION

As assistive robots are becoming more popular, artificial hearing capabilities, i.e. robot audition [START_REF] Okuno | Robot audition: Its rise and perspectives[END_REF], are widely recognized as essential for robot perception. Embedded microphones grant robots the ability to estimate and track the spatial location of sound sources over time [START_REF] Marković | Active speaker localization with circular likelihoods and bootstrap filtering[END_REF], [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF], [START_REF] Vincent | Audio source localization by optimal control of a mobile robot[END_REF]. Having a better knowledge about location of the sources can help the robot to localize itself with respect to known sound sources [START_REF] Nakadai | Real-time sound source localization and separation for robot audition[END_REF], [START_REF] Valin | Robust recognition of simultaneous speech by a mobile robot[END_REF], [START_REF] Martinson | Discovery of sound sources by an autonomous mobile robot[END_REF] as well as to improve the performance of source separation and speech recognition and, consequently, to efficiently interact with humans and the environment. For these reasons, source localization from multi-microphone recordings plays a central role in robot audition.

Broadly speaking, source localization techniques can be classified into three families [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF]. One approach consists in computing the time delay of arrival (TDOA) between every pair of microphones using generalized cross-correlation with phase transform (GCC-PHAT) [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF] and to derive the source position by triangulation. This approach is typically outperformed by steered response power (SRP) [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF] or multiple signal classification (MUSIC) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] techniques that compute the pseudo-likelihood of each candidate position on a grid, and pick the maxima on that grid (see [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF], [START_REF] Badali | Evaluating realtime audio localization algorithms for artificial audition in robotics[END_REF], [START_REF] Blandin | Multi-source TDOA estimation in reverberant audio using angular spectra and clustering[END_REF] for experimental comparisons). Binaural variants of these techniques have been developed for situations when the array is mounted on the robot head [START_REF] Nakadai | Active audition for humanoid[END_REF]. In the typical situation when the distance to the source is larger than the array size, all of these techniques can only estimate the source angle of arrival (AoA) but not its distance.

Most of these approaches have been implemented on robots [START_REF] Nakadai | Active audition for humanoid[END_REF], [START_REF] Kim | Speaker localization using the TDOA-based feature matrix for a humanoid robot[END_REF], [START_REF] Nakamura | Real-time super-resolution sound source localization for robots[END_REF] and used to estimate the source AoA in a similar way as with a static microphone array. Robots can actually provide more information by exploiting motion: this is known as active audition. The absolute position of a source can be tracked by fusing the motion of the robot and the auditory perception. This is done using probabilistic filters such as nonlinear extensions of Kalman filtering [START_REF] Portello | Acoustic models and Kalman filtering strategies for active binaural sound localization[END_REF], [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF], particle filters [START_REF] Marković | Active speaker localization with circular likelihoods and bootstrap filtering[END_REF], [START_REF] Lehmann | Particle filter design using importance sampling for acoustic source localization and tracking in reverberant environment[END_REF], [START_REF] Valin | Robust localization and tracking of simultaneous moving sound sources using beamforming and particle filtering[END_REF], [START_REF] Fallon | Acoustic source localization and tracking of a time-varying number of speakers[END_REF], or occupancy grids [START_REF] Martinson | Auditory evidence grids[END_REF], [START_REF] Dejong | Auditory occupancy grids with a mobile robot[END_REF], [START_REF] Vincent | Audio source localization by optimal control of a mobile robot[END_REF]. Compared to a static microphone array, the head or whole body movements performed in an active audition setting allow the robot to avoid front-back confusion and to estimate the distance to the source. Most of these techniques consider the situation when the sound source to be tracked is continuously active [START_REF] Portello | Acoustic models and Kalman filtering strategies for active binaural sound localization[END_REF], [START_REF] Martinson | Auditory evidence grids[END_REF], [START_REF] Vincent | Audio source localization by optimal control of a mobile robot[END_REF].

However, the sound sources in real life, e.g. speech, are not always active. Silence exists between successive words, successive sentences, and speaker turns. The problem of tracking one intermittent and moving source by using a nonlinear mixture Kalman filter (MKF) is presented in [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF]. In this work, the source activity (i.e., whether the speaker is active or inactive in a given time frame) is supposed to be known a priori. The experimental evaluation in anechoic conditions showed promising results [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF]. However, in the real world, reverberation (i.e., sound echoes on the walls of the room) and noise are present that severely degrade the performance of source localization [START_REF] Vermaak | Nonlinear filtering for speaker tracking in noisy and reverberant environments[END_REF], [START_REF] Blandin | Multi-source TDOA estimation in reverberant audio using angular spectra and clustering[END_REF] and Source Activity Detection (SAD) [START_REF] Ramírez | Voice activity detection. fundamentals and speech recognition system robustness[END_REF], [START_REF] Zhang | Deep belief networks based voice activity detection[END_REF]. Joint estimation of the activity and false measurement from AoA measurements alone is introduced in [START_REF] Portello | Active binaural localization of intermittent moving sources in the presence of false measurements[END_REF]. The accuracy of this activity detection which depends on the AoA measurement is presumably poor compared to an independent SAD. We claim that joint estimation of the source location and activity from AoA measurement and SAD is required to achieve robust estimation in such real-world conditions.

The contributions of this paper are twofold. First, we formulate a framework that applies to any microphone array geometry for explicitly tracking both the location and activity of one intermittent moving sound source in a reverberant environment, with a nonlinear MKF. This choice of filter is motivated by the fact that the state vector consists both of continuous and discrete variables, as well as the fact that the observation model is nongaussian. Also, it facilitates comparison with the algorithm in [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF]. Second, we provide an experimental validation of our algorithm in several realistic scenarios.

The remainder of the paper is organized as follows. Section II presents the mathematical model of MKF. The evaluation protocol and the results are respectively described in Sections III and IV. Section V concludes the paper.

II. MIXTURE KALMAN FILTER

We consider the problem of estimating the absolute spatial position of a single intermittent, moving audio source by a moving robot. We represent the belief about the source position, motion and activity by a mixture of Gaussians. At the beginning, this mixture expresses all possible locations and activities of the source in the room. The goal of the MKF is to update the belief over time given a sequence of measurements. In this section, we first define the state vector and the observation vector, and then detail the proposed MKF algorithm.

A. State vector

By contrast with [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF], we take the source activity into account in the state vector. We define the state vector as:

X a =   X r X s a   =               x r y r θ r x s y s θ s v s w s a               (1)
where X r is the pose of the robot, i.e., its absolute position [x r , y r ] and its orientation θ r w.r.t. the x-axis; X s is the state of the sound source, i.e., its absolute position [x s , y s ] , its orientation θ s w.r.t. the x-axis, and its linear and angular velocities [v s , w s ]; a is source activity, where a = 1 indicates that the source is active, otherwise a = 0. We assume that the pose of the robot X r is known and we need to estimate state of the source X s and its activity a. With the state of the robot in the state vector, this model have potential to deal with non-zero process noise of the robot motion when we have additional observation for the robot.

B. Observation vector

As mentioned above, audio source localization techniques can estimate the source AoA but not its distance. Therefore, we assume that the observation vector Z k in a given time frame k consists of one AoA measurement Z l k (obtained via a localization technique) and one source activity measurement Z a k (obtained via an SAD technique). The likelihood of the state vector w.r.t. this observation can be expressed as

P (Z k |X k , a k ) = P sn (Z l k |X k )P (Z a k |a k ) for a k = 1 P n (Z l k )P (Z a k |a k ) for a k = 0
(2) with P sn and P n denoting the distribution of the measured AoA when the source is active or inactive, respectively. In the latter case, it is supposed that the recorded signal consists of spatially diffuse noise, so P n does not depend on X k .

An example of P sn is shown in Fig. 1 for our linear microphone array and the localization technique considered in Section III-A. The probability density concentrates around the true AoA and its symmetric w.r.t the microphone axis: this phenomenon is known as front-back confusion.

The probability density for other AoAs is nonzero, but much smaller. Therefore we can approximate the observation model by a mixture of 2 Gaussians:

P sn (Z l k |X k ) = M k j=1 ω j N Z [h j (X k ), R i,j k ] (3)
with

M k j=1 ω j = 1 and M k = 2.
For a nonlinear microphone array, the observation model could be represented by only a single Gaussian.

C. Recursive Bayesian estimation

The transition probability between time steps is given by

P (X k , a k |X k-1 , a k-1 ) = P (a k |X k-1 , a k-1 )P (X k |X k-1 , a k-1 , a k ). (4)
The source activity a k and the state X k are conditionally independent so we can rewrite the above equation as

P (X k , a k |X k-1 , a k-1 ) = P (a k |a k-1 )P (X k |X k-1 ). (5)
The state transition probability P (X k |X k-1 ) is defined by the dynamic model

X k = f (X k-1 , u) + d (6)
where u are the robot commands and d is the process noise of the robot and the sound source with covariance matrix Q.

Note that the control input u for the robot is given during the estimation, so for simplicity it is later omitted from the equations. The source activity transition probability P (a k |a k-1 ) is defined by P appear = P (a k = 1|a k-1 = 0) which is the source appearance probability or P disappear = P (a k = 0|a k-1 = 1) which is the source disappearance probability.

The posterior probability of the state vector can be recursively computed by alternating these two steps:

• prediction: compute P (X k , a k |Z 1:k-1 ) given the previ- ous belief P (X k-1 , a k-1 |Z 1:k-1
) and the state transition model, as shown in [START_REF] Martinson | Discovery of sound sources by an autonomous mobile robot[END_REF],

P (X k , a k |Z 1:k-1 ) = a k-1 P (X k , a k |X k-1 , a k-1 )P (X k-1 , a k-1 |Z 1:k-1 )d(X k-1 ) = a k-1 P (a k |a k-1 )P (X k |X k-1 )P (X k-1 , a k-1 |Z 1:k-1 )d(X k-1 ) (7)
P (X k , a k |Z 1:k-1 ) = N k|k-1 i=1 ω i k|k-1 N ( Xi k|k-1 , P i k|k-1 )δ(a k = 0) + N k|k-1 (M k-1 +1) i=N k|k-1 +1 ω i k|k-1 N ( Xi k|k-1 , P i k|k-1 )δ(a k = 1) (8)
ω i k|k-1 =          ω i k-1|k-1 P (a k = 0|a k-1 = 0), if i ∈ 1, N k-1|k-1 ω i k-1|k-1 P (a k = 0|a k-1 = 1), if i ∈ N k-1|k-1 + 1, N k-1|k-1 (M k-1 + 1) ω i-N k-1|k-1 (M k-1 +1) k-1|k-1 P (a k = 1|a k-1 = 0), if i ∈ N k-1|k-1 (M k-1 + 1) + 1, N k-1|k-1 (M k-1 + 2) ω i-N k-1|k-1 (M k-1 +1) k-1|k-1 P (a k = 1|a k-1 = 1), if i ∈ N k-1|k-1 (M k-1 + 2) + 1, N k-1|k-1 (2M k-1 + 2) (9)
P (X k , a k |Z 1:k ) =η N k|k-1 i=1 ω i k|k-1 P (Z a k |a k )Pn(Z l k )N ( Xi k|k-1 , P i k|k-1 )δ(a k = 0)+ η 2N k|k-1 i=N k|k-1 +1 ω i k|k-1 P (Z a k |a k )Psn(Z l k |X k )N ( Xi k|k-1 , P i k|k-1 )δ(a k = 1) (10)
P (X k , a k |Z 1:k ) =η N k|k-1 i=1 ω i k|k-1 P (Z a k |a k )Pn(Z l k )N ( Xi k|k-1 , P i k|k-1 )δ(a k = 0)+ η 2N k|k-1 i=N k|k-1 +1 M k j=1 ω i k|k-1 ω j P (Z a k |a k )NZ [h j (X k ), R i,j k ]N ( Xi k|k-1 , P i k|k-1 )δ(a k = 1) (11)
λ i,j = 1 2π(HP i k|k-1 H T + R i,j k ) e -1 2 Z j k -h( Xi k|k-1 )) T HP i k|k-1 H T +R i,j k -1 Z j k -h( Xi k|k-1 )) (12)
P (X k , a k |Z 1:k ) =η N k|k-1 i=1 ω i k|k-1 P (Z a k |a k )Pn(Z l k )N ( Xi k|k , P i k|k )δ(a k = 0) + η 2N k|k-1 i=N k|k-1 +1 M k j=1 ω i k|k-1 ω j P (Z a k |a k )NZ [h j ( Xk ), R i,j k ]N ( Xi k|k-1 , P i k|k-1 )δ(a k = 1) = N k|k i=1 ω i k|k N ( Xi k|k , P i k|k )δ(a k = 0) + N k|k (M k +1) i=N k|k +1 ω i k|k N ( Xi k|k , P i k|k )δ(a k = 1) (13)
P (X k-1 , a k-1 |Z 1:k-1 ) = N k-1|k-1 i=1 ω i k-1|k-1 N ( Xi k-1|k-1 , P i k-1|k-1 )δ(a k-1 = 0)+ N k-1|k-1 (M k-1 +1) i=N k-1|k-1 +1 ω i k-1|k-1 N ( Xi k-1|k-1 , P i k-1|k-1 )δ(a k-1 = 1) (14)
• update: recompute the belief P (X k , a k |Z 1:k ) given the prediction and the new measurement

Z k P (X k , a k |Z 1:k ) = ηP (Z k |X k , a k )P (X k , a k |Z 1:k-1 ) (15
) where η is a normalizing constant.

D. Derivation of the MKF

Since the state vector includes both continuous and discrete variables and the observation model is a mixture of Gaussians, we propose a MKF to address these two issues.

1) Prediction step: Assume that at the previous time step k -1 the belief about the state (X k-1 , a k-1 ) is given by the mixture of Gaussians in [START_REF] Kim | Speaker localization using the TDOA-based feature matrix for a humanoid robot[END_REF] with weights ω i k-1|k-1 such that i ω i k-1|k-1 = 1. Applying the prediction rule [START_REF] Martinson | Discovery of sound sources by an autonomous mobile robot[END_REF] to this density yields the predicted density [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF]. This is a mixture of Gaussians, whose weights ω i k|k-1 are expressed in ( 9) and means and variances are given by

Xi k|k-1 = f ( Xi k-1|k-1 , u) (16)
F i k-1 = ∂f k-1 (X, u) ∂X | X= Xk-1|k-1 (17)
P i k|k-1 = F i k-1 P i k-1|k-1 F iT k-1 + Q k-1 . ( 18
)
2) Update step: By applying the update rule (15) to the predicted density and replacing the observation model by the mixture of 2 Gaussians in (3), we obtain the new belief in [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] and [START_REF] Badali | Evaluating realtime audio localization algorithms for artificial audition in robotics[END_REF]. The product of every 2 Gaussians

N Z [h j (X k ), R i,j k ] N ( Xi k|k-1 , P i k|k-1
) can be computed in closed form as λ i,j N ( Xi k|k , P i k|k ) where λ i,j is a constant defined in [START_REF] Blandin | Multi-source TDOA estimation in reverberant audio using angular spectra and clustering[END_REF]. The new belief can therefore be expressed as in [START_REF] Nakadai | Active audition for humanoid[END_REF], where

ω i k|k = ω i k|k-1 P n (Z l k )P (Z a k |a k )η if a k = 0 ω i k|k-1 ω j P (Z a k |a k )λ i,j η if a k = 1 (19)
If

a k = 0, then Xi k|k = Xi k|k-1 . If a k = 1, then Xi k|k = Xi k|k-1 + G i k Z j k -h( Xi k|k-1 ) (20)
H i k = ∂h(X) ∂X | X= Xk|k-1 (21)
P i k|k = P i k|k-1 -G i k H i k P i k|k-1 (22)
S i k = H i k P i k|k-1 H iT k + R i,j k ( 23
)
G i k = P i k|k-1 H iT k (S i k ) -1 . ( 24
)
3) Pruning: From (11), we can realize that the number of hypotheses in the MKF increases over time, which will consume a lot of memory. To deal with this problem, when the number of hypotheses is larger than N max we keep only the N min hypotheses with the highest weights and prune the other hypotheses which have much lower weights.

III. EXPERIMENTAL PROTOCOL

We conducted numerical experiments to evaluate our MKF algorithm for tracking one intermittent and moving speech source with room reverberation and noise. Our experimental settings mimic the smart room at Inria Nancy, where the robot is a Turtlebot equipped with a Kinect sensor. In this work, we are only using the linear array of 4 microphones included in the Kinect.

A. Data

Due to the statistical nature of false measurements, a large number of experiments is needed to obtain statistically meaningful results. Such a large number of experiments can hardly be conducted with a real robot. Therefore, we resort to simulation of the robot movements, the source movements, and the resulting location and activity measurements. We employ state-of-the-art techniques for the simulation of reverberation and acoustic noise, whose parameters are fixed as in [START_REF] Vincent | Audio source localization by optimal control of a mobile robot[END_REF] and closely match the real conditions in that room. More specifically, the reverberation time (250 ms), the intensity of speech and noise, and the noise spectrum match those of the real environment. The source AoA is estimated by MUSIC with generalized singular value decomposition (GSVD) [START_REF] Nakamura | Real-time super-resolution sound source localization for robots[END_REF] as implemented in HARK [START_REF] Nakadai | Design and implementation of robot audition system 'HARK' -open source software for listening to three simultaneous speakers[END_REF]. The probability distribution of AoAs estimated by MUSIC-GSVD for each of 360 true AoAs (from 0 • to 359 • ) and 5 distances (from 0.5 to 3 m) was constructed and used to simulate the observed source location. We consider an SAD error rate of 5%. We considered four different scenarios, depending whether the sound source is static or mobile (v s = 0.07 m/s, w s = 8 • /s) and inactive for many short time intervals (0.5 s) or a long time interval (2 s). For each scenario, we randomly generated 100 source trajectories for a duration of 10 s. The robot trajectory was fixed in all experiments with a maximum speed of 0.38 m/s.

B. Algorithm settings

We set the parameter values of the MKF as follows. The time step is 0.1 s. The covariance matrix Q was set as Q = diag(0, 0, 0, 0.00095 m 2 , 0.00062 m 2 , 6.2 • 2 , 0, 0), the initial covariance of the pose of the robot is P R 0|0 = diag(0, 0, 0), the variance R i,j varied as a function of the source distance between 0.8 • 2 at 0.3 m and 4.5 • 2 at 3 m, the source appearance/disappearance probabilities were set to P appear = 0.5 and P disappear = 0.5, and the number of hypotheses in the MKF to N max = 300 and N min = 24. By testing our model with different values of the hyperparameters (e.g., P appear , P disappear and SAD error rate) in separated experiments, we see that our model is robust to those hyper-parameters.

We evaluate the performance of our MKF method with vs. without tracking of the source activity. Our method without

IV. RESULTS

A. Example scenario -Visualization

Fig. 2 shows the first few seconds of tracking one intermittent, moving source. At time t = 0 s, the mixture is initialized with several components evenly distributed over the room in order to approximate a uniform prior. After 1 s, half of the hypotheses for the source position are distributed along the direction from the source to the robot and the rest are symmetric w.r.t. the microphone axis. This symmetrical uncertainty is due to the front-back confusion phenomenon illustrated in Fig. 1. These symmetrical hypotheses become smaller and disappear after 3 s, thanks to the robot motion. More precisely, the motion of the source for the symmetric components is bigger, less coherent and therefore less probable than of the correct components. For a nonlinear microphone array, the transitory phase with two directions would not exist.

B. Example scenario -Estimated trajectories

Fig. 3 compares the source trajectory with the estimations of our MKF and Portello MKF. As both models are mixture models, the posterior distribution is usually not unimodal. In order to generate a single point estimate, we simply compute the mean of the distribution. After the first few seconds discussed above, both trajectories follow the sound source. We can also observe that, there are some moments when the estimated source location of Portello MKF is far away from actual source location, however, our MKF still can track the source location with lower estimation error.

Fig. 4 shows the estimation error, that is the distance between the estimated source position and the true position. During the first 3 s, both MKFs have high estimation error due to the front-back ambiguity. Between 3.5 and 3.9 s, the source is inactive and the SAD is correct. During this period, the uncertainty about the source position increases because of the lack of measurements. The uncertainty gets smaller when source becomes active again.

At time t = 2.8 s, a false measurement of the source activity occurs: the source is active but SAD detects it as inactive. The estimation error of Portello MKF becomes larger than ours but only for one time step. Conversely, from t = 3.2 s to 7 s, the estimation error of Portello MKF is lower than ours but by 2 cm only.

A false measurement of AoA occurs at time t = 4.8 s. The AoA difference between observation and ground truth is 9 • , this is not a big value. As a result, both the estimated error of our MKF and Portello MKF have a small rise.

At time t = 7 s, a false AoA measurement occurs: the ground truth AoA is 81 • , but the measured AoA is 62 • . Although such a false measurement can occur with very low probability, it can have a major impact. Indeed, the estimation error of Portello MKF increases drastically and remains large. By contrast, the estimation error of our MKF does not change much. This is an unexpected benefit of the proposed approach: when a false AoA measurement occurs, the weight of the hypotheses corresponding to an inactive source increases, so that the belief is little affected.

At time t = 8.2 s, a false measurement of the source activity occurs: the source is inactive but SAD detects it as active. Again, it appears that our MKF can handle this situation but Portello MKF severely suffers from it. This was expected, since the tracking of the source activity implemented by our method was designed precisely to address this issue.

C. Statistical analysis

In order to conduct a statistical comparison of our method and Portello's, we compare the distributions of errors between the estimate and the ground truth. These distributions, shown in Fig. 5, accumulate the error for all time steps and all experimental runs. We can see that the error for our MKF is mostly concentrated below 1 m while Portello MKF shows a bigger error distribution with a heavy tail. This is confirmed with a Wilcoxon signed-rank test that assesses that our MKF has a significantly smaller error than Portello with p < 0.01. Interestingly, our MKF outperforms Portello MKF in all four scenarios, whether the source is static or moving and includes short or long silences.

V. CONCLUSIONS

We presented an MKF that applies to any microphone array geometry for tracking one intermittent, moving sound source in a reverberant environment using a mobile robot. The main theoretical contribution of our method is the explicit estimation of the source activity, which allows us to cope with imperfect SAD algorithms.

Experimental results have demonstrated the significantly better performance of our algorithm compared to recent work in the literature [START_REF] Portello | Active localization of an intermittent sound source from a moving binaural sensor[END_REF]. An additional advantage of our algorithm is its ability to track the location of the source in the presence of false AoA measurements.

This work is for now done in a realistic simulation but evaluation on a real robot should highlight the advantage of using a realistic sensor model in the estimation. Another perspective of this work is its extension to multiple sound sources, integrating signal characteristics to help disambiguate the data association issue.

Fig. 1 .

 1 Fig. 1. Distribution of the measured AoA when the actual source is at 96 • and 0.5 m from the microphone array.

Fig. 2 .

 2 Fig. 2. Visualization of our MKF in the example scenario. Robot positions are shown as red squares, and the actual source position as a green circle. Blue ellipses represent 95% confidence regions of source location estimation of various hypotheses in the mixture with a transparency proportional to the weight of the components.

Fig. 3 .

 3 Fig. 3. Estimated trajectories in the example scenario. Dashed black line: trajectory of the robot, Solid black line: ground truth trajectory of the sound source, Green line: estimate of Portello MKF, Red line: estimate of our MKF.

Fig. 4 .

 4 Fig. 4. Top: Estimation error of our MKF vs Portello MKF. Bottom: Ground truth of source activity.

Fig. 5 .

 5 Fig. 5. Estimation error distribution of our MKF method in red and Portello MKF method in blue.

Inria, Villers-lès-Nancy, F-54600, France

CNRS, Loria, UMR n o 7503, Vandoeuvre-lès-Nancy, F-54500, France

Université de Lorraine, Loria, UMR n o 7503, Vandoeuvre-lès-Nancy, F-54500, France firstname.lastname@inria.fr

The difference with their method lies in the number of Gaussians used to approximate the observation likelihood.


Holger Heumann

Lukas Drescher

Kersten Schmidt

A High Order Method for Contour Integrals with an Application to Plasma Modeling in Nuclear Fusion

Keywords: contour integrals, high order quadrature, curved interfaces

We introduce a novel method to compute approximations of contour integrals. The new method is based on the coarea formula in combination with a Galerkin projection. As such it fits seamlessly into the spirit of hp/spectral finite element methods and circumvents the expensive and technical computation of contours. We provide convergence estimates showing that high order convergence can be achieved provided the data is sufficiently smooth. The theoretical results are supplemented by extensive numerical experiments for an example application from plasma modeling in nuclear fusion.

f (x) ds(x) , y ∈ (0, 1),

where ψ : Ω ⊂ R n → [0, 1], n = 2, 3 is a Lipschitz continuous, scalar function, and f : Ω → R ∪ {-∞, ∞} is another integrable scalar function. The values g f,ψ (y) are well defined for all y for which the sets {x ∈ R n , ψ(x) = y} are (n -1)-dimensional hypersurfaces and f is bounded on these sets. The functions g f,ψ represent the average of f over each hypersurface of ψ. In the particular case f = 1 we have that g f,ψ is the area (or length) of the hypersurface. We will assume throughout this work that the area (length) g 1,ψ (y) is bounded for almost all y ∈ [0, 1].

If the hypersurfaces of ψ have an explicit parametrization, standard quadrature formulas provide an adequate mean to approximate [START_REF] Atkinson | Numerical evaluation of line integrals[END_REF]. It is even possible to devise higher order methods that avoid the computation of the derivatives of the parametrization [START_REF] Atkinson | Numerical evaluation of line integrals[END_REF].

In here we want to focus on the case where such a parametrization is a priori not known. The standard approach is the following: The scalar functions ψ are approximated by piecewise linear finite element functions ψ h and the values g f,ψ h (y) can be assembled from contributions in each element. This approach is a variation of a family of algorithms in computer graphics [START_REF] Newman | A survey of the marching cubes algorithm[END_REF] that are used to visualize hypersurfaces. In the case of linear finite elements this is the method of marching triangles or tetrahedrons [START_REF] Payne | Surface mapping brain function on 3d models[END_REF][START_REF] Shirley | A polygonal approximation to direct scalar volume rendering[END_REF], while projecting onto the bilinear finite element space is the even older method of marching squares or cubes [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]. There does not exist, at least to our knowledge, rigorous convergence analysis. But clearly the projection of ψ onto low order spaces restricts the order of the method. Very high accuracy can be only achieved through meshes with a huge number of elements.

It is a well known fact in spectral and higher order methods, that one can find way more economical finite dimensional approximations of ψ, if one exploits properties such as smoothness. Nevertheless we would like to stress that the direct elementwise assembling of g f,ψ h with ψ h piecewise polynomials of degree p will be more difficult and expensive. In the case of linear finite elements the intersections of hypersurfaces with the mesh elements are planar, while in the case of higher order finite elements the intersection are curved. Hence, it is necessary to determine numerically a local parametrization if the accuracy of ψ h shall be inherited to the approximation of g f,ψ h . Finding such local parametrizations requires similar techniques as finite element methods on curved domains [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF]Chapter 12] and amounts to solve repeatedly non-linear problems. For higher order polynomial approximations of ψ this approach is fairly expensive and technical. Such efforts might be reasonable if one is interested in values of (1) for one specific value y as it is the case for Date: August 17, 2016, 2010 Mathematics Subject Classification: 65D30, 65N15, 65N30, 76W05 keywords: contour integrals, high order quadrature, curved interfaces. hypersurface integrals that appear in level set methods and fictitious domain methods or extended finite element methods. A non-exhaustive list of references about this topic is [START_REF] Engquist | Discretization of Dirac delta functions in level set methods[END_REF][START_REF] Min | Geometric integration over irregular domains with application to level-set methods[END_REF][START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF] and the many references therein.

In this article we are interested in computing efficiently approximation of (1) on the whole interval y ∈ (0, 1). We propose a novel approach fitting seamlessly into the spirit of higher order finite element methods. We will see that the hypersurface averaging functions (1) are unique solutions of the variational formulation: Seek g f,ψ ∈ L 2 (0, 1) such that

1 0 g f,ψ (y)λ(y)dy = Ω f (x)|∇ψ(x)|λ(ψ(x)) dx ∀λ ∈ L 2 (0, 1). (2)
Moreover, we will introduce associated Galerkin methods, where only integrals over the domain Ω have to be evaluated.

The main motivation for this work is the magnetohydrodynamic equilibrium equation, an elliptic semi-linear equation for the so-called poloidal flux ψ. Computing magnetohydrodynamic equilibria is a central problem for modeling of plasma in nuclear fusion reactors. Very important functionals of the solution of such equilibrium problems are the hypersurface averaged quantities (1) that are the focus of this article. In the subsequent section we will adapt the parlance from plasma physics and call the functions (1) the geometric coefficients.

The outline of this paper is the following. In the next section, we will recall the coarea formula and show (Lemma 2) that the geometric coefficients verify the weak formulation [START_REF] Beretta | An inverse problem originating from magnetohydrodynamics. II. The case of the Grad-Shafranov equation[END_REF]. With this it is straight forward to introduce in Section 3 the Galerkin projection methods to compute approximations of g f,ψ (y) in some finite dimensional space, e.g., a polynomial space. We provide also a convergence estimate that balances a best approximation result and the consistency error. Section 4 presents the plasma equilibrium problem, elaborates on details for computing the geometric coefficients and states a couple of regularity results. Extensive numerical experiments highlight the practicability of the new method.

GEOMETRIC COEFFICIENTS AND THE COAREA FORMULA

We recall that ψ : Ω ⊂ R n → (0, 1), n = 2, 3 is a Lipschitz continuous, scalar function, not necessary monotone and that the geometric coefficients are the integrals

g f,ψ (y) := {x∈Ω, ψ(x)=y} f (x) ds(x) , with f : Ω → R ∪ {-∞,
∞} another integrable scalar function. We have the following result, the coarea formula, highlighting the relationship between integrated geometric coefficients g f,ψ (y) and volume integrals.

Theorem 1 (Coarea formula). Let f : Ω → R ∪ {-∞, ∞} be integrable and ψ : Ω → [0, 1] Lipschitz, then g f,ψ ∈ L 1 (0, 1) and 1 0 g f,ψ (y)dy = Ω f (x)|∇ψ(x)|dx .
Proof. See [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF]Theorem 3.2.12].

An immediate corollary of the coarea formula is the following identity.

Corollary 2. Let the assumptions of Theorem 1 be fulfilled. If moreover f |∇ψ| 1 2 ∈ L 2 (Ω) and g 1,ψ ∈ L ∞ (0, 1) then g f,ψ ∈ L 2 (0, 1) and

(3)

1 0 g f,ψ (y) λ(y)dy = Ω f (x) λ(ψ(x)) |∇ψ(x)| dx ∀λ ∈ L 2 (0, 1) .

Inria

Proof. The assertion is an immediate consequence of the obvious identity g f,ψ (y)λ(y) = g f •λ(ψ),ψ (y) and Theorem 1, and so

1 0 g f,ψ (y) λ(y)dy = 1 0 g f •λ(ψ),ψ (y)dy = Ω f (x)λ(ψ(x))|∇ψ(x)| dx.
The Cauchy-Schwarz inequality implies

Ω f (x)λ(ψ(x))|∇ψ(x)| dx ≤ f |∇ψ| 1 2 L 2 (Ω) λ • ψ |∇ψ| 1 2 L 2 (Ω) ,
and so the right-hand side of (3) is bounded in L 2 (0, 1) due to

λ • ψ |∇ψ| 1 2 2 L 2 (Ω) = Ω λ 2 (ψ(x))|∇ψ(x)| dx = 1 0 g 1,ψ (y)λ 2 (y)dy ≤ λ 2 L 2 (0,1) g 1,ψ L ∞ (0,1) .
Hence, the variational formulation (3) for g f,ψ ∈ L 2 (0, 1) is well-posed and the proof is complete.

We are going to use this formula to compute approximations of the geometric coefficients via a Galerkin projection, rather then introducing directly discretizations of (1).

We introduce P P (0, 1) the space of all polynomials of degree less or equal to P on [0, 1] and look for an approximation g P f,ψ ∈ P P (0, 1) of g f,ψ such that (4) where we have as a standard approximation result [28, Theorem 3.17]

(5) g P -g f,ψ L 2 (0,1) ≤ inf

v∈P P (0,1) v -g f,ψ L 2 (0,1) ≤ CP -k g f,ψ H k (0,1) ,
with C > 0, C ∈ R independent of P but dependent of k. The strength of this idea is that we can use the coarea formula (Corollary 2) to rewrite the integration over [0, 1] in the right hand side of (4) in an integration over Ω which, when working with the hp/spectral FEM approximation rather than with ψ, is our new approach.

Remark 3 (Continuity of geometric coefficients). To obtain in general geometric coefficients g f,ψ that are continuous in (0, 1) and so H k (0, 1) regular with k ≥ 1 2 that can be approximated in P P (0, 1) with a convergence rate larger than 1 2 we have to assume that ψ takes almost all values in [0, 1]. If ψ does not take values in a subinterval I ⊂ (0, 1) then the geometric coefficients are zero for all y ∈ I and might be discontinuous on ∂I. In this case the convergence rate would be below 1 2 .

DISCRETIZATION AND CONVERGENCE

Let V h,p ⊂ C 0 (Ω) be a finite element space defined over a triangulation of Ω, where h > 0 corresponds to the largest diameter of an element of the triangulation and p to the polynomial degree on the reference element. In the following we think of ψ h,p ∈ V h,p as approximation of ψ with ψ h,p : Ω → [0, 1] and we assume standard approximation results

ψ -ψ h,p L 2 (Ω) ≤ Ch p+1 and ψ -ψ h,p H 1 (Ω) ≤ Ch p .
So ψ h,p can be for example the L 2 or another kind of Galerkin projection of a function ψ. We refer to Sec. 4, where ψ h,p is the scaled to (0, 1) Galerkin projection of the solution ψ of a partial differential equation. We introduce the Galerkin projection of geometric coefficients: Find g P f h,p ,ψ h,p ∈ P P (0, 1) such that (7)

1 0 g P f h,p ,ψ h,p (y) λ(y)dy = Ω f h,p (x) λ(ψ h,p (x))|∇ψ h,p (x)|dx ∀λ ∈ P P (0, 1) .
We introduced here also f h,p (x) to cover cases where f depends on ψ and, hence, f h,p on its approximation such as f (x) = |∇ψ(x)| 2 , where f h,p (x) will be |∇ψ h,p (x)| 2 . The implementation of ( 7) compromises two steps/ingredients. Firstly, a standard (P + 1) × (P + 1) mass matrix containing the L 2 inner products of pairs for basis functions of P P (0, 1) has to be assembled. This mass matrix becomes diagonal when Legendre polynomialsmapped to (0, 1) -are used as basis functions. Second, P + 1 entries of the right hand side vector has to be assembled. As λ(ψ) is piecewise smooth with respect to the mesh underlying the definition of V h,p , these P + 1 entries can be approximated by the normal quadrature rules in hp/spectral FEM. In particular there is no need to determine explicitly the hypersurfaces of ψ h,p .

We have the following well-posedness and convergence results, the latter balancing the best approximation error inf v∈P P (0,1) v -g f,ψ L 2 (0,1) with the L 2 -error of ψ h,p .

Theorem 4 (Well-posedness and error of discretized geometric coefficients). Let the assumptions of Corollary 2 be fulfilled and let the sets

{x ∈ R n , ψ(x) = y} be (n -1)- dimensional hypersurfaces for all y ∈ [0, 1]. If moreover for f h,p : Ω → R ∪ {-∞, ∞} and ψ h,p ∈ V h,p it holds f h,p |∇ψ h,p | 1 2 ∈ L 2 (Ω) then the solution g P f h,p ,ψ h,p
of (7) is well-defined in L 2 (0, 1) and for any k > 0 there exits a constant C > 0 such that

g P f h,p ,ψ h,p -g f,ψ L 2 (0,1) ≤ C P -k g f,ψ H k (0,1) + (f |∇ψ| -f h,p |∇ψ h,p |)|∇ψ| -1 2 L 2 (Ω) + 8 √ 3 P 4 ψ -ψ h,p L 2 (Ω) f h,p |∇ψ h,p | L 2 (Ω) .
Proof. The well-posedness follows from the Lax-Milgram lemma, where the assumption on f h,p and ψ h,p leads similarly to the proof of Corollary 2 to g P f h,p ,ψ h,p ∈ L 2 (0, 1).

To prove the estimate on the discretization error we introduce first the notation

(f, λ • ψ) Ω = Ω f (x)λ(ψ(x))dx .
Then, the Strang lemma [5, Chap. III] yields

g P f h,p ,ψ h,p -g f,ψ L 2 (0,1) ≤ inf λ∈P P (0,1) λ -g f,ψ L 2 (0,1)
+ sup

λ∈P P (0,1) (f |∇ψ|, λ • ψ) Ω -(f h,p |∇ψ h,p |, λ • ψ h,p ) Ω λ L 2 (Ω)
and we estimate the numerator in the second term on the right hand side:

(f |∇ψ|, λ • ψ) Ω -(f h,p |∇ψ h,p |, λ • ψ h,p ) Ω ≤ (f |∇ψ| -f h,p |∇ψ h,p |, λ • ψ) Ω + (f h,p |∇ψ h,p |, λ • ψ -λ • ψ h,p ) Ω ≤ (f |∇ψ| -f h,p |∇ψ h,p |)|∇ψ| -1 2 L 2 (Ω) |∇ψ| 1 2 (λ • ψ) L 2 (Ω) + (f h,p |∇ψ h,p | • (ψ -ψ h,p ), λ • η) Ω (8)
where we used the Cauchy-Schwarz inequality and the mean value theorem, stating that there is some function η:

Ω → (0, 1) such that for all x ∈ Ω λ(ψ(x)) -λ(ψ h,p (x)) = λ (η(x))(ψ(x) -ψ h,p (x))
Inria and the fact that C ∞ (0, 1) is dense in L 2 (0, 1).

Using the Cauchy-Schwarz inequality we can bound the last term in ( 8) as

(f h,p |∇ψ h,p | • (ψ -ψ h,p ), λ • η) Ω ≤ sup y∈(0,1) |λ (y)| f h,p |∇ψ h,p | L 2 (Ω) ψ -ψ h,p L 2 (Ω) ,
where Sobolev embedding estimates and inverse estimates for polynomials lead to

sup y∈(0,1) |λ (y)| ≤ 8 √ 3P 4 λ L 2 (0,1) .
More precisely, it is easy to verify that the Sobolev embedding estimate [START_REF] Schwab | p-and hp-finite element methods : theory and applications in solid and fluid mechanics[END_REF]Theorem A.21]

λ L ∞ (0,1) ≤ 2 √ 3 λ H 1 (0,1)
holds, and with inverse estimates [28, Theorem 3.91] for λ ∈ P P (0, 1) we find

λ L 2 (0,1) ≤ 2 √ 3P 2 λ L 2 (0,1) , λ L 2 (0,1) ≤ 2 √ 3(P -1) 2 λ L 2 (0,1) ≤ 12P 2 (P -1) 2 λ L 2 (0,1) .
Hence, we can bound

λ 2 L ∞ (0,1) ≤ 4 3 ( λ 2 L 2 (0,1) + λ 2 L 2 (0,1) ) ≤ 192 P 8 λ 2 L 2 (0,1) .
Moreover, the assumptions of the theorem on f h,p and ψ h,p and V h,p ⊂ C 0 (Ω) imply f h,p |∇ψ h,p | ∈ L 2 (Ω). Finally, by the coarea formula (3)

|∇ψ| 1 2 (λ • ψ) 2 L 2 (Ω) = (|∇ψ|λ • ψ, λ • ψ) Ω = 1 0 λ 2 (y)g 1,ψ (y)dy.
and the assertion follows.

∈ R n , ψ(x) = y} are (n -1)- dimensional hypersurfaces for all y ∈ [0, 1] is not very critical. If for y 0 ∈ [0, 1] the set {x ∈ R n , ψ(x) = y 0 } is not a (n -1)-dimensional hypersurface, then the geomet- ric coefficient g f,ψ is infinite at y = y 0 . However, with the assumptions of Corollary 2 g f,ψ is bounded in L 2 (0, 1) even so if ψ is constant in each connected component of a n- dimensional open set Ω 0 ⊂ Ω. Note, that the term (f |∇ψ| -f h,p |∇ψ h,p |)|∇ψ| -1 2 L 2 (Ω)
in Theorem 4 can be replaced in this case by

(f |∇ψ| -f h,p |∇ψ h,p |)|∇ψ| -1 2 L 2 (Ω\Ω0) + 4P 2 f h,p |∇ψ h,p | L 1 (Ω0) .
As a simple corollary, inverse inequalities yield a convergence result for the error of the Galerkin approximation g P f h,p ,ψ h,p measured in the H 1 -seminorm and the L ∞ -norm.

Theorem 7 (H 1 -and L ∞ -error of the discretized geometric coefficients). Let the assumptions of Theorem 4 be fulfilled. Then, for any k > 0, there exists a constant C > 0 such that

|g P f h,p ,ψ h,p -g f,ψ | H 1 (0,1) ≤ CP 2 g P f h,p ,ψ h,p -g f,ψ L 2 (0,1) + P -k g f,ψ H k (0,1)
and the same bound holds for the error in the L ∞ (0, 1)-norm.

Proof. According to [START_REF] Schwab | p-and hp-finite element methods : theory and applications in solid and fluid mechanics[END_REF]Theorem 3.17] we have a projection operator I P : H 1 (0, 1) → P P (0, 1) such that for all v ∈ H k (0, 1)

|v -I P v| H 1 (0,1) ≤ CP 1-k v H k (0,1) , v -I P v L 2 (0,1) ≤ CP -k v H k (0,1) .
Combining these results with triangle inequalities and inverse inequality [START_REF] Schwab | p-and hp-finite element methods : theory and applications in solid and fluid mechanics[END_REF]Theorem 3.91] we obtain the estimate in the H 1 -seminorm as

|g P f h,p ,ψ h,p -g f,ψ | H 1 (0,1) ≤|g P f h,p ,ψ h,p -I P g f,ψ | H 1 (0,1) + |I P g f,ψ -g f,ψ | H 1 (0,1) ≤CP 2 g P f h,p ,ψ h,p -I P g f,ψ L 2 (0,1) + CP 1-k g f,ψ H k (0,1) ≤CP 2 g P f h,p ,ψ h,p -g f,ψ L 2 (0,1) + CP 1-k (1 + P ) g f,ψ H k (0,1) .
The L ∞ (0, 1)-estimate follows from the Sobolev embedding theorem, e.g. [START_REF] Schwab | p-and hp-finite element methods : theory and applications in solid and fluid mechanics[END_REF]Theorem A.21].

Remark 8. In the setting of (7) we know for λ ∈ P P (0, 1) and ψ h,p ∈ V h,p that also the restriction λ • ψ h,p|K to an element K of τ h is the image of a polynomial defined on the reference element, and hence the consistency error due to quadrature falls into the usual setting of variational crimes (see e.g. [19, pp. 155-166] and [6, pp. 201-203]).

Remark 9. It is straightforward to replace P P (0, 1) in [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF] with some other finite dimensional approximation space, such as spline spaces, for which we can prove a similar convergence results as in the preceding theorem. On the other hand it is a priori not obvious, whether standard quadrature on elements of the mesh τ h underlying the definition of V h,p is sufficient.

Remark 10. If one is only interested in the approximation in specific values of g f,ψ (y), e.g. at y = y c , y c ∈ (0, 1), as it would be the case for level set methods, fictitious domain method or extended finite element methods, it might be worth to introduce a smooth window function w ε : R → R with w ε (0) = 1 and compact support (-ε, ε). Clearly g f •wε(ψ-yc),ψ (y c ) = g f,ψ (y c ) and

1 0 g f •wε(ψ-yc),ψ (y)λ(y)dy = Ω f (x)w ε (ψ(x) -y c )λ(ψ(x))|∇ψ(x)|dx .
The integrand in the right hand side of the last equation vanishes in large parts of the domain Ω which reduces the amount of work required to assemble the right hand side of the Galerkin projection [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF]. Theorem 7 provides a convergence result. Alternatively, depending on the application and the concrete setting, it is also possible to pursue a more localized approach, where the global integral g f,ψ (y) is split into a sum of integrals

g f,ψ (y c ) = T g T f,ψ (y c ) = T {x∈Ω∩T,ψ(x)=yc} f (x)ds(x)
where the summation runs over the elements T of a non-overlapping decomposition of Ω.

Exemplary for this are the fictitious domain methods, where one needs to evaluate integrals g f,ψ for a large number of different f with local support. Applying our method to the approximation of g T f,ψ (y c ) replaces the technical integration on exact, but non-polynomial boundaries with quadrature on T , that is in general readily available. The convergence assertion of Theorems 4 and 7 control the approximation error and in the simplest case of exact f and ψ, this reduces to the standard error (5) of polynomial approximation in 1D. In the following we will consider a continuous, bounded scalar function ψ : Ω ⊂⊂ (R + × R) → R and introduce the normalization

ψ(r, z) = ψ(r, z) -ψ 0 (ψ) ψ 1 (ψ) -ψ 0 (ψ) ∈ (0, 1) .
with ψ 1 (ψ) = inf (r,z)∈Ω ψ(r, z), ψ 0 (ψ) := sup (r,z)∈Ω ψ(r, z). In our application, the axisymmetric magnetohydrodynamic equilibrium problem, ψ is the so-called poloidal flux and solution to the following elliptic quasi-linear partial differential equation:

(9) -∇ • 1 µ 0 r ∇ψ(r, z) = j p (r, ψ) in Ω ; ψ(r, z) = ψ bd at ∂Ω ;
with ψ bd ∈ R, ∇ the gradient in two dimensions and µ 0 the constant magnetic permeability of vacuum. The right hand side j p (r, ψ) is supposed to be known and positive, but in general non-linear in ψ. Hence, by the maximums principle we have ψ 1 = inf x∈Ω ψ(x) = ψ bd . The equation ( 9) is the celebrated Grad-Shafranov-Schlüter equation [START_REF] Grad | Hydromagnetic equilibria and force-free fields[END_REF][START_REF] Shafranov | On magnetohydrodynamical equilibrium configurations[END_REF][START_REF] Lüst | Axialsymmetrische magnetohydrodynamische Gleichgewichtskonfigurationen[END_REF], that is one of the central equations for modeling of magnetically confined plasma in tokamaks.

The domain Ω corresponds to a section through the torus-shaped domain covered by the plasma in the tokamak and it can be easily shown that ( 9) is the equlibrium problem of magnetohydrodynamics (MHD) in systems with axial symmetry. The MHD equlibrium is the combination of magnetostatics curl B = µ 0 J, div B = 0 and force balance grad p = J × B with B the magnetic induction, J the electric current density and p the hydrodynamic pressure. The source term j p and ψ/r in (9) are the axial components of J and A, where A is the electric vector potential with B = curl A. We refer to standard text books [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF][START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF][START_REF] Goedbloed | Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF][START_REF] Jardin | Computational methods in plasma physics[END_REF] for the details of the derivation. While ( 9) is usually referred to as the fixedboundary equilibrium problem as Ω is fixed, there is also the even more important freeboundary equilibrium problem that is formulated in an unbounded domain with boundary conditions at infinity and where the domain covered by the plasma is an unknown. Clearly, the numerical solution of free-boundary equilibrium problems faces more difficulties [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF] than the fixed-boundary equilibrium problem. But as the focus of this work is the approximation of the geometric coefficients we restrict the presentation to the fixed-boundary equilibrium problem while remarking that our new approach to compute geometric coefficients generalizes with minor modifications to the free-boundary equilibrium problem.

In the famous Grad/Hogan-approach [START_REF] Grad | Hydromagnetic equilibria and force-free fields[END_REF][START_REF] Grad | Classical diffusion in a tokamak[END_REF][START_REF] Hinton | Theory of plasma transport in toroidal confinement systems[END_REF] the transient MHD equations are reformulated as the so-called coupled problem of equilibrium and resistive diffusion and transport equations. While the equilibrium problem is the one we introduced in (9), the resistive diffusion and transport equations are basically Faraday's law and the hydrodynamic equations written in the curvi-linear coordinate system that is induced by the contour lines of the poloidal flux function ψ. The variation along the flux lines is neglected and the geometric coefficients [START_REF] Ph | Concepts -An Object-Oriented Software Package for Partial Differential Equations[END_REF] c f,ψ (y) = {(r,z)∈Ω,ψ(r,z)=y} f (r, z) |∇ψ(r, z)| rds that appear due to this coordinate transformation are related to the coefficients (1) for n = 2, ds(x) = rds and dx = rdrdz via

g f,ψ (y) = c f •|∇ψ|,ψ (y) .

RR n°8948

The unknowns of the resistive diffusion and transport equation on the other hand determine in some non-linear fashion the profile of the current density j p in the equilibrium problem [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF]. Due to the unusual non-linear coupling the equations in the Grad/Hogan-approach were dubbed queer differential equations. The Grad/Hogan-approach is the method of choice, when it comes to the simulation of plasma evolution on the longest timescale, the timescale of resistive diffusion, where entire discharge scenarios can be modeled. Devising stable algorithms is very challenging and an active area of research with many open problems. Numerical methods for computing accurately not only the geometric coefficients but also sensitivities are an important step towards stable and efficient methods for the coupled problem.

4.1. Galerkin formulation. We introduce the spaces

V := ψ : Ω → R, ψ L 2 (Ω) < ∞, |∇ψ| L 2 (Ω) < ∞ ∩ C 0 (Ω) V 0 := ψ ∈ V, ψ |∂Ω = 0 with ψ 2 L 2 (Ω) = Ω ψ 2 r -1 drdz and the mappings A : V × V → R and J p : V × V → R with (11) A(ψ, ξ) := Ω 1 µ 0 r ∇ψ • ∇ξ drdz and J p (ψ, ξ) := Ω j p (r, ψ)ξ drdz .
Then, we have the weak formulation of ( 9): find ψ-ψ 1 ∈ V 0 such that (12) A(ψ, ξ) = J p (ψ, ξ), ∀ξ ∈ V 0 .

In the general non-linear case it is convenient to use Newton's method to solve the nonlinear problem: the Newton update ψ n+1 is the solution of the infinite dimensional linear system

A(ψ n+1 , ξ) -D ψ J p (ψ n , ξ)(ψ n+1 ) = J p (ψ n , ξ) + D ψ J p (ψ n , ξ)(ψ n ), ∀ξ ∈ V 0 .
The derivative of J p follows from [4, Lemma I.4]: (13)

D ψ J p (ψ, ξ)( ψ) = Ω ∂j p (r, ψ) ∂ψ N ψ ξ ψ 1 -ψ 0 (ψ) drdz, + ψ(r 0 (ψ), z 0 (ψ)) Ω ∂j p (r, ψ) ∂ψ 1 -ψ ψ 1 -ψ 0 (ψ) ξ drdz .
We have implemented the Galerkin method for [START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF] in CONCEPTS [START_REF] Ph | Concepts -An Object-Oriented Software Package for Partial Differential Equations[END_REF][START_REF] Schmidt | Computation of the band structure of twodimensional photonic crystals with hp finite elements[END_REF][START_REF] Wang | On highorder FEM applied to canonical scattering problems in plasmonics[END_REF] (www. concepts.math.ethz.ch), using finite dimensional hp-FEM approximation spaces V h,p for V . CONCEPTS is a hp-FEM library that is based on basis functions defined on rectangular reference domains. Transfinite interpolation techniques ensure that boundaries of the mapped reference elements coincide with the (curved) boundaries of the domain.

For the most general case rigorous theoretical convergence results are not yet available, mainly due the j p depending on the normalized poloidal flux. If j p is independent of the normalized poloidal flux ψ we end up with a linear problem and have the usual convergence assertion.

Theorem 11. Let ψ be the solution of (12) and ψ h,p ∈ ψ h,p ∈ V h,p ∩ V 0 the solution of its Galerkin discretization satisfying

(14) A(ψ h,p , ξ) = J p (ψ h,p , ξ) + A(ψ 1 , ξ), ∀ξ ∈ V h,p ∩ V 0 .
For j p (r, ψ) independent of ψ we have

ψ -ψ h,p L 2 (Ω) ≤ Ch p+1 and ψ -ψ h,p H 1 (Ω) ≤ Ch p .
with constants C > 0 independent of h and p.

Inria

A High Order Method for Contour Integrals with an Application to Plasma Modeling in Nuclear Fusion11

Moreover we have the pointwise estimates for any x ∈ Ω

(15) (ψ -ψ h,p )(x) ≤ C h p+1 log(h -1 ) p = 1 h p+1 p > 1 .
Proof. Convergence estimates in L 2 and H 1 are standard. The pointwise estimate follow as in [START_REF] Wahlbin | Local behavior in finite element methods[END_REF]Theorem 11.1] The Galerkin projection of geometric coefficients for this application is: Find c P f h,p ,ψ h,p ∈ P P (0, 1) such that ( 16)

1 0 c P f h,p ,ψ h,p (
y) λ(y)dy = Ω f h,p (r, z) λ(ψ h,p (r, z))r dr ds ∀λ ∈ P P (0, 1) .

Similar as in Theorem 4 we can deduce a convergence result, that balances in this case the best approximation error inf v∈P P (0,1) v -c f,ψ L 2 (0,1) with the L 2 -error of ψ h,p and the error of the maximal ψ h,p -value.

Theorem 12. Let the assumptions of Corollary 2 be fulfilled, let the sets {x ∈ R n , ψ(x) = y} be (n -1)-dimensional hypersurfaces for all y ∈ [0, 1] and ψ h,p be the solution of (14) that converges pointwise to ψ the solution of [START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF].

If moreover for f h,p : Ω → R ∪ {-∞, ∞} and ψ h,p ∈ V h,p it holds f h,p |∇ψ h,p | -1 2 ∈ L 2 (Ω) then the solution c P f h,p ,ψ h,p
of ( 16) is well-defined in L 2 (0, 1) and for any k > 0 there exists a constant C > 0 such that, asymptotically,

c P f h,p ,ψ h,p -c f,ψ L 2 (0,1) ≤ C(P -k c f,ψ H k (0,1) + f -f h,p L 2 (Ω) ) + CP 4 ψ -ψ h,p L ∞ (Ω) .
Proof. From Theorem 4 we have

c P f h,p ,ψ h,p -c f,ψ L 2 (0,1) ≤ C(P -k c f,ψ H k (0,1) + f -f h,p L 2 (Ω) ) + CP 4 ψ -ψ h,p L 2 (Ω) .
We note

ψ -ψ h,p L 2 (Ω) ≤ ψ -ψ h,p + ψ 0 (ψ) -ψ 0 (ψ h,p ) ψ 1 -ψ 0 (ψ) L 2 (Ω) + ψ 1 (ψ h,p ) -ψ 0 (ψ h,p ) ψ 1 -ψ 0 (ψ) -1 ψ h,p L 2 (Ω) ,
from which we have

1 - ψ 1 (ψ h,p ) -ψ 0 (ψ h,p ) ψ 1 -ψ 0 (ψ) -1 ψ -ψ h,p L 2 (Ω) ≤ ψ -ψ h,p + ψ 0 (ψ) -ψ 0 (ψ h,p ) ψ 1 -ψ 0 (ψ) L 2 (Ω) + ψ 1 (ψ h,p ) -ψ 0 (ψ h,p ) ψ 1 -ψ 0 (ψ) -1 ψ L 2 (Ω) .
As for h small enough or p large enough

ψ 1 (ψ h,p ) -ψ 0 (ψ h,p ) ψ 1 -ψ 0 (ψ) -1 < 1, we conclude ψ -ψ h,p L 2 (Ω) ≤ C( ψ -ψ h,p L 2 (Ω) + |ψ 0 (ψ) -ψ 0 (ψ h,p )|) ≤ C ψ -ψ h,p L ∞ (Ω) ,
and the assertions follow. Again we can deduce a convergence result for the error measured in the H 1 -norm.

Theorem 14. Let the assumptions of Theorem 12 be fulfilled. Then, for any k > 0, there exists a constant C > 0 such that

|c P f h,p ,ψ h,p -c f,ψ | H 1 (0,1) ≤ CP 1-k (1 + P ) c f,ψ H k (0,1) + CP 2 c P f h,p ,ψ h,p -c f,ψ L 2 (0,1) .

4.2.

Sensitivities of geometric coefficients. In many practical application it is useful to compute not only the geometric coefficients but also the sensitivities with respect to the poloidal flux map ψ in some direction δψ:

(17) D ψ c f,ψ (y), δψ = lim ε→0 c f,ψ+εδψ (y) -c f,ψ (y) ε .
Examples are Newton-type methods for the coupled problem of equilibrium and resistive diffusion and transport. Similarly, in optimization problems for the coupled problem such sensitivities are required. By the coarea formula we find as in Corollary 2 the following characterization of the sensitivities of the geometric coefficients:

(18) 1 0 D ψ c f,ψ (y), δψ λ(y)dy = L f,ψ (λ, δψ) ∀λ ∈ L 2 (0, 1) , with L f,ψ (λ, δψ) = Ω f (r, z) λ (ψ(r, z)) δψ(r, z) ψ 1 (ψ) -ψ 0 (ψ) r drdz - Ω f (r, z) λ (ψ(r, z)) 1 -ψ(r, z) ψ 1 (ψ) -ψ 0 (ψ) δψ(r 0 (ψ), z 0 (ψ)) r drdz ,
where we used the non-trivial derivative

D ψ ψ 0 (ψ)(δψ) = δψ(r 0 (ψ), z 0 (ψ)) .
The functional (r 0 (ψ), z 0 (ψ)) = arg sup (r,z)∈Ω ψ(r, z) is the location of the supremum. Likewise as for the geometric coefficients we have a simple Galerkin method for computing approximations of the sensitivities. Galerkin projection of sensitivities of geometric coefficients Find Dc P f h,p ,ψ h,p (δψ h,p ), ∈ P P (0, 1) such that (19)

1 0 Dc P f h,p ,ψ h,p (δψ h,p )(y) λ(y)dy = L f ,ψ h,p (λ, δψ h,p ) ∀λ ∈ L 2 (0, 1) .
And again the method can be implement with standard quadrature rules on elements of the mesh underlying the definition of V h,p . If the same quadrature rules are used on ( 7) and [START_REF] Knabner | Numerical methods for elliptic and parabolic partial differential equations[END_REF], then the discretization of the continuous sensitivity is identical to the sensitivity of the discretization of the geometric coefficient. This will be different, when free-boundary problems are considered and we refer to [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF] for the details on quadrature and sensitivities for intersection of triangles and level sets. The assumption c f,ψ ∈ H k (0, 1) in Theorem 4 is fairly strong, seen the definition [START_REF] Ph | Concepts -An Object-Oriented Software Package for Partial Differential Equations[END_REF]. First of all it seems to be obvious that c f,ψ looses regularity in the neighborhood of vaules y ∈ [0, 1] where ψ attains a critical point. Second, we know that solutions of elliptic problems in non-smooth domains loose regularity in the vicinity of corners of the boundary of the domain and it can be expected that this loose of regularity influences the regularity of the geometric coefficients. While precise characterization of the regularity of solutions of the non-linear elliptic problem ( 9) in domains with corners can be found in [START_REF] Beretta | An inverse problem originating from magnetohydrodynamics. II. The case of the Grad-Shafranov equation[END_REF][START_REF] Beretta | An inverse problem originating from magnetohydrodynamics. III. Domains with corners of arbitrary angles[END_REF], all the test cases and examples in fusion engineering literature are arbitrarily smooth. Moreover it is an open problem how the regularity of ψ influences the regularity of c f,ψ . There is also very few theory available explaining how critical points of ψ influence the regularity of c f,ψ . The only result, we are aware of, is due to G. Vigfusson, stated without proof in [32, Theorem 2].

Theorem 15 (G. Vigfusson 1979). Assume that ψ has only one critical point, that is the maximum of ψ in Ω and

ψ |∂Ω ≡ ψ 1 ∈ R, then c f,ψ ∈ C n ([0, 1]), if f ∈ C 2n (Ω) and ψ ∈ C 2n+2 (Ω) .
Proof. The proof for n = 0 and n = 1 follows from various results in [START_REF] Vigfússon | The averaged Green function with applications to quasi-static plasma equilibrium[END_REF] and is based on a technique called elliptic expansions that compares the contour lines of ψ with ellipses and quantifies the deviation. The Theorems 9 and 10 in [31, p. 62] provide

(20) C 1,ψ (•) := • 0 c 1,ψ (y )dy ∈ C n+1 ([0, 1]) if ψ ∈ C 2n+2 (Ω) ,
with n = 0, n = 1 or n = 2. Furthermore, another theorem in [31, p. 75] states

(21) c f,ψ c 1,ψ ∈ C n ([0, 1]), if f ∈ C 2n (Ω) and ψ ∈ C 2n+1 (Ω) ,
with n = 0 or n = 1 and the assertion of Theorem 15 follows.

Such a smoothness result is not true if ψ has saddle points in Ω as we see for the example in the following remark.

Remark 16. We consider Ω = [0, 1] 2 , ψ(r, z) = 1 -rz and f (r, z) = i,j≥0 a ij r i z j with a ij ∈ R. ψ(r, z) has a saddle point at (r X , z X ) = (0, 0) with y X = 1. We compute for

f 1 (r, z) = r n z m , n, m ∈ N, n = m: C f1,ψ (y) = 1 1-y 1 (1-y)/z r n z m dr dz = (m -n) -(m + 1)(1 -y) n+1 + (n + 1)(1 -y) m+1 (m + 1)(n + 1)(m -n)
and for f 2 (r, z) = r n z n , n ∈ N:

C f2,ψ (y) = 1 1-y 1 (1-y)/z r n z n dr dz = 1 -(1 -y) n+1 (n + 1) 2 + (1 -y) n+1 log(1 -y) n + 1 , which yields C f1,ψ (y) = c f1,ψ (y) = (1 -y) n -(1 -y) m m -n . and C f2,ψ (y) = c f2,ψ (y) = -log(1 -y)(1 -y) n .
Hence we find

c f,ψ (y) = i,j≥0 i =j a ij (1 -y) i -(1 -y) j j -i -log(1 -y) i≥0 a ii (1 -y) i .
In the neighborhood of the saddle point y = y X (= 1) the geometric coefficients c f,ψ (y) behave like log(y X -y) whenever a 00 = 0. In particular, c f,ψ (y X ) is unbounded, while the quotient c f,ψ (y X )/c 1,ψ (y X ) is bounded:

lim y→y X c f,ψ (y) c 1,ψ (y) = f (r X , z X ) = a 00 .
4.4. Numerical Experiment: Elliptic Equilibrium. We consider an elliptic equilibrium described by ( 22)

ψ E (r, z) = ψ 0 - (r -r 0 ) 2 a 2 + (z -z 0 ) 2 b 2 ,
with ψ 1 = 0. We have an explicit parametrization for the iso-surfaces:

{ r z ; ψ(r, z) = s} = { r E (t) z E (t) = r 0 + a ψ 0 -s cos(t) z 0 + b ψ 0 -s sin(t) ; 0 ≤ t ≤ 2π} .
and the right hand side for the equation ( 9) evaluates to:

-∇ • 1 µ 0 r ∇ψ E (r, z) = 1 µ 0 2 a 2 r + 2 b 2 r - 2(r -r 0 ) a 2 r 2 .
In the subsequent calculations we are using always ψ 0 = 2, (r 0 , z 0 ) = (2, 0) and a = 1, b = 3.

4.4.1.

Convergence of ψ E . To begin with we are using the Galerkin method [START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF] to determine an approximation ψ h,p of ψ E in V h,p . As this example case is linear the standard analysis for elliptic problems asserts convergence estimates [START_REF] Lütjens | The CHEASE code for toroidal MHD equilibria[END_REF] ψ 4.4.2. Convergence of the the geometric coefficients. Moving on to the geometric coefficients we find using

E -ψ h,p L 2 (Ω) = O(h p+1 ) and ψ E -ψ h,p H 1 (Ω) = O(h p ) ,
|∇ψ E | 2 = (ψ 0 -s)(b 2 cos(t) 2 + a 2 sin(t) 2 ) 4 a 2 b 2 , ṙE (t) 2 + żE (t) 2 = (ψ 0 -s)(b 2 cos(t) 2 + a 2 sin(t) 2 ) , the following analytic values c 1 r 2 ,ψ E (y) = ψ E =y r -1 |∇ψ E | ds = 2π 0 ṙ(t) 2 + ż(t) 2 r(t)|∇ψ E | dt = |ψ 0 |abπ r 2 0 + a 2 ψ 0 y c |∇ψ E | 2 r 2 ,ψ E (y) =|ψ 0 | ψ E =y |∇ψ E | r ds = 2π 0 ṙ(t) 2 + ż(t) 2 |∇ψ E | |∇ψ E | 2 r dt = 4π|ψ 0 | (b 2 -a 2 )r 0 (r 0 -r 2 0 + a 2 ψ 0 y) -a 4 ψ 0 y a 3 b r 2 0 + a 2 ψ 0 y .
Fixing the polynomial degree P of approximation of the geometric coefficients the Theorems 12 and 14 assert 10 -1

c P f h,p ,ψ h,p -c f,ψ L 2 (0,1) c P f h,p ,ψ h,p -c f,ψ H 1 (0,1) = O(h p+1 log(h -1 )) p = 1 O(h p+1 ) p > 1 ,
refinement level l (hl = 2 -l h0) p = 1 p = 2 p = 3 p = 4 p = 5 p = 6 O(h p l ) FIGURE 1. Elliptic equilibrium: Relative L 2 -error E 0 (l, p) := ψ h l ,p - ψ E L 2 (Ω) / ψ E L 2 (Ω) (top) and relative H 1 -error E 1 (l, p) := ψ h l ,p - ψ E H 1 (Ω) / ψ E H 1 (Ω) (bottom)
f = f 1 = 1/r 2 and f = f 2 = |∇ψ| 2 /r 2 .
for the Galerkin projection c P f h,p ,ψ h,p with f h,p = f = 1/r 2 and

c P f h,p ,ψ h,p -c f,ψ L 2 (0,1) c P f h,p ,ψ h,p -c f,ψ H 1 (0,1) = O(h p ) , with f = |∇ψ| 2 /r 2 and f h,p = |∇ψ h,p | 2 /r 2 .
In the following experiments we fix P = 15. The convergence for c 1 r 2 ,ψ E is in general better than predicted by Theorem (12) (see Figure 3), when the polynomial degree is larger than 1. The experimental convergence rates seem to indicate that it might be possible to FIGURE 3.

Elliptic equilibrium: Relative L 2 -error e 0 (l, p) := c P f,ψ h l ,p -c f,ψE L 2 (0,1) / c f,ψE L 2 (0,1) and H 1 -error e 1 (l, p) := c P f,ψ h l ,p -c f,ψE H 1 (0,1) / c f,ψE H 1 (0,1) for f = 1/r 2 for refinement in polynomial degree p (left) and refinement level l of the mesh (right) with P = 15.

improve the assertion of Theorem ( 12) by a duality argument: The experimental rate of convergence for refinement in the mesh for fixed polynomial degree p > 2 seems to be twice as large as predicted by our theory. Moreover (see Figure 5), also for

c |∇ψ E | 2 r 2
,ψ E the convergence is better than predicted and except for the lowest order case the additional consistency error due to f h,p does not seem to spoil the rate of convergence.

The experimental results for the error measured in the H 1 -norm (see Figures 3 and5) highlights even more that our new method for computing the geometric coefficients is most powerful, when combined with high order polynomial approximation of ψ. In this case the approximation error in ψ decays very quickly and the positive powers of P in front of this approximation error in Theorems ( 12) and ( 14) do not harm too much.

Finally, our experiments suggest (see Figures 4 and6) also that the high powers of the discretization parameter P in the estimates of Theorem 12 and 14 are too pessimistic. We observe in these examples a convergence like P 1 2 rather than P 4 for the error in L 2 (0, 1). For the error in H 1 (0, 1) we see a convergence like P 5 2 rather than P 6 . Additionally we see in the experiments the saturation levels for small P which highlight that the discretization parameter P needs to be balanced accordingly with h and p, to achieve a desired accuracy level with minimal effort. 4.5. Numerical Experiment: Solov'ev Equilibrium. The next example, taken from [START_REF] Lütjens | The CHEASE code for toroidal MHD equilibria[END_REF], is a very popular benchmarking case in computational plasma physics. We take Elliptic equilibrium: Relative L 2 -error e 0 (l, 1) := c P f,ψ h l ,1 -c f,ψE L 2 (0,1) / c f,ψE L 2 (0,1) (top left) and H 1 -error e 1 (l, 1) := c P f,ψ h l ,1 -c f,ψE H 1 (0,1) / c f,ψE H 1 (0,1) (top right) for f = 1/r 2 for refinement in polynomial degree P and refinement level l of the mesh and relative L 2 -error e 0 (l, 1) and H 1 -error e 1 (l, 1) at fixed level l = 5 for refinement in polynomial degree P (bottom). such that it is the solution to [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF] with right hand side j p (r, ψ S ) = 1 + κ 2 µ 0 κr 3 0 q 0 r where r 0 and q 0 are the major radius and the safety factor at the magnetic axis, and a and κ are the effective minor radius and elongation of the last closed flux surface, given by ψ = 0. A parametrization of the boundary is ∂Ω = {(r, z) : r 2 = r 2 0 + 2ar 0 cos(t), z = κa r 0 r sin(t), 0 ≤ t ≤ 2π}

In the following example we will chose the same parameter values r 0 = 1, a = 0.32, κ = 1.7 and q 0 = 1.0 as in [START_REF] Lee | Ecom: A fast and accurate solver for toroidal axisymmetric MHD equilibria[END_REF].

4.5.1. Convergence of ψ S . As in the case of the elliptic equilibrium our experiments (see Figure 7) reproduce the theoretical convergence results [START_REF] Lütjens | The CHEASE code for toroidal MHD equilibria[END_REF]. The numerical solutions for the first 3 mesh refinements are shown in Figure 8.

4.5.2. Convergence of the geometric coefficients. For this example we do not have analytic values of the various geometric coefficients c f,ψ with f = 1/r 2 and f = |∇ψ| 2 /r 2 . Hence, we explore the convergence behavior using a reference c ref f,ψS := c 15 f,ψ h 2 ,16 computed at a discretization with P = 15 for the 1D space and p = 16 for the 2D space on the mesh of level l = 2.

Here again we observe spectral convergence for geometric coefficients c f,ψ that involve weighting functions f independent of ψ (see Figure 10) or weighting functions that are dependent of ψ (see Figure 11). Elliptic equilibrium: Relative L 2 -error e 0 (l, 1) := c P f,ψ h l ,1 -c f,ψE L 2 (0,1) / c f,ψE L 2 (0,1) (top left) and H 1 -error e 1 (l, 1) := c P f,ψ h l ,1 -c f,ψE H 1 (0,1) / c f,ψE H 1 (0,1) (top right) for f = |∇ψ E | 2 /r 2 and f h l ,p = |∇ψ h l ,p | 2 /r 2 for refinement in polynomial degree P and refinement level l of the mesh and relative L 2 -error e 0 (l, 1) and H 1 -error e 1 (l, 1) at fixed level l = 5 for refinement in polynomial degree P (bottom).

  ψ (y) λ(y)dy ∀λ ∈ P P (0, 1) ,

1 0g

 1 P (y) λ(y)dy = Ω f (x) λ(ψ(x)) |∇ψ| dx ∀λ ∈ P P (0, 1) ,

Remark 5 .

 5 Our numerical experiments in Section 4 indicate that the high powers of P in the bound of Theorem 4 are by far too pessimistic. Remark 6. The assumption in Theorem 4, that the sets {x

Inria

4 .

 4 APPLICATION: THE AXISYMMETRIC PLASMA EQUILIBRIUM PROBLEM

Remark 13 .

 13 The high powers of P in the bound of Theorem 12 are a consequence of the estimates of Theorem 4. Once again, our numerical experiments in Section 4 indicate they are too pessimistic.

Inria 4 . 3 .

 43 A few regularity results for geometric coefficients.

  which are clearly reproduced by our experiments (see Figure 4.4.1).

κ 2 r 2 z 2 -a 2 r 2 0 1

 221 FIGURE 4.Elliptic equilibrium: Relative L 2 -error e 0 (l, 1) := c P f,ψ h l ,1 -c f,ψE L 2 (0,1) / c f,ψE L 2 (0,1) (top left) and H 1 -error e 1 (l, 1) := c P f,ψ h l ,1 -c f,ψE H 1 (0,1) / c f,ψE H 1 (0,1) (top right) for f = 1/r 2 for refinement in polynomial degree P and refinement level l of the mesh and relative L 2 -error e 0 (l, 1) and H 1 -error e 1 (l, 1) at fixed level l = 5 for refinement in polynomial degree P (bottom).

  FIGURE 6. Elliptic equilibrium: Relative L 2 -error e 0 (l, 1) := c P f,ψ h l ,1 -c f,ψE L 2 (0,1) / c f,ψE L 2 (0,1) (top left) and H 1 -error e 1 (l, 1) := c P f,ψ h l ,1 -c f,ψE H 1 (0,1) / c f,ψE H 1 (0,1) (top right) for f = |∇ψ E | 2 /r 2 and f h l ,p = |∇ψ h l ,p | 2 /r 2 for refinement in polynomial degree P and refinement level l of the mesh and relative L 2 -error e 0 (l, 1) and H 1 -error e 1 (l, 1) at fixed level l = 5 for refinement in polynomial degree P (bottom).

  Une Nouvelle Méthode d'ordre élevé pour des Intégrales de surface et curviligne Résumé : Nous présentons une nouvelle méthode pour calculer des approximations des intégrales de surface et des intégrales curvilignes. Cette méthode est basée sur la formule de co-aire combinée avec une projection de type Galerkin. Elle entre directement dans le cadre des méthodes d'éléments finis d'ordre élevé. Nous démontrons des théorèmes de convergence et les résultats théoriques sont illustrés par des exemples numériques venant d'une application en fusion nucléaire.

	A HIGH ORDER METHOD FOR CONTOUR INTEGRALS WITH AN
	APPLICATION TO PLASMA MODELING IN NUCLEAR FUSION
	HOLGER HEUMANN, LUKAS DRESCHER, AND KERSTEN SCHMIDT
	Mots-clés : intégrale de surface, intégrale de curviligne, intégration numérique, éléments
	finis d'ordre élevé
		1. INTRODUCTION
	This article is concerned with the approximation of the integrals
	(1)	g f,ψ (y) :=

{x∈Ω, ψ(x)=y}

RR n°8948

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

refinement level l (hl = 2 -l h0)

Elliptic equilibrium: Relative L 2 -error e 0 (l, p) := c P f h l ,p ,ψ h l ,p -c f,ψE L 2 (0,1) / c f,ψE L 2 (0,1) and H 1 -error e 1 (l, p)

r 2 for refinement in polynomial degree p (left) and refinement level l of the mesh (right) with P = 15.