The Magneto-hydrodynamic Richtmyer-Meshkov Instability

1
The Magneto-hydrodynamic Richtmyer-Meshkov
Instability

• Ravi Samtaney
• Computational Plasma Physics Group
• Princeton Plasma Physics LaboratoryPrinceton
  University
• GALCIT Fluid Mechanics SeminarCaltech,
  04/18/2003

2
AMR Simulations of the Magneto-hydrodynamic
Richtmyer-Meshkov Instability

• Ravi Samtaney
• Computational Plasma Physics Group
• Princeton Plasma Physics LaboratoryPrinceton
  University
• CPPG SeminarPPPL, 04/09/2003

3
Acknowledgement

• Phillip Colella and Applied Numerical Algorithms
  Group, LBNL
• Steve Jardin, PPPL
• DOE SciDAC program

4
Outline

• MHD primer (MHD with shocks)
• Richtmyer-Meshkov Instability (RMI)
• RMI suppression in MHD
• Numerical Method
• Adaptive Mesh Refinement (AMR)
• Unsplit upwinding method
• div(B) issues
• Conclusion and future work

5
Princeton Plasma Physics Laboratory

• DOE Collaborative National Center for plasma and
  fusion science
• Primary mission is to develop the scientific
  understanding and key innovations which will lead
  to an attractive fusion energy source
• Magnetic fusion research at Princeton began in
  1951 under the code name Project Matterhorn
• PPPL has been a leader in magnetic confinement
  fusion experiments utilizing the tokamak approach
• Experimental work at PPPL includes TFTR, NSTX,
  NCSX, MHD Hall Thrusters etc.
• Many research groups including theory,
  computational plasma physics group, plasma
  science and technology etc.

100 physicists 100 engineers 250 support staff
6
Plasma Confinement

• Definition (bland) A plasma is a hot ionized
  gas whose motions both generate electromagnetic
  fields and are influenced by them. It is the
  fourth state of matter.

7
Typical Time Scales in a next step experiment
with B 10 T, R 2 m, ne 1014 cm-3, T 10 keV
ENERGY CONFINEMENT
SAWTOOTH CRASH
ELECTRON TRANSIT
TURBULENCE
CURRENT DIFFUSION
?LH-1
??A
ISLAND GROWTH
?ci-1
?ce-1
10-10
10-2
104
100
10-8
10-6
10-4
102
SEC.
Neglect displacment current, integrate over
velocity space, average over surfaces, neglect
ion electron inertia Transport Codes discharge
time-scale
Single frequency and prescribed plasma
background RF Codes wave-heating and
current-drive
Neglect displacement current, average over
gyroangle,
(some) with electrons Gyro-kinetics
Codes turbulent transport
Neglect displacement current, integrate over
velocity space, neglect electron
inertia Extended MHD Codes macro stability
8
Challenge to Theory Simulations
atomic mfp
electron-ion mfp
system size
skin depth
tearing length
Huge range of spatial and temporal scales
Overlap in scales often means strong (simplified)
ordering not possible
ion gyroradius
debye length
electron gyroradius
Spatial Scales (m)
10-2
10-6
10-4
100
102
pulse length
current diffusion
Inverse ion plasma frequency
inverse electron plasma frequency
confinement
ion gyroperiod
Ion collision
electron gyroperiod
electron collision
10-10
10-5
100
105
Temporal Scales (s)
9
MACROSCOPIC (MHD) SIMULATIONS
DIFFERENT
LEVELS OF ANALYSIS CAPABILITY
Less complex model, valid for high-collisionality,
strong fields, long times
More computationally demanding. Required to
describe many important but subtle phenomena.
Two Fluid MHD plus energetic gyro-particles
Gyro-particle ions and fluid electrons
Full orbit particle ions and fluid electrons
Single Fluid Resistive MHD
Two Fluid MHD (electrons and ions)
MHD modes destabilized by wave-particle resonance
with energetic species
External kink modes
Neoclassical tearing mode (including rotation)
Collisionless reconnection
Tilting and interchange modes with large
gyro-orbits
Kinetic stabilizationof internal MHD modes by ions
10
MHD Primer

