Terascale Computations of Multiscale Magnetohydrodynamics for Fusion Plasmas: ORNL LDRD Project

1. Terascale Computations of Multiscale Magnetohydrodynamics for Fusion Plasmas:ORNL LDRD Project D. A. Spong*, E. F. D’Azevedo†, D. B. Batchelor*, D. del-Castillo-Negrete*, M. Fahey‡, S. P. Hirshman*, R. T. Mills‡ *Fusion Energy Division, †Computer Science and Mathematics Division, ‡National Center for Computational Sciences Oak Ridge National Laboratory S. Jardin, W. Park, G.Y.Fu, J. Breslau, J. Chen, S. Ethier, S. Klasky Princeton University Plasma Physics Laboratory

2. Our project identified six tasks to work on over two years. Our focus during the past year has been on the first three. Adapt/optimize M3D for the Cray Adapt/optimize DELTA5D for the Cray Develop particle closure relations Couple M3D and particle closures Demonstrate code on realistic examples Final report, documentation • Our focus will shift to tasks 3 - 6 for next year • Access to ORNL X1E/XT3 computers required -more significant as project progresses

3. Significant progress on tasks that are critical to this LDRD project: • DELTA5D porting and optimization • Vectorized particle stepping + rapid field evaluation (3D spline, SVD compression) • M3D porting and optimization • New vectorized sparse matrix-vector multiply algorithm • Particle-based closure relations • Improved compatibility with MHD (flows/n/T/ isolated from particle model), benchmarked with other codes (DKES) • PPPL Subcontract added (May - September, 2005) • Addresses issues for fusion scientific end station applications

4. DELTA5D and M3D optimization

5. Particle simulation performance has been significantly improved by optimizing the magnetic field evaluation routines: DELTA5D • DELTA5D converted to cylindrical geometry for compatibility with M3D • 3D B-spline routine optimized by vectorization • For larger problems, spline memory requirements will limit number of particles per processor • New data compression techniques developed

6. DELTA5D High performance + small memory footprint SVD* fits of magnetic/electric field data have been developed -rank approximation N x N matrix Strategy: combine 2-D SVD* fit (R,Z) with 1-D Fourier series () *SVD: Singular value decomposition

7. New sparse matrix-vector multiply routine (vectorized) developed: 10 times faster M3D • M3D spends much time in: Sparse elliptic equation solvers • PETSC matrix storage • CSR (Compressed Sparse Row) • Unit stride good for scalar/poor for vector processors • CSRP (CSR with permutation) • Reorders/groups rows to match vector register size • “strip-mining” breaks up longer loops • CSRP algorithm encapsulated into a new PETSC matrix object Matrix-vector multiplies ILU preconditioner (forward/backward solves)

8. This substantial improvement in matrix-vector operations has exposed the next layer of opportunity for M3D performance gains: M3D • For small test problems up to 15% improvement is observed with the vectorized CSRP routine • Incomplete LU factorization uses sparse triangular solves that are not easily vectorizable • Other preconditioning techniques offer better vectorization - tradeoff between effectiveness vs. cost of iterations • Multi-color ordering has been used on Earth Simulator to vectorize sparse triangular solves • Forces diagonal blocks in the triangular matrices to be diagonal • Other novel vectorizable preconditioners • Sparse approximate inverse (SPAI) • Sparse direct solvers

9. Neoclassical Closure Relations

10. Our goal is to couple kinetic transport effects with an MHD model - important for long collisional path length plasmas such as ITER • Closure relations: enter through the momentum balance equation and Ohm’s law: • Moments hierarchy closed by= function of n, T,V, B, E • Requires solution of Boltzmann equation: f = f(x,v,t) • High dimensionality: 3 coordinate + 2 velocity + time

11. Neoclassical transport closures introduce new challenges: • Collisions introduce new timescales • lengthy evolution to steady state, especially at low collisionalities • Our particle model (DELTA5D) assigns particles to processors, assumes global data can be gathered at each MHD step • Multi-layer time-stepping • orbit time step (v||t1 << L||) - determined by LSODE • collisional time step (t2 << 1) • For low collisionality t1 << t2 • Want to avoid calculating quantities (flows, macroscopic gradients) that are already evolved by the MHD model • New f partitioning

12. Main components of DELTA5D Monte Carlo code: Set initialconditions andassign groupsof particlesto processors N x (time step followed by collisional step) Merge diagnostic data from different processors   • process, depending on • physics problem • - write out from PE 0 • - Collect diagnostic data • - Check for particles that • exit the solution domain • (i.e., last closed surface) • reseed or remove • Check for collisions that • scatter particles outside • of |v||/v| ≤ 1 or energy ≥ 0 • For beams or alphas, • check for thermalized • particles • Neutral beam heating/slowing down • Global energy confinement • Orbit dispersion/diffusion • PDF data for nonlocal transport • Bootstrap current • ICRF heating • Alpha particle/runaway elect. confinement

13. DELTA5D equations were converted from magnetic to cylindrical coordinates Uses bspllib 3D cubic B-spline fit to data from VMEC

14. A model perturbed field has been added to mock up tearing modes: B = BVMEC + (BVMEC):   

15. Coulomb collision operator for collisions of test particles (species a) with a background plasma (species b):

16. Monte Carlo (Langevin) Equivalent of the Fokker-Planck Operator[A. Boozer, G. Kuo-Petravic, Phys. Fl. 24 (1981)]

17. Local Monte-Carlo equivalent quasilinear ICRF operator (developed by J. Carlsson)

18. A f method is used that follows the distribution function partitioning used inH. Sugama, S. Nishimura, Phys. Plasmas 9, 4637 (2002).

19. New MHD viscosity-based closure relations are more consistent with the MHD model DELTA5D || benchmark with DKES || from DELTA5D: running time averages - 1 flux surface 3200 electrons * DKES: Drift Kinetic Equation Solver

20. DELTA5D In the next phase, kinetic closure relations will be further developed and coupled with the MHD model: • Closure relations • Collisionality scans, convergence studies • Examine/check 2D/3D variation of stress tensor • Accelerate slow collisional time evolution of viscosity coefficients • Test pre-converged restarts • Equation-free projective integration extrapolation methods • Extension to more general magnetic field models: • Green-Kubo molecular dynamics methods - direct viscosity calculation • DELTA5D/M3D coupling • Interface, numerical stability, data locality issues

21. For next step: need data from several saturated M3D tearing mode test cases at different helilicities (2/1, 3/2, 4/3, etc.)

22. ’ and island width evolution from my earlier calculations with the FAR/KITE reduced MHD codes

23. Neoclassical island width evolution from the FAR/KITE reduced MHD codes

24. M3D positive ’ runs (based on case from W. Park) Nonlinear M3D 2/1 tearing mode

25. Port and optimize M3D for Cray X1 Port and optimize DELTA5D for Cray X1 Strategy/tasks for Terascale MHD LDRD Adapt particle model to cylindrical geometry, M3D fields, parallelize (by processor or region?) Develop efficient representation for E, B broadcast to all processors Develop ’ > 0 tearing mode test cases optimize/test • Develop closure relations: • current or pressure tensor • neo. ion polarization drift effects • generalized Sugama • Ohm’s law • delta-f or full-f • new methods: • viscosity from dispersion vs.time of momentum (?) • track parallel friction between particles/background Couple currents, Pressure tensor,… back into M3D - need equations from flows/viscosities to JBS optimize/test numerical stability • Monte Carlo issues: • Collision operator conservation properties • Test particle/background consistency • Smoothed particles or bin and then smooth • Resolution of trapped/passing boundary layers