# Presentation: 475643

## Adaptive Mesh Simulations of Pellet Injection in Tokamaks

- Ravi Samtaney
- Computational Plasma Physics Group
- Princeton Plasma Physics LaboratoryPrinceton University
- SIAM Annual Conference
- July 10-14 2006, Boston, MA
- Acknowledgement: DOE SciDAC Program

## Collaborators

- P. Colella and Applied Numerical Algorithms Group (LBNL)
- S. C. Jardin (PPPL)
- P. Parks (GA)
- D. Reynolds (UCSD)
- C. Woodward (LLNL)
- Funded through the TOPS CEMM and APDEC SciDAC projects.  RS supported by US DOE Contract No. DE-AC020-76-CH03073

## Outline

- Introduction and motivation
- Description of physical phenomenon
*Spatial and temporal scales *
- Equations and models
- Adaptive mesh refinement (AMR) for shaped plasma in flux-surface coordinates
- Results
*HFS vs. LFS Pellet injection*
- Newton-Krylov fully implicit method
- Future directions and conclusion

## Pellet Injection & Edge Localized Modes

- Motivation
*Injection of frozen hydrogen pellets is a viable method of fueling a tokamak *
*Presently there is no satisfactory simulation or comprehensive predictive model for pellet injection (esp. for  ITER )*
*H-mode operation of ITER will be accompanied by edge localized modes (ELMS) **(ITER Physics Experts Group,Nucl. Fusion 1999)*
*Pellet injection related to ELMS **(Gohill et al. PRL, 2001; Lang et al. Nucl. Fusion 2000)** *
- Objectives
*Develop a comprehensive simulation capability for pellet injection and ELMs in tokamaks (esp. ITER) with modern technologies such as adaptive mesh refinement for spatial resolution and fully implicit Newton-Krylov approach for temporal stiffness*

***Pellet injection in TFTR***

***HFS***

***LFS***

## Physical Processes: Description

- _Non-local_ electron transport along field lines rapidly heats the pellet cloud (_****__**e**_).
*Frozen pellet encounters hot plasma and ablates rapidly*
*Neutral gas surrounding the solid pellet is ionized *
*Ionized, but cool plasma, continues to get heated by electrons*
*A high **** “plasmoid” is created*
- Ionized plasmoid expands
*Fast magnetosonic time scale *_******__**f.**_
- Pellet mass moves across flux surfaces _****__**a.**_
*So-called “anomalous” transport across flux surfaces is accompanied by reconnection*
- Pellet cloud expands along field lines _****__**c**_.
*Pellet mass distribution continues along field lines until pressure equilibration*
- Pellet lifetime p

**Figure from M****ü****ller et al., Nuclear Fusion 42 (2002)**

## Scales and Resolution Requirements

- Time Scales e < f  < a < c   < p
- Spatial scales: Pellet radius rp << Device size L ~O(10-3)
- Presence of magnetic reconnection further complicates things
*Thickness of resistive layer scales with ~ ****1/2*
*Time scale for reconnection is ~ ****-*1/2
- Pellet cloud density ~ O(104) times ambient plasma density
- Electron heat flux is non-local
- Large pressure and density gradients in the vicinity of cloud
- Pellet lifetime ~ O(10-3) s long time integrations
_**Resolution estimates**_

| Tokamak | MajorRadius | N |     Nsteps | Spacetime Points |
| --- | --- | --- | --- | --- |
| CDXU (Small) | 0.3  | 2 x 107  |   2 x 105 | 4 x 1012 |
| DIIID (Medium) | 1.75 | 3.3 x 109 |   7 x 106 | 2.3 x 1017 |
| ITER (Large) | 6.2 | 1.5 x 1011 | 9 x 107  | 1.4 x 1019 |

## Related Work - Local vs. Global Simulations

- Earliest ablation model by Parks (Phys. Fluids 1978)
- Detailed multi-phase calculations in 2D of pellet ablation (MacAulay, PhD thesis, Princeton Univ 1993, Nuclear Fusion 1994)
- Detailed 2D Simulations of pellet ablation by Ishizaki, Parks et al. (Phys. Plasmas 2004)
*Included atomic processes – ablation, dissociation, ionization, pellet fluidization  and distortion; semi-analytical model for electron heat flux from background plasma *
- In above studies, the domain of investigation was restricted to only a few cm around the pellet
*Also, in these studies the magnetic field was *_*static*_
- 3D Simulations by Strauss and Park (Phys. Plasmas, 1998)
*Solve an initial value problem. Initial condition consisted of a density “blob” to mimic a *_*fully ablated*_* pellet cloud which, compared with device scales, was relatively large due to resolution restrictions*
*No motion of pellet modeled*
- 3D Adaptive Mesh Simulation of pellet injection by Samtaney et al. (Comput. Phys. Comm, 2004)