• A continuum or fluid description of a plasma
• MHD equations are mathematical models describing
  the flow of a conducting fluid coupled with
  electromagnetism
• Governing are the equations of hydrodynamics
  coupled with Maxwells equationsFaraday ? B/? t
  -r E (implies r B0)Ohm E-v B ?
  JAmpere Jr B
• Momentum equations include J B forceEnergy
  equation includes J E (Ohmic heating)
• Several ways to combine hydrodynamics with
  Maxwell
• Weakly coupled, Strongly coupled, vector
  potential formulations etc.
• Mathematically, the result is a system of coupled
  nonlinear partial differential equations which
  usually must be solved numerically
• For ideal MHD the equations are hyperbolic PDEs
• Major problem is dealing with solenoidal property
  of the magnetic field

11
Electromagnetic Coupling

• Weakly coupled formulation
• Hydrodynamic quantities in conservative form,
  electrodynamic terms in source term
• Hydrodynamic conservation jump conditions
• One characteristic wave speed (ion-acoustic)
• Tightly coupled formulation
• Fully conservative form
• MHD conservation and jump conditions
• Three characteristic wave speeds (slow, Alfvén,
  fast)
• One degenerate eigenvalue/eigenvector

j E
12
Single-fluid resistive MHD Equations

• Equations in conservation form

Parabolic
Hyperbolic
Reynolds no.
Lundquist no.
Peclet no.
13
MHD Waves

