Computational MHD: A Model Problem for Widely Separated Time and Space Scales

1. Computational MHD:A Model Problem for Widely Separated Time and Space Scales Dalton D. Schnack Center for Energy and Space Science Center for Magnetic Self-Organization in Laboratory and Space Plasmas Science Applications International Corp. San Diego, CA 92121 USA

2. Computational MHD • Computational MHD is challenging because of: • Widely separated spatial scales • Widely separated time scales • Extreme anisotopy • Closure uncertainties • Different types of MHD problems • Kinetic energy dominant (convection, dynamo) • Magnetic energy dominant (corona, lab plasmas) • Require different computational approaches • Will review approaches that have been successful of simulating highly magnetized plasmas

3. Strongly Nonlinear KE ~ ME Flow is important Shocks Turbulence Dynamos Convection Jets Properties of “statistical” state with all wavelengths represented Relatively simple BCs and geometry Important stiffness comes from nonlinearity Advection is important Weakly Nonlinear KE << ME Dominated by magnetic field Coronal loops Magnetic fusion Slow departures from equilibrium (instabilities) Culmination in disruptive behavior Complex BCs and geometry Important stiffness comes from linear properties Resistive layers, current sheets Slowly growing instabilities Parasitic modes Advection not so important Types of MHD Problems

4. Different Problems Require Different Approaches • Magneticallydominated plasmas • Coronal loops, laboratory fusion plasmas • Magnetic forces are strong (b ~ 1 in corona, b << 1 in fusion plasmas) • Low collisionality • Fluid equations not valid at small scales • Closures: characterize non-local kinetic effects (collective particle dynamics) Need to develop new algorithms • Flow dominated plasmas • Stellar interiors, convection, dynamo • Magnetic forces are weak (b >> 1) • Collisional • Fluid equations good at all scales • Closures: characterize sub-grid scale turbulence Can adapt hydrodynamic algorithms

5. “Sub-grid Scale” Models • Computational resources limit the number of degrees of freedom in discrete model • Important physical processes occur below the scale of the grid • Large Reynolds’ number turbulence • Kinetic effects in low collisionality plasmas • Must be captured in a sub-grid scale model • Turbulence • Averaging - characterize effect of small scales on large scales (new dependent variables) u = <u> + w, <<u>> = <u>, <w> = 0, <wiwj> ≠ 0 • Closure - expressions or equations for new variables, <ui∂iuj> => -n∂2<uj> • Maintain consistency - Stokes’ theorem, flux conservation,…. • Kinetic effects • Characterize non-local kinetic effects as local stress tensor, heat flux, ….. • Chapman-Enskog-like closures (integration along field lines) • Subcycling kinetic/particle • We still don’t know exactly what equations to solve! • Extensive validation against experiment and observation is required

6. To come…… • Building discrete models of physical systems • Physical basis for MHD • MHD boundary conditions • Spatial approximations • Temporal approximations (mostly implicit) • Examples • Extensions to MHD • Where can we go from here?

7. Finite Dimensional Models

8. Discrete Models • Physics described by PDEs on the continuum • “Infinite” degrees of freedom (eigenmodes) • Build discrete model with finite (N) degrees of freedom: PDEs ==> N algebraic equations • Need convergence: Discrete solution => PDE solution as N  • Requirements for convergence: • Consistency: discrete equations => differential equations as N  (including boundary conditions!) • Stability: small errors cannot grow without bound • Lax’s Theorem: Stability implies convergence for well posed consistent approximations

9. Guides for Building Discrete Models • No unique, consistent, stable discrete analog • Guide: Discrete system should have as many properties of physical system as possible • Conservation laws • Normal modes (eigenfunctions and eigenvalues) • Self-adjointness, symmetries, anti-symmetries • Boundary conditions • Not a purely mathematical exercise • Requires experience and intuition

10. The MHD Model

