Response to the FSP questionnaire on Energetic Particle physics component based on NOVA-K code

Generic questions

Physics focus and programmatic scope of the component candidate (CC)

1. Give a brief, high-level description of CC's functionality
   
   NOVA-K suite of codes are linear hybrid MHD/kinetic codes for EP driven ideal and kinetic MHD eigenmode instabilities. It is able to predict growth and damping rates. It also implements the theoretical model for AE saturation amplitude. It is being worked on to include the quasilinear diffusion model for fast ions in the presence of multiple instabilities.
   
   
2. Describe CC's user base and application scope (who uses CC and for what).
   
   NOVA-K codes are widely used for AE (Alfven eigenmode) structure and stability calculations. Among its users are GA, MIT, JAERI, Zhejing University, China. Mostly used for AE structure and comparison with the observed modes. Recently it is being used more often for predictions of the future experiments stability against various energetic particle driven modes. Several publications exist especially for ITER basing in part on the simulations with NOVA [ITER publications - see item 6].
   

Physical and mathematical models

1. What are the equations solved in CC?
   
   NOVA solves ideal MHD equations and finds eigenmodes, such as TAEs. NOVA-K evaluates fixed mode TAE kinetic growth rates by evaluating the quadratic form with the perturbed distribution function coming from the drift kinetic equation.
   
   
2. What are the limitations of CC imposed by orderings or by neglected terms?
   
   Main limitations are caused by neglecting thermal ion FLR, toroidal rotation, and drift effects in the eigenmode equations. Thus NOVA can not reproduce some important modes, such as kTAE, kinetic RSAE modes. It can not well describe some of the  dampings, such as radiative damping.
   
   
3. What feedback loops are included in the model?
   
   None yet.
   

Numerical approaches

1. What are the discretization approaches for time and space?
   
   Finite element methods are used in radial direction and Fourier harmonics are used in poloidal and toroidal directions. Emloyment of the theoretical model for the mode saturation assumes that the system is near the threshold all the time.
   
   
2. What are the linear and nonlinear solvers involved?
   
   The choice of linear solvers is available within NOVA, such as sparse, SUPER-LU, and old homemade solver from UT at Austin 1976 year.
   
   
3. Describe algorithmic scalability.
   
   NOVA-K code can potentially scale well with the choice of super-LU linear solver.
   
   
4. List performance-engineering tools used to enhance performance.
   
   None.
   

Software engineering issues

1. Give a complete list of CC's inputs, e.g. the set of input parameters, the range of valid values for each, and their dependence on each other
   
   Plasma equilibrium is computed separately and is used in various parts of the simulations. To compute mode structures additional information is provided for the plasma density. Finally for the stability run of the TAE modes a detailed set of plasma parameters is given. Several options are available. First, the data is read from TRANSP output of a certain shot run and used, which is the most convenient way because it has all the required information. Second way to interface with the code is via the configuration file, which contains all the parameters required for the code. Example (default) case with the set of plasma parameters follows.
   
      Toroidal mode number                                                   
    ntor
      3
      Key=1=>TRANSP.dat,=0 =>anal.prfs,=2=>anal+ T(r)<=p(r)/n(r)             
    itransp
      0
      Major rad. of geom. center [cm] if(itransp|=1)                         
    rmaj
    0.262E+03
      Minor rad. of last surface [cm]                                        
    amin
    0.950E+02
      Vacuum mag. field at geom. center [Gauss]                              
    B0
    0.445E+05
      Electron density [cm^-3] at magnetic axis                              
    dn_e
    0.500E+14
      Electron temperature [eV] at magnetic axis                             
    T_e
    0.683E+04
      ICTs central density [cm^-3] for D,T,H,C                               
    dn_i
    0.500E+14 0.500E+08 0.500E+08 0.500E+08
      ICTs central temperature [eV] for D,T,H,C                              
    T_i
    0.143E+04 0.143E+04 0.143E+04 0.143E+04
      Central fast particle beta: D,T,alphas                                 
    betah0
    0.800E-01 0.200E-09 0.300E-09
      Fast particle energy [eV]: D,T,alphas                                  
    energyh
    0.100E+07 0.800E+05 0.352E+07
      1st parameter for plasma density                                       
    alphar
    0.000E+00
      2nd parameter for plasma density                                       
    prho
    0.144E+01
      3rd parameter for plasma density                                       
    arho
   -0.164E+01
      Fast particle mass                                                     
    rmhp
    0.200E+01 0.300E+01 0.400E+01
      Fast particle charge                                                   
    zh
    0.100E+01 0.100E+01 0.200E+01
      Fast particle index 1 to 3                                             
    ihsps
      1
      alphart,prhot,arhot:Tthermal~(1-alphart*rsq**prhot)**arhot             
    alphart
    0.000E+00 0.200E+01 0.100E+01
      alpharh,prhoh,arhoh:bet_h~exp(-|(rsq^1/2-alpharh)/arhoh|^prhoh)        
    alpharh
    0.000E+00 0.200E+01 0.100E+01
      3D distr. funct. params                                                
    chi0
    0.500E+00 0.000E+00 0.400E+00 0.200E+00 0.130E+00 0.200E+00
      Finite Orbit Width: 0< xfow <1                                         
    xfow
    0.100E+01
      Key to m1 stability calculations im1                                   
    im1
      0
      Distr.function type s-slow.down,m -mxwll.(see taem.f)                  
    _dtype
   s
   
   
