ADVANCED MHD ALGORITHM FOR SOLAR AND SPACE SCIENCE

1
ADVANCED MHD ALGORITHM FOR SOLAR AND SPACE SCIENCE
2
GOALS

• Develop an efficient 3-D representation of the
  resistive MHD model on an unstructured grid of
  tetrahedra
• Truly arbitrary geometry
• Use cartesian coordinates
• Avoids coordinate singularities and complicated
  metrics
• Apply to a variety of problems
• Solar physics
• Structure and dynamics of active regions
• Coronal mass ejections
• Modeling of inner heliosphere
• Fusion
• Stellarators
• Incorporate adaptive mesh refinement

3
CHALLENGES

• Discrete representation of differential operators
• Compactness (coupling of nearest neighbor points
  only)
• Self-adjointness
• Annihilation properties (e.g., )
• Solution of implicit system
• Grid generation
• Implementation
• Code and data structure
• Parallelism
• Grid decomposition

4
RESISTIVE MHD MODEL
5
CURRENT AND MAGNETIC FIELD
Vector potential, magnetic field, and current
density
Both J and B are solenoidal Current density
operator is self-adjoint
Seek discrete operators that satisfy these same
conditions
6
TETRAHEDRAL GRID
7
FINITE VOLUME METHOD
8
APPLY TO MAGNETIC FIELD
9
DIVERGENCE OF B

• Apply Gauss theorem to dual median volume
  element surrounding vertex n

• After some algebra, contributions from common
  sides of adjoining tetrahedra cancel!

10
ALTERNATE DERIVATION OF B

• A varies linearly within tetrahedron

• Take the curl of this function

• Identical with finite volume result.

11
CURRENT DENSITY
12
CURL-CURL IS SELF-ADJOINT
13
VARIATIONAL PRINCIPLE
14
DISCRETE VARIATION
15
BOUNDARY CONSTRAINT FOR B

• Discrete minimization makes no reference to
  boundary conditions
• Discrete expression for curl curl operator is
  3Nn equations in 3Nn unknowns
• Could solve for all unknowns, including values at
  the Mn boundary vertices
• Absence of surface term implies that solution
  will satisfy the natural boundary conditions,
  i.e.,

• Since At is not fixed, this implies that

• Constraint on surface field and volume current

In general, we must specify At on the boundary
16
PLACEMENT OF VARIABLES ON GRID

• Vertices
• A, J, v
• Centroids
• r, p, B
• Velocity averaged to faces or centroids, as
  required
• Apply conservation laws to control
  volume
• Equation of motion not in conservation form
• Use anisotropic semi-implicit operator

17
ADVECTION

• Control surfaces for upwind advection

Cell centered quantity
Vertex centered quantity
18
TIME SCALES IN RESISTIVE MHD
Explicit time step impractical
Require implicit methods
19
PARTIALLY IMPLICIT TIME DIFFERENCING

• MHD operator contains widely separated time
  scales (eigenvalues)

• Treat only fast part of operator implicitly to
  avoid time step restriction

• Precise decomposition for complex nonlinear
  system is often difficult or impractical to
  achieve

20
OPERATOR SPLITTING

• In MHD, F and W are known, but an expression for
  S is difficult to achieve
• Use operator splitting

• Explicit expression for S is not required

21
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 weather modeling
• Has proven to be very useful for resistive and
  extended MHD

22
SEMI-IMPLICIT OPERATOR FOR MHD

• Linearized, ideal MHD wave equation
• 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 1000

23
SEMI-IMPLICIT OPERATOR
24
DISCRETE SEMI-IMPLICIT OPERATOR
25
COMPUTING ISSUES

• F90 implementation
• Use object-oriented features
• Facilitate code modification and maintenance
• Use existing software implementations
• MPI for parallelism
• LaGrit (LaGrit Team, 1999) for mesh generation
• METIS (Karypis Kumar, 1999) for partitioning
  grid among processors
• PETSc (Baley, et al., 2000) for preconditioned CG
  solver on unstructured grid