## Current Work

- Combine global MHD simulations in a tokamak geometry with detailed local physics including ablation, ionization and electron heating  in the neighborhood of the pellet
- AMR techniques to mitigate the complexity of the multiple scales in the problem
- Newton-Krylov approach for wide range of temporal scales

## Mathematical Model 

- Mathematical Model

- Additional constraint r¢ **B **=0*** ***

***Hyperbolic terms***

***Diffusive terms***

***Density: Ablation******Energy :Electron heat flux***

## Ablation Model 

*Above equation uses cgs units*
- Abalation occurs on the pellet surface
*Regularized as a truncated Gaussian of width 10 r**p*
*Pellet shape is spherical for all t*
*Pellet trajectory is specified as either HFS or LFS*
*Monte Carlo integration to determine average source in each finite volume *

- Ablation Model

## Electron Heat Flux Model 

*Assumes Maxwellian electrons and neglects pitch angle scattering*
- Solve for opacities as a “steady-state” solution to an advection-reaction equation
*Solve by using an upwind**method*
*Advection velocity is ****b***
- Ansatz for energy conservation
*    **Sink term on **     flux surface**     outside cloud*

- Electron Heat Flux Model

## Curvilinear coordinates for shaped plasma

- Adopt a flux-tube coordinate system (flux surfaces  are determined from a separate equilibrium calculation)
*R**** ́**** R (****, ****), and Z**** ́**** Z (****, ****)*
-   ** ́*** **** (R,Z), and **** **** ́**** ****(R,Z)*
-  Flux surfaces:  *= ***0* ***
-   *coordinate is retained as before*
- Equations in transformed coordinates

## Numerical method

- Finite volume approach
- Explicit second order or third order TVD Runge-Kutta time stepping
- The hyperbolic fluxes are evaluated using upwinding methods
*seven-wave Riemann solver:  **F**** ́**** F(U**L**, U**R**) = 1⁄2(F(U**L**) + F(U**R**) - ****k** ****k ** | ****k** | r**k **) where ****k** = l**k** (U**R** – U**L**)*
*Harten-Lee-vanLeer (HLL) Method **(SIAM Review 1983)** **F**** ́**** F(U**L**, U**R**) = ****min**F(U**L**) + ****max**F(U**R**) +****min** ***max*(U**R** – U**L**) /(***max*-***min*)*
- Diffusive fluxes computed using standard second order central differences
- The solenoidal condition on B
*imposed using the *Central Difference* version of *_*Constrained Transport*_* **(Toth JCP 161, 2000)*
*Including the non-conservative source term in the equation to advect ****r¢**** B errors **(Powell et al. , JCP 1999)*
*By projection at n+1/2 time step **(Samtaney et al., SciDAC 2005)*
***r**** ****¢**** ****B ********** 0 on coarse mesh cells adjacent to coarse-fine interfaces*
_**Initial Conditions:**_  Express B=1/R( **£** **r**  + g() )  fnc(). Initial state is an MHD equilibrium obtained from a Grad-Shafranov solver.
_**Boundary Conditions**_: Perfectly conducting for =o,  zero flux (due to zero area) at =i, and periodic in  and .

## 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](http://www.seesar.lbl.gov/ANAG/chombo)_)
*(Chombo is an AMR developer’s toolkit)*
- Adaptivity in both space and time
- Mesh generation: necessary to ensure volume preservation and areas of faces upon refinement
- Flux-refluxing step at end of time step ensures conservation

******

******

## Pellet Injection: AMR

- Meshes clustered around pellet
- Computational space mesh structure shown on right
- Mesh stats
*32**3** – base mesh with 5 levels, and refinement factor 2*
*Effective resolution: 1024**3*
*Total number of finite volume cells:113408*
*Finest mesh covers 0.015 % of the total volume*
*Time adaptivity: **	**1 (**** t)**base**=32 (**** t)*finest

## Pellet Injection: Zoom into Pellet Region

## Pellet Injection: Zoom in

## Pellet Injection: Pellet in Finest Mesh

## Pellet Injection: Pellet Cloud Density 

******

******

******

