INRIA, Evaluation of Theme 1 Project-team CALVI Eric Sonnendr - - PowerPoint PPT Presentation

INRIA, Evaluation of Theme 1 Project-team CALVI

Eric Sonnendr¨ ucker

IRMA, Universit´ e de Strasbourg and CNRS projet CALVI INRIA Lorraine

Paris, 17-19 March 2009

Outline

1

Presentation of the project-team

2

Mathematical modeling of charged particles

3

Numerical methods for Vlasov equations Semi-Lagrangian methods on a fixed grid Adaptive semi-Lagrangian methods PIC solvers

4

Other results

5

Objectives for next 4 years

The CALVI Project

Scientific topics : Mathematical modeling and analysis, development and analysis of numerical schemes, Development of simulation codes on modern computers. Difficulties include complicated geometries, multiple scales and large size of discretized problems. Interdisciplinary project between mathematics and computer science involving strong links with physicists. Participating laboratories:

IECN, UMR 7502, Univ. Henri Poincar´ e-CNRS-INRIA, Nancy. IRMA, UMR 7501, University of Strasbourg and CNRS, Strasbourg. LSIIT, UMR 7005, University of Strasbourg and CNRS, Strasbourg.

Aims of the project

Modeling of plasmas and charged particle beams involving multiple scales through the development and justification of approached models using asymptotic analysis ; Analysis and development of efficient numerical methods adapted to the problem at hand. Efficient parallelization of complex codes for use on massively parallel high performance computers One major application: ITER

Controlled fusion

Magnetic confinement (ITER, Cadarache) Inertial confinement (LMJ, Aquitaine)

Kinetic models for plasmas and particle beams

Our reference model is based on the collisionless relativistic Vlasov-Maxwell equations ∂fs ∂t + p msγs · ∇xfs + qs(E(x, t) + p msγs × B) · ∇pfs = 0, ∂tE − c2∇ × B = − J ǫ0 , ∇ · E = ρ ǫ0 , ∂tB + ∇ × E = 0, ∇ · B = 0, where γ2

s = 1 + |p|2 m2

sc2 and the source terms are computed by

ρ =

• s

qs

• fs dp,

J =

• s

qs ms

• fs

p γs dp. In some cases Maxwell’s equations can be replaced by a reduced model like Poisson’s equation.

Reduced models: e.g. Tokamak plasma

Identification of physically relevant parameters v2

th = kT

m , ωc = qB m , ωp = q2n0 ε0m , rL = vth ωc , λD = vth ωp Write dimensionless (Vlasov and Poisson) equations ∂f ∂t + ρ ν v · ∇xf + (ηρ ν E + 1 ν v × B) · ∇vf = 0. −ηλ2∆φ = ni − ne. Relate different scales and pass to the limit.

New results involving mathematical analysis

Existence, uniqueness and stability result for a laser plasma interaction model (Bostan, Labrunie) Strong magnetic field limit in different configurations:

Guiding center (Bostan) Finite Larmor radius (Bostan) New derivation of drift-kinetic limit (Degond - Hirstoaga - ES

• Vignal)

The gyrokinetic Vlasov-Poisson model

In the frame of our work and ITER and the ANR EGYPT we are interested in simulations based on the following reduced model ∂f ∂t + dX dt · ∇xf + dV dt ∂f ∂v = 0, with B∗ dX dt = b × ∇J(φ) + 1 q (mV 2

∇ × b + µb × ∇B) + VB

B∗ dV dt = −(B + m q V∇ × b) · ( µ m∇B + q m∇J(φ)) and B∗ = B + m

q V∇ × b · b. J is the gyroaverage operator.

Numerical simulation of kinetic equations

Difficulties: Model defined in phase space (4D, 5D or 6D). Appearance of small scales. Particle methods: more efficient in high dimensions.

Good qualitative results at relatively low cost Numerical noise and slow convergence → difficult to have accurate results in some situations.

Methods using a grid of phase space:

Huge amount of computations in more than 3 dimensions. No numerical noise, but diffusion.

Spectral methods: Fourier - Hermite.

The backward semi-Lagrangian Method

f conserved along characteristics Find the origin of the characteristics ending at the grid points Interpolate old value at origin of characteristics from known grid values → High order interpolation needed Typical interpolation schemes.

Cubic spline (Cheng-Knorr) Cubic Hermite with derivative transport (Nakamura-Yabe)

Often used with directional splitting.

Problem with non conservative Vlasov solver