2. Give a compete list of CC's outputs.
   
   On the output ideal part of the code produces mode structures depending on the options used in the run. It can output several modes at one run within the given frequency range. This is ideal MHD part of NOVA.
   Kinetic NOVA extension, called NOVA-K, produces various outputs:
   
   - growth rates due to fast ions
   - damping rates due to various mechanisms: Landau, collisional, radiative, continuum
   - saturated mode amplitude for a single TAE  mode based on quasilinear theory
   - various particle drift motion frequencies (optional)
   - partial contributions to the growth rates from different particles (optional)
   - resonant curves in the velocity space with the width of the resonance layer (optional)
   
   
3. Give a list of CC's software dependencies.
   
   ncar graphics, nag, pgplot, superlu
   
   
4. Give a list of smaller components contained in CC; for example, CC1 and CC2.
   
   - equilibrium q-solver
   - equilibrium mapper for different coordinate choices
   - package for matrix preparation, novain
   - package to solve ideal MHD equation, novast
   - stability computations, nova-k
   - nonperturtabive stability, nova-kn (partially ready)
   
   
5. List supported platforms and describe portability
   
   Linux, unix systems. Easily portable to various architectures: PC, suns, cray.
   

Verification

1. Give a list of verification and validation tests performed (with preference to community-wide benchmark tests); highlight disagreements to identify problem areas.
   
   A list of publications exists for the V&V using NOVA-K code. We have assembled the website on this, which is a good source of the verification information. It URL is https://w3.pppl.gov/~ngorelen/NOVA_ref.html. Recently a set of benchmark cases was developed with the goal to verify various codes within ITPA group. These cases are documented here https://w3.pppl.gov/~ngorelen/NOVA_vv.html
   
   
2. Give the appropriate/valid subsets of the equations/models/parameters that can be used in an independent way.
3. Illustrate convergence to analytic or asymptotic solutions in special cases.
   
   Various reduced models for shear Alfven wave can be used to recover simple dispersion relation of TAE modes, which are typically computed by the NOVA code. Another model to compare with is fishbone instabilities driven by EPs. Theory reduced model agrees well with the code results [C. Z. Cheng, "Kinetic extensions of magnetohydrodynamics for axisymmetric toroidal plasmas", Phys. Reports, v.211, p.1, 1992]
   
   
4. Rate of convergence studies to show the numerical methods are behaving as expected.
   
   Excellent convergence in the mode structure, frequency illustrate that the level of the understanding is high. The stability analysis need improvement. Althouth main driving and damping mechanisms are understood and implemented in the code, many aspects of the EP physics require detailed studies. In particular the radiative and continuum dampings need are not converged well in the code with the theoretical expectations.
   
   
5. Can CC be instrumented to provide RHS source terms? This is to facilitate the use of the Method of Manufactured solutions to demonstrate convergence for a sufficiently rich test problem to showcase the physics of interest.
   
   I believe so and we are working on the incorporating the quasilinear theory into the code so that the result of NOVAK run will be the relaxed fast ion distribution function.
   
   
6. For what problems and in which parameter regimes does CC perform well (badly), in terms of numerical stability and valid physics result production?
7. Can CC be applied to ITER physics studies, and with what limitations?

Performance

1. Document processor scaling of time-to-solution on topical verification or other physically-relevant problems.
   
   The code is fast. It takes several minutes to computed one eigenmode, but to do the full stability analysis for a given plasma equilibrium it may take up to 1 day of user time. Depending on the case it may also take up to a few days of the computer time for the complete stability including the stable modes.
   
   
2. Which machines in the US is CC mostly running on? How many processors (cores) does it typically use for production runs?
   
   It runs on a local clusters because of the convenience with 1 core per run.
   
   
3. Describe performance variation with complexity of physics.
   
   If only the mode structures are required, NOVA can be run fast. Stability analysis requires more time. Typically order of magnitude longer.
   
   
4. List the major serial and parallel bottlenecks (e.g., I/O, message-passing).

Developmental issues

1. What are your plans, and what problems would you like to solve with more development? How would this change the equations, discretization, or numerical methods that you use?
   
   More kinetic physics are planned to be incorporated into the code. The plans are to include the toroidal rotation, Hall physics, two fluid extension, and drift effects. This will require carefull consideration of the mapping of the velocity space of fast ions to the real space with the new technique for the descretization of the phase space. Numerical method will remain the same.
   
   
2. What tools do you wish you had available to you in your code development processes?
   
   Better graphics. Better debugger.
   

Topical-area-specific questions

Energetic Particles

1. Are there codes, theories and/or reduced fast ion transport models, which can be integrated into the plasma codes, such as TRANSP (or FSP0).
2. What physics is required to develop such models?
3. What boundary conditions are required for successful modeling of fast ion driven instabilities?
4. Are the wall heat loads and blistering effects from escaping high energy alphas tolerable?
5. Does alpha stabilization of core MHD (e.g., sawteeth) allow these modes to reach high amplitudes, followed by strong relaxation oscillations?