Spectral Elements for Anisotropic Diffusion and Incompressible MHD

1. Spectral Elements for Anisotropic Diffusion and Incompressible MHD Paul Fischer Argonne National Laboratory

2. Terascale Simulation Tools and Technologies • Goal: Enable high-fidelity calculations based on multiple coupled physical processes and multiple physical scales • Adaptive methods • Advanced meshing strategies • High-order discretization strategies • Technical Approach: • Develop interchangeable and interoperable software components for meshing and discretization • Push state-of-the-art in discretizations

3. Terascale Simulation Tools and Technologies • Goal: Enable high-fidelity calculations based on multiple coupled physical processes and multiple physical scales • Adaptive methods • Advanced meshing strategies • High-order discretization strategies • Technical Approach: • Develop interchangeable and interoperable software components for meshing and discretization • Push state-of-the-art in discretizations

4. Outline • Driving physics • SEM overview • Costs • MHD formulation • Convective transport properties of SEM • Anisotropic Diffusion • Brief summary of current results on MHD project • Conclusions

5. Driving Physics • Anisotropic Diffusion: ( S. Jardin, PPPL) • Central to sustaining high plasma temperatures • Want to avoid radial leakage • Need to understand significant instabilities • Incompressible MHD:( F. Cattaneo UC, H. Ji PPPL, … ) • Several liquid metal experiments are under development to understand momentum transport in accretion disks • Magneto rotational instability (MRI) is proposed as a mechanism to initiate turbulence capable of generating transport • The magnetic Prandtl number for liquid metals is ~ 10-5  Re=106 & Rm=10 for experiments • Numerically, we can achieve Re=104 & Rm=103 ( 2005 INCITE award ) • Free-surface MHD ( H. Ji, PPPL ) • Proposed as a plasma-facing material for fusion • Desire to understand the effect of B on the free surface

6. 2D basis function, N=10 Spectral Element Overview • Spectral elements can be viewed as a finite element subset • Performance gains realized by using: • tensor-product bases (quadrilateral or brick elements) • Lagrangian bases collocated with GLL quadrature points • Operator Costs: Standard FEM* SEM # memory accesses: O( EN 6 ) O( EN 3 ) # operations: O( EN 6 ) O( EN 4 ) N= order, E= number of elements, EN 3 = number of gridpoints • Dramatic cost reductions for large N ( > 5 )

7. Computational Advantages of Spectral Elements • Exponential convergence with N • minimal numerical dispersion / diffusion • for anisotropic diffusion, lements of ker(A||) well-resolved – sharp decoupling of isotropic and anisotropic modes • Matrix-free form • matrix-vector products cast as efficient matrix-matrix products • number of memory accesses identical to 7-pt. finite difference • no additional work or storage for anisotropic diffusion tensor

8. A*B A*B A+B A+B mxm vs m+m SEM Computational Kernel on Cached-Based Architectures • matrix-matrix products, C=AB: 2N 3 ops for 2N 2 memory references • much of the additional work of the SEM is covered by efficient use of cache – e.g., time for A*B vs A+B for N=10 is: 2.0 x on DEC Alpha, 1.5 x on IBM SP

9. t — plus appropriate boundary conditions on u and B Typically, Re >> Rm >> 1 Semi-implicit formulation yields independent Stokes problems for u and B Incompressible MHD

10. t Convection: dominates transport dominates accuracy requirements often the challenging part of the discretization treated explicitly in time Diffusion: “easy” ( ?? ) Projection: div u = 0 div B = 0 dominates work isotropic SPD operator multiple right-hand side information scalable multilevel Schwarz methods ( 1999 GB award ) SE multigrid ( Lottes & F 05 ) Incompressible MHD in a Nutshell

11. High-Order Methods for Convection-Dominated Flows Phase Error for h vs. p Refinement: ut + ux = 0 h-refinementp-refinement

12. High-Order Methods for Convection-Dominated Flows • Fraction of accurately resolved modes (per space direction) is increased only through increased order • Savings cubed inR3 • Rate of convergence is extremely rapid for high N • Important for multiscale / multiphysics problems ( Q: Why do we want 10 9 gridpoints? ) • Stability issues are now largely understood • stabilization via DG, filtering, etc. • dealiasing • Still, mustresolve structures (no free lunch…) • Computational costs are somewhat higher • Data access costs are equivalent to finite differences

