Keywords: Anisotropic elliptic problem, Variational formulation, Asymptotic-Preserving Scheme, Plasmas, Euler equations, Lorentz force, Large Magnetic field, Low-Mach number, Drift-fluid limit 1 AMS subject classification: 65N06, 65N12, 65M06, 65M12, 76W05, 76X05, 76N17

This paper is devoted to the numerical approximation of a degenerate anisotropic elliptic problem. The numerical method is designed for arbitrary space-dependent anisotropy directions and does not require any specially adapted coordinate system. It is also designed to be equally accurate in the strongly and the mildly anisotropic cases. The method is applied to the Euler-Lorentz system, in the drift-fluid limit. This system provides a model for magnetized plasmas.

Degenerate anisotropic elliptic problems and magnetized plasma simulations.

Stéphane Brull (1) , Pierre Degond (2,3) , Fabrice Deluzet (2,3) 1-Institut de Mathématiques de Bordeaux UMR 5251 Equipe Mathématiques Appliquées de Bordeaux (MAB)

Université Bordeaux 1 351, cours de la Libération -33405 TALENCE cedex FRANCE email: Stephane.Brull@math.u-bordeaux1.fr 2-Université de Toulouse; UPS, INSA, UT1, UTM ; Institut de Mathématiques de Toulouse ; F-31062 Toulouse, France. email: fabrice.deluzet@math.univ-toulouse.fr

Introduction

This paper discusses the numerical resolution of degenerate anisotropic elliptic problems of the form:

-(b • ∇) (∇ • (b φ ε )) + εφ ε = f ε , in Ω, (1.1) (b • ν) ∇ • (b φ ε ) = 0 on ∂Ω, (1.2)
where Ω ⊂ R 2 or R 3 , f ε is a given function, b is a normalized vector field defining the anisotropy direction and ε measures the strength of this anisotropy. In this expression ∇ and ∇• are respectively the gradient and divergence operators. The unit outward normal at x ∈ ∂Ω is denoted by ν. In the context of plasmas, ε is related to the gyro period (i.e. the period of the gyration motion of the particles about the magnetic field lines), and the anisotropy direction b satisfies b = B/|B| with the magnetic field B verifying ∇ • B = 0. Eq. (1.1) may also arise in other contexts, such as rapidly rotating flows, shell theory and may also be found when special types of semi-implicit time discretization of diffusion equations are used. The elliptic equation is not in the usual divergence form due to an exchange between the gradient and divergence operators. However, the methodology would apply equally well to the operator ∇ • ((b ⊗ b) • ∇φ)), up to some simple changes. The expression considered here is motivated by the application to the Euler-Lorentz system of plasmas. This application has already been considered in a previous study [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF] but we introduce two important developments. First the present numerical method does not request the development of a special coordinate system adapted to b. In [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF], b was assumed aligned with one coordinate direction. Second, the present paper considers Neumann boundary conditions instead of Dirichlet ones as in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF]. Although seemingly innocuous, this change brings in a considerable difficulty, linked with the degeneracy of the limit problem, as explained below.

A classical discretization of problem (1.1), (1.2) leads to an ill-conditioned linear system as ε → 0. Indeed setting formally ε = 0 in (1.1), (1.2), we get:

-(b • ∇)∇ • (b ψ) = f (0) , in Ω, (1.3) (b • ν)∇ • (b ψ) = 0 , on ∂Ω, (1.4)
with f (0) = lim ε→0 f ε . The homogeneous system associated to (1.3), (1.4) admits an infinite number of solutions, namely all functions ψ satisfying ∇•(bψ) = 0. This degeneracy results from the Neumann boundary conditions (1.4) and would also occur if periodic boundary conditions were used. On the other hand, (1.3) is not degenerate if supplemented with Dirichlet or Robin conditions, which was the case considered in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF]. A standard numerical approximation of (1.3), (1.4) generates a matrix whose condition number blows up as ε → 0, leading to very time consuming and/or poorly accurate solution algorithms.

To bypass these limitations, we follow the idea introduced in [START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF] and use a decomposition of the solution in its average along the b-field lines and a fluctuation about this average. This decomposition ensures an accurate computation of the solution for all values of ε. In [START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF], this decomposition approach was developed for a uniform b and a coordinate system with one coordinate direction aligned with b. To extend this approach to arbitrary anisotropy fields b, a possible way is to use an adapted curvilinear coordinate system with one coordinate curve tangent to b. This is the route followed by [START_REF] Besse | Three dimensional simulation of ionospheric plasma disturbances[END_REF], which proposes an extension of [START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF] in the context of ionospheric plasma physics, where the anisotropy direction is known analytically (given by the earth dipolar magnetic field). The approach developed here is different and aims at a method which does not request the generation of special curvilinear coordinates. Indeed, in the general case, computing such coordinates can be complex and costly, especially for time-dependent problems where b evolves in time.

For this purpose, we solve a variational problem for each of the terms of the decomposition. The main difficulty lies in the discretization of the functional spaces in which each component of the solution is searched. In the present paper, this difficulty is solved by introducing two kinds of variational systems, one corresponding to a second-order elliptic problem (for the average) and one, to a fourth order system (for the fluctuation). An alternative to this method is proposed in [START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF]. It avoids the resolution of a fourth-order problem at the price of the introdution of Lagrange multipliers which lead to a larger system. In the present paper, we design a method which breaks the complexity of the problem in smaller pieces and requires less computer ressources.

As an application of the method and a motivation for studying problem (1.3), (1.4), the drift-fluid limit of the isothermal Euler-Lorentz system is considered. These equations model the evolution of a magnetized plasma. In this case, the anisotropy direction is that of the magnetic field and the parameter ε is the reciprocal of the non dimensional cyclotron frequency. The drift-fluid limit ε → 0 of the Euler-Lorentz system is singular because the momentum equation becomes degenerate. In this paper, we propose a scheme able to handle both the ε ∼ 1 and ε ≪ 1 regimes, giving rise to consistent approximations of both the Euler-Lorentz model and its drift-fluid limit, without any constraint on the space and time steps related to the possible small value of ε. Schemes having such properties are referred to as Asymptotic-Preserving (AP) schemes. These schemes are particularly efficient in situations in which part of the simulation domain is in the asymptotic regime and part of it is not. Indeed, in most practical cases, the parameter ε assumes a local value which may change from one location to the next or which may evolve with time.

The usual approach for dealing with such occurences is through domain decomposition: the full Euler-Lorentz model is used in the region where ε = O(1) and the drift-fluid limit model is used where ε ≪ 1. There are several drawbacks in using this approach. The first one is the choice of the position of the interface (or cross-talk region), which can influence the outcome of the simulation. If the interface evolves in time, an algorithm for interface motion has to be devised and some remeshing must be used to ensure compatibility between the mesh and the interface, which requires heavy code developments and can be quite CPU time consuming. Determining the right coupling strategy between the two models can also be quite challenging and the outcome of the numerical simulations may also depend on this choice. Because these questions do not have straigthforward answers, domain decomposition strategies often lack robustness and reliability. Here, using the original model with an AP discretization method everywhere prevents from these artefacts and permits to use the same code everywhere for both regimes.

