NIMROD Team Meeting Minutes

August 27-28, 2008

San Diego, California

Model Development

Carl Sovinec started this session reporting on two development projects.  The first adds a numerically implicit hyper-diffusive operator to the continuity equation.  The associated particle flux is proportional to third derivatives of number density, and like hyper-viscosity in fluid turbulence computations, it allows strong smoothing near the node spacing with little effect on scales of interest.  This has been helpful for 3D two-fluid simulations.  The second project is the development of a Grad-Shafranov solver using the NIMROD computational framework—work being done by U-WI grad student Eric Howell.  The solver is working and getting benchmarked on simple equilibria.  Further tests are needed for cases with pressure profiles and a separatrix.  Development work for writing dump files based on the GS solution is underway.

Eric Held gave an overview of Jeong-Young Ji’s work on his generalized velocity-moment method as a kinetic description of plasmas.  It is based on complete orthogonal polynomials and the system of equations is derived with the exact linear Coulomb collision operator.  Ji is publishing separate manuscripts on derivations and results for heat flux and stress.  One significant physical effect from the results is how the ion heat flux depends on collisions with electrons when Te<Ti.  Details of a system truncated at 21 moments were presented, and needs for implementation in the NIMROD framework were described.

Held also presented recent work and new ideas related to drift kinetic equation (DKE) calculations.  He and his group are trying to make better use of data from integral computations (to deposit accumulated information at more locations than just the launch point for a characteristic), but new work for solution of the DKE in differential form shows promise.  The differential DKE computations are motivated by success with the mixed finite element method (MFEM), where the parallel component of the diffusive heat flux density is solved simultaneously with temperature.  This results in impressive efficiency improvements over standard high-order FE with large conductivity ratios.  The differential DKE computations solve for coefficients of a polynomial expansion of the kinetic distortion simultaneously with T, with some simplifications wrt to collisions and speed dependence for feasibility testing.  Held showed the results of heat flux over an island in the collisional limit, and they already show reasonable agreement with the fluid closure.  He will be considering temporal centering and how to incorporate a finite differencing grid or FE approximation for the speed independent variable.

Charlson Kim reported on recent progress and plans for the energetic particle model.  The full-kinetic model is an option, along with drift kinetics, and can use orbit averaging.  More work is going into the use of high-order fields for the particle push and using high-order shape functions for the particle deposition.  The particle shapes are defined with respect to physical position and deposited onto all nearby quadrature points where the shape is nonzero.  This will affect domain decomposition for parallel computation—more needs to be done, and the decomposition for particles will likely need to be separate from that of the fluid algorithm.  Kim has also modified the lagr_1D routine to use a table for the coefficients, which improves performance.  Schnack and Kruger expressed interested in the updates for application to the giant sawtooth problem.

Scott Kruger reported on code coupling for the slow-MHD part of the SWIM project and on diagnostics created by Dylan Brennan for tearing-mode analysis.  Kruger reviewed the status of the coupling project at the time of their review earlier this summer.  [Phase 0: symm. ad hoc drive—done, Phase 1: asymm. ad hoc drive—in progress, Phase 2: pass NIMROD fields to GENRAY—in progress, Phase 3: use quasi-linear diffusivity in FRF computation—formulated.]  The coupling is going through EQDSK files, and Kruger has modified NIMROD to write EFIT-format files (called ‘fake’ because they may not represent an accurate GS solution).  Here, the challenging part is going through different spatial meshes.  This part is working.  Tom Jenkins is working to understand GENRAY so that the quasi-linear computation can be returned to NIMROD.  This is also nontrivial, because the data is generated along rays, which deposit unevenly across the NIMROD mesh, an issue similar to making better use of the integral closure information, which is mention above.  Through numerical trials, it has been found that the standard QSHEP2 routine is inadequate for this application.  The tearing diagnostic determines the poloidal Fourier components of NIMROD data using a straight field-line angle and flux-coordinate components.  It is now working reliably and is showing phase differences between different rational surfaces, in addition to the amount of reconnected field.  Kruger also showed equilibria where packing nodes during initialization helps field-line quality.

