# Presentation: 475642

## Tokamak Pellet Fueling Simulations

- Ravi Samtaney
- Computational Plasma Physics Group
- Princeton Plasma Physics LaboratoryPrinceton University
- SIAM Conference on Parallel Processing
- February 22-24 2006, San Francisco, CA
- 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)
- Special thanks to Mark Adams (Columbia)
- 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: Objective and Motivation

- 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 )*
- Objectives
*Develop a comprehensive simulation capability for pellet injection into tokamaks*
*Quantify the differences between “inside launch”  (HFS) and “outside launch” (LFS)*
- Approach
*Adaptive mesh refinement for large range of spatial scales*
*Implicit Newton-Krylov for large range of temporal scales*

***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

## Equations and Models 

- Equations and Models

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

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

***Hyperbolic terms***

***Diffusive terms***

- Mass source is given using the ablation model by Parks and Turnbull (Phy. Plasmas 1978) and Kuteev (Nuclear Fusion 1995)
*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 *

## Electron Heat Flux Model 

*Assumes Maxwellian electrons and neglects pitch angle scattering*
*Where                                        ,                  and*
- 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*
*Harten-Lee-vanLeer (HLL) Method (SIAM Review 1983) *
- Diffusive fluxes computed using standard second order central differences
- The solenoidal condition on B is imposed using the *Central Difference* version of _Constrained Transport_ (Toth JCP 161, 2000)
***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***

## Limitations of Explicit Approach

- Time step set using explicit CFL condition of fastest wave:
- 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 for a mesh with ***x* = 10**-4** m, ***t* = 3.3 x 10**-11** s:    *_*need O(10*__*7*__*) time steps*_
- Longer time steps (implicit methods) are a practical necessity!

## Adopt a Jacobian-Free Newton-Krylov Approach
The choice of implicit method affects stability, accuracy, extensibility

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)]
Pushes difficulty of larger time step sizes onto solver

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

- The choice of implicit method affects stability, accuracy, extensibility
- 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**)]*
*Pushes difficulty of larger time step sizes onto solver*
- 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*

- Fully Implicit Approach (Reynolds, Woodward)

## Choose a model problem with the similar separation of time scales
Single fluid resisitive MHD
Instantaneous heating by electrons
Initial conditions: Taylor state 
g0 and 0 chosen to give a strong “toroidal” field

*Single fluid resisitive MHD*
*Instantaneous heating by electrons*
*Initial conditions: Taylor state** *
*g**0** and ****0** chosen to give a **strong “toroidal” field*

- Application to Pellet Injection

- Verification and CPU Timings

- Good agreement between explicit and implicit methods

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

**1**

**8**

**64**

**256**

## Conclusion 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
- Fully implicit Newton-Krylov, with no preconditioning, applied to a 3D model pellet injection problem
- Future work
*Implicit NK method for mapped grids*
*Physics-based preconditoners *
*Long term plan: combine implicit with adaptive mesh refinement*