11. Physics Problem: Modeling Magnetized Plasmas • Dynamics of electrical conducting fluids in the presence of a magnetic field • Must solve material dynamics and Maxwell simultaneously • Simplest model is MHD, but….. Not just hydrodynamics with Lorentz force!

12. What Equations to Solve? • Complete description => kinetic equation • 6 dimensions + time • Useless in practice • Take successive velocity moments of the kinetic equation ==> fluid model (3 dimensions + time) • Each moment equation contains the next higher order moment • Must truncate, or close, the system • Closure assumption: highest order moment can be written in terms of lower order moments • Assume V2/c2 << 1 • Ignore displacement current • Implies quasi-neutrality

13. “Single-fluid” form (me=0, ne=ni=n) • Combine equations for different species • Equivalent to two-fluid model

14. Possible Fluid Models

15. Hydrodynamics Cannot support shearing displacements Linear: Sound wave Nonlinear (Riemann problem): 3 waves Shock Expansion Contact discontinuity 4 possible combinations What a Difference B Makes!

16. Hydrodynamics Cannot support shearing displacements Linear: Sound wave Nonlinear (Riemann problem): 3 waves Shock Expansion Contact discontinuity 4 possible combinations MHD Can support shearing displacements Linear: Sound wave Compressional Alfvén wave Shear Alfvén wave Nonlinear (Riemann problem): 7 waves 3 shocks 3 contact discontinuities Expansion 648 possible combinations! What a Difference B Makes!

17. The Challenges of Computational MHD

18. Challenges for Computational MHD • Strongly magnetized plasmas exhibit: • Extreme separation of time scales • Extreme separation of spatial scales • Extreme anisotropy • Each has implications for algorithms

19. 1. Extreme Separation of Time Scales • Lundquist number: • Explicit time step impractical: Requires implicit methods

20. 2. Extreme Separation of Spatial Scales • Important dynamics occurs in internal boundary layers • Structure is determined by plasma resistivity or other dissipation • Small dissipation cannot be ignored • Long wavelength along magnetic field • Extremely localized across magnetic field: d/L ~ S-a << 1 for S >> 1 • It is these long, thin structures that evolve nonlinearly on the slow time scale • Requires specialized gridding

21. 3. Extreme Anisotropy • Magnetic field locally defines special direction in space • Important dynamics are extended along field direction, very narrow across it • Propagation of normal modes (waves) depends strongly on local field direction • Transport (heat and momentum flux) is also highly anisotropic • Requires accurate treatment of Inaccuracies can lead to “spectral pollution” and anomalous perpendicular transport

22. MHD Boundary Conditions • Hydrodynamics: • Conservation laws in flux-divergence form: • Specify normal flux at boundary: • Depends on inflow/outflow, super/sub-sonic • Magnetic flux in MHD: • Must also specify Etan, and nothing more! If algorithm requires more than this, it is inconsistent!!

23. Building a Discrete Model: Spatial Discretization

24. Types of Spatial Discretization • 2 general types of approximation • Local (minimize error locally) • Finite differences • Taylor series expansion • Finite volumes • Local integral formulation • Global/Galerkin (minimize error globally) • Based on expansion in basis functions • Spectral methods • Basis functions defined globally • Finite element methods • Basis functions defined locally • All are used in computational MHD

25. Spatial Approximation • Replace: • With: • Cannot separate boundary conditions from approximation and/or grid • Discrete problem should not require more BCs

26. Finite Volumes • An “algorithm” for generating finite difference formulas • Physically motivated • Integrate conservation laws of small “finite volume” defined by grid • Exactly conservative • Consistent BCs • MHD also requires “finite area” • Magnetic flux conservation • Consistent BCs

27. y Centroid: r, p, rV Vertex: A, E x Face: B, V z Finite Volumes for MHD • Primary and dual grids (one possibility) Conserves mass, momentum, magnetic flux • Faces and edges on physical boundary to preserve BCs

28. Galerkin (Non-local) Methods • Finite differences and finite volumes minimize error locally • Based on Taylor series expansion • Galerkin methods minimize error globally • Based on expansion in basis functions • Solve “weak form” of problem Minimize global error by expansion in basis functions