We conclude this introductory section by some bibliographical remarks. In magnetized plasma simulations, many works are based on the use of curvilinear coordinate systems where one of the coordinate curves is tangent to the magnetic field (see e.g. [START_REF] Miyamoto | Controlled fusion and plasma physics[END_REF], the gyro-kinetic and gyro-fluid developments [START_REF] Beer | Toroidal gyrofluid equations for simulations of tokamak turbulence[END_REF][START_REF] Dorland | Gyrofluid turbulence models with kinetic effects[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-lagrangian code[END_REF][START_REF] Hammett | Developments in the gyrofluid approach to tokamak turbulence simulations[END_REF] and the many attempts for generating specialized coordinate systems [START_REF] Beer | Field-aligned coordinates for nonlinear simulations of tokamak turbulence[END_REF][START_REF] Boozer | Establishment of magnetic coordinates for a given magnetic field[END_REF][START_REF] Hitchon | Flux coordinates and magnetic field structure. A guide to a fundamental tool of plasma theory[END_REF][START_REF] Dimits | Fluid simulations of tokamak turbulence in quasiballooning coordinates[END_REF][START_REF] Hamada | Hydromagnetic equilibria and their proper coordinates[END_REF][START_REF] Kaveeva | Poloidal and toroidal flows in tokamak plasma near magnetic islands[END_REF][START_REF] Stern | Geomagnetic Euler potentials[END_REF]). The present work, together with [START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF] is one of the very few attempts to design numerical methods free of the use of special coordinate systems (see also [START_REF] Ottaviani | An alternative approach to field-aligned coordinates for plasma turbulence simulations[END_REF]). The key idea behind this method is the concept of Asymptotic Preserving (AP) schemes as described above. APschemes have first been introduced by S. Jin [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF] in the context of diffusive limits of transport models. They have recently found numerous applications to plasma physics in relation e.g. to quasineutrality [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Crispel | An Asymptotic Preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | Asymptotic-Preserving Particle-In-Cell method for the Vlasov-Poisson system near quasineutrality[END_REF][START_REF] Degond | Numerical approximation of the Euler-Poisson-Boltzmann model in the quasineutral limit[END_REF][START_REF] Degond | Analysis of an Asymptotic Preserving scheme for the Euler-Poisson system in the quasineutral limit[END_REF] and strong magnetic fields [START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF] as well as to fluid-mechanical problems such as the small Mach-number limit of compressible fluids [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equation[END_REF]. Other applications of AP-schemes can be found in [START_REF] Buet | Diffusion limit of the Lorentz model: Asymptotic Preserving schemes[END_REF][START_REF] Buet | Asymptotic Preserving and positive schemes for radiation hydrodynamics[END_REF][START_REF] Carrillo | Simulation of fluid and particles flows: Asymptotic Preserving schemes for bubbling and flowing regimes[END_REF][START_REF] Filbet | A class of Asymptotic Preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Klar | An Asymptotic Preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF][START_REF] Lemou | A new Asymptotic Preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Mcclarren | The effects of slope limiting on Asymptotic-Preserving numerical methods for hyperbolic conservation laws[END_REF]. Numerical methods for anisotropic problems have been extensively studied in the literature using numerous techniques such as domain decomposition techniques [START_REF] Giraud | Schur complement preconditioners for anisotropic problems[END_REF][START_REF] Khoromskij | Robust Schur complement method for strongly anisotropic elliptic equations[END_REF], Multigrid methods, smoothers [START_REF] Llorente | Robust multigrid smoothers for three dimensional elliptic equations with strong anisotropies[END_REF], the hp-finite element method [START_REF] Melenk | hp-finite element methods for singular perturbations[END_REF]. However, these methods are based on a discretization of the anisotropic PDE as it is written. The method presented here as well as in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF][START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF][START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF] relies on a totally different concept, namely viewing the anisotropy as a singular perturbation and using Asymptotic-Preserving techniques. This paper is organized as follows. In section 2 the solution methodology for the degenerate anisotropic elliptic problem (1.1), (1.2) is detailed. Section 3 is devoted to the discretization strategy. In section 4 the drift-fluid limit of the isothermal Euler-Lorentz system is introduced. The AP-scheme is derived, giving rise to the anisotropic elliptic problem (1.1), (1.2). The numerical method for the anisotropic elliptic problem is validated in section 5. Finally a numerical application to the Euler-Lorentz system is given in section 6.

A decomposition method for degenerate anisotropic elliptic problems

We first present the methodology in the simpler case of a uniform b-field. The method will then be extended to an arbitrary b-field.

Overview of the method in the uniform b-field case

A two dimensional configuration is considered in this section, with the position variable (x, y) belonging to a square domain (x, y)

∈ Ω = [0, 1] × [0, 1] ⊂ R 2
. The b field is assumed uniform, equal to the unit vector pointing in the y direction. In this case, the singular perturbation problem (1.1), (1.2) reads:

εφ ε (x, y) - ∂ 2 ∂y 2 φ ε (x, y) = f ε (x, y), in ]0, 1[×]0, 1[, (2.1)
∂ ∂y φ ε (x, y) = 0, for y = 0 or y = 1. (2.2)
We assume that:

lim ε→0 1 ε 1 0 f ε (x, y) dy exists and is finite, ∀x ∈ [0, 1]. (2.3)
This framework is similar to [START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF]. Here, we recall the bases of the methodology. The problem is well posed for all ε > 0 but a standard discretization may lead to ill-conditionned matrices when ε ≪ 1. Indeed if ε is formally set to zero, we get the following degenerate problem

- ∂ 2 ∂y 2 ψ(x, y) = f (0) (x, y), in ]0, 1[×]0, 1[, ∂ ∂y ψ(x, y) = 0, for y = 0 or y = 1, (2.4)
assuming that f ε has the following expansion f ε = f (0) + εf (1) + o(ε). This system admits a solution under the compatibility condition On the other hand, the limit φ (0) = lim ε→0 φ ε is unique. Indeed, it is easy to see that the solution ψ of (2.4) such that 1 0 ψ(x, y) dy = 0 for all x ∈ [0, 1] is unique. Since φ (0) is a particular solution of (2.4), it can be written

φ (0) = ψ + ζ(x).
(2.5)

In order to determine ζ , we integrate (2.1) with respect to y and get

1 0 φ ε (x, y) dy = 1 ε 1 0 f ε (x, y) dy , (2.6)
Taking the limit ε → 0 in this equation and inserting (2.5), we get ζ(x) = 1 0 f (1) (x, y) dy, which determines φ (0) uniquely. Now, if a standard numerical method is applied to (2.1), (2.2), the resulting matrix will be close, when ε ≪ 1, to the singular matrix obtained from the discretization of (2.4). Therefore, its condition number will blow up as ε → 0, resulting in either low accuracy, or high computational cost. To overcome this problem, we decompose φ ε according to

φ ε = p ε + q ε , p ε (x) = 1 0 φ ε (x, y) dy, (2.7)
i.e. p ε is the average of φ ε along straight lines parallel to b and q ε is the fluctuation of the solution with respect to this average. p ε and q ε satisfy:

∂p ε ∂y (x, y) = 0, ∀(x, y) ∈ Ω, (2.8
)

1 0 q ε (x, y) dy = 0, ∀x ∈ [0, 1]. (2.9)
They are orthogonal for the scalar product of L 2 , i.e. Ω p ε q ε dx dy = 0. Inserting this decomposition into (2.6) yields

p ε (x) = 1 ε 1 0 f ε (x, y) dy, ∀x ∈ [0, 1]. (2.10) Moreover, p ε satisfies lim ε→0 p ε (x) = lim ε→0 1 ε 1 0 f ε (x, y) dy = 1 0 f (1) (x, y) dy = ζ(x),
where ζ is defined by (2.5). Now, q ε is the solution of the following problem:

- ∂ 2 ∂y 2 q ε (x, y) + εq ε (x, y) = ξ ε (x, y), ∀(x, y) ∈ [0, 1]×]0, 1[, (2.11) 1 0 q ε (x, y) dy = 0, for x ∈ [0, 1], (2.12)
∂ ∂y q ε (x, y) = 0, for y = 0 or y = 1, (

where

ξ ε = f ε - 1 0 f ε dy = f ε -εp ε ,
is the projection of f ε on the space of functions satisfying (2.9). Compared to (2.1), (2.2), system (2.11)-(2.13) involves the additional condition (2.12). This condition is important: it makes the system uniformly well-posed when ε → 0. Additionally, the limit system is

- ∂ 2 ∂y 2 q (0) (x, y) = f (0) , ∀(x, y) ∈ [0, 1]×]0, 1[, (2.14)
1 0 q (0) (x, y) dy = 0, for x ∈ [0, 1], (2.15)
∂ ∂y q (0) (x, y) = 0, for y = 0 or y = 1, (

and has a unique solution equal to ψ. Consequently, as ε → 0

φ ε = p ε + q ε → ζ + ψ = φ (0) .
Therefore, the proposed decomposition leads to two uniformly well-posed problems when ε → 0, which allows to reconstruct the limit solution φ (0) of the original problem.

The numerical approximations of conditions (2.10) or (2.12) is delicate if the mesh is not aligned with the y coordinate axis. In order to overcome this problem, a weak formulation is introduced. Define

V = H 1 (0, 1), K = {v ∈ V | ∂ y v = 0}. Then, φ ε is the solution of the variational formulation Find φ ε ∈ V such that Ω ∂φ ε ∂y ∂ψ ∂y dx dy + ε Ω φ ε ψ dx dy = Ω f ε ψ dx dy, ∀ψ ∈ V. (2.17)
Let K ⊥ be the orthogonal space to K in L 2 (0, 1). Now, the decomposition (2.7), corresponds to the decomposition of φ ε on K and K ⊥ . Indeed, it is easily checked that p ε ∈ K and q ε ∈ K ⊥ and they are orthogonal, as already noticed. Now, inserting ψ ∈ K in (2.17), we get that p ε is the solution of

Find p ε ∈ K such that Ω (p ε - 1 ε f ε ) ψ dx dy = 0, ∀ψ ∈ K, (2.18)
which means that p ε is the orthogonal projection of

ε -1 f ε onto K. Now, inserting ψ ∈ K ⊥ in (2.17) leads to Find q ε ∈ K ⊥ such that Ω ∂q ε ∂y ∂ψ ∂y dx dy + ε Ω q ε ψ dx dy = Ω (f ε -εp ε ) ψ dx dy, ∀ψ ∈ K ⊥ , (2.19)
which is the variational formulation of (2.11)-(2.13). The use of these variational formulations allows for the discretization of (2.1), (2.2) on arbitrary meshes compared to the anisotropy direction. This is an important advantadge over the strong formulations (2.10) or (2.11)-(2.13). These formulations are now generalized to arbitrary anisotropy fields b in the next section.

Presentation of the method for a general anisotropy field 2.2.1 Preliminaries

This subsection is devoted to the resolution of degenerate elliptic problems (1.1), (1.2) for general anisotropy fields b. we first introduce the space

V = {φ ∈ L 2 (Ω) /∇ • (bφ) ∈ L 2 (Ω)}, K = {φ ∈ V / ∇ • (bφ) = 0 on Ω}, W = {h ∈ L 2 (Ω) /(b • ∇)h ∈ L 2 (Ω)}, W 0 = {h ∈ W / (b • ν)h = 0 on ∂Ω}.
The projection of a function on K is the generalization of the average operation (2.10), while the projection on K ⊥ corresponds to computing its fluctuation. The space W 0 is used to characterize K ⊥ . The projections on K and K ⊥ are well-defined thanks to the: Theorem 2.1. We have the following properties

1. K is closed in L 2 (Ω). 2. W 0 equiped with the norm h W 0 = (b • ∇)h L 2 (Ω) is a Hilbert space and (b • ∇)W 0 is a closed space of L 2 (Ω). 3. K ⊥ = (b • ∇)W 0 . Proof. 1) Let φ n ∈ V such that φ n → φ in L 2 (Ω).
Then, φ n → φ in the distributional sense and the operation φ → ∇ • (bφ) is continuous for the topology of distributions. Therfore, ∇ • (bφ) = 0, which shows that φ ∈ V.

2) W 0 is a Hilbert space for the norm

h = h L 2 (Ω) + (b • ∇)h L 2 (Ω)
. According to the Poincaré inequality, the norms and W 0 are equivalent. The closedness of W 0 for the L 2 topology follows from 3).

3) The inclusion (b • ∇)W 0 ⊆ K ⊥ is obvious. We sketch the proof of the converse inclusion and leave the details to the reader. We make the hypothesis that all b-field lines are either tangent to a non-zero measure set of ∂Ω or intersect ∂Ω at two points x -and

x + such that ±(b • ν)(x ± ) > 0.
The points x -and x + are called the conjugate points of the b-field line and are respectively the incoming and outgoing points of this field line to the domain. These assumptions can certainly be weekened at the expense of technical difficulties which are outside the scope of this paper. Let ψ ∈ K ⊥ . By taking the primitive of ψ along the b-field lines, there exists φ ∈ W such that ψ = (b • ∇)φ. We can additionally impose that φ = 0 on ∂Ω -where

∂Ω ± = {x ∈ ∂Ω | ± (b • ν)(x) > 0}. Let θ ∈ K. We have 0 = Ω ψ θ dx = Ω (b • ∇)φ θ dx = ∂Ω (b • ν)φ θ dS(x), (2.20)
where dS(x) is the superficial measure on ∂Ω. Since, θ ∈ K its values at conjugate points are related by a linear relation. In particular, they can be taken simultaneously non-zero. Then, since the values of φ on ∂Ω -vanish, (2.20) implies that the values of φ on ∂Ω + vanish as well. Consequently, (b • ν)φ = 0 on ∂Ω, which shows that φ ∈ W 0 . This proves the result.

Therefore, we can decompose φ ε uniquely as

φ ε = p ε + q ε , p ε ∈ K, q ε ∈ K ⊥ , (2.21)
and state problem (1.1), (1.2) as

-(b • ∇) (∇ • (b q ε )) + ε(p ε + q ε ) = f ε , in Ω, (2.22) (b • ν) ∇ • (b q ε ) = 0, in ∂Ω, (2.23) p ε ∈ K and q ε ∈ K ⊥ .
(2.24)

Next, we introduce the variational approach. We multiply (2.22) by a test function ψ ∈ V, and integrate it on Ω. Using a Green formula together with the boundary condition (2.23), we find that

