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

1
Computational MHDA 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
Types of MHD Problems

• 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 ltlt 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

4
Different Problems Require Different Approaches

• Flow dominated plasmas
• Stellar interiors, convection, dynamo
• Magnetic forces are weak (b gtgt 1)
• Collisional
• Fluid equations good at all scales
• Closures characterize sub-grid scale turbulence
• Can adapt hydrodynamic algorithms

• Magnetically dominated plasmas
• Coronal loops, laboratory fusion plasmas
• Magnetic forces are strong (b 1 in corona, b ltlt
  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

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 ltugt
  w, ltltugtgt ltugt, ltwgt 0, ltwiwjgt ? 0
• Closure - expressions or equations for new
  variables, ltui?iujgt gt -n?2ltujgt
• 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 dont 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 gt N algebraic equations
• Need convergence Discrete solution gt PDE
  solution as N ? ?
• Requirements for convergence
• Consistency discrete equations gt differential
  equations as N ? ? (including boundary
  conditions!)
• Stability small errors cannot grow without bound
• Laxs 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 gt kinetic equation
• 6 dimensions time
• Useless in practice
• Take successive velocity moments of the kinetic
  equation gt 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 ltlt 1
• Ignore displacement current
• Implies quasi-neutrality

13
Single-fluid form (me0, nenin)

• Combine equations for different species

• Equivalent to two-fluid model

14
Possible Fluid Models
15
What a Difference B Makes!

• Hydrodynamics
• Cannot support shearing displacements
• Linear
• Sound wave
• Nonlinear (Riemann problem)
• 3 waves
• Shock
• Expansion
• Contact discontinuity
• 4 possible combinations

16
What a Difference B Makes!

• 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!

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 ltlt 1 for S gtgt 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
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
• aiexp(ikx) gt Fourier spectral methods
• ailocalized polynomial gt 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 hp1
• 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 n1 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 n1 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
Old explicit code
Semi-implicit term

• Stabilizes by dispersion
• Slows down unstable waves
• k-dependent inertia

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

• Newtons method
• Evaluation of Jacobian
• Jacobian-free methods
• Newton-Krylov methods

48

• Example
• MHD Relaxation and the Laboratory Dynamo

49
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

51
Lab Dynamo Mechanism

• Laboratory dynamo is a nonlinear driven system
• Plasma is resistive
• Poloidal flux (toroidal current) sustained by
  applied voltage
• Toroidal flux is sustained by a dynamo
• Mediating dynamics
• Nonlinear behavior of long wavelength MHD
  instabilities
• Driven by resistive evolution of mean fields
• Nonlinear evolution converts poloidal flux gt
  toroidal flux
• Like differential rotation
• Finite resistivity makes flux conversion
  irreversible
• No conversion of toroidal flux gt poloidal flux
• Poloidal flux supplied by external circuit
• Is it really a dynamo??

One-way flux conversion
52
Examples of Dynamo Results

• Taylor relaxation during sustainment

Self-reversal

• Perturb by 1 part in 106
• Dynamo is nonlinear and chaotic

• Begins far from relaxed state
• Time dependence of mode energy and field at wall
• Dynamo sustains discharge after t 0.07

53

• Example
• Tokamak Disruption

54
Toroidal Fusion Plasmas Tokamak

• Confine hot plasma to extract fusion energy
• Most promising concepts are toroidal - tokamak
• Magnetic field lines trace out flux surfaces
• Geometry is axisymmetric, dynamics are 3-D

55
Experimental Tokamak Disruption

• Time dependence at disruption onset
• Growing 3-D magnetic perturbation (n  1)
• Nonlinear evolution?
• Effect on confinement?
• Can this be predicted?

• Increase in neutral beam power
• Plasma pressure increases
• Sudden termination (disruption)

56
Simulation of Disruption

• Resistive MHD simulation of disruption
• Experimentally measured initial conditions
• Anisotropic heat conduction
• Vacuum region
• Ideal modes grows with finite resistivity
  (S  105)
• Magnetic field becomes stochastic
• Non-symmetric heat load on wall and divertor
• Time for crash 200 msec.
• Power 5 GW

57
Disruption Dynamics
58

• Example
• Dynamics of the Solar Corona

59
Solar Coronal Dynamics

• Coronal magnetic field energized by stresses
  (motions) at the photosphere
• Approximately force-free current slowly builds
• Sudden disruption
• Violent relaxation?
• Coronal mass ejection (CME)
• Structure of heliosphere
• Use observed photospheric fields
• Space weather

60
Solar Corona, Solar Wind, and Heliosphere
61
Launch and Propagation of CME

• Disruptions of corona result in CMEs propagates
  into heliosphere
• Coronal dynamics (lt 20 Rs) is an implicit problem
• Large variation of Alfvén and sound speeds
• Subsonic to transonic flows
• Wave dominated
• Inner heliosphere (gt 20 Rs) is an explicit
  problem
• Supersonic and super-Alfvénic flows
• Advection dominated
• Steep gradients
• Shocks
• Can couple implicit (Mikic Linker) and explicit
  (Odstrcil) codes beyond sonic and Alfvénic points
  to model launch and propagation of CME into
  heliosphere

62
CME Launch and Propagation Using Coupled Codes
63

• The Next Steps

64
Beyond MHD

• MHD is not a great model for hot magnetized
  plasmas
• Collisionality is low but not zero
• Current layers width comparable to ion Larmor
  radius
• Separate ion and electron dynamics important
• Need lowest order FLR corrections
• Full Ohms law (Hall diamagnetic)
• Gyro-viscosity (non-dissipative momentum
  transport)
• Extended MHD
• Implications for algorithms

65
Extended MHD Algorithms

• FLR terms admit new waves
• Hall term (electrons, Ohms law)
• Whistlers (correction to shear branch)
• Gyro-viscosity (ions, momentum equation)
• Corrections to compressional branch
• These new waves are dispersive
• Cannot be treated explicitly
• Require new SI operators (4th order)
• Under development

66
Limits on Modeling
Balance of algorithm performance and problem
requirements with available cycles

• Algorithms
• N - of meshpoints for each dimension
• a - of dimensions
• 1.5 - transport
• 3 (spatial) fluid
• 5-6 kinetic (spatial velocity)
• Q - code-algorithm requirements (Tflop /
  meshpoint / timestep)
• Dt - time step (seconds)

• Constraints
• P - peak hardware performance (Tflop/sec)
• e - hardware efficiency
• eP - delivered sustained performance
• T - problem time duration (seconds)
• C - of cases / year
• 1 case / week gt C 50

67
Implications Need for Better Closures

• Assumptions
• Performance is delivered
• Implicit algorithm
• Q ind. of Dt (!!)
• Requirements
• At least 3-D physics required
• Required problem time 1 msec - 1 sec
• Conclusions
• 3-D (i.e., fluid) calculations for times of 10
  msec within reach
• Longer times require next generation computers
  (or better algorithms)
• Higher dimensional (kinetic) long time
  calculations unrealistic
• Integrated kinetic effects must come through low
  dimensionality fluid closures

68
Outlook

• Interesting fluid systems have many degrees of
  freedom
• Kolmogorov gt Re per spatial dimension
• Realistic discrete model gt N Re3
• Small scale (kinetic or turbulent) physics
  affects large scale dynamics
• The computers will never be big enough or fast
  enough for realistic direct simulation!
• Need more work on closures partnerships with
  theory validation with experiment
• Need a better idea?
• Abandon traditional initial value/time marching
  approach?
• Equation-free multiscale methods (Kevrikidis
  Gear)
• Project microscale dynamics onto macroscales
• Perhaps these will lead to economical, realistic
  (predictive?) modeling