When non conservative splitting is used for the numerical solver, the solver is not exactly conservative. Does generally not matter when solution is smooth and well resolved by the grid. The solver is still second order and yields good results. However: Fine structures develop in non linear simulations and are at some point locally not well resolved by the phase space grid. In this case a non conservative solvers can exhibit a large numerical gain or loss of particles which is totally unphysical. Lack of robustness.

Vortex in Kelvin-Helmholtz instability

• 0.02

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 200 400 600 800 1000 ’diagt1kh2.plot’ u 1:3 ’diagt16kh2.plot’ u 1:3 ’diagt17kh2.plot’ u 1:3

• 0.35
• 0.3
• 0.25
• 0.2
• 0.15
• 0.1
• 0.05

0.05 200 400 600 800 1000 ’diagt1kh2.plot’ u 1:4 ’diagt16kh2.plot’ u 1:4 ’diagt17kh2.plot’ u 1:4

Conservative semi-Lagrangian method

Start from conservative form of Vlasov equation ∂f ∂t + ∇ · (fa) = 0.

• V f dx dv conserved along characteristics

Three steps:

High order polynomial reconstruction. Compute origin of cells Project (integrate) on transported cell.

Splitting on conservative form: always conservative. Implementation with splitting in 1D conservative equations easy as cells are then defined by their 2 endpoints. Filters can be designed to impose positivity.

Link between classical and conservative semi-Lagrangian methods

For constant coefficient advections it can be shown that C-Lag(2d) ⇐ ⇒ SL-Lag(2d+1) PSM ⇐ ⇒ SPL Consequences :

1 Classical and conservative semi-Lagrangian methods

equivalent for constant coefficients split equations.

2 The PFC method (Filbet-ES-Bertrand, JCP 2001)

corresponds for the Vlasov-Poisson (or Vlasov-Maxwell) systems to a classical semi-Lagrangian method with cubic Lagrange interpolation.

The forward semi-Lagrangian method

f conserved along characteristics Characteristics advanced with same time schemes as in PIC method. Leap-Frog Vlasov-Poisson Runge-Kutta for guiding-center or gyrokinetic Values of f deposited on grid of phase space using convolution kernel. Identical to PIC deposition scheme but in whole phase space instead of configuration space only. Similar to PIC method with reconstruction introduced by Denavit (JCP 1972).

Massive parallelism

Two supercomputers dedicated for magnetic fusion in near future (J¨ ulich 2009, Japan 2012). More than 10000 processors usable for gyrokinetic codes. New programming constraint (with respect to < 100 processors) for efficient use:

No transfer from one processor to all. No global data redistribution. Sophisticated adaptive methods probably less competitive due to overhead. Charge balance problems for particle methods. Advantage to local methods with static charge balance.

GYSELA 5D

Strong involvement in the development of the 5D gyrokinetic code GYSELA 5D developed at CEA Cadarache by V. Grandgirard. Backward semi-Lagrangian method. Local spline algorithm. Hybrid MPI + OpenMP parallelization.

Simulation of Zonal Fluxes with GYSELA 5D

32 64 128 256 512 1024 2048 4096

• Nb. of processors

20 40 60 80 100 120 advection 1D v|| advection 1D ϕ advection 2D field solver diagnostics total for one run

Relative efficiency (BULL/Itanium2)

16 64 256 1024 4096

• Nb. of processors

50 100 150 200 GYSELA runs ideal speedup

Relative speedup for GYSELA 5D code

New results on semi-Lagrangian schemes for the Vlasov equation on fixed grids

Besse - Segr´ e - ES (TTSP 2005): Semi-Lagrangian methods on unstructured meshes. Grandgirard et al (JCP 2006): Semi-Lagrangian solver for drift-kinetic model. Besse (SINUM 2008): Convergence of semi-Lagrangian method with cubic Hermite interpolation. Besse - Mehrenberger (Math of Comp 2008): Convergence SL method for different classes of high-order interpolators. Crouseilles - Latu - ES (JCP 2009) : Spline interpolation on patches for efficient parallelisation on massively parallel computers.

Adaptive semi-Lagrangian method

Semi-Lagrangian method consists of two stages : advection and interpolation Interpolation can be made adaptive : approximate f n with as few points as possible for a given numerical error using non linear approximation. Construct approximation layer by layer, starting from coarse approximation and adding pieces to improve precision where needed, using nested grids. It is possible to modify hierarchical decomposition so as to exactly conserve mass and any given number of moments even when grid points are removed.

Construction of a hierarchical approximation