Ω ∇ • (b q ε ) ∇ • (b ψ) dx + ε Ω (p ε + q ε )ψdx = Ω f ε ψdx. (2.25)
The aim now is to decompose problem (2.25) into a problem for p ε and a problem for q ε . Hence in the following two subsections the test function ψ is chosen successively in K and in K ⊥ .

Equation for p

ε ∈ K Chosing ψ = r ∈ K in (2.25), we obtain the problem Find p ε ∈ K such that Ω (εp ε -f ε ) r dx = 0, ∀r ∈ K. (2.26)
This problem admits a solution in K which is uniformly bounded in L 2 (Ω) as ε → 0 under the compatibility condition lim ε→0 1 ε Ω f ε r dx exists and is finite, ∀r ∈ K.

(2.27)

Assuming that f ε has the following decomposition

f ε = f (0) + εf (1) + o(ε) . in L 2 (Ω), this condition implies that f (0) ∈ K ⊥ . Next, since εp ε -f ε ∈ K ⊥ , according to Theorem 2.1 there exists g ε ∈ W 0 such that εp ε -f ε = (b • ∇)g ε . (2.28)
Taking the product with b and the divergence of the result, we obtain the following Proposition 2.1. p ε is given by

p ε = 1 ε (f ε + b • ∇g ε ) in Ω. (2.29)
where g ε satisfies the problem:

-∇ • ((b ⊗ b)∇g ε ) = ∇ • (f ε b) in Ω, (2.30) (b • ν)g ε = 0 on ∂Ω, (2.31)
or, in variational form

Find g ε ∈ W 0 such that Ω (b • ∇g ε )(b • ∇θ) = Ω f ε b • ∇θ dx, ∀θ ∈ W 0 . (2.32) 2.2.3 Equation for q ε ∈ K ⊥ Taking ψ = s ∈ K ⊥ in (2.25) gives: Ω ∇ • (bq ε )∇ • (bs) dx + ε Ω q ε s dx = Ω f ε s dx. (2.33)
But since q ε and s ∈ K ⊥ , theorem 2.1 implies that there exists h ε and θ ∈ W 0 such that q ε = b • ∇h ε and s = b • ∇θ. Therefore, we get the following Proposition 2.2. q ε is given by:

q ε = b • ∇h ε , (2.

34)

where h ε satisfies the following fourth-order problem: Remark 1. In [START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF], the characterization of K ⊥ as (b • ∇)W 0 is not used. Instead, the constraint that q ∈ K ⊥ is taken into account through a mixed formulation. The number of unknowns and the size of the problem are therefore larger in [START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF] than in the present work. In practice, the resolution of the fourth order problem (2.35), (2.36), (2.37) can be reduced by solving two second-order problem, as shown below. Therefore, the introduction of a fourth order problem does not bring specific difficulties.

-∇ • [(b ⊗ b)∇(∇ • (b ⊗ b)∇h ε )] + ε∇ • ((b ⊗ b)∇h ε ) = ∇ • (bf ε ), in Ω, (2.35) (b • ν)∇ • ((b ⊗ b)∇h ε ) = 0, on ∂Ω, (2.36) (b • ν)h ε = 0, on ∂Ω, (2.37) or, in variational form Find h ε ∈ W 0 such that Ω ∇ • ((b ⊗ b)∇h ε ) ∇ • ((b ⊗ b)∇θ) dx + ε Ω (b • ∇h ε ) (b • ∇θ) dx = = Ω f ε (b • ∇θ) dx, (2.

Extension to non-homogeneous Neumann boundary conditions

The application targeted in this paper, and detailed in section 4, requires the handling of non-homogeneous Neumann boundary conditions. In this subsection φ ε is solution to the following inhomogeneous Neuman problem:

εφ ε -(b • ∇)(∇ • (b φ ε )) = b • ∇κ + f ε 2 , on Ω, (2.39) (b • ν)∇ • (b φ ε ) = -(b • ν)κ, on ∂Ω. (2.40)
where κ is a given function in W. We denote by f 1 = b • ∇κ and by

f ε = f 1 + f ε 2 .
Using the same decomposition (2.21) as before, we find that p ε satisfies (2.29) and g ε is the solution of (2.30), (2.31) or (2.32) with f ε replaced by f ε 2 (and satisfying (2.27)). Similarly, q ε satisfies (2.34) and h ε is the solution of (2.35), (2.36), (2.37), or of (2.38) with '0' at the right-hand side of (2.36) replaced by (b • ν)κ, the other terms being unchanged. The details are left to the reader.

Space discretization

The problem is discretized using a finite volume method. The domain is decomposed into a familly R of rectangles M i-1/2,j-1/2 =]x i-1 , x i [×]y j-1 , y j [ with x i = i∆x and y j = j∆y. We look for a piecewise constant approximation p ε R of p ε on each M i-1/2,j-1/2 and denote by p i-1/2,j-1/2 its constant value on this rectangle. The function g ε is approximated by a constant function on a dual mesh D, consisting of rectangles D i,j =]x i-1/2 , x i+1/2 [×]y j-1/2 , y j+1/2 [ where x i-1/2 = (i -1/2)∆x, y i-1/2 = (i -1/2)∆y. Then g ε is approximated by a piecewise constant function g ε D with its constant values denoted by g ε i,j . We approximate (2.29) by

p i-1 2 ,j-1 2 = 1 ε f ε (x i-1 2 , y j-1 2 ) + b(x i-1 2 , y j-1 2 ) • (∇g ε ) i-1 2 ,j- 1 2
.

We now define approximations (b • ∇) app and ∇ • ( • b) app of operators Ψ → (b • ∇Ψ) and Φ → ∇ • (b Φ) such that they are discrete dual operators to each other. For this purpose, we define L R and L D the space of piecewise constant functions on meshes of types R and D respectively.

Definition 3.1. The operator (b • ∇) app : L D → L R is defined by ((b • ∇Ψ) app ) i-1 2 ,j-1 2 = = b(x i-1 2 , y j-1 2 ) • Ψ i,j -Ψ i-1,j 2∆x + Ψ i,j-1 -Ψ i-1,j-1 2∆x , Ψ i,j -Ψ i,j-1 2∆y + Ψ i-1,j -Ψ i-1,j-1 2∆y . (3.1)
The operator ∇ • ( • b) app : L R → L D is defined by Definition 3.2. We define:

(∇ • (bΦ) app ) i,j = = 1 2 ∆x b x (x i+ 1 2 , y j-1 2 ) - 1 2 ∆y b y (x i+ 1 2 , y j-1 2 ) Φ i+ 1 2 ,j-1 2 + 1 2 ∆x b x (x i+ 1 2 , y j+ 1 2 ) + 1 2 ∆y b y (x i+ 1 2 , y j+ 1 2 ) Φ i+ 1 2 ,j+ 1 2 - 1 2 ∆x b x (x i-1 2 , y j-1 2 ) + 1 2 ∆y b y (x i-1 2 , y j-1 2 ) Φ i-1 2 ,j-1 2 - 1 2 ∆x b x (x i-1 2 , y j+ 1 2 ) - 1 2 ∆y b y (x i-1 2 , y j+ 1 2 ) Φ i-1 2 ,j+ 1 2 . ( 3
(∇ • (b ⊗ b • ∇)) app = (∇ • (• b)) app • (b • ∇) app , (3.3)
where • is the composition operation.

Finally, the approximation of problem (2.30), (2.31) is by solving the discrete problem for the piecewise constant function g on D:

(∇ • (b(⊗)b • ∇)) app g = (∇ • (bf )) app , (3.4)
together with Dirichlet boundary conditions on g, where f is a piecewise constant function on R. Now, problem (2.35), (2.36), (2.37) for q ε can be decomposed in two decoupled second-order elliptic problems of the type (2.30), (2.31) and can be solved by a similar method. Indeed by setting u = -∇ • ((b ⊗ b)∇h), we get that (2.35), (2.36), (2.37) is equivalent to the following two elliptic problems:

∇ • ((b ⊗ b)∇u) -εu = ∇ • (bf ) in Ω (3.5) (b • ν)u = 0 on ∂Ω (3.6)
and

-∇ • ((b ⊗ b)∇h) = u in Ω (3.7) (b • ν)h = 0 on ∂Ω. (3.8)
To summarize, the resolution of problem (1.1), (1.2) reduces to three independent resolutions of problems similar to (3.4).

4 Application to the Euler-Lorentz system in the drift limit

Introduction

In this section the drift-fluid limit of the isothermal Euler-Lorentz is investigated. This regime is representative of strongly magnetized plasma, for which the pressure term equilibrates the Lorentz force. It is obtained by letting a dimensionless parameter ε, representing the non-dimensional gyro-period as well as the square Mach number, go to zero. This limit is singular because the momentum equation in the direction of the magnetic field degenerates. Since the field may not be uniformly large, we wish to derive an Asymptotic-Preserving (AP) scheme which guarantees accurate discretizations of both the limit regime for strongly magnetized plasma (ε ≪ 1) and the standard Euler-Lorentz system when the field strength is mild (ε ∼ 1). With this aim, the Euler-Lorentz system is discretized in time by a semi-implicit scheme. This scheme has already been studied in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF] for a uniform and constant magnetic field aligned with one coordinate and for physically less meaningful Dirichlet boundary conditions. The present methodology allows us to investigate the case of non-uniform magnetic fields and Neumann boundary conditions. Indeed, the anisotropic elliptic equation (1.1), (1.2) appears as the central building block of the scheme, which allows for the computation of the field-aligned momentum component. In this presentation, we will mainly focus on this aspect, the other ones being unchanged compared to [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF].

The Euler-Lorentz model and the AP scheme

The scaled isothermal Euler-Lorentz model takes the form:

∂ t n ε + ∇ • (n ε u ε ) = 0 , (4.1)
ε ∂ t (n ε u ε ) + ∇ • (n ε u ε ⊗ u ε ) + T ∇n ε = n ε (E + u ε × B) , (4.2)
where n ε , u ε and T are the density, the velocity and the temperature of the ions, respectively. Here, the electric field E and the magnetic field B are assumed to be given functions. The parameter ε is related to the gyroperiod of the particles about the magnetic field lines, and simultaneously to the squared Mach number. We refer to [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF] for more details on the model, the scaling and the drift-fluid limit ε → 0. Now we introduce the time dicretization of the model. Let B m be the magnetic field at time t m , |B| m its magnitude and b m = B m /|B| m its direction. For a given vector field v, denote by (v) m and (v) m ⊥ its parallel and perpendicular components with respect to b

m ie v = (v) m b m + (v) m ⊥ , (v) m = v • b m , (v) m ⊥ = b m × (v × b m ).
Similarly, we denote by ∇ m and ∇ m • the parallel gradient and divergence operators respective to this field. The time semi-discrete scheme proposed in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF] is as follows:

Definition 4.3. The AP scheme is the scheme defined by:

n m+1 -n m ∆t + ∇ • (nu) m+1 = 0 , (4.3)
ε (nu) m+1 -(nu) m ∆t + ∇ • (nu ⊗ u) m + T ∇n # m+1 = n m E m+1 + (nu) m+1 × B m+1 , (4.4)
where ∇n # m+1 is given by,

∇n # m+1 = (∇n m ) m+1 ⊥ + ∇n m+1 m+1 b m+1 . (4.5)
By considering the scalar product of (4.4) with b m+1 , we get

ε( (nu) m+1 -(nu) m ∆t + ∇ • (n ε u ε ⊗ u ε ) m ) • b m+1 = -T ∇ m+1 n m • b m+1 + n m E m+1 • b m+1
and after easy computations [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF], we find that (nu) m+1 satisfies the following anisotropic elliptic problem:

ε ∆t (nu) m+1 -T ∆t ∇ m+1 ∇ m+1 • (nu) m+1 m+1 = T ∆t ∇ m+1 ∇ • (nu) m+1 m+1 ⊥ -T ∇ m+1 n m + ε ∆t (nu) m -ε ∇ • (nu ⊗ u) m + n m E m+1 m+1 . (4.6)
By setting (nu) m+1 = φ ε and by taking f = f 1 + f 2 with

f 1 = 1 ∆t b • ∇(∇ • (nu m+1 ⊥ )), (4.7
)

f 2 = - ε T (∆t) 2 (n u) m - ε T ∆t ∇ • (n u ⊗ u) m + n m E m+1 m+1 (4.8) - 1 ∆ t (b • ∇n m ). (4.9)
this problem can be put in the framework of (1.1). In [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF], because b was chosen parallel to one of the coordinate axes, a direct discretization of (4.6) using finite differences could be performed. Here, for an arbitrary anisotropy direction b, we use the method developed in the previous sections. We do not detail the description of the discretization of the other equations, since it follows [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF].

The right-hand side (4.9) can be decomposed as

f ε 2 = f (0) 2 +εf (1)
2 with f (0) 2 corresponding to the first two terms and f

(1)

2 , to the last two one. Moreover if we suppose that

n m E m+1 m+1 - 1 ∆ t (b • ∇n m ) ∈ K ⊥ , (4.10)
the compatibility condition (2.27) is satisfied. This property amounts to saying that the integrated force along a magnetic field line is zero. If the property is not satisfied, parallel velocities of order O(1/ε) are generated, which is physically unrealistic (because collisions will ultimately slow down the plasma ions). Therefore, assuming (4.10) is physically justified.

As in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF], we will compare the AP scheme with the classical semi-discrete scheme for the Euler-Lorentz model, given by: Definition 4.4. The 'classical' semi-discrete scheme is defined by:

n m+1 -n m ∆t + ∇ • (nu) m = 0 , (4.11)
ε (nu) m+1 -(nu) m ∆t + ∇ • (nu ⊗ u) m + T (∇n) m = n m E m+1 + (nu) m+1 × B m+1 . (4.12)
In [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF], it is shown that this scheme is not uniformly stable with respect to ε and so that it cannot be AP.

Except from the parallel momentum equation, which has just been discussed, the other equations of the model are discretized following [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF]. For the sake of brevity, we will not reproduce their presentation here.

Boundary conditions

The following boundary conditions are set up for test purposes only. We impose Dirichlet boundary conditions on the density n m+1 = n B with n B independent of time. For the perpendicular momentum, we impose the relation obtained after taking the limit when ε → 0 in (4.2),

nu m+1 ⊥ = - 1 |B| m+1 b × T ∇n m + n m E m+1 .
By considering the mass conservation equation at the domain boundary, we have

n m+1 -n m ∆t + ∇ • (b n m+1 u m+1 ) + ∇ • (nu) m+1 ⊥ = 0, on ∂Ω.
Therefore, as the density satisfies Dirichlet boundary conditions with timeindependent Dirichlet values, we get ). Then, we can apply the framework of section 2.2.4. When nu has been calculated, an approximation is employed in order to provide values of nu in a layer of fictious cells surrounding the boundary, by using homogeneous Neumann boundary conditions. The values in the fictitious cells are then useful to compute gradient terms which occur in the other equations of the Euler-Lorentz model.

(b • ν)∇ • (b nu m+1 ) = -(b • ν)∇ • (n u m+1 ⊥ ), on ∂Ω.

Numerical results for the elliptic problem 5.1 Introduction

In this section the efficiency of the numerical method introduced in sections 2 and 3 for the singular perturbation problem (1.1), (1.2) is investigated through numerical experiments. These experiments are carried out on a two dimensional uniform Cartesian mesh. Two sets of test cases are presented. In the first one, the anisotropy, or magnetic field, is oblique, which means that it is assumed uniform in space, but not necessarily aligned with any coordinate axis. In the second set, the field direction is non uniform. In both cases, the strength of the anisotropy is assumed uniform and is given by the value of ε. An analytical solution φ a is constructed for the singular perturbation problem (1.1), (1.2) and is compared with its approximation φ h computed on the mesh. For the test cases, the following L 1 , L 2 and L ∞ norms are used to estimate the errors between the numerical approximation φ h and the analytical solution φa :

e 1 = φa -φ h L 1 φa L 1 = i,j |φ a (x i , y j ) -φ h (i, j)| i,j |φ a (x i , y j )| , e 2 = φa -φ h L 2 φa L 2 = ( i,j |φ a (x i , y j ) -φ h i,j | 2 ) 1 2 ( i,j |φ a (x i , y j )| 2 ) 1 2 , e ∞ = φa -φ h L ∞ φa L ∞ = max i,j |φ a (x i , y j ) -φ h i,j | max i,j |φ a (x i , y j )| .
(5.1)

5.2 Numerical results for an oblique magnetic field

Introduction and test case settings

For these numerical experiments the simulation domain is the square Ω = [0, 1] × [0, 1]. The magnetic field is defined by B = (sin α, cos α, 0), with α the angle of the b-field with the x-axis ranging from 0 to π/2. In order to validate the numerical method an analytical solution denoted φ a for problem (1.1), (1.2) is constructed. It is written

φ a (x, y) = sin x sin(α) -y cos(α) + b • ∇H(x, y) , f ε a (x, y) = -b • ∇ ∇ • ((b ⊗ b)∇H(x, y)) + ε sin x sin(α) -y cos(α) + b • ∇H(x, y) , H(x, y) = (x -1)(y -1)xy 3 .
(5.

2)

The function φ a is the solution of problem (1.1), (1.2) with the right-hand side f ε a . φ a presents itself as decomposed into p ε (first terms) and q ε (second term). Note also that f ε a can be decomposed as

f ε a = f (0) a + εf (1)
a with f ∈ K ⊥ and the compatibility condition (2.27) is satisfied. However, for the simulations carried out below, the construction of the right-hand side f ε a is performed using the discrete operators (b • ∇) app and ∇ • ( • b) app in order to ensure that the compatibility condition (2.27) is satisfied by the discrete operators, namely f

(0) a ∈ K ⊥ app ,
where

K app = {φ / ∇ • (bφ) app = 0} , K ⊥ app = (b • ∇) app (W 0 ) .

Homogeneous Neumann boundary conditions

This simulation is run with α = π/3. On figure 1, we represent the relative errors as functions of the mesh sizes for different values of ε ranging from 10 -3 to 10 -9 . The curves of figure 1 are plotted using logarithmic decimal scales. We observe a linear decrease of the errors with vanishing mesh sizes, with a slope equal to 2, which proves that the global scheme is second order accurate. More importantly, we observe from figures 1(a) and 1(b), that the precision remains the same while ε is decreased by three orders of magnitude. However, for the more refined grids using the smallest value of ε of this simulation set (10 -9 , see figure 1(c)), a slight degradation of the convergence is observed for small mesh sizes. This slight degradation can be explained. Indeed, p ε is given by a stiff problem, since εp ε is obtained as the difference of two quantities scaling as ε 0 = O(1) (see (2.29), (2.30)). To investigate the influence of ε on the accuracy of the approximation of p ε , the L ∞ norm of the relative error made on p ε and on ∇ • (bp ε ) as functions of ε are plotted on figure 2.

Figure 2(a) shows a linear behavior of ∇ • (bp ε ) with vanishing ε (in log scale). To explain this feature, we note that the discretization of the second order operator in (2.30) provides a computation of ε (∇ • (• b)) app (p ε ) with the precision of the linear system solver used for the computation of g ε , which is limited by round-off errors. This error is amplified after multiplication by the factor 1/ε. This analysis still holds for the accuracy of p ε as a function of ε represented on figure 2(b) with slight differences. For the largest values of ε, we observe a plateau (red dashed line) explained by the discretization error of the discrete operators. The space discretization introduced here is second order accurate, i.e. is O(h 2 ) where h = max(∆x, ∆y). Since the right-hand side is well prepared this error only applies to the εf (1) part of f ε and is then -20 -18

-16 -

ε (decimal log scale).

∇ • (bp ε

) norm (decimal log scale).

(a) Infinity norm for ∇•(bp ε ) as a function of ε in decimal log. scales.

-20 -18 -16 - scale as ε, providing then an approximation of p ε independent of ε. The numerical methods introduced in [START_REF] Degond | Duality-based Asymptotic-Preserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF] have been developed under this assumption that f (0) = 0. The present paper is developed under a weaker hypothesis, required by the application to the Euler-Lorentz model in the drift-limit. This explains why a comparable accuracy cannot be reached. Therefore, strictly speaking, our scheme is AP for the computation of p ε only when f (0) = 0, or, when f (0) = 0, only if the round-off errors brought by the linear system solver are smaller than the discretization error. Still, it is AP without any restriction for the computation of q ε (i.e. even when

f (0) = 0).
The next simulation is aimed at investigating whether the accuracy depends on the angle between b and the coordinate axes. For this purpose, simulations are carried out on a mesh composed of 40 × 40 cells and for α ranging form 0 to π/2. When α = 0 the b field is aligned with the x-axis and when α = π/2, it is aligned with the y-axis. The relative errors are displayed as functions of α on figure 3. We observe that the variations of the errors are small on the whole range of angles. This confirms that the method provides accurate results, even when the mesh is far from consistent with the b-field direction. Direction of the magnetic field.

Error in function of the direction of the magnetic field.

error in L1 norm error in L2 norm error in infinite norm

Inhomogeneous Neumann boundary conditions

We remark that φ ε (x, y) = 2x 2 + y 2 is an analytical solution of system (2.39), (2.40) for f 2 (x, y) = ε(2x 2 +y 2 ) and κ = -∇•(bf ). For this analytical solution and ε = 10 -9 , we have checked that the relative error does not exceed 10 -13 .

5.3

