Integrated Modeling of Wave-Plasma Interactions and Wave Interactions with MHD

1. Integrated Modeling of Wave-Plasma Interactions and Wave Interactions with MHD D. B. Batchelor Oak Ridge National Laboratory Workshop on ITER Simulation Beijing May 15-19, 2006 • Program and recent results from Center for Simulation of Wave-Plasma Interactions – CSWPI • Focused Integration Initiatives – toward a comprehensive Fusion Simulation Project • Program and computational architecture of center for Simulation of Wave Interactions with MHD – SWIM DBB

2. DBB

3. Antenna model 3D launch structure Plasma impedance (RANT3D,TOPICA) Wave solvers Propagation and heating (PICES, TORIC, AORSA) Plasma response Fokker-Planck solver (CQL3D, ORBIT-RF) Plasma dynamics Transport, stability, turbulence Antenna-edge interactions (sheath models) A complete simulation of RF/plasma effects requires coupling of several elements, previously studied more or less independently The focus of our SciDAC project is to integrate these elements to treat the full, nonlinear interactions between them. However, for this talk, we concentrate only on the wave solution and response of the distribution function (blue box). DBB

4. Objectives of CSWPI: understand heating to ignition, plasma control through localized heat, current and flow drive, t < 10-7 sec plasma wave current: an integral operator on E Plasma response is highly non-local  solve integral equation, large dense linear system Quasi-linear – time average distribution function f0 evolves slowly, described by time averaged Boltzmann equation particle sources quasi-linear RF flow along field lines radial transport collisional relaxation DBB

5. We are advancing two main wave solver codes within our project Blowup region • All Orders Spectral Algorithm (AORSA) – 1D, 2D & 3D (Jaeger) • Spectral in all 3 dimensions • Cartesian/toroidal coordinates • Includes all cyclotron harmonics • No approximation of small particle gyro radius r compared to wavelength l • Produces huge, dense, non-symmetric, indefinite, complex matrices • TORIC – 2D (Brambilla/Bonoli/Wright) • Mixed spectral (toroidal, poloidal), finite element (radial) • Flux coordinates • Up 2nd cyclotron harmonic • Expanded to 2nd order in r/l • Sparse banded matrices Slow ion cyclotron wave Electrostatic ion Bernstein wave DBB

6. Details of electromagnetic field structure are compared using two very different codes – TORIC and AORSA2D TORIC at 240Nr x 255 NmAORSA at 230Nx x 230 Ny • Both codes are using the same equilibrium from an Alcator C-Mod discharge with mixture of D-3He-H in (21%-23%-33%) of ne proportion. • The previous state of the art was qualitative comparison of integrated quantities such as power deposition profiles ICW IBW FW DBB

7. Calculation of wave fields self-consistently with the plasma distribution requires closed loop coupling of four significant physics models METS non-Maxwellian s Quasilinear operator f0(v||, v, y, qp) AORSA2D wave solver CQL3D Fokker Planck solver f0(v||, v, y, qp = 0) • High power waves can drive the distribution far from Maxwellian • Significant non-Maxwellian components can be produced by neutral injection or fusion alpha particles Non-Maxwellian distributions can: • Affect local damping rate and wavelength • Modify heating and current drive profiles • Change partition of power deposition among plasma species • Affect mode-conversion DBB

8. A conductivity module allowing general plasma distributions has been developed – all orders, all cyclotron harmonics, shared by all wave solvers where For Maxwellians, these 2D velocity space integrals can be expressed analytically in terms of Bessel functions and plasma dispersion function Z. For non-Maxwellians, they must be calculated for every mode in the wave spectrum at each point in space (1010) – R. Dumont, C. Phillips, D. Smithe Collaboration with Math/CS – accelerated so that evaluation of non-Maxwellian takes only about 4 times as long as evaluation of analytic Bessel/Z functions for Maxwellian DBB

9. Power absorption PRF for non-Maxwellian plasmas The local power absorption can be expressed as: (McVey, Sund, and Scharer, 1985; M.Brambilla,1988; D. Smithe 1989) where Wl is the local energy absorption kernel, The sums over k1 and k2 can be very costly to evaluate in 2D (4 nested loops) and 3D (six nested loops) with 2D velocity space integrals inside. Calculation of quasilinear operator is similar Physics analysis + Code restructuring in collaboration with Math/CS resulted in factor of 100 speedup DBB