Val Izzo presented development work that uses NIMFL to assess the acceleration of superthermal electrons into runaways that could be damaging to a device with the stored magnetic energy of ITER.  She uses existing NIMROD solutions to determine E and B for each of a set of particles.  The present computations use drift kinetic orbits and a drag term, but they can be generalized to include FLR effects.  The modifications are mostly in int_segment; though, it was suggested that the integration can be done by modifying the derivative-function routine that is called by LSODE.  Other planned development includes drag from synchrotron radiation and to automate the computations over multiple dump files.

Performance Development and Scaling

Ping Zhu presented information on porting NIMROD to the Cell Broadband Engine processors that are used in IBM blade systems, such as Roadrunner, and in PS3 computers.  The environment is inhomogeneous with the power processor (PPE) being similar to other IBM power chips but the SPEs being accelerators that have no memory.  To take advantage of the SPEs, one has to call libraries, such as a special CBE BLAS library, that use them.  In an effort to test NIMROD, Zhu has focused on the matrix-vector multiplication routine, which is used intensively if SuperLU is not called for preconditioning.  So far, there is no systematic improvement in speed, and this may be due to the fact that the SPEs have been developed for single-precision computation, while the NIMROD computations use double-precision.

Sovinec presented recent work on improving preconditioning for the two-fluid magnetic field advance.  He reviewed the fact that the Fourier representation for the toroidal coordinate provides a natural physics basis for preconditioning magnetic-confinement computations; perturbations are small, so the diagonal-in-n block is naturally dominant.  This has been the basis for preconditioning all of NIMROD’s 3D solves.  The use of a preconditioner based on planes at different toroidal angles was expected to be complementary, but it suffered from lack of diagonal dominance and was not effective even in combination with what is normally done.  Using selected off-diagonal couplings has proven to be more effective; tests show that the iteration count does not grow with increasing toroidal resolution.  The importance of using some electron inertia was also discussed.  Finally, a new weak scaling study on Franklin shows that the new approach works well with SuperLU_DIST for the diagonal-in-n inversions up to more than 1000 processors.  Some data and looping rearrangement is expected to further improve the Fourier scaling part.

Held reviewed processor scaling for the integral closures computations.  He and Kruger had developed code to separate closure-calculation processors from fluid-calculation processors, and this approach scales reasonably well to a few thousand processors.  Synchronization between the separate groups is the primary issue for further scaling.  Running an entire parameter scan simultaneously is another approach to using large numbers of processors, and this was used for heat flux computations that were presented at the RMP workshop.

Kim discussed scaling issues for the energetic particle model.  The present bottleneck is load balancing with NIMROD’s standard domain decomposition being based on grid blocks and Fourier layers.  The particles can be decomposed separately but may need duplication of the field data that is used for the particle push.  An action item from this discussion is to do a literature search on different domain decomposition strategies that are used in PIC computations to bring in as many ideas as possible before further development.


Bonita Squires presented computations of peeling/ballooning based on the toroidal circular cross section TOQ equilibrium developed by Scott Kruger and Phil Snyder for benchmarking.  For this case, ideal plasma behavior is reproduced for S³5ę107.  Results also asymptote when S for the edge region is below unity.  [There was discussion concerning length scales for the edge.]  To see a growth-rate maxima at low n (associated with peeling), the resistivity gradient needs to be very narrow and fit between the pressure pedestal and the rational surface of the externally resonant mode.  Squires also described cases that show unusual linear numerical behavior.  For example, when the resistivity gradient is broad and has some extent across the pressure gradient, growth rates computed from magnetic and kinetic energies do not remain at the same value.  A nonlinear demonstration computation that uses a fixed resistivity profile was also presented.  An action item related to this work is to send the TOQ equilibrium and specification to Steve Jardin, who is doing similar testing with M3D-C1.