Hierarchical approximation is constructed by defining an interpolation method enabling to go from coarse grid to fine grid. Two methods have been tried:

1 Interpolating wavelets based on Lagrange polynomial

• interpolation. Classical wavelet compression technique.

Addressed moment conservation issues (Gutnic, Haefele, ES)

2 Hierarchical approximation based on finite element

• interpolants. More local, cell based → simpler and

potentially more efficient parallelization. (Campos-Mehrenberger, Hoenen-Violard)

Extended to relativistic Vlasov-Maxwell (Besse et al).

New results for adaptive methods

Convergence result for linear elements (Campos Pinto, Mehrenberger) Extension to relativistic Vlasov-Maxwell. More complex computation of current densities. No splitting (Besse, Ghizzo, Latu, ES). Work on data structures and code optimization and paralelization

Optimization and parallelization of Hierarchical FE version in 4D phase-space. Load balancing using Hilbert curve (Hoenen-Violard). Development of fast data structure based on two levels of dense arrays each containing two grid levels. OpenMP parallelization (Latu).

Parametric Vlasov-Maxwell instability (1/3)

Parametric Vlasov-Maxwell instability (2/3)

Parametric Vlasov-Maxwell instability (3/3)

Major problem of electromagnetic PIC solver

Maxwell’s equation ∂tE − c2∇ × B = − J ǫ0 , ∇ · E = ρ ǫ0 , ∂tB + ∇ × E = 0, ∇ · B = 0, depend on continuity equation for well posedness

∂ρ ∂t + ∇ · J = 0

Not necessarily satisfied at the discrete level by numerical Vlasov solver. How to overcome ill-posedness ?

1 Find a more robust Maxwell solver. 2 Numerical methods that enforce exactly divh curlh = 0 and ρn+1

h

−ρn

h

∆t

+ divhJh = 0.

Work both on Maxwell solver and on charge and current deposition schemes.

A high order electromagnetic PIC solver : ANR HOUPIC

Development of a high order conforming Maxwell solver. Needs only solution time dependent equations. Others are satisfied automatically. Current density needs to be computed consistently with charge density to ensure discrete charge formulation. Solution found for arbitrary grid and order, provided conforming Maxwell solver is used. For exact charge conservation current density must satisfy J

n+ 1

2

i

= 1 ∆t

• k

v

n+ 1

2

k

· tn+1

tn

ϕ1

i (xk(t)) dt.

Trajectory xk(t) is affine. Hence ϕ1

i (xk(t)) polynomial of

degree 2p − 1 on quads and 2p on triangles if particle stays within cell. Use Gauss quadrature exact for these polynomials.

And also ...

Solution of Maxwell’s equations in singular domains (Ciarlet Jr, Labrunie) Development of two time scale solvers based on homogenization techniques (1D Euler, 2-scale PIC method and semi-Lagrangian solver for particle accelerators, gyrokinetic limit) (Mouton, Fr´ enod, ES) Develoment of Asymptotic Preserving semi-Lagrangian solver in the quasi-neutral limit (Belaouar, Crouseilles, Degond, ES) Development of steering and visualization tools for plasma simulations (Haefele, Dischler, Coulaud : ANR MASSIM). Development of a 1D simulation platform embedding most

• f 1D Vlasov codes.

Objectives for the next 4 years

A large scale initiative from INRIA on magnetic fusion was just started involving several INRIA project-teams as well as others. Calvi will play a leading role in this project, so that our

• bjectives will be closely linked to the work proposed in the

large scale initiative.

Mathematical analysis

Theoretical investigation of gyrokinetic and gyrofluid limits. Hierarchy of models, different levels of approximation. Rigorous foundations for Hamiltonian techniques used by physicists Study of the gyroaverage operator. Inclusion of electromagnetic effects.

Development of simulation tools

Improve semi-Lagrangian method, in particular for gyrokinetic equations

Numerical framework for efficient simulation of different dynamics across and along magnetic field lines: curvilinear coordinates, field aligned meshes. Electromagnetic effects. Drift-kinetic electrons. Optimization of codes for massively parallel computers.

Asymptotic Preserving methods, in particular for ITER injector. Continue study of high order PIC methods.

Integrated simulation framework and visualization

Involvement in EUFORIA project (7 FP of EU, 2008-2010)

• n the development of an infrastructure for integrated

tokamak simulation Involvement in ITM Taks Force of EFDA, in particular development of visualization tools. Code coupling and steering using the Kepler workflow.