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 Presentation of the project-team 1 Mathematical modeling of charged particles 2 Numerical methods for Vlasov equations 3 Semi-Lagrangian methods on a fixed grid Adaptive semi-Lagrangian methods PIC solvers Other results 4 Objectives for next 4 years 5
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 ∂ f s p p ∂ t + · ∇ x f s + q s ( E ( x , t ) + × B ) · ∇ p f s = 0 , m s γ s m s γ s ∂ t E − c 2 ∇ × B = − J ∇ · E = ρ , , ǫ 0 ǫ 0 ∂ t B + ∇ × E = 0 , ∇ · B = 0 , s = 1 + | p | 2 where γ 2 s c 2 and the source terms are computed by m 2 q s p � � � � ρ = q s f s d p , J = f s d p . m s γ s s s 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 ω p = q 2 n 0 th = kT ω c = qB r L = v th λ D = v th v 2 m , m , ε 0 m , , ω c ω p Write dimensionless (Vlasov and Poisson) equations ∂ f ν E + 1 ∂ t + ρ ν v · ∇ x f + ( ηρ ν v × B ) · ∇ v f = 0 . − ηλ 2 ∆ φ = n i − n e . 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 dt · ∇ x f + dV � ∂ f ∂ t + d X ∂ f = 0 , dt ∂ v � with b × ∇ J ( φ ) + 1 B ∗ d X q ( mV 2 = � ∇ × b + µ b × ∇ B ) + V � B dt B ∗ dV � − ( B + m q V � ∇ × b ) · ( µ m ∇ B + q = m ∇ J ( φ )) dt 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.2 0.05 ’diagt1kh2.plot’ u 1:3 ’diagt1kh2.plot’ u 1:4 ’diagt16kh2.plot’ u 1:3 ’diagt16kh2.plot’ u 1:4 0.18 ’diagt17kh2.plot’ u 1:3 ’diagt17kh2.plot’ u 1:4 0 0.16 -0.05 0.14 0.12 -0.1 0.1 -0.15 0.08 -0.2 0.06 0.04 -0.25 0.02 -0.3 0 -0.02 -0.35 0 200 400 600 800 1000 0 200 400 600 800 1000
Conservative semi-Lagrangian method Start from conservative form of Vlasov equation ∂ f ∂ t + ∇ · ( f a ) = 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
Relative efficiency (BULL/Itanium2) 120 100 80 60 advection 1D v || advection 1D ϕ 40 advection 2D field solver 20 diagnostics total for one run 0 32 64 128 256 512 1024 2048 4096 Nb. of processors
Relative speedup for GYSELA 5D code 200 GYSELA runs ideal speedup 150 100 50 0 16 64 256 1024 4096 Nb. of processors
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).
Recommend
More recommend
