# Presentation: 403858

## Spectral Elements for Anisotropic Diffusion and Incompressible MHD

- Paul Fischer
- Argonne National Laboratory

## Terascale Simulation Tools and Technologies

*Goal:* Enable high-fidelity calculations based on multiple coupled physical processes and multiple physical scales
*Adaptive methods*
*Advanced meshing strategies*
*High-order discretization strategies*
*Technical Approach:*
- Develop interchangeable and interoperable software components for meshing and discretization
- Push state-of-the-art in discretizations
*   *

## Terascale Simulation Tools and Technologies

*Goal:* Enable high-fidelity calculations based on multiple coupled physical processes and multiple physical scales
*Adaptive methods*
*Advanced meshing strategies*
*High-order discretization strategies*
*Technical Approach:*
- Develop interchangeable and interoperable software components for meshing and discretization
- Push state-of-the-art in discretizations
*   *

## Outline

- Driving physics
- SEM overview
- Costs
- MHD formulation
- Convective transport properties of SEM
- Anisotropic Diffusion
- Brief summary of current results on MHD project
- Conclusions

## Driving Physics

*Anisotropic Diffusion:                                   *                   ( S. Jardin,  PPPL)
*Central to sustaining high plasma temperatures*
*Want to avoid radial leakage*
*Need to understand significant instabilities*
*Incompressible MHD:** **	**	**           *( F. Cattaneo UC, H. Ji PPPL, ... )
*Several liquid metal experiments are under development to understand momentum transport in accretion disks*
*Magneto rotational instability (MRI) is proposed as a mechanism *
*    **to initiate turbulence capable of generating transport*
*The magnetic Prandtl number for liquid metals is ~ 10**-5** *
*	**** **Re=10**6 **& Rm=10 for experiments*
*Numerically, we can achieve Re=10**4 **& Rm=10**3                  *( 2005 INCITE award )
*Free-surface MHD**	**	**	**	**                      *( H. Ji, PPPL )
*Proposed as a plasma-facing material for fusion*
*Desire to understand the effect of ****B**** on the  free surface*

## Spectral Element Overview 

*Spectral elements can be viewed as a finite element subset*
*Performance gains realized by using:*
*tensor-product bases** (quadrilateral or brick elements)*
*Lagrangian bases collocated with GLL quadrature points*
*Operator Costs:**	**	**Standard  FEM***** **	**	**  SEM*
*	**	**# memory accesses:**	**    **O( **EN** 6 **) **	**	**O( **EN** 3 **)**	*
*	**	**# operations:**	**	**    **O( **EN** 6 **) **	**	**O( **EN** 4 **)*
*	**	**N** **=** order, **E** **=** number of elements, ** EN **3 **= **number of gridpoints*
*Dramatic cost reductions for large **N** ( > 5 )*

## Computational Advantages of Spectral Elements

- Exponential convergence with *N*
- minimal numerical dispersion / diffusion
- for anisotropic diffusion, lements of ker(A||) well-resolved         – sharp decoupling of isotropic and anisotropic modes
- Matrix-free form
- matrix-vector products cast as efficient matrix-matrix products
- number of memory accesses identical to 7-pt. finite difference
- no additional work or storage for anisotropic diffusion tensor

##  SEM Computational Kernel on Cached-Based Architectures

- matrix-matrix products, *C=AB**: **2N** 3* ops for *2N** 2* memory references
*much of the additional work of the SEM is covered by efficient use of cache  – e.g., time for A*B vs A+B for N=10 is:*
*     *
*	**2.0 x on DEC Alpha,                                                 1.5 x on IBM SP*

## Incompressible MHD

- t

*— **plus appropriate boundary conditions on ****u ****and**** B***
- Typically,* Re >> Rm >> 1*
- Semi-implicit formulation yields independent Stokes problems for **u **and** B**

## Incompressible MHD in a Nutshell

- t

- Convection:
*dominates transport *
*dominates accuracy requirements*
*often the challenging part of the discretization*
*treated explicitly in time*
- Diffusion:
*“**easy”     *( ?? )
- Projection:  *div ****u**** = 0   div ****B**** = 0*
*dominates work*
*isotropic SPD operator*
*multiple right-hand side information*
*scalable multilevel Schwarz methods**	**	*( 1999 GB award )
*SE multigrid**	**	**	**	**	*( Lottes & F 05 )

## High-Order Methods for Convection-Dominated Flows

- 	   Phase Error for *h* vs. *p* Refinement:  *u**t** + u**x** = 0** *
*	**	**h-refinement*                                     *p-refinement*

## High-Order Methods for Convection-Dominated Flows

*Fraction *of accurately resolved modes  (*per space direction*) is increased only through increased order
*Savings cubed in** **R**3*
*Rate *of convergence is extremely rapid for high *N*
*Important for multiscale / multiphysics problems*
*	**( Q: Why do we want 10** 9** gridpoints? )*
*Stability *issues are now largely understood
*stabilization via DG, filtering, etc.*
*dealiasing*
*Still*, *must* *resolve* structures (no free lunch...)
*Computational costs are somewhat higher*
*Data access costs are equivalent to finite differences*

