# Presentation: 301549

## An Implicit Method for Magnetic Fusion Extended MHD using Adaptive, High-Order, High-Continuity Finite Elements

- SIAM Annual Meeting
- Boston, MA
- July 13 2006

- S. C. Jardin, J. Breslau, N. Ferraro
*Princeton Plasma Physics Laboratory*
- Bauer, M. Shephard
*Renssalier Polytechnic Institute*

## The Extended MHD equations for a magnetized (fusion) plasma are a high-order system of 8 scalar variables that are characterized by a wide range of space and timescales.  

Our approach is as follows:

Multiple space scales   	unstructured adaptive elements
Multiple time scales 	  	implicit time differencing
High order derivatives 	C1 continuity elements (up to 4th order)
8 scalar variables	  	split implicit time advance & compact rep.
Strong magnetic field	  	stream function/potential representation

- Our approach is as follows:
- Multiple space scales   	unstructured adaptive elements
- Multiple time scales 	  	implicit time differencing
- High order derivatives 	*C**1 *continuity elements (up to 4th order)
- 8 scalar variables	  	split implicit time advance & compact rep.
- Strong magnetic field	  	stream function/potential representation

## Extended MHD Equations:

## Scalar data is represented using 18 degree of freedom quintic triangular finite elements  Q18

-   All data is at nodes:  function + first 5 derivatives (6 dof)
-   Complete quintic polynomial has 21 coefficients
-   18 values come from the 3 nodes (3 x 6)
-    3 values come from requirement that the normal derivative along each edge be only a (univariate) cubic....leads to *C**1* continuity
-   Contains a complete Taylor series through 4th order...error ~ *h**5*
*  *Compact representation ... only 3 dof/triangle
-  *C**1* continuity allows up to 4th derivatives in space without introducing auxiliary variables
-   Unstructured triangular mesh allows adaptive zoning

## Implicit velocity time-advance substitutes in from field equations to contain all Ideal MHD wave phenomena

- Use SuperLU to invert this linear operator to get from time n to (n+1)

## A similar technique is used on the magnetic field equations.  Fully implicit Extended MHD (2-fluid) equations-- time step determined by accuracy only:

- Alfven Wave physics

- Whistler, KAW, field diffusion physics

- density

- pressure

-  4 sequential matrix solves per time step
-  3 non-trivial subsets with 6,4,2 variables

## GEM Nonlinear Benchmark

- Poloidal Flux 

- Current Density J

- GEM Reconnection Problem

- Resistive MHD
- High and Low Viscosity
-  ( = 10 ,  = 0.1 )
- Two-Fluid

-  Provides a non-trivial, convenient test problem for code verification and validation and cross-code comparison
-  Also, extending this by adding an equilibrium magnetic field into the plane (guide field)

## Resistive MHD gives convergent results.  Time step depends on accuracy requirement only

-   Dependence of kinetic energy on viscosity (varies by 100)
-   Lower viscosity cases require smaller timestep for accuracy
-   60 x 60 grid gives adequate spatial resolution

## Current Density contours for resistive MHD case

- Resolution requirements are modest for resistive MHD;   = .005, =10

## 2-fluid reconnection qualitatively different,  requires high resolution for convergent results

-  Note sudden transition where velocity abruptly increases
- These calculations used a hyperviscosity term in Ohm’s law proportional to (x)2 . . .     required for a stable calculation

## Reconnection velocity field contains both incompressible and compressible components

-    Without a strong (guide) magnetic field,  the velocity field is not close to incompressible during reconnection

## Current Density contours for 2-fluid MHD 

- Current Density contours for 2-fluid MHD

-   Starts like resistive MHD
-   Dramatic change in configuration for t > 20

## Close-up of 2-fluid current density at  t=32

-  Note very localized region of high current density in center

- midplane

- These calculations did not assume any symmetry, except for initial and boundary conditions

## Midplane Current density collapses to the width of 1-3 triangular elements

-                t=32 time of previous contour plot
- ( note sudden collapse at t=24)

## Reconnection rate and maximum velocity for resistive and 2-fluid cases

-   note transition at     t ~ 24 for 2-fluid case
-   After transition,     2-fluid reconnection rate is over 50 times larger than resistive MHD rate

## Midplane electric field before and after transition

- Reconnection rate:

## Midplane electric field before and after transition

- one triangle

- Hyper-resistivity coefficient must be large enough that current density collapse is limited to 1-2 triangles:  *reason for factor (****x)**2*

## Test of sensitivity to hyperresistivity coefficient H

-  These calculations had (120)2 grid, everything else fixed
-  In caption is value of of ***H**(****x)**2*
-  Appears to be converging to a unique answer for ***H **(****x)**2** * 0
-  Need for small value of ***H** (****x)**2* implies need for small x                                                                                                                            (to avoid current layer collapsing to less than one zone width)

## Effect of adding a guide field on 2-fluid reconnection

-   Adding a guide field delays transition time and reduces maximum reconnection rate
-   Small effect for B0 < 0.2
-   To be studied further

## Change in velocity field with toroidal field strength

-   Velocity field becomes more like incompressible flow as toroidal field strength increases

## Adaptive Meshing

- Andy Bauer (RPI) has implemented an arbitrary Adapted Mesh in the M3D-*C**1** *code and is exploring different adaptive strategies.  This should greatly improve the efficiency of the 2-fluid reconnection problem.

## Summary and Conclusions

- Full 8-field E-MHD equations solved in 2D slab geometry with stream function/potential form
- Q18 elements allow high accuracy, compact representation of solution
- Implicit solution technique implies time step can be set by accuracy requirements, not stability
- 2-fluid reconnection problems require localized regions with high resolution...natural for adaptive refinement
- Extensions to toroidal geometry and 3D underway

## Supplementary vgs if time permits

## Tests of the 2-fluid and Gyroviscous terms:Gravitational Instability

**x**

**y**

## Braginskii gyro-viscosity in M3D-C1:

- Ramos

- Breslau

- Braginskii gyro-viscosity in M3D-*C**1*:

## M3D-C1  Gravitational Instability stabilized by Gyroviscosity

- M3D-*C**1*  Gravitational Instability stabilized by Gyroviscosity

-  Ferraro

- We have calculated the stabilization of the gravitational instability by Gyroviscosity:   low beta (left) and high beta (right)

## M3D-C1...Gravitational Instability:  nonlinear 

- Ferraro

## Anisotropic Diffusion

- Solve to steady state

- Shows greater than N-5 convergence

## Non-linear evolution of tilting cylinder in full 6-field 2-fluid model

- J:    t=0.8

- J:    t=3.2

- J:    t=4.0

- J:    t=4.8

- :    t=0.8

- :    t=3.8

- :    t=4.0

- :    t=4.8

## Adaptive Meshing

- Andy Bauer (RPI) has implemented an arbitrary Adapted Mesh in the M3D-C1 code and is exploring different adaptive strategies