29. Galerkin Discrete Approximation • Solution generally requires inverting the mass matrix • Different basis functions give different methods • Usually: bi = ai • ai=exp(ikx) => Fourier spectral methods • ai=localized polynomial => finite element methods

30. Finite Elements • Project onto basis of locally defined polynomials of degree p, e.g.: p = 1 • Integrate by parts: • Polynomials of degree p converge as hp+1 • Natural implementation of boundary conditions • Automatically preserves self-adjointness • Excellent for smooth functions • Works well with arbitrary grid shapes • Now widely used in MHD

31. Solenoidal Constraint • Faraday: • Depends on • Cannot be guaranteed in discrete model • Modified wave system • Projection • Diffusion • Grid properties

32. Building a Discrete Model: Temporal Discretization

33. Temporal Discretization • Explicit methods • Solve directly for values at n+1 in terms of values at time step n • Computationally efficient, but….. Time step limited by condition for numerical stability:

34. Temporal Discretization • Implicit methods • Values at n+1 in defined implicitly terms of values at time step n • Requires inversion of operator: • More work than explicit method, but… Unconditionally stable for any time step…..BUT: Must follow time scale of interest!

35. Multiple Time Scales(Parasitic Waves) • MHD operator contains widely separated time scales (eigenvalues) • Treat only “fast” part of operator implicitly to avoid time step restriction • Precise decomposition of W for complex nonlinear system is often difficult or impractical to achieve

36. Dealing with Parasitic Waves • Original idea from André Robert (1971) • In MHD, F and W are known, but an expression for S is difficult to achieve • W : full MHD operator • F: linearized MHD operator • Use operator splitting: • Expression for S not needed

37. Semi-Implicit Method • Recognize that the operator F is completely arbitrary!! • G can be chosen for accuracy and ease of inversion • G should be easier to invert than F (or W!) • G should approximate F for “modes of interest” • Some choices are better than others! • The semi-implicit method originated decades ago in climate modeling • Has proven to be very useful for resistive and extended MHD

38. How SI Works

39. How SI Works Old explicit code

40. How SI Works Old explicit code

41. How SI Works Old explicit code

42. How SI Works Old explicit code Semi-implicit term

43. How SI Works Old explicit code Semi-implicit term

44. How SI Works Old explicit code Semi-implicit term

45. How SI Works • Stabilizes by dispersion • Slows down unstable waves • k-dependent inertia Old explicit code Semi-implicit term

46. Semi-Implicit Operator for MHD • Linearized, ideal MHD wave operator • Wide spectrum of normal modes • Highly anisotropic spatial operator • Basis of many implicit formulations • Not a simple Laplacian • Requires specialized pre-conditioners • Challenge: find optimum algorithm for inverting this operator with CFL ~ 104

47. Fully Implicit Approach • Semi-implicit method is efficient when stiffness comes from linearities • Can split linear terms from full operator • Fusion and coronal plasmas • Fully implicit approach required when stiffness comes from non-linearities • Turbulence, etc. • Must solve non-linear discrete equations • Newton’s method • Evaluation of Jacobian • “Jacobian-free” methods • “Newton-Krylov” methods

48. Example: MHD Relaxation and the Laboratory “Dynamo”

49. Relaxed State Nonlinear Relaxation Dynamo Diffusion Instabilities Plasma Relaxation • System attempts to minimize energy subject to constraints: • Leads to preferred (“relaxed”) configurations (Taylor) • Occurs in a variety of situations • Reversed-field pinch (“dynamo”) • Tokamak (disruption) • Solar corona (CMEs?) • Underlying dynamics are MHD-like • Turbulence? • Subset of long wavelength modes? • Cyclic process

50. Laboratory Dynamo • Toroidal Z-pinch (RFP) • Positive toroidal field in center • Negative toroidal field at edge • Cyclic relaxation