Zhu presented ballooning computations that go into the late nonlinear stage (defined by perturbation amplitude scaling with respect to the k-parallel/k-perp ratio).  They are helped by the new particle hyper-diffusivity and are being used in comparison with recent analytical work.  The equilibrium has circular cross section and is generated by the ESC code.  There is no dissipation except for the particle hyper-diffusivity, and the nonlinear computations go a factor of 10 further in amplitude than they did with the Fick’s law particle flux.  Development to track displacement is underway.

Chris Hegna presented work from graduate student Mark Schlutt that applies NIMROD to 3D cylindrical stellarator equilibria.  The motivation is the existence of experimental results beyond the conditions where interchange-type and low-n instability is expected.  The computations start with helically asymmetric magnetic perturbations that are frozen in the wall.  Computational results show that when perturbed with shearing motions, the 3D equilibrium remains intact.  A source of heat is then added with anisotropic thermal conduction to develop pressure gradients.  A Chapman-Enskog-like closure will eventually be added to the study.

Hegna also described a study that graduate student Andrea Montgomery is starting.  The general topic is the resistive wall mode (RWM) and field-error penetration in a cylinder.  It is motivated by a ’95 publication from J. Finn that considered how resonant modes can provide the damping that is needed for rotational stabilization.

Kruger presented tearing-mode applications by Brennan (with some help from Kruger).  He reviewed the different categories of neoclassical tearing modes (NTMs) outlined in a paper by Brennan.  To assess different simulations, the DĘ-matrix theory is being applied.  The importance of rotation profiles is also being considered, and the new diagnostics described earlier are being used to examine relative phases of perturbations at different resonant surfaces.  In the case of an n=2 mode, it appears as a secondary instability (being made linearly unstable during the course of a nonlinear computation) but is not driven directly by the primary mode.  There is no evidence of secondary islands at this point, but the response when viscosity is changed is not understood.

Sovinec presented two-fluid reconnection and tearing work by graduate student Jacob King and himself.  He reviewed some history of two-fluid 1/1 internal kink computations.  Recent NIMROD simulation results use the new preconditioning approach and have been spatially converged with S=106, Pm=0.1, di=0.11, and de£0.01.  A comparison with a similar MHD computation clearly shows X-point reconnection geometry in the two-fluid case, but the accelerating growth-rate first reported by Aydemir has not been obvious.  The NIMROD computations are at poloidal-beta <1, which may be related.  Results on a tearing mode in a non-reversed pinch are a start on multihelicity tearing with two-fluid modeling.  The new preconditioner and physically realistic electron inertia keep the time-step from crashing through saturation for a case with moderate rs (3% of minor radius), but a computation with larger rs still runs into numerical difficulties.  The cases run so far do not produce significant multihelicity coupling, as expected in reversed profiles.

Izzo presented her recent study of runaway electrons in DIII-D and comparing with conditions for ITER.  The runaway electrons form when electric field exceeds the Rosenbluth Ecrit, and relative to this criterion, ITER will have safer conditions than DIII-D.  Massive gas injection is too slow to penetrate, so the new technique is dilution cooling with D2.  Izzo starts nonlinear computations after the cooling has penetrated, and MHD activity is quickly triggered.  The simulations reproduce the key issue of thermal quench occurring prior to current quench (hence large E with T-dependent resistivity).  Unlike DIII-D, the dilution cooling appears to prevent runaways in ITER conditions.

Dalton Schnack presented Tom Jenkins’ recent work for RF stabilization of slow-MHD; Jenkins participated remotely.  Schnack reviewed the cartoon model for asymmetric ad hoc current drive, and Jenkins has worked out the analytics for the dissipation of shear Alfven waves excited by a local source.  The symmetric n=0 mode is the only one that is not damped, which seems reasonable, but there was debate on the spatial scale that would appear in the decay rates.  When using the ad hoc drive in toroidal tearing computations, cases with varying integrated drive suggest that a small spot size (less than island width) for the drive is not very effective.  However in a series with fixed ad hoc-induced current, the smaller spot size is actually more effective in reducing gamma (measured by resetting and running the linear phase of the mode with the ad hoc current-modified profile).  Heating associated with ECRH can also be modeled.