Numerical results for a non uniform magnetic field

Introduction and test case settings

In this subsection Ω =]1, 2[×]1, 2[ and the magnetic field is given by:

B = |B| b, b = (sin(θ), -cos(θ)) , tan(θ) = y x . (5.3)
For this case, an analytical solution of (1.1), (1.2) can be found. We consider

H var defined on [1, 2] × [1, 2] by H var (x, y) = (1 -x) 3 (1 -y) 3 (2 -x) 3 (2 -y) 3 . According to Theorem 2.1, b • ∇H var ∈ K ⊥ . So φ = 1 + b • ∇H var is the solution of (1.1), (1.
2) when the right-hand f ε of (1.1) has the expression

f ε = -b • ∇ (∇ • (b ⊗ b)∇H var ) + ε (1 + b • ∇H var ) .

Homogeneous Neumann boundary conditions

On figures 4(a), 4(b) and 4(c), we have represented the relative errors as functions of the mesh size when ε goes from 10 -3 to 10 -9 . We observe that all the three norms decrease when the mesh sizes decrease, in a similar fashion as in the oblique uniform b-field.

Inhomogeneous Neumann boundary conditions

We take the test case of subsubsection 5.2.3 again, and we find a similar conclusion: with ε = 10 -9 , the relative error in L ∞ norm does not exceed 10 -11 .

6 Numerical results for the Euler-Lorentz system in the drift limit

Introduction and test case settings

This part is devoted to the validation of the AP-scheme (4.3), (4.4), (4.5) for the Euler-Lorentz system. Due to the lack of analytical solutions, the validation procedure consists in comparisons of the AP-scheme with the classical discretization (4.12). The classical discretization is subject to a CFL stability condition that imposes the time step to resolve (i.e. to be smaller than) the fastest time scales involved in the system. These time-resolved simulations require a time step which scales like √ ε (because the CFL condition involves the acoustic wave speed which scales like 1/ √ ε). The AP-scheme is designed to be stable independently of ε when ε → 0. In these situations, the time step cannot resolve the fastest time scales involved in the system, which leads to under-resolved simulations. The stability of the AP-scheme in under-resolved situations has be demonstrated in [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF]. In this case, the requested CFL condition only involves the fluid velocity, which is an O(1) quantity, and not the acoustic speed [START_REF] Degond | An Asymptotic-Preserving scheme for the Euler equations in a strong magnetic field[END_REF] and explains the possibility of using large time steps, independent of ε. We want to check this feature again when the scheme is equipped with our new elliptic solver. Two test cases are presented, one for an oblique uniform magnetic field, another one for a non uniform magnetic field with the same expressions as in section 5. In both cases, the electric field is chosen as E = (0, 0, B x + B y ), where B x and B y are the components of the magnetic field. The initial condition is defined by the following uniform data: n = 1, (nu) x = 1, (nu) y = -1 and (nu) z = 0 which defines a stationary solution of the Euler-Lorentz system. A local perturbation of order ε in then applied to this stationary state and the evolution of the system is observed for both the AP and the classical schemes.

Numerical results for an oblique uniform magnetic field

The results for the AP and the classical schemes are compared on figure 5 in a resolved case. Both schemes provide comparable results. However we observe the formation of a thin boundary layer on the domain frontiers for the AP-scheme but it is not responsible for the development of an instability. Next we consider the same test case with an under-resolved time step ∆t which is 10 times larger than the time step provided by the CFL condition of the classical scheme. These simulation results are collected on figure 6. In this case, the conventional scheme leads to unstable results contrary to the AP scheme and proves the capability of the AP-scheme to provide stable computations for time steps that resolve neither the acoustic wave-speed nor the gyration period.

Numerical results for a non uniform magnetic field

For the non uniform case, n = 1, (nu) x = 1, (nu) y = -1 and (nu) z = 0 are not stationary solutions to the Euler-Lorentz system. In particular, with the chosen initial condition, sharp boundary layers are generated. But the the AP scheme can still be compared with the classical scheme in the resolved case for a validation procedure. Then we take the same initial conditions as for the oblique magnetic field case. Figures (7(b), 7(a), 7(d), 7(c), 7(f), 7(e)) show that the two schemes provide similar results.

Next we consider the under-resolved time step 10∆t. In this situation Fig. 8(a), 8(c), 8(e) show that the classical scheme is unstable. By contrast, Fig. 8(b), 8(d), 8(f) demonstrate that the AP-scheme provides stable results. The increased numerical diffusion generated by the large time step gives rise to a widening of the boundary layer. Keeping the boundary layer accurate would require some mesh refinment in the vicinity of the boundary. This point is deferred to future work.

Moreover as the initial conditions of the present test case are not stationary solutions of the Euler-Lorentz model, it is important to check if the results obtained in the non resolved case by the AP scheme correspond to the proper limit regime. So we compare the results obtained with and without the local perturbation on the initial conditions. The difference between the results obtained with the two simulations remain of the same order as the perturbation of the initial condition. Fig. (9(a), 9(b), 9(c)) present the difference between the solutions obtained with the perturbed and non-perturbed initial condition, for n, (nu) x , (nu) y . The figures show that this difference is actually of 10 -10 for the density and 10 -6 for the momenta. The difference with the value of ε = 10 -9 , can be explained by the accumulation of the truncation error over the simulation time.

Conclusion and perspectives

A numerical method for degenerate anisotropic elliptic problems has been investigated. This method is based on a variational formulation together with a decomposition of the solution. This problem has been applied to the resolution of an Asymptotic-Preserving scheme for the isothermal Euler-Lorentz system. Numerical simulations demonstrate the ability of the scheme to handle under-resolved situations where the time-step exceeds the CFL stability condition of the classical scheme.

Forthcoming works will be devoted to the generalization of this approach for the full Euler system with a non linear pressure law. In this case non linear anisotropic elliptic problem have to be handled. Moreover we can also deal with the more physical situation of a plasma constituted by a mixture of ions and electrons. In this situation the model can be described by the two-fluid Euler-Lorentz system coupled with quasi-neutrality equation.

. 2 )

 2 Proposition 3.4. b • ∇ app and ∇ • b • app are adjoint operators to each other . Proof. Easy and left to the reader, thanks to a discrete Green formula. Next, we define (∇ • ((b ⊗ b) • ∇)) app by the composition of the two operators b • ∇ app and ∇ • b • app :