10. CQL3D: Bounce-averaged Fokker-Planck Code Solves for bounce-averaged distribution at torus equatorial plane (qP = 0), f0(r, v||, v, t) r = generalized radial coordinate labeling (non-circular) flux surface l = field line connection length GE = velocity space flux due to toroidal electric field (Ohmic) GRF = = velocity space flux due to full, bounce average, RF Quasi-linear operator (all harmonics, Bessel functions, all wave modes) Gcoll = full, nonlinear, 2D, relativistic collisional operator R(f) = Radial diffusion and pinch operator with v dependent coefficients SNB = Monte Carlo neutral beam source (NFREYA) L(v) = velocity dependent prompt loss term DBB

11. Iteration of AORSA2d and CQL3D for C-mod minority H case [128x128 grid, 256 processors for 3 hrs on Cray XT3 – ORNL] DBB

12. Energetic ion tails measured during ICRH in Alcator C-Mod H+ ~0.1% B • Diagnostic used is a compact neutral particle analyzer (CNPA) • The spectra from 75keV-~175keV is dominated by fast ion CX w/ DNB Hydrogen. The spectra from ~250-400keV is dominated by Boron CX. • Using a weighted Maxwellian fit, we get Tion~70keV from our simple model. This corresponds to a RF power of ~5W/cc. Courtesy of Vincent Tang, MIT (APS, 2005) DBB

13. Bounce averaged distribution function for H minority fast wave heating on Alcator C-Mod • AORSA2D and CQL3D iterated to convergence • Tail ‘temperature” peaks at ~ 72 keV near r/a = 0.46 • Experimental fit gave ~ 70keV, but may not be comparable because of volume and pitch-angle weighting DBB

14. At 116 MHz, about 95% of the power is absorbed by the deuterium beam in ICRH – NBI on DIII-D 0th iteration PRF = 0 1st iteration PRF = 1.6 MW 7th iteration PRF = 1.6 MW 8th iteration PRF = 1.6 MW CQL3D: Neutral beam D 8Ω(D) 8Ω(D) 8Ω(D) 8Ω(D) 4Ω(H) e e e e P(D) = .55 MW = 34.5 % P(D) = 1.54MW = 96.2 % P(D) = 1.47 MW = 92.3 % P(D) = 1.50 MW = 94.6 % DIII-D shot 122080 (B. Pinsker) at f = 116 MHz (8th harmonic D; non-Maxwellian) and 2% minority H (4th harmonic H, Maxwellian) This is in disagreement with the experiment which shows little power absorbed at the 8th harmonic. One possible explanation is that radial diffusion and finite orbit effects are becoming important at the higher frequency DBB

15. Monte-Carlo Orbit Solver (ORBIT-RF) and Full Wave Solver (TORIC4) • ORBIT-RF solves Hamiltonian guiding center drift orbit equations • Interaction of resonant ion with wave is modeled as stochastic RF resonant kicks • Allows treatment of finite orbit width effects and radial transport of energetic particles. • Stay tuned for next talk by Vincent Chan DBB

16. Wave-Plasma interactions: Summary • We have begun to explore full-wave phenomena in 2D and 3D with much improved resolution and physics • Mode conversion in 2D on ITER scale devices • Fast wave heating in 3D • Effects of non-Maxwellian distribution on waves • We have achieved self consistent 2D wave solutions with distribution function solutions • Applications to DIII-D, Alcator C-Mod, ITER • To make feasible the integration of advanced RF models with other disciplines, such as transport or stability, or to carry out extensive studies within RF, careful attention to code speedup and optimation • It’s amazing how much you can optimize • Our partnership with Computer Science and Applied math has been essential Support your local applied math and computer science departments DBB

17. The ISOFS subcommittee (2002) recommended that the US undertake a major initiative – Fusion Simulation Project (FSP) “The purpose is to make a significant advance toward the ultimate objective of fusion simulation – to predict in detail the behavior of any discharge in a toroidal magnetic fusion device on all important time and space scales.” DBB