• GMV (Ortega, 2000) for visualization of data on
  tetrahedral grid
• Expedited code development

26
EXAMPLE GRID DECOMPOSITION

• Decomposition of cubic, cylindrical, and
  spherical domains for parallel processing using
  METIS

27
EXAMPLE POTENTIAL CORONAL FIELD
28
EXAMPLE CORONAL POTENTIAL FIELD
29
EXAMPLE LINEAR SOUND WAVES
30
SOUND WAVES IN A BOX
X-Component of velocity
Pressure
31
SOUND WAVES IN A SPHERE
32
NONLINEAR SHOCK PROBLEM
G. A. Sod, J. Comp. Phys. 27,1 (1978)

• Diaphragm separating left and right states of
  fluid
• Diaphragm is broken at t 0
• Expansion fan moves to left
• Shock and contact discontinuity move to right
• Well documented nonlinear solution of
  hydrodynamic equations

33
NONLINEAR SHOCK PROBLEM
Temporal evolution of the density
34
MHD TORSIONAL ALFVEN WAVES
35
MHD TORSIONAL ALFVEN WAVES
Magnetic Energy
Perturbed magnetic field vectors
36
MHD NON LINEAR KINK MODE
68441 nodes, 398948 cells, 16 processors
37
MHD NON LINEAR KINK MODE
Kinetic energy
Initial conditions unstable Gold-Hoyle
equilibrium
Magnetic energy
At t0 a random perturbation Is applied and the
m1 kink instability is triggered
38
MHD NON LINEAR KINK MODE
39
MHD SHOCK IN CYLINDRICAL COORDINATES
Modified from Brio, M. and Wu, C. C., J. Comp.
Phys. 75, 400 (1988), and adapted to cylindrical
geometry

• Diaphragm separating left and right states of
  fluid
• Diaphragm is broken at t 0
• Fast rarefaction and slow compound waves move to
  left
• Slow shock, contact discontinuity, and fast
  rarefaction wave move to right.

40
MHD SHOCK IN CYLINDRICAL COORDINATES
482007 nodes, 2717151 cells, 16 processors
41
MHD SHOCK IN CYLINDRICAL COORDINATES
Cutlines at t1
42
MHD SHOCK IN CYLINDRICAL COORDINATES
Cutplane of density at t1
43
THE SOLAR WIND FROM 30R? TO 5 A.U.

• We simulate the propagation of the hydrodynamic
  solar wind in the heliosphere.
• The mesh consists of 148596 nodes and 875520
  cells and extends from 30R? to 5 A.U.
• At 30R? we specify the boundary conditions a
  30?-degree-wide belt of dense and slow solar wind
  inclined of 20? degrees in respect to the
  rotation axis, surrounded by the fast solar wind.
  
• The angular rotation speed is 14 ? degrees per
  day.
• We advance the hydrodynamic equations for 30 days.

44
THE SOLAR WIND FROM 30R? TO 5 A.U.
A cut of the mesh and an enlargement showing the
inner boundary
45
THE SOLAR WIND FROM 30R? TO 5 A.U.
Cutplanes of the flow speed
Cutplanes of density times r2
46
THE SOLAR WIND FROM 30R? TO 5 A.U.
Enhanced density regions near the ecliptic plane
47
MH4D STATUS

• Formulated discrete algorithm for resistive MHD
  on a tetrahedral grid
• Based on variational principle
• Compact, self-adjoint, etc.
• Implicit viscosity and resistivity
• Used available tools for implementation (F90,
  LaGrit, METIS, PETSc, GMV
• Expedited development schedule
• Validation
• Potential coronal field computed from boundary
  data
• Linear sound waves in cubic and spherical domains
• Nonlinear shock tube problem
• Linear torsional Alfvén waves in a cylinder
• Nonlinear MHD shock problem
• Propagation of the supersonic solar wind in the
  heliosphere
• Next steps
• Optimize preconditioners
• Apply to solar and heliospheric problems
• Adaptive mesh refinement (AMR)
• Implement web page
• Goal Distribute code to user community as open
  source project