• Disturbances propagate with finite speed
• Three sound speeds
• Fast magneto-acoustic speedcf21/2(a2A2((a2A2
  )2-4a2An2)1/2
• Slow magneto-acoustic speedcs21/2(a2A2-((a2A2
  )2-4a2An2)1/2
• Alfvén speed (also occurs in incompressible
  fluids)An2 B n/?
• cs lt a lt cfcs lt An lt cf

14
MHD Discontinuities

• MHD shock waves
• Fast shocks
• Slow shocks
• Intermediate shocks (unstable if ?0)
• Compound shocks (shock-rarefaction structures do
  not occur in ideal MHD)
• Alfvén shocks

• MHD contact discontinuities
• No tangential jump in velocity if normal
  component of B is present
• Shear flow discontinuities allowed if B.n0

Figures from PhD thesis of Hans De Sterck,
Katholieke Universiteit Leuven
15
Richtmyer-Meshkov Instability (RMI)

• RMI
• occurs at a fluid interfaceaccelerated by a
  shock-wave
• impulsive Rayleigh-Taylor
• Linear stability analysis by Richtmyer (1960)
• Experimental confirmation by Meshkov (1970)

16
Richtmyer-Meshkov Instability (RMI)

• Occurs
• in inertial confinement fusion where it is the
  main inhibiting hydrodynamic mechanism
• in astrophysical situations
• Approx. 200 peer-reviewed publications numerous
  proceedings in last 15 years.

X-ray images from Nova laser expts.Barnes et al.
(LANL Progress Report 1997)
17
Parameters and Initial/Boundary Conditions

• Parameters which completely define the posed
  Mathematical problem
• Mach number of the incident shock, M
• Density ratio (or Atwood number at the fluid
  interface), ?, At(?-1)/(?1)
• Perturbation ?
• Magnetic field strength ?-1B02/2p0

18
Simulation Results

• Parameters which completely define the posed M2,
  ?45o, ?3, ?-10 (B00) or?-10.5 (B01)
  t0.385

19
Simulation Results

• t1.86

20
Simulation Results

• t3.33

Results! Why, man, I have gotten a lot of
results. I know severalthousand things that
wont work- Thomas Edison
21
Simulation Results Early Refraction Phase

• t3.33

22
Vortex Dynamical Interpretation

• During incident shock refraction on the density
  interface vorticity is generated due to
  misalignment of density and pressure gradients
  (Hawley and Zabusky, PRL 1989)

Light
Heavy
r p
r p
Heavy
Light
r ?
r ?
23
Why does B suppress the instability? A
mathematical explanation

• Consider jump conditions across a stationary
  discontinuity? un 0? un2 p ½ B2 ½
  Bn2 0? un ut Bn Bt 0Bn 0Bt un
  Bn ut 0(e p ½ B2) un Bn (u B) 0
• For un0 (no flow through discontinuity) implies
• ut 0 if Bn!0 or
• pBt2 0 and ut and Bt can jump if Bn0
• For un!0, ut jumps across the discontinuity.
• In MHD shocks can support shear
• Jump conditions originally derived by Teller and
  do Hoffman (1952) Other sources Friedrichs and
  Kranzer (1958), Anderson (1963)

Young man, in mathematics you dont understand
things, you just get used to them J. Von
Neumann
24
Vortex sheets Early Refraction Phase

• t0.385
• The shocked interface is a vortex sheet (?-10)
• In MHD shocks can support shear

25
Vortex Dynamical Interpretation

• Baroclinic vorticity is taken away by slow shocks
  in the presence of a magnetic field. These shocks
  are stable (at least in our simulations)
• Net circulation in the domain

26
Current Sheets Early Refraction Phase

• t0.385

27
I dont believe this! What if B is really small?

• In case the magnetic field is very small, the
  slow shocks move closer to the interface. The
  interface amplitude grows due to an entrainment
  effect.

?-10.05
28
Results for ? lt 1

• The same suppression effect is observed

t1.91
t3.31
29
Adaptive Mesh Refinement

• Adaptive mesh refinement provides a numerical
  microscope
• Provides resolution where we need it
• Optimally should be based on a local truncation
  error analysis
• Mesh refined where local error exceeds a user
  defined threshold
• In practice, we refine the mesh where a threshold
  is exceeded (usually depends on gradients of flow
  quantities)
• LTE analysis is not very meaningful if
  discontinuities are present
• Usually leads to more efficient compuations
• Various AMR frameworks/packages
• GrACE
• SAMRAI
• Paramesh
• CHOMBO
• AMRITA

30
Why AMR?

• For the RMI case
• Base mesh 256x32
• Plus 3 levels of refinement (each with a factor
  of 4)
• Effective uniform mesh of 16384x 2048
• Estimated computational saving about a factor of
  25
• Each run took about 80 hours on NERSC
• Uniform mesh run would have taken 25x802000
  hours
• Area covered
• Finest mesh 2.1
• Level 2 6.5
• Level 1 15.6

It is unworthy of excellent mento lose hours,
like slaves, in thelabors of calculation
-Gottfried Wilhelm Leibnitz
31
Why AMR?

• Area covered
• Finest mesh (321) 7.6
• Level 2 (38) 18.5
• Level 1 (5) 50
• Level 0 (2) 100

32
AMR Basics

• Flag cells which need to be refined
• Using a clustering algorithm make patches
• Copy from previous overlapping fine mesh and
  interpolate from coarse mesh where new refinement
  occurs
• Surround each fine patch with a layer of ghost
  cells
• Interpolation must be consistent to the order of
  the method

33
AMR Basics

• Update coarse mesh solution
• Update fine mesh r times
• Coarse mesh solution interpolated in space and
  time to provide boundary conditions at CFI
• Synchronize at the end of each time step by flux
  refluxing to maintain strong conservation
• Average fine mesh solution replaces coarse
  solution where coarse/fine meshes overlap

Adaptive time stepping
Synchronizeby replacingcoarse meshfluxes by
fine mesh fluxes at CFI
34
Adaptive Mesh Refinement with Chombo

• Chombo is a collection of C libraries for
  implementing block-structured adaptive mesh
  refinement (AMR) finite difference calculations
  (http//www.seesar.lbl.gov/ANAG/chombo)
• (Chombo is an AMR developers toolkit)
• Mixed language model
• C for higher-level data structures
• FORTRAN for regular single grid calculations
• C abstractions map to high-level mathematical
  description of AMR algorithm components
• Reusable components. Component design based on
  mathematical abstractions to classes
• Based on public-domain standards
• MPI, HDF5
• Chombovis visualization package based on VTK,
  HDF5
• Layered hierarchical properly nested meshes
• Adaptivity in both space and time

35
Chombo Layered Design

• Chombo layers correspond to different levels of
  functionality in the AMR algorithm space
• Layer 1 basic data layout
• Multidimensional arrays and set calculus
• data on unions of rectangles mapped onto
  distributed memory
• Layer 2 operators that couple different levels
• conservative prolongation and restriction
• averaging between AMR levels
• interpolation of boundary conditions at
  coarse-fine interfaces
• refluxing to maintain conservation at coarse-fine
  interfaces
• Layer 3 implementation of multilevel control
  structures
• Berger-Oliger time stepping
• multigrid iteration
• Layer 4 complete PDE solvers
• Godunov methods for gas dynamics
• Ideal and single-fluid resistive MHD
• elliptic AMR solvers

36
Numerical Method Upwind Differencing

• The one-way wave equation propagating to the
  right
• When the wave equation is discretized upwind
  (i.e. using data at the old time level that comes
  from the left the wave equations becomes
• Advantages
• Physical The numerical scheme knows where the
  information is coming from
• Robustness The new value is a linear
  interpolation between two old values and
  therefore no new extrema are introduced

37
Symmetrizable MHD Equations

• The symmetrizable MHD equations can be written in
  a near-conservative form (Godunov, Numerical
  Methods for Mechanics of Continuum Media, 1,
  1972, Powell et al., J. Comput. Phys., vol 154,
  1999)
• Deviation from total conservative form is of the
  order of ??B truncation errors
• The symmetrizable MHD equations lead to the
  8-wave method. The eigenvalues are
• The fluid velocity advects both the entropy and
  div(B)in the 8-wave formulation

38
Numerical Method Finite Volume Approach

• Conservative (divergence) form of conservation
  laws
• Volume integral for computational cell
• Fluxes of mass, momentum, energy and magnetic
  field entering from one cell to another through
  cell interfaces are the essence of finite volume
  schemes. This is a Riemann problem.

39
Numerical method Riemann Solver

• The eigenvalues and eigenvectors of the Jacobian,
  dF/dU are at the heart of the Riemann solver
• Each wave is treated in an upwind manner
• The interface flux function is constructed from
  the individual upwind waves
• For each wave the artificial dissipation
  (necessary for stability) is proportional to the
  corresponding wave speed

• Discontinuous initial condition
• Interaction between two states
• Transport of mass, momentum, energy and magnetic
  flux through the interface due to waves
  propagating in the two media
• Riemann solver calculates interface fluxes from
  left and right states

40
Unsplit method Basic concept

• Original idea by P. Colella (Colella, J. Comput.
  Phys., Vol 87, 1990)
• Consider a two dimensional scalar advection
  equation
• Tracing back characteristics at t? t
• Expressed as predictor-corrector
• Second order in space and time
• Accounts for information propagating across
  corners of zone

Corner coupling
I
II
41
Unsplit method Hyperbolic conservation laws

• Hyperbolic conservation laws
• Expressed in primitive variables
• Require a second order estimate of fluxes

42
Unsplit method Hyperbolic conservation laws

• Compute the effect of normal derivative terms and
  source term on the extrapolation in space and
  time from cell centers to cell faces
• Compute estimates of Fd for computing 1D Flux
  derivatives ? Fd / ? xd

- I
I
43
Unsplit method Hyperbolic conservation laws

• Compute final correction to WI,,d due to final
  transverse derivatives
• Compute final estimate of fluxes
• Update the conserved quantities
• Procedure described for D2. For D3, we need
  additional corrections to account for (1,1,1)
  diagonal couplingD2 requires 4 Rieman solves
  per time stepD3 requires 12 Riemann solves per
  time step

II
44
The r B0 Problem

• Conservation of ??B 0
• Analytically if ??B 0 at t0 than it remains
  zero at all times
• Numerically In upwinding schemes the curl and
  div operators do not commute
• Approaches
• Purist Maxwells equations demand ??B 0
  exactly, so ??B must be zero numerically
• Modeler There is truncation error in components
  of B, so what is special in a particular
  discretized form of ??B?
• Purposes to control ??B numerically
• To improve accuracy
• To improve robustness
• To avoid unphysical effects (Parallel Lorentz
  force)

45
Approaches to address the r B0 constraint

• 8-wave formulation r B O(h?) (Powell et al,
  JCP 1999 Brackbill and Barnes, JCP 1980)
• Constrained Transport (Balsara Spicer JCP 1999,
  Dai Woodward JCP 1998, Evans Hawley Astro. J.
  1988)
• Field Interpolated/Flux Interpolated Constrained
  Transport
• Require a staggered representation of B
• Satisfy r B0 at cell centers using face values
  of B
• Constrained Transport/Central Difference (Toth
  JCP 2000)
• Flux Interpolated/Field Interpolated
• Satisfy r B0 at cell centers using cell
  centered B
• Projection Method
• Vector Potential (Claim CT/CD schemes can be
  cast as an underlying vector potential. Evans
  and Hawley, Astro. J. 1988)
• Require ad-hoc corrections to total energy
• May lead to numerical instability (e.g. negative
  pressure ad-hoc fix based on switching between
  total energy and entropy formulation by Balsara)

46
r B0 by Projection

• Compute the estimates to the fluxes Fn1/2i1/2,j
  using the unsplit formulation
• Use face-centered values of B to compute r B.
  Solve the Poisson equation r2? r B
• Correct B at faces BB-r?
• Correct the fluxes Fn1/2i1/2,j with projected
  values of B
• Update conservative variables using the fluxes
• The non-conservative source term S(U) ? r B has
  been algebraically removed
• On uniform Cartesian grids, projection provides
  the smallest correction to remove the divergence
  of B. (Toth, JCP 2000)
• Does the nature of the equations change?
• Hyperbolicity implies finite signal speed
• Do corrections to B via r2?r B violate
  hyperbolicity?
• Conservation implies that single isolated
  monopoles cannot occur. Numerical evidence
  suggests these occur in pairs which are spatially
  close.
• Corrections to B behave as ? 1/r2 in 2D and 1/r3
  in 3D
• Projection does not alter the order of accuracy
  of the upwinding scheme and is consistent

47
Unsplit Projection AMR Implementation

• Implemented the Unsplit method using CHOMBO
• Solenoidal B is achieved via projection, solving
  the elliptic equation r2?r B
• Solved using Multgrid on each level (union of
  rectangular meshes)
• Coarser level provides Dirichlet boundary
  condition for ?
• Requires O(h3) interpolation of coarser mesh ? on
  boundary of fine level
• a bottom smoother (conjugate gradient solver)
  is invoked when mesh cannot be coarsened
• Physical boundary conditions are Neumann d?/dn0
  (Reflecting) or Dirichlet
• Multigrid convergence is sensitive to block size
• Flux corrections at coarse-fine boundaries to
  maintain conservation
• A consequence of this step r B0 is violated on
  coarse meshes in cells adjacent to fine meshes.
• Code is parallel
• Second order accurate in space and time

48
Conclusion
(These, gentlemen, are the opinions upon which I
will base my facts Winston Churchill)

• Numerical evidence presented that the RMI is
  suppressed in the presence of a magnetic field
• Mathematical explanation
• Physical explanation still lacking
• Conjectures Main result will remain essentially
  same
• For non-ideal MHD
• For Hall MHD
• In three dimensions
• A conservative solenoidal B AMR MHD code was
  developed
• Unsplit upwinding method for hyperbolic fluxes
• r B0 achieved via projection

(If the facts dont fit the theory, change the
facts- Albert Einstein)
49
Future Work

• Linear stability analysis (a la Richtmyer)
• Self-similar solutions
• Triple point analysis
• Etc.

50
Related AMR Work

• High resolution parallel 2D magnetic reconnection
  runs.
• Implicit treatment of viscous/conductivity terms
• Pellet injection AMR simulations