18. Fundamental challenges of fusion simulation – can we design the ultimate fusion simulation code starting right now? Basic description of plasma is 7D f(x, v, t), evolution determined by non-linear Boltzmann equation + Maxwell’s equations • High dimensionality • Extreme range of time scales – wall equilibration/electron cyclotron O(1014) • Extreme range of spatial scales – machine radius/electron gyroradius O(104) • Extreme anisotropy – mean free path in magnetic field parallel/perp O(108) • Non-linearity • Sensitivity to geometric details New physics will have to be incorporated over a long time scale Cannot predict progress in computer hardware, systems, or mathematical algorithm Particle sources convection in velocity space convection in space Collisional relaxation toward Maxwellian in velocity space Self -consistent charge and current source densities for Maxwell’s eqs obtained from velocity moments DBB

19. FSP Approach – start simple, pair-wise physics couplings “Focused Integration Initiatives” (FIIs) • Address one or a few of the fundamental issues for physics integration (wide range of time scales, range of space scales, extreme anisotropy, non-linearity, models of different dimensionality) • Address problems of immediate programmatic importance • Test different approaches to computational framework/code integration that might ultimately lead to a basis for complete simulation Focused Integration Initiative  center for Simulation of Wave Interactions with MHD (SWIM) DBB

20. SWIM brings together two mature sub-disciplines of fusion plasma physics, each with a demonstrated code base High power wave-plasma interactions – CSWPI Extended MHD – CEMM Why couple these particular two disciplines? • Macroscopic instabilities can limit plasma performance • RF waves can mitigate and control instabilities DBB

21. There are several experimentally demonstrated mechanisms by which RF waves can control sawtooth behavior L.-G. Eriksson et al., PRL 92, 235004 (2004) Sawtooth control on JET with Minority Current Drive on JET ICRF stabilization on JET • Sawteeth can limit plasma performance themselves, or can trigger other instabilities – disruptions, neoclassical tearing modes • Many physics processes interact – qualitative understanding exists but quantitative verification and prediction is lacking • ICRF minority current drive can either increase or decrease period and amplitude • Likely stabilization/destabilization mechanism – RF modification of current profile • ICRF heating can produce “monster” sawteeth – period and amplitude increased • Likely stabilization mechanism – energetic particle production by RF DBB

22. It has been demonstrated experimentally that suppression of NTM by RF leads to improvement in confinement • Empirical scaling of NTM pressure limits in ITER leave no margin in performance • “Understanding the physics of neoclassical island modes and finding means for their avoidance or for limiting their impact on plasma performance are therefore important issues for reactor tokamks and ITER” – ITER Physics Basis (1999) • Electron cyclotron current drive drives down mode amplitude • keeps mode rotating (no drop in frequency) • improves energy confinement R. Prater APS 2003 DBB

23. SWIM has two sets of physics goals distinguished by the time scale of unstable MHD motion MHD << HEATING MHD ~ HEATING Fast MHD phenomena – separation of time scales • Response of plasma to RF much slower than fast MHD motion • RF drives slow plasma evolution, sets initial conditions for fast MHD event • Example: sawtooth crash Slow MHD phenomena – no separation of time scales • RF affects dynamics of MHD events  MHD modifications affect RF drive plasma evolution • Deals with multi-scale issue of parallel kinetic closure including RF – a new, cutting edge field of research • Example: Neoclassical Tearing Mode Te0 time Te0 time We are approaching these regimes in two campaigns of architecture development and physics analysis and validation DBB

24. Simulation of plasma evolution requires complete model – Integrated Plasma Simulator (IPS) MHD << HEATING MHD ~ HEATING • Heating and current drive sources • Particle sources • Transport • Magnetic field evolution Te0 time Te0 time Integrated Plasma Simulator will allow coupling of virtually any fusion fusion code, not just RF and MHD, and should provide the framework for a full fusion simulation DBB

25. IPS design – Component based architecture, Plasma State component plays a central role, components implemented using existing fusion codes XPLASMA NUBEAM FRANTIC TORIC GENRAY AORSA CQL3D NUBEAM NCLASS TSC GCNM TEQ TSC DCON ELITE NOVA-K PEST-2 BALLOON M3D NIMROD GLF23 MMM95 DBB

26. Overall Project Strategy • Leverage previous investments to maximum possible extent • Fusion SciDAC projects, predictive TRANSP, NTTC • Computer Science – SciDAC ISICs, Fusion Grid, user interaction system from LEAD weather prediction system, data management methods from Common Instrument Management Architecture CIMA, authentication technology from myProxy • Component-oriented architecture (but not formal component system like CCA) • Design for flexibility and extensibility by abstracting interfaces at a high level – include multiple codes in each component from the beginning • Start with a strictlyfile-mediated exchange of data between application components • No changes needed to the application component code; SWIM overall interactions handled by scripting • Allows isolating of any bugs or problems • Eliminates issues of re-entrancy and overall system state consistency • Highest end computer systems handled from the start, not as an afterthought • Porting codes to Cray XT3 at ORNL DBB