13. c = (-x,y) c = (-y,x) Stabilizing convective problems: Models of straining and rotating flows: • Rotational case is skew-symmetric. • Filtering attacks the leading-order unstable mode. • Dealiasing (high-order quadrature) yields imaginary eigenvalues– vital for MHD N=19, M=19 N=19, M=20 l straining field rotational field

14. Diffusion – easy ??

15. CEMM Challenge Problems S. Jardin, PPPL • Anisotropic diffusion in a toroidal geometry • Two-dimensional tilt mode instability • Magnetic reconnection in 2D • Provides a problem suite that • captures essential physics of fusion simulation • stresses traditional numerical approaches • identifies pathways for next generation fusion codes • Excellent vehicle for initiating SciDAC interatctions.

16. Anisotropic Diffusion in Toroidal Domains • b – normalized B-field, helically wrapped on toroidal surfaces • thermal flux follows b.

17. High degree of anisotropy creates significant challenges • For this problem is more like a (difficult) hyperbolic problem than straightforward diffusion. • This is reflected in the variational statement for the steady case with

18. High degree of anisotropy creates significant challenges • Numerical challenges: • radial diffusion, • nearly singular, with large (but finite) null space • avoid grid imprinting • unsteady case constitutes a stiff relaxation problem • preconditioning nearly singular systems • CEMM challenge: • establish spatial convergence for steady state case • check unsteady energy conservation when • investigate the behavior of the tearing mode instability High degree of anisotropy creates significant challenges

19. Steady State Error – 3D, k|| = 10 8 centerpoint error total number of gridpoints It is advantageous to use few elements of high order • Fewer gridpoints are required • CPU time proportional to number of gridpoints (N odd, 3-7) (N even, 2-6)

20. Steady-State L2 – error over a range of discretizations E N High Anisotropy Demands High Accuracy • A = k|| A|| +AI • AI controls radial diffusion • A|| must be accurately represented when k|| >> 1 • Error must scale as ~ 1 / k||

21. 18th mode in circular geometry, k|| = 108 lk N=2 N=16 k k High Anisotropy Demands High Accuracy • Difficulty stems from high-frequency content in null space of A|| • High-order discretizations are able to accurately represent these functions.

22. Error vs. t T vs. r,t N=14 Error vs. r, N=12 Evolution of Gaussian Pulse for • minimal radial diffusion • no grid imprinting • careful time integration required (e.g., adaptive DIRK4) f-averaged temperature vs. time Evolution of Gaussian pulse for

23. W SEM is able to identify critical physics – • Tearing mode instability • radial perturbation: • b = b0 + e cos(mq-nf) r • field lines do not close on m-n rational surface • magnetic island results, with significant increase in radial conductivity. • island width scales as: • W ~ in accord with asymptotic theory Note: this is a subtle effect! Island width vs. k|| at onset. W W

24. max dT/dr a32 a32 Outstanding Challenges for Anisotropic Diffusion Simulation • Preconditioning • need null-space control • condition number scales as k|| • Non-aligned grids predict early island formation f-averaged temperature vs. time

25. Incompressible MHD Results

26. Axisymmetric Hydro Simulations of Taylor-Couette w/ Rings Re=6200 • Re=620 steady • Re=6200 unsteady • Axisymmetric MHD simulations are being carried out now. • Starting point for 3D simulations, which are being compared with experiments at PPPL. Computation by Obabko, Fischer, & Cattaneo inner cylinder outer cylinder Normalized Torque Vorticity

27. Computational MRI: preliminary results Computations Fischer, Obabko & Cattaneo • Nonlinear development of Magneto-Rotational Instability • Cylindrical geometry similar to Goodman-Ji experiment • Hydrodynamically stable rotation profile • Weak vertical field • Use newly developed spectral element MHD code • Try to understand differences between experiments and simulations • Simulations ReRm (moderate). Experiments Re>>Rm (Rm smallish)

28. Summary & Conclusions • Block-structured SEM provides an efficient path to high-order • accurate treatment of challenging physics • fast cache-friendly operator evaluation ( N 4 vs. N 6) • For anisotropic diffusion • effects of grid imprinting are minimized • able to capture physics of high anisotropy • MRI experiment • MRI has been observed with axial periodicity at Re=Rm=1000. • preliminary axisymmetric results indicate hydrodynamic unsteadiness at Re=6000 for two-ring boundary configuration, may be mitigated in 3D… • Free-surface MHD • Free-surface NS is working • Coupling with full MHD is underway