Therefore, nu m+1

  is a solution to the anisotropioc elliptic problem with inhomogeneous Neumann boundary conditions (2.39), (2.40), with κ = -(b•ν)∇• (n u m+1 ⊥

  and h = ∇ • ((b ⊗ b)∇H(x, y)). The function h verifies homogeneous Dirichlet boundary conditions on the domain boundaries, which implies, according to theorem 2.1, that f (0) a

ε

  (decimal log scale). Relative error infinity norm (decimal log scale).(b) Relative error in infinity norm for p ε as a function of ε in decimal log. scales.

Figure 2 :

 2 Figure 2: Oblique magnetic field test case for α = π/3 and ∆x = ∆y = 1/60. Approximation of the p ε part of the solution.

Figure 3 :

 3 Figure 3: Oblique magnetic field test case for ε = 10 -9 and ∆x = ∆y = 1/40. Norms of the relative error (5.1) as a function of the angle of the magnetic field with the x-axis α.

  log scale Relative error in log scale Error w.r.t the step size in space. log scale Relative error in log scale Error w.r.t the step size in space. log scale Relative error in log scale Error w.r.t the step size in space. L1 norm L2 norm L infinity (c) ε = 10 -9

Figure 4 :

 4 Figure 4: Non uniform magnetic field test case: error norms, defined by (5.1), for the solution φ ε as a function of the mesh size, in decimal logathimic scales, and for different values of ε.

  nu x (AP-scheme).

  nu x (classical scheme).

  nu y (classical scheme).

Figure 5 :

 5 Figure5: Euler-Lorentz test case for an oblique magnetic field in the resolved case at time t = 3.95 10 -6 s: density (n) and momentum (nu x , nu y ) computed by the AP-scheme (left) and the classical scheme (right) for ε = 10 -9 and ∆x = ∆y = 1/40. The angle of the magnetic field with the x-axis is α = π/3.

  nu y (classical-scheme)

Figure 6 :

 6 Figure6: Euler-Lorentz test case for an oblique magnetic field in the underresolved case at time 3.95 10 -5 s: density (n) and momentum (nu x , nu y ) computed by the AP-scheme (left) and the classical scheme (right) for ε = 10 -9 and ∆x = ∆y = 1/40. The angle of the magnetic field with the x-axis is α = π/3.

  nu y (AP-scheme)

Figure 7 :

 7 Figure 7: Euler-Lorentz test case for a non uniform magnetic field in the resolved case at time t = 3.95 10 -6 s: density (n) and momentum (nu x , nu y ) computed by the AP-scheme (right) and the classical scheme (left) for ε = 10 -9 and ∆x = ∆y = 1/40.

  nu y (AP-scheme)

Figure 8 :

 8 Figure 8: Euler-Lorentz test case for a non uniform magnetic field in the under-resolved case at time t = 3.95 10 -5 s: density (n) and momentum(nu x , nu y ) computed by the AP-scheme (right) and the classical scheme (left) for ε = 10 -9 and ∆x = ∆y = 1/40.

  Difference for nu y

Figure 9 :

 9 Figure 9: Difference between the solutions obtained with an initial perturbation of order ε = 10 -9 and the solution without any perturbation for the variable magnetic field after 1.58 s of simulation in the non resolved case.

0 f (0) (x, y) dy = 0 for all x ∈ [0, 1], which is satisfied thanks to hypothesis (2.3). However the solution is not unique. Indeed, if ψ verifies (2.4) then ψ + ζ is also a solution for all functions ζ = ζ(x) which depend on the x-coordinate only.

This work has been partially supported by the Marie Curie Actions of the European Commission in the frame of the DEASE project (MEST-CT-2005-021122), by the 'Fédération de recherche CNRS sur la fusion par confinement magnétique', by the Association Euratom-CEA in the framework of the contract 'Gyro-AP' (contract # V3629.001 avenant 1) and by the University Paul Sabatier in the frame of the contract 'MOSITER'. This work was performed while the first author held a postdoctoral position funded by the Fondation 'Sciences et Technologies pour l'Aéronautique et l'Espace', in the frame of the project 'Plasmax' (contract # RTRA-STAE/2007/PF/002). The authors wish to express their gratitude to G. Falchetto, X. Garbet and M. Ottaviani from CEA-Cadarache and G. Gallice and C. Tessieras from CEA-Cesta for fruitful discussions and encouragements.

proportional to εO(h 2 ) in b • ∇g ε , giving rise to a O(h 2 ) consistency error for p ε . The value of the plateau is thus only dependent of the mesh sizes and does not depend on the values of ε. With vanishing values of ε the round-off errors due to the linear system solver grow linearly (in log scale) until they reach the consistency error (O(h 2 )). This occurs for a value of ε which, for this test case, can be estimated as approximately ε = 10 -9 . For smaller ε, the discretization error is negligible compared to the round-off errors amplified by the factor 1/ε and the accuracy of p ε deteriorates linearly with vanishing ε.

The accuracy of the approximation of p ε can be made totally independent of ε under the assumption that f (0) = 0. In this case, both b • ∇g ε and f ε


Moncef Hidane

Jean-François Aujol

Yannick Berthoumieu

Charles-Alban Deledalle

Jean-François Aujol

Super

J.-F Aujol

SUPER-RESOLUTION FROM A LOW-AND PARTIAL HIGH-RESOLUTION IMAGE PAIR

Keywords: super-resolution, inpainting, nonlocal patch-based methods, graph-regularization, total variation, Douglas-Rachford algorithm

INTRODUCTION

We consider in this paper the problem of recovering a highresolution (HR) image from a pair consisting of a complete low-resolution (LR) image and an incomplete HR one. Before entering technical details, we first describe the context of our work as well as some related works.

Motivations and Context

The application that motivates our work concerns the estimation of petrophysical parameters of petroleum reservoirs through image processing techniques. In this context, a cylindrical sample of a rock is extracted from a reservoir and then imaged through computed tomography (CT) scanning in This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux -CPU (ANR-10-IDEX-03-02). J-F. Aujol is a member of Institut Universitaire de France. order to produce a three-dimensional digital representation. This 3D acquisition is further analyzed by an image analysis software in order to automatically segment and classify different regions according to visual cues related to some underlying petrophysical parameters. A typical example concerns the automatic classification of different textured regions, which in turn relates to different porosity levels of the underlying material.

In a practical setting, the limited resolution provided by the non invasive 3D reconstruction from CT acquisition does not allow for automatic inspection. Thus a need for higher resolution image data emerges. In the context we are interested in, this need is alleviated by performing a second acquisition, this time focusing on a specific part of the volume. This second acquisition thus provides incomplete but higher resolution slices of the extracted cylindrical object. Now, the goal is to combine the image data coming from both acquisitions in order to obtain a complete high-resolution volume.

Modeling

In this work, we put aside the 3D aspect and focus on the case of 2D images. Furthermore, we adopt a simplified framework where we assume that the forward model for the HR/LR pair is completely known. In particular, we do not discuss possible registration issues. Letting y 1 ∈ R p and y 2 ∈ R n respectively denote the lexicographical ordering of the complete LR and partial HR images, we adopt the following discrete forward model:

y 1 = SHf 0 + η 1 ∈ R p , y 2 = M (f 0 + η 2 ) ∈ R n , (1)
where 0 < p < n; f 0 ∈ R n denotes the unknown HR image; S ∈ R p×n stands for spatial downsampling by a dyadic factor r = 2 k in each direction; H ∈ R n×n accounts for spatial blurring of the image, modelled in our case by a circular convolution with a known point spread function h: Hx = h ⊛ x; M = diag(m 1 , . . . , m n ) ∈ R n×n accounts for some sort of occlusion in the image in the form of a binary mask: m i = 1, if pixel i is observed, and m i = 0 otherwise; the vectors η 1 ∈ R p and η 2 ∈ R n are white Gaussian noise.

Related Work

Considering the reconstruction of f 0 from y 1 alone is a zooming problem, often also called magnification, interpolation, upsampling or single-image super-resolution. This problem has been the subject of many investigations and it would be difficult within this note to sketch the different trends developed. We instead refer the interested reader to the recent works [START_REF] Bourquard | Anisotropic interpolation of sparse generalized image samples[END_REF][START_REF] Mallat | Super-resolution with sparse mixing estimators[END_REF][START_REF] Glasner | Superresolution from a single image[END_REF] and to references therein. The reconstruction of f 0 from y 2 alone is an inpainting problem. As for the zooming problem, inpainting has been extensively studied; the interested reader can consult [START_REF] Shen | Mathematical models for local nontexture inpaintings[END_REF] as a starting point for this large research theme. Broadly speaking, inpainting algorithms can be categorized as either geometryoriented [START_REF] Masnou | Disocclusion: a variational approach using level lines[END_REF][START_REF] Igehy | Image replacement through texture synthesis[END_REF] or texture-oriented [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF]; the latter class being also termed examplar-based (see [START_REF] Aujol | Exemplar-based inpainting from a variational point of view[END_REF][START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF] and references therein).

Considering the reconstruction of f 0 from both y 1 and y 2 is a super-resolution problem that can be seen either as an inpainting problem with additional LR consistency constraints or as a zooming problem with HR (soft) constraints. While the traditional super-resolution setting consists in the reconstruction of f 0 from multiple LR images related by subpixel shifts, we note that this setting becomes appropriate when we consider the full 3D problem sketched earlier. We refer the interested reader to [START_REF] Milanfar | Super-resolution imaging[END_REF] for a review article and to [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] for a recent monograph.

Plan of the paper

Our contributions are presented in Section 2 where we introduce a strategy for exploiting both the available HR and LR data for the super-resolution task. Therein, we instantiate a nonlocal variational problem whose solution gives an estimation of the HR image. Links to total variation regularization [START_REF] Buades | Image denoising methods. a new nonlocal principle[END_REF] are discussed and a minimization using Douglas-Rachford algorithm is presented. We perform a set of numerical experiments in Section 3, and show that it allows to recover fine detail structures. We conclude in Section 4.

Notations

In all this paper, a standard LR image is denoted y 1 ∈ R p while y 2 ∈ R n denotes a standard HR image, p < n. The matrices M, S, H are the ones defined above. Let Ω = {1, . . . , n} denote the set of HR pixels, D ⊂ Ω denote the set of pixels where the HR information is missing, so that D c denotes the set of pixels where the HR information is available (see Figure1). For an image f and a pixel i, p i (f ) denotes an image patch extracted from f and centered at i.

For a matrix A and a vector y, we write {A . = y} for the set {x, Ax = y}. For a set C, ι C denotes the indicator function: ι C (x) = 1 if x ∈ C, and ι C (x) = +∞ otherwise. We denote by Γ 0 (R n ) the set of all convex lower semi-continuous and proper functions on R n . Further notations will be introduced as needed.

NONLOCAL VARIATIONAL FRAMEWORK

The recovery of an estimate f of f 0 from the measurements y 1 and y 2 according to the forward model ( 1) is an ill-posed inverse problem. As is usual for such inverse problems in imaging, we formulate the estimation task in a variational setting leading to the minimization an energy function of the form

F (f ) = F 1 (f ) + F 2 (f ).
The function F 1 is a regularization term forcing the solutions to have pre-specified properties, while the term F 2 (f ) penalizes the discrepancy between f and y 1 , y 2 , according to the forward model [START_REF] Bourquard | Anisotropic interpolation of sparse generalized image samples[END_REF]. In the presence of Gaussian noise, the latter term is usually taken, in its penalized form, as the squared ℓ 2 distance and can be written in our case as

F 2 (f ) = λ 1 SHf -y 1 2 + λ 2 M f -y 1 2
. If the noise level in the HR image can be neglected, we will instead consider the following form:

F 2 (f ) = λ SHf -y 1 2 + ι {M . =y2} (f ).

LR/HR-Driven Graph Regularizer

To obtain an estimate f of f 0 we need to define a sensible prior F 1 . The regularizer we propose in this paper tries to enforce the coherence of patches in a way similar to the NL-Means algorithm [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF]. The first step is to interpolate the LR image y 1 ∈ R p to match the domain of the incomplete HR image y 2 ∈ R n . To this end, we use bicubic interpolation and get ỹ1 ∈ R n . The image ỹ1 is used to build similarities between patches whose centers lie in D and patches whose centers lie outside D. In practice, we form a directed graph G whose vertices are the (locations of the) pixels of y 2 and whose weighted adjacency matrix W ∈ R n×n satisfies w i,j = 0 if and only if (i, j) ∈ D × D c and the patch p j (ỹ 1 ) belongs the k-nearest neighbors of p i (ỹ 1 ) according to the Euclidean distance. In the experiments we performed, we opted for binary weights w i,j ∈ {0, 1}.

Equipped with the graph encoded by the weighted adjacency matrix W , we propose the following regularizer, in-spired by [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Franc | Edge direction preserving image zooming: a mathematical and numerical analysis[END_REF]:

F 1 (f ) = n i=1 ( n j=1 w i,j (f i -f j ) 2 ) 1/2 .
Note that, due to the discussed properties of W , the outer sum actually runs over D, while for i ∈ D, the inner sum runs only over a subset of D c consisting of the neighbors of i. This is in contrast with the works in [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Franc | Edge direction preserving image zooming: a mathematical and numerical analysis[END_REF], where the matrix W is assumed symmetric.

To summarize, we use an interpolated version ỹ1 of the LR image y 2 to infer a non-symmetric binary relation between the occluded patches and the visible HR patches and use this relation to infer the values of missing pixels by adopting F 1 as a prior and F 2 as a data fidelity measure.

Let us note that if we allow to define edges in Ω × Ω instead of restricting to D × D c , then one can recover the isotropic total variation (TV) regularizer [START_REF] Buades | Image denoising methods. a new nonlocal principle[END_REF] by connecting each pixel to its right and bottom neighbors. Works using the TV prior for image interpolation and super-resolution can be found in [START_REF] Tony F Chan | Superresolution image reconstruction using fast inpainting algorithms[END_REF][START_REF] Marquina | Image superresolution by tv-regularization and bregman iteration[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

Douglas-Rachford Algorithm in a Product-Space

From the previous section we are led to consider the following convex optimization problem:

minimize f ∈R n F 1 (f ) + λ 2 SHf -y 1 2 + ι {M . =y2} (f ). (2)
To proceed with the derivation of the algorithm, we rewrite the penalizer F 1 in the form F 1 (f ) = ∇ w (f ) 1,2 , where ∇ w : R n → R n×n is given, for i, {1, . . . , n}, by (∇ w f ) i,j = √ w i,j (f jf i ) and . 1,2 is the mixed-norm given, for p ∈ R n×n , by p 1,2 = n i=1 n j=1 p 2 i,j . As the energy in (2) is not differentiable, we use proximal splitting algorithms [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] in order to tackle the optimization. The base operation in this class of algorithms is the computation of the proximity operator of a function

F ∈ Γ 0 (R n ): prox F z = argmin x∈R n 1 2 x -z 2 + F (x).
The energy in (2) is composite (mixing the operators M, S, H) and contains two nonsmooth terms. The direct application of the forward-backward or the primal-dual algorithm of Chambolle and Pock [START_REF] Patrick | A proximal decomposition method for solving convex variational inverse problems[END_REF] would lead to iterative algorithms with difficult optimization sub-problems at each iteration. In order to take advantage of the properties of the involved operators, and to allow partial splitting, we adopt the strategy suggested in [START_REF] Hidane | Nonlinear multilayered representation of graph-signals[END_REF] by transcripting the Douglas-Rachford algorithm in a product space. More precisely, letting

X = R n×n × R n×n × R n×n and C = {(∇ w f ′ , Hf ′ , f ′ ), f ′ ∈ R n }, we rewrite (2) in the form minimize (p,g,f )∈X p 1,2 + λ 2 Sg-y 1 2 +ι {M . =y2} (f )+ι C (p, g, f ).
(3) We solve problem (3) by applying the Douglas-Rachford algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

Letting F 1 (p, g, f ) = p 1,2 + λ 2 Sg - y 1 2 + ι {M . =y2} (f ) and F 2 (p, g, f ) = ι C (p, g, f )
, we recall that in each iteration of the Douglas-Rachford algorithm we need to compute prox F1 and prox F2 . Let us start with the evaluation of prox F1 . Since F 1 is separable, the evaluation of prox F1 (p, g, f ) amounts to evaluating each part separately. The first component, prox . 1,2 , corresponds to the matrix form of soft-thresholding [24]. The matrix S being diagonal, the evaluation of prox S.-y1 2 amounts to solving a diagonal linear system of equations. Finally, the evaluation of the third proximity operator is immediate since its corresponds to projecting on the set {M . = y 2 }. The evaluation of prox F2 leads to a linear system whose matrix is A = ∇ t w ∇ w + H t H + I, where ∇ t w denotes the adjoint of ∇ w , with respect to the standard Euclidean inner product. We solve this inner system with a conjugate gradient method, using a 'warm start' strategy, meaning that we start the conjugate gradient solver at iteration k + 1 with the result it yielded at iteration k. The algorithm we use can now be summarized as follows:

1: γ > 0, β ∈]0, 2[ 2: p 1,0 ∈ R n×n , g 1,0 ∈ R n , f 1,0 ∈ R n 3: p 2,n = ST α (p 1,n ) (matrix soft-thresholding) 4: g 2,n = (γλS t S + I) -1 (g 1,n + γλS t y 1 ) 5: f 2,n = proj M . =y2 (f 1,n ) (Euclidean projection) 6: a n = A -1 (∇ t (2p 2,n -p 1,n )+ H t (2g 2,n -g 1,n )+2f 2,n - f 1,n ) (using conjugate gradient) 7: p 1,n+1 = p 1,n + β(∇a n -p 2,n ) 8: g 1,n+1 = g 1,n + β(Ha n -g 2,n ) 9: f 1,n+1 = f 1,n + β(a n -f 2,n )
As remarked at the end of the previous section, if we allow to define edges in {1, . . . , N } × {1, . . . , N } instead of D ×D c only, then one can recover the isotropic total variation regularizer. Then, by adopting circular boundary conditions, the operator ∇ t w ∇ w corresponds to the 4-stencil discretization of the Euclidean Laplacian operator. Hence, the matrix A in step 6 becomes diagonal in the discrete Fourier domain and thus the system can be solved by applying the discrete Fourier transform, modulating, and applying the inverse discrete Fourier transform. This is the version we use when we compare in the next section our regularization scheme with TV-based regularization.

EXPERIMENTS

We illustrate in this section the behaviour of our method on a set of 3 images. The first and second images are part of two standard test images; the last image is an image of a slice of a rock captured by an electronic microscope 1 .

In each case we start with the ground truth full HR image f 0 and simulate the incomplete image y 2 by masking part of f 0 . Similarly, we simulate the complete LR image y 1 by applying a Gaussian convolution kernel of bandwidth σ followed by spatial downsampling of factor r in each direction. The parameters used are given in the caption of corresponding figure. We then interpolate y 1 to the size of y 2 and get ỹ1 . Then, we construct a weighted graph as explained in Section 2, based on 5 × 5 patches extracted from ỹ1 . We choose 15 nearest neighbors in each case. Once the corresponding graph has been estimated, we minimize the energy (2) using the algorithm of Section 2. We finally compare the solution we get with the one obtained by adopting the total variation as a prior. The corresponding energy is minimized again using the algorithm of Section 2, the difference being that the inner system of linear equation can be solved by Fourier transform as explained above. We also compare those two solutions with a simple recopy of ỹ1 (see, e.g., image (d) in Figure 2). Regarding the regularization parameter λ, it has been chosen by trial and error, each time seeking to optimize the peak signal to noise ration between the solution and the ground truth.

From the Figures 2, 3 and 4 we see that our approach invariably performs better than the bicubic interpolation and the TV-based super-resolution both in term of PSNR and visual inspection. In particular, for the textured area of the slice rock image, our method allows to recreate the textured part in a faithful manner.

CONCLUSION

We have studied in this paper an inverse problem which combines the single-image super-resolution and the inpainting/completion problems. We have motivated the importance of this problem and one of its possible applications. Then we developed a method capable of taking advantage of both the HR and LR image data in order to perform super-resolution. The proposed method consists in minimizing a convex but

Fig. 1 :

 1 Fig. 1: Example of a partial HR (a) and complete LR (b) image pair. Our goal is to complete (a) using information from (b), or equivalently, to super-resolve (b) using information from (a).

Fig. 2 :

 2 Fig. 2: (a) ground-truth image 128 × 128; (b) LR image (σ = 1 and r = 2); (c) masked HR (50%percent missing); (d) bicubic interpolation of missing part (psnr=31.3);(e) TV (psnr=34.69); (f) our result (psnr=37.66 ).

Fig. 3 :Fig. 4 :

 34 Fig. 3: (a) ground-truth image 128 × 128; (b) LR image (σ = 2 and r = 2); (c) masked HR (50%percent missing); (d) bicubic interpolation of missing part (psnr=29.27);(e) TV (psnr=33.39); (f) our result (psnr=34.91).

The authors would like to thank the team "Sismage" from the Group TOTAL for providing CT data.


M Hidane

O Lézoray

A Elmoataz

GRAPH SIGNAL DECOMPOSITION FOR MULTI-SCALE DETAIL MANIPULATION

Keywords: Graph signals, multi-scale decomposition, detail manipulation

In this paper we introduce a new unified framework for multiscale detail manipulation of graph signals. The key to this unification is to model any kind of data as signals defined on appropriate weighted graphs. Graph signals are represented as the sum of successive layers, each capturing a given scale of detail. Detail layers are obtained through a series of regularization procedures based on total variation penalization over graphs. Layers are then processed separately before being recombined, thus achieving detail manipulation. The benefit of the approach is shown on images, 3D meshes and 3D colored point clouds.

INTRODUCTION

Edge-preserving image filtering has recently emerged as a key tool in image processing and computational photography. Most of the existing methods decompose an image into a piecewise smooth base layer and a detail layer. The obtained decomposition can then be used in a large variety of applications including sharpening, HDR tone-mapping, image editing [START_REF] Choudhury | Perceptually motivated automatic sharpness enhancement using hierarchy of non-local means[END_REF][START_REF] Durand | Fast bilateral filtering for the display of high-dynamic-range images[END_REF][START_REF] Fattal | Edge-avoiding wavelets and their applications[END_REF]. Since it is essential to have control over details at a variety of scales rather than on a single scale, multi-scale decompositions have received much interest [START_REF] Fattal | Multiscale shape and detail enhancement from multi-light image collections[END_REF]. Multi-scale image representations usually take the form of decompositions of the form f = s k + d k + . . . d 1 where s k is a coarse approximation of the original image f and the d i 's are images capturing details at specific spatial scales. Modern approaches for obtaining multi-scale decompositions retain the principle of the Laplacian pyramid but use nonlinear filters [START_REF] Paris | A fast approximation of the bilateral filter using a signal processing approach[END_REF][START_REF] Xu | Image smoothing via l 0 gradient minimization[END_REF][START_REF] Subr | Edgepreserving multiscale image decomposition based on local extrema[END_REF] in order to avoid edge blurring. Images are not the only data for which detail manipulation is of interest. For instance, with the recent availability of cheap 3D sensors, the acquisition of 3D colored point cloud is now commonplace. As for images, disposing of a hierarchical decomposition of 3D data can be very useful for their editing. If such models have recently been proposed for 3D point clouds [START_REF] Ohtake | A multi-scale approach to 3d scattered data interpolation with compactly supported basis functions[END_REF], there exists actually no unified way of obtaining multiscale decompositions whatever the considered data (images, 3D meshes, 3D colored points clouds, etc.). Therefore, there is a huge need towards the conception of methods that enable to process, under the same framework, all these different types of data. In the same time there is actually much interest in the development of algorithms that enable to process high-dimensional data that reside on the vertices or edges of a graph, referred to as graph signals [START_REF] Shuman | The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains[END_REF]. Indeed, graphs are a very natural and convenient way to represent any data [START_REF] Lézoray | Image Processing and Analysis with Graphs: Theory and Practice, Digital Imaging and Computer Vision[END_REF].

In this paper we propose a common framework allowing to obtain multi-scale decompositions of graphs signals (e.g., images and 3D points clouds). To this aim, we propose multiscale decompositions of any graph signal by iterative application of regularization procedures, in the spirit of [START_REF] Tadmor | A multiscale image representation using hierarchical (bv, l 2) decompositions[END_REF][START_REF] Hidane | Nonlinear multilayered representation of graph-signals[END_REF].

MULTILAYERED REPRESENTATION OF SIGNALS ON GRAPHS

Graphs signals

We begin by premising the definitions and notations we use.

A weighted graph G is a triple (V, E, w) where V is a finite nonempty set of vertices, E ⊂ V × V is the set of edges, and w : E → R is a weight function. We restrict ourselves to undirected graphs, meaning that for all x, y ∈ V , (x, y) ∈ E if and only if (y, x) ∈ E. We further exclude self-loops, that is edges of the type (x, x), and consider only positive symmetric weight functions w. In all this paper, we interpret w(x, y) as a measure of similarity between vertices x and y. This measure of similarity is usually obtained as a decreasing transformation of a prior distance measure. We will come back to this point in subsequent sections. Without loss of generality, we let V = {1, . . . , N } for N ≥ 1. The graph G can now be characterized by its weighted adjacency matrix W ∈ R N ×N given by

w i,j = w(i, j) if (i, j) ∈ E, 0 otherwise.
We denote by X = R N the set of real-valued functions defined on V (i.e., graph signals). Each f ∈ X assigns a real value f i to each vertex i. The set X is endowed with the usual inner product: u, v X = N i=1 u i v i for u, v ∈ X. We denote the associated norm by . X . We will also consider the space Y = R N ×N endowed with the following inner product: p, q Y = N i,j=1 p i,j q i,j for p, q ∈ Y . We denote the associated norm by . Y . We further consider on Y the following dual mixed norms defined for all p ∈ Y by:

p 1,2 = N i=1 p ⊺ i,. X = N i=1   N j=1 p 2 i,j   1/2 , ( 1
)
p ∞,2 = max 1≤i≤N p ⊺ i,. X = max 1≤i≤N   N j=1 p 2 i,j   1/2 . (2)
We consider the following first order weighted difference operator ∇ w : X → Y defined for all u ∈ X and i, j ∈ V by

(∇ w u) i,j = √ w i,j (u j -u i ). (3)
Being linear on finite-dimensional spaces, ∇ w admits a unique adjoint ∇ * w : Y → X with respect to the inner products ., . X and ., . Y . The divergence operator div w : Y → X is defined as div w = -∇ * w and its expression is given for all p ∈ Y by (div w p) i = N j=1 √ w i,j (p i,j -p j,i ).

Denoising graph signals

We consider the following degradation model f = u 0 + η where f, u 0 ∈ R N and η is random noise perturbation. In order to construct an estimator for the true signal u 0 , classical methods in signal and image processing rely on energy minimization methods [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]. Among these methods, total variation (TV) minimization [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] has played a prominent role. In the graph setting of the present paper, the isotropic TV functional is given for all u in X by [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing[END_REF]

J w (u) := ∇ w u 1,2 = N i=1   N j=1 w i,j (u j -u i ) 2   1/2 (4)
The functional J w can be regarded as a measure of "smoothness" of data living on general graphs. It is the sum of the local variations around all the vertices of the graph. The choice of the ℓ 1 norm can be seen as a regularizing prior promoting the sparsity of the local variations vector. If we assume the noise term η to be white and Gaussian, the TV-based estimator is obtained by solving the following convex optimization problem

minimize u∈X E w (u; f, λ) = λJ w (u) + 1 2 u -f 2 X , (5)
where the parameter λ is related to the noise level. Though convex, the energy function E w (u; f, λ) is not differentiable. Thus, gradient-based methods can not be employed to solve problem [START_REF] Paris | A fast approximation of the bilateral filter using a signal processing approach[END_REF]. By now, many algorithms based on splitting techniques are available to tackle this problem [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]. Due to its excellent performance, we choose in this paper to use the recent primal-dual algorithm proposed by Chambolle and Pock in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. We show in the sequel how this can be used with weighted graphs. Let us consider the following general optimization problem:

min x∈X F (Kx) + G(x), (6)
where F ∈ Γ 0 (Y ), G ∈ Γ 0 (X) and K : X → Y a linear operator. Recently, Chambolle and Pock have proposed the following iterative algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] to solve efficiently (6):

           x 0 = x0 = f, y 0 = 0 y n+1 = prox σF * (y n + σK xn ), x n+1 = prox τ G (x n -τ K * y n+1 ), xn+1 = x n+1 + θ(x n+1 -x n ), (7)
where F * is the conjugate of F , K * is the adjoint operator of K, and prox the proximity operator defined as:

prox f (x) = arg min y∈Y f (y) + 1 2 ||y -x|| 2 . ( 8
)
The convergence of algorithm ( 7) is guaranteed if θ = 1 and 0 < τ σL 2 < 1 where L = ||K|| = max ||x||≤1 ||Kx||. The optimization problem ( 5) can be formulated on graphs with

F = ||.|| 1,2 , K = ∇ w , K * = -div w and G = 1 2λ ||. -f || 2 X .
By replacing F , K, K * and G, in [START_REF] Subr | Edgepreserving multiscale image decomposition based on local extrema[END_REF], we can simplify the algorithm. The norms • 1,2 and • ∞,2 being dual, we obtain an algorithm that produces three sequences

(p n ) n ∈ Y N , (u n ) n ∈ X N and (u n ) n ∈ X N given by      p n+1 = F(p n + σ∇ w u n ), u n+1 = G(u n + τ div w p n+1 ), u n+1 = 2u n+1 -u n , (9)
where (∀p ∈ Y ), (F(p)) i,j = pi,j max(1, pi,. 2 ) and (∀u ∈ X), G(u) = 1 τ +λ (λu + τ f ). The interested reader can refer to [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Hidane | Nonlinear multilayered representation of graph-signals[END_REF] for details concerning the derivation of ( 9). This algorithm is parametrized by the structure of the graph (topology and weight function w), the functions u, u 0 , and several parameters (λ, τ , θ and σ).

Multilayered Representation of Signals on Graphs

The denoising task described in the previous subsection yields a decomposition of a noisy signal f into a clean component û and a residual v: f = û+v. In this context, the parameter λ in (5) is related to the noise level. If we relax the requirement of relating λ to the variance of η, we obtain more general decompositions of the form f = u + v. A prominent family of such decompositions is the structure-texture decomposition introduced in [START_REF] Meyer | Oscillating Patterns in Image Processing and Nonlinear Evolution Equations: The Fifteenth Dean Jacqueline B[END_REF]. In this context, u represents an image containing the contrasted shapes initially present in the visual scene, while v contains the texture and repetitive patterns. Such decompositions have been employed for image inpainting [START_REF] Bertalmio | Image inpainting[END_REF], contour detection [START_REF] Buades | Fast cartoon+ texture image filters[END_REF] and image segmentation [START_REF] Casaca | Spectral image segmentation using image decomposition and inner product-based metric[END_REF]. The structure-texture decompositions are dependent on a unique scale. In order to be able to manipulate graph signals at multiple scales, it is important to turn those decompositions into multi-scale ones [START_REF] Talebi | Global image editing using the spectrum of affinity matrices[END_REF]. For images, Tadmor et al. [START_REF] Tadmor | A multiscale image representation using hierarchical (bv, l 2) decompositions[END_REF] proposed an algorithm to turn such single-scale decompositions into multiscale ones. When implemented for signals defined on graphs with the energy E w , this algorithm is the following

       v -1 = f , u i = argmin u∈X E(u; v i-1 ; λ i ), i ≥ 0, v i = v i-1 -u i , i ≥ 0, (10)
where (λ i ) i is a sequence of scales. Starting with an initial parameter λ 0 , a first decomposition of f is obtained by applying ( 5) with λ = λ 0 , yielding f = u 0 + v 0 . The layer u 0 should be interpreted as a first sketch of the f , while the residual v 0 is to be understood as a detail coefficient. Modifying the scale parameter from λ 0 to λ 1 and applying ( 5) to v 0 with λ = λ 1 yields the subsequent decomposition v 0 = u 1 + v 1 . Now the term u 1 can be interpreted as a second layer extracted from f by means of v 0 . Iterating the same process n times leads to the following representation

(∀n ≥ 0) f = n i=0 u i + v n . (11)
The u i 's thus represent different layers of f captured at different scales. These layers are parametrized by three variables: the graph topology through the adjacency matrix W , the energy function E , and the sequence λ 0 , . . . , λ n involved in the successive minimizations. It is clear that in order to extract the successive layers in a coherent manner, the sequence of scales (λ i ) i≥0 should be decreasing. In terms of image decomposition, this assumption has the following simple interpretation: as the process (10) evolves, the successive minimizers extract more details from the original graph signal [START_REF] Tadmor | A multiscale image representation using hierarchical (bv, l 2) decompositions[END_REF][START_REF] Hidane | Nonlinear multilayered representation of graph-signals[END_REF].

EXAMPLES

First, to illustrate how the proposed decomposition works on different graph signals, we consider images and 3D meshes.

For each graph signal, a specific graph topology is considered. For images, a 8-adjacency grid graph with Gaussian

weights w i,j = e -f i -f j 2 2 2σ 2
is considered. For 3D meshes, the graph corresponds to the triangulation topology defined by the 3D mesh and is weighted by

w i,j = 1 ǫ+ fi-fj 2 2
. The graphs signal f i is a color vector for images and a vector of 3D coordinates for 3D meshes. The weight function depends on the input image and allows to introduce adaptability into the decomposition process. The initial scale λ 0 is chosen by hand so that the first layer u 0 corresponds to a first sketch of the original image. We adopt a dyadic progression of scales: λ i+1 = λ i /2 for the extraction of the next layers. Top row of Figure 1 shows sums of layers extracted from the image of a woman having several acne marks and freckles, that can be considered as very subtle detail layers of the skin. One can see clearly the reconstruction being performed as the successive layers are summed up. This example shows how well our method is able to separate details according to their scale. Bottom row of Figure 1 shows the sum of successive extracted layers for a 3D mesh of a head. One can see the reconstruction being performed as the successive layers are extracted and summed up. This example also shows the ability of our method to separate details according to their geometrical scale. Second, to manipulate details of graphs signals, we first decompose the original graph signal f into n layers. Let u i , i ∈ {0, . . . , n -1} denote the obtained layers. Then, we edit the graph signal by weighting each layer and adding the layers back together. We consider three levels of detail manipulation (coarse, intermediate and fine, denoted as g i with i ∈ [START_REF] Choudhury | Perceptually motivated automatic sharpness enhancement using hierarchy of non-local means[END_REF][START_REF] Fattal | Edge-avoiding wavelets and their applications[END_REF]) depending on which weighted layers are summed up (specified by parameters l i ). They are obtained with

g 1 = i=l1 i=0 (1 + iδ 1 )u i , g 2 = g 1 + i=l2 i=l1+1 (δ 2 + (i -l 1 -1)δ 1 δ 2 )u i , g 3 = g 2 + i=n-1 i=l2+1 (δ 2 2 + (i -l 2 -1)δ 1 δ 2 )u i .
We consider three different types of graphs signals: images, 3D meshes and 3D point clouds. The last type of graph signal associates a color vector to vertices and the graph is constructed from the 3D coordinates with a k nearest neighbor graph. For images, we set δ 1 = 2.5 and δ 2 = 0.25, δ 1 = 0.15 and δ 2 = 1 for 3D meshes and point clouds. Figure 2 shows the results. As it can be seen each level of detail manipulation enables to boost details on multiple scales. For 3D meshes, this enables to manipulate the mesh to sharpen it without any computation of the normals of the faces. Finally, it is important to note that, to the best of our knowledge, it is the first time that a method is proposed to sharpen 3D colored point clouds and this shows the innovation of our proposal. The potential of this method is huge since many 3D scanners now also acquire a color per vertex, but this information is never used to manipulate the details of the point cloud.

CONCLUSION

We have proposed a new framework for detail manipulation of graph signals. Our method represents an image or a 3D model as a graph signal and generates successive layers, each capturing a given level of detail. The obtained layers are then processed separately before being recombined, achieving detail sharpening or attenuation. The proposed method is very general and can be used for any type of graph signals and we have presented an innovative detail manipulation of colored 3D meshes.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Sum of successive layers from an image (top row) and a 3D mesh (bottom row) hierarchical decomposition.