27. SWIM sample execution model Typical component Simulation Control File(s) Select specific component codes (e.g. TORIC vs AORSA) Code-specific initialization data Fusion machine specific data Simulation control data: Computer environment settings Plasma initial conditions Plasma source and control, time sequences (events) Simulation Controller Computer environment initialization Initialize plasma state Initialize component Advance Sequence: [t  t+Dt] Step RF Step Particle Source Step Fokker Planck Solver Solve j||(E) Step Profiles Step Magnetics & Equilibrium Test linear stability Apply Reduced MHD if unstable Evaluate Simulation Control finalize, take another step, or adjust and restart current step RF Component RF_step: Get_plasma_state Convert plasma state to code-specific form and write code-specific initialization data Launch Implementing code: (i.e. AORSA, TORIC, or GENRAY Convert code specific output to plasma_state form Save internal state of component Write output to plasma_state DBB

28. Focused Integration Initiative – Center for Simulation of Wave Interactions with MHD (SWIM) D. B. Batchelor (PI),L. A. Berry, S. P. Hirshman, W. A. Houlberg, E. F. Jaeger, R. Sanchez – ORNL Fusion Energy D. E. Bernholdt, E. D’Azevedo, W. Elwasif, S. Klasky– ORNL Computer Science and Mathematics S. C. Jardin, G-Y Fu, D. McCune, J. Chen, L. P Ku, M. Chance, J. Breslau – PPPL R. Bramley – Indiana University,D. Keyes – Columbia University, D. P. Schissel – General Atomics, R. W. Harvey – CompX, D. Schnack –SAIC, J. Ramos, P. T. Bonoli, J.C. Wright– MIT Unfunded participants: L. Sugiyama – MIT, C. C. Hegna – University of Wisconsin, H. Strauss – New York University, P. Collela –LBNL H. St. John – General Atomics, G. Bateman, A. Kritz – Lehigh Univ. See our fun website at:www.cswim.org DBB

29. Summary • Our objective is to develop the SWIM architecture such that that the community will adopt the backbone for computational fusion research • The component-oriented architecture, with a flexible control level, will allow coupling of virtually any fusion code, not just RF and MHD, and should provide the framework for a full fusion simulation • Addressing ultra-scale computing from the start • Three physics goals • Unique Integrated Plasma Simulator for comprehensive long-time simulation of tokamak discharges • Sawtooth instability stabilization and destabilization • Neoclassical tearing mode stabilization DBB

30. Backup slides DBB

31. Mode conversion of plasma wavesA beautiful story of science – theory  simulation  experiment Plasma waves have an unpleasant habit of changing their character in the middle of a non-uniform plasma – mode conversion n|| = S Ion Bernstein Wave (IBW) conversion in 1D Fast wave Fast wave • On the right (low magnetic field) the fast wave has long wave length and the IBW has short, imaginary wavelength (evanescent) • In the center (near the ion-ion hybrid resonance) the modes interact • On the left (high magnetic field) the fast wave has long wave length, the IBW has short wavelength, which must be resolved, but is well separated from the fast wave. DBB

32. Surprise – fast waves can be converted to slow electromagnetic ion-cyclotron waves, as well as the expected Bernstein waves Quantitative explanation of PCI measurements on C-mod • This could be important – ICW much better at flow drive than fast wave or IBW. • Could this be used to control transport barriers? – using simulations to design experiments to find out • Could this work on ITER? – higher resolution simulations may show Simulation with TORIC code IBW/ICW DBB

33. 2D RF simulation gives complete, quantitative picture –understood by comparison with theory, compared in detail with experiment • Approximate, 1D, analytic theory (F.W. Perkins, 1977) • Provided valuable paradigms for mode conversion • Indicated several conversions were possible • Did not give quantitative information for real 2D situations • High-resolution simulations across the full plasma cross section • Includes arbitrary cyclotron harmonics • Very short wavelength structures – limited by computer size and speed, not formulation DBB