Charlson Kim reviewed computations being performed at PSI-Center.  The center is returning to the nimdevel version of the code.  Their NIMROD applications include FLR in RFP tearing, HIT-SI, LDX, Swarthmore SSX modifications, FRC work for thrusters, and ZAP.  The LDX computations are using a heat source to generate equilibria from the vacuum fields.  They have progressed to getting a global interchange linearly unstable.

Sovinec presented graduate student Chris Carey’s MHD study related to astrophysical jets.  The 3D simulations retain symmetry in regions where the current profile is very collimated, and linear initial-value computations indicate that viscosity and axial flow are not the stabilizing effect.  The azimuthal rotation is key.  Carey has done eigenvalue analysis in both Lagrangian and Eulerian references frames and found that rigid rotation stabilizes pinches where either scalar pressure or magnetic pressure provide force-balance with centrifugal acceleration.  With the perturbed mass term added to NIMROD’s linear advective force, the initial-value computations agree quatitatively.  [This change makes the linear behavior consistent with ‘full’ continuity in nonlinear computations instead of ‘fix profile’ continuity.]  A paper describing this linear behavior is in preparation.

Sovinec provided an overview on helicity injection computations at Wisconsin.  He reviewed results from Adam Bayliss on simulation of HIT-II in weakly relaxing conditions and direct comparison with experimental results on plasma current and injector current as a function of toroidal field strength.  There is good quantitative agreement, but experiment and computation deviate somewhat and in different directions from a simple equilibrium estimate.  For flux compression on Pegasus, graduate student John O’Bryan has developed the thermal conduction model to include demagnetization in Braginskii’s fluid model.  In this application, the low-temperature plasma outside the current channel is very collisional, so the modeling reverts to what is usually parallel conduction, which is smaller than perpendicular in these conditions.  Undergraduate student Tom Bird is using a toroidally localized ad hoc current drive to induce a helix of current, as in Pegasus miniature gun injection.  He will be investigating relaxation as more current is induced; however, greater resolution is needed than in the current helix result that was shown.

Schnack presented Bayliss’ recent work to simulate Cary Forest’s new plasma dynamo experiment, which has rings of permanent magnets and electrodes that spin plasma from the surface.  A small-scale experiment uses cylindrical geometry, and this series of computations model that experiment.  The effect of the electrodes is modeled with a boundary condition on tangential electric field.  While the dynamics near the surface are complicated, the flow transitions to rigid rotor at the interior of the midplane.  The flows excite an n=1 perturbation in 3D computations.

Project Improvement

Related to the website, the following action items were discussed:

Š      Provide updates to the publication list to either Sovinec or Kruger.

Š      Modify the opening page so that the links go to text on the same page.

Š      The team member list needs updating.

Š      Send pdfs of this meeting’s presentations to Sovinec.

Related to discussions with Steve Jardin:

Š      Send him the peeling/ballooning benchmark equilibrium.

Š      ITER wants computational information on structural forces during disruption.

Š      We should consider modeling the Zakharov wall mode.

Rostom Dagazian told us that OFES is looking for a new computational milestone for FY 2010.  We discussed possibilities related to sawteeth and RMP.

The final segment included a discussion of using workflow software to automate jobs and manage data.  There are different opinions on how useful this would be for NIMROD, but it may provide some political benefits.  The Doxygen software can be useful for documentation.  It automatically extracts comments from code and puts them in linked html pages.  A minor modification in the writing of nimrod.out will allow it to be used for generating a file.  Finally, there was a request to establish a new release version of nimdevel.