## Results - HFS vs. LFS

- BT = 0.375T
- n0=1.5**£** 1019/m3
- Te1=1.3Kev
- =0.05
- R0=1m, a=0.3 m
- Pellet: rp=1mm, vp=1000m/s

******

***t=100***

***t=7***

***t=256***

## HFS vs. LFS - Average Density Profiles 

***HFS Pellet injection shows better core fueling than LFS***
*Arrows indicate average pellet location*

***Core***

***Edge***

## HFS vs. LFS: Instantaneous Density Profiles

***Radially outward shift in ******both cases indicates higher ******fueling effectiveness for HFS***

*********=0***

*********=************/4***

*********=0***

*********=************/4***

*********=0***

*********=************/4***

## Pellet Injection: LFS/HFS Launch

***Density******Instantaneous temp equilibration on flux surfaces***

## Time step set using explicit CFL condition of fastest wave:
Pellet Injection:  pellet radius rp = 0.3 mm, injection velocity vp = 450 m/s, fast magneto-acoustic speed cf 1⁄4 106 m/s:
To resolve pellet need O(107) time steps
Longer time steps (implicit methods) are a practical necessity
Fixed time step, two-level -scheme using a Jacobian-Free Newton-Krylov nonlinear solver [KINSOL]: 
f(Un) = Un – Un-1 – t [ g(Un) + (1- g(Un-1)],     g(U) = r¢(Fp(U) – Fh(U))
 = 1  )  Backward Euler [O(t)];  =0.5  )  Cranck-Nicholson [O(t2)]
Adaptive time step, adaptive order, BDF method for an up to 5th order accurate implicit scheme [CVODE]:
f(Un) = Un – i=1:qn,i Un-i – tn n,1 g(Un-1) – tn 0 g(Un)
Time step size and order adaptively chosen based on heuristics balancing accuracy, nonlinear & linear convergence, stability
Joint work with Dan Reynolds (UCSD) and Carol Woodward (LLNL)

- Pellet Injection:  pellet radius *r**p* = 0.3 mm, injection velocity *v**p* = 450 m/s, fast magneto-acoustic speed *c**f* 1⁄4 106 m/s:
*To resolve pellet *_*need O(10*__*7*__*) time steps*_
- Longer time steps (implicit methods) are a practical necessity
- Fixed time step, two-level **-scheme using a Jacobian-Free Newton-Krylov nonlinear solver [KINSOL]:
- f*(*Un*)** **=** *Un* **–** *Un-1* **–** ***t* **[** *g*(*Un*)** **+** **(1-*** g*(*Un-1*)**]**,     *g*(*U*)** **=** **r¢**(*Fp*(*U)* **–** *Fh*(*U*)**)*
- * **=** **1  **)**  **Backward Euler **[O(****t)];**  ****=****0.5  **)**  **Cranck-Nicholson** [O(****t**2**)]*
- Adaptive time step, adaptive order, BDF method for an up to 5th order accurate implicit scheme [CVODE]:
- f*(*Un*)** **=** *Un* **–** ***i*=*1*:*qn,i* *Un-i* **–** ***tn* *n,1* *g*(*Un-1*)** **–** ***tn* *0* *g*(*Un*)*
*Time step size and order adaptively chosen based on heuristics balancing accuracy, nonlinear & linear convergence, stability*
- Joint work with Dan Reynolds (UCSD) and Carol Woodward (LLNL)

- JFNK Fully Implicit Approach for Resistive MHD

## Choose a model problem with asimilar separation of time scales(Reyolds et al. JCP 2006)

- Pellet Injection - Implicit Simulations

- Good agreement between explicit and implicit methods

- Implicit (no preconditioners) overtakesexplicit method as problem size getslarger.

-  Implicit simulations in a toroidalgeometry.  t = 100  texplicit

## Summary and Future Plan

- Preliminary results presented from an AMR MHD code
*Physics of non-local electron heat flux included*
*HFS vs. LFS pellet launches *
*HFS core fueling is more effective than LFS*
*Numerical method is upwind, conservative and preserves the solenoidal property of the magnetic field*
-  AMR is a practical necessity to simulate pellet injection in a tokamak with detailed local physics
- Preliminary results from a fully implicit Newton-Krylov method for pellet injection in tokamaks
- Future work
*Physics-based pre-conditoners for fully implicit JFNK method*
_*Proposed work under SciDAC-2*_*: Combine adaptive and fully implicit methods to manage the wide range of spatial and temporal scales*