## Stabilizing convective problems:			

- Stabilizing convective problems:

*Models of straining and rotating flows*:
- Rotational case is skew-symmetric.
*Filtering* attacks the leading-order unstable mode.
*Dealiasing*  (*high-order quadrature*) yields imaginary eigenvalues – *vital for MHD*
* **	**	**	**	**                     **N=19, M=19              N=19, M=20*
*                                                                              ***

*straining*
*field *
*rotational*
*field *

- Diffusion – easy ??

## CEMM Challenge Problems 		S. Jardin, PPPL

*Anisotropic diffusion in a toroidal geometry*
*Two-dimensional tilt mode instability*
*Magnetic reconnection in 2D*
- Provides a problem suite that
- captures essential physics of fusion simulation
- stresses traditional numerical approaches
- identifies pathways for next generation fusion codes
- Excellent vehicle for initiating SciDAC interatctions.

## Anisotropic Diffusion in Toroidal Domains

**b** – normalized ***B***-field, helically       wrapped on toroidal surfaces
- thermal flux follows **b**.

## High degree of anisotropy creates significant challenges 

- For                            this problem is more like a (difficult) hyperbolic problem than straightforward diffusion.
- This is reflected in the variational
- 	statement for the steady case with

## Numerical challenges:
radial diffusion,
nearly singular, with large (but finite) null space
avoid grid imprinting
unsteady case constitutes a stiff relaxation problem
preconditioning nearly singular systems

CEMM challenge:
establish spatial convergence for steady state case
check unsteady energy conservation when 
investigate the behavior of the tearing mode instability

- radial diffusion,
- nearly singular, with large (but finite) null space
- avoid grid imprinting
- unsteady case constitutes a stiff relaxation problem
- preconditioning nearly singular systems
- CEMM challenge:
- establish spatial convergence for steady state case
- check unsteady energy conservation when
- investigate the behavior of the tearing mode instability

- High degree of anisotropy creates significant challenges

- High degree of anisotropy creates significant challenges

## It is advantageous to use few elements of high order

*Fewer gridpoints are required*
*CPU time proportional to number of gridpoints*

- (*N* odd, 3-7)

- (*N* even, 2-6)

## High Anisotropy Demands High Accuracy

*A = ****|| **A**||** +** **A**I** *
*A**I ** *controls *radial* diffusion
*A**||** must be accurately represented when ****|| **>> 1*
*Error must scale as ~ 1 / ****||** *

## High Anisotropy Demands High Accuracy

- Difficulty stems from high-frequency content in null space of *A*||
- High-order discretizations are able to accurately represent these functions.

## Evolution of Gaussian Pulse for 

- Evolution of Gaussian Pulse for

*minimal radial diffusion*
*no grid imprinting *
*careful time integration required (e.g., adaptive **DIRK4**)**  *

***-averaged temperature vs. time*

- Evolution of Gaussian pulse for

## SEM is able to identify critical physics – 

*Tearing mode instability*
- radial perturbation:
***b**** = ****b******0**** + **** cos(m****-n****) ****r***
- field lines do not close     on *m-n* rational surface
- magnetic island results, with significant increase    in radial conductivity.
- island width scales as:
- W ~                   in accord with asymptotic theory
- 	Note: this is a subtle effect!

- W

*Island width vs. ****|| **at onset*.

- W

## Outstanding Challenges for Anisotropic Diffusion Simulation

- Preconditioning
- need null-space control
- condition number scales as **||
- Non-aligned grids predict early island formation

***-averaged temperature vs. time*

- Incompressible MHD Results

## Axisymmetric Hydro Simulations of Taylor-Couette w/ Rings

- Re=620 steady
- Re=6200 unsteady
- Axisymmetric MHD simulations are being carried out now.
- Starting point for 3D simulations, which are being compared with experiments at PPPL.

*Computation by Obabko, Fischer, & Cattaneo*

*Normalized Torque                   Vorticity*

*inner cylinder*
*outer cylinder*

*Re=6200*

## Computational MRI: preliminary results

**Computations Fischer, Obabko & Cattaneo**

- Nonlinear development of Magneto-Rotational Instability
- Cylindrical geometry similar to Goodman-Ji experiment
- Hydrodynamically stable rotation profile
- Weak vertical field
- Use newly developed spectral element MHD code
- Try to understand differences between experiments and simulations
- Simulations ReRm (moderate). Experiments Re>>Rm (Rm smallish)

## Summary & Conclusions

- Block-structured SEM provides an efficient path to high-order
- accurate treatment of challenging physics
- fast cache-friendly operator evaluation ( *N** 4* vs. *N** 6* )
- For anisotropic diffusion
- effects of grid imprinting are minimized
- able to capture physics of high anisotropy
- MRI experiment
- MRI has been observed with axial periodicity at Re=Rm=1000.
- preliminary axisymmetric results indicate hydrodynamic unsteadiness at Re=6000 for two-ring boundary configuration, may be mitigated in 3D...
- Free-surface MHD
- Free-surface NS is working
- Coupling with full MHD is underway