PIES

1. Code Name: PIES (Princeton Iterative Equilibrium Solver)
2. Code Category: Equilibrium
3. Primary Developer: D. Monticello
4. Other Developers and Users: A. Reiman, S. Lazerson, M. Drevlak (Germany). Other occasional users at CIEMAT (Spain) and Kyoto University (Japan). Has also had several additional users within the lab.
5. Short description (one line if possible): Calculates free-boundary 3D equilibria with islands and stochastic regions.
6. Computer Language (Fortran77, Fortran90, C, C++, etc) and approx # of lines: Fortran 90. 112,000 lines, but this includes multiple alternative algorithms for virtually every section of the code.
7. Type of input required (including files and/or output from other codes). Is there any special input preparation system (eg, GUI): Namelist input file for choosing algorithms and controlling parameters. Usually initialized with a VMEC solution.
8. Type of output produced (including files that are read by other codes and sizes of large files and synthetic diagnostics): NETCDF, gmeta plot file, ASCII file
9. Describe any postprocessors which read the output files: Scripts used for plotting.
10. Status and location of code input/output documentation: Extensive documentation in code.
11. Code web site? No.
12. Is code under version control? What system? Is automated regression testing performed? CVS. No regression testing.
13. One to two paragraph description of equations solved and functionality including what discretizations are used in space and time:
The PIES code solves the equation curl(B) = j[B], where j[B] is a complicated, nonlinear function of B determined by j_\perp = B \times \nabla p / B^2, B \cdot \nabla (j_\|/B) = -\nabla \cdot j_\perp. The equation for j_\perp follows from the MHD force-balance equation, while that for j_\| follows from \nabla \cdot j = 0. This formulation has an advantage that it is not necessary to evaluate curl(B) to calculate j_\perp. (Near a rational surface, it can be the case that |j_\|| >> |j_\perp|, requiring a cancelation between two terms large compared to j_\perp in the evaluation of curl(B). There can also be localized currents near the rational surfaces, further complicating the evaluation of curl(B).) Note that force-balance parallel to the magnetic field does not enter into these equations, allowing for the solution of the equilibrium in stochastic regions having finite pressure gradient. For this purpose, the magnetic differential equation along the chaotic field line trajectories which determines j_\| is solved using methods of statistical turbulence theory, more specifically resonance broadening theory. PIES uses finite differencing in the radial coordinate and Fourier harmonics in the poloidal and toroidal directions. For the free boundary calculation, the code utilizes a user-specified "reference surface" that serves as the boundary of the PIES grid, with the Green's function solution used outside the reference surface, and all components of the field continuous at the reference surface. Realistic coils can be specified, with the boundary condition that the field goes to zero at infinity. Three different algorithms are available to solve the nonlinear equilibrium equation: 1) Picard iteration with underrelaxation; 2) A Jacobian-free Newton-Krylov scheme with adaptive preconditioning, globalized by subspace restricted backtracking; 3) a full Newton algorithm, which parallelizes more readily than JFNK, and can be preferable in some cases when a large number of CPUs are available.
14. What modes of operation of code are there (eg: linear, nonlinear, reduced models, etc): Code can run in mode where it flattens pressure and current profiles in islands and stochastic regions, or it can specify current profiles in island interiors (for simulating neoclassical tearing modes), or it can keep profiles fixed if these are determined experimentally (for analyzing experimental data).
15. Journal references describing code:

"Calculation of Three-Dimensional MHD Equilibria With Islands and Stochastic Regions", A. Reiman and H. Greenside, Compt. Phys. Commun. 43, 157-167 (1986).
"Numerical Solution of Three-Dimensional Magnetic Differential Equations,'' A. H. Reiman and H. S. Greenside, J. Comput. Phys. 75, 423-443 (April 1988).
"Convergence Properties of a Nonvariational 3D MHD Equilibrium Code,'' H. S. Greenside, A. H. Reiman, and A. Salas, J. Comput. Phys. 81, 102 (1989).
"Three-Dimensional Plasma Equilibrium Near a Separatrix,'' A. H. Reiman, N. Pomphrey, and A. H. Boozer, Phys. Fluids B 1, 555-562 (March 1989).
"Tools for Three-Dimensional MHD Equilibrium Configuration Studies,'' J. L. Johnson, D. Monticello, N. Pomphrey, and A. H. Reiman, Proceedings of the Thirteenth Conference on the Numerical Simulation of Plasmas, (Santa Fe, NM, September 1989).
"Computation of Zero Beta Three-Dimensional Equilibria with Magnetic Islands,'' A. H. Reiman and H. S. Greenside, J. Comput. Phys. 87, 349-365 (1990).
"Computation of Magnetic Coordinates and Action-Angle Variables,'' A. H. Reiman and N. Pomphrey, J. Comput. Phys. 94, 225-249 (May 1991).
"The PIES Code: Three-Dimensional Equilibria with Pressure, Islands and Stochastic Regions,'' D. A. Monticello, A. H. Reiman, and C. Y. Wang, Proceedings of the Ninth International Workshop on Stellarators, (International Atomic Energy Agency, Vienna, Austria, 1993), pp. 114-120.
"Comparison of ATF and TJ-II Stellarator Equilibria as Computed by the 3D VMEC and PIES Codes,'' J. L. Johnson, D. A. Monticello, A. H. Reiman, A. Salas, A. L. Fraguas and S. P. Hirshman, Comp. Phys. Commun. 77, 1-10 (September 1993).
"Stellarator Equilibrium Studies with the PIES Code,'' P. Merkel, J. L. Johnson, D. A. Monticello, A. H. Reiman, A. Salas, and A. L. Fraguas, Proceedings of the Fifteenth International Conference on Plasma Physics and Controlled Nuclear Fusion Research, (International Atomic Energy Agency, Vienna, Austria, September 1994).
"Recent Developments and Applications of the PIES Code,'' D. A. Monticello, J. L. Johnson, A. Reiman, and P. Merkel, Proceedings of the Tenth International Workshop on Stellarators}, EUR-CIEMAT 30 (1995) p. 45-48.
"The PIES algorithm and its application to stellarators, heliacs and tokamaks,'' D. A. Monticello, A. H. Reiman, J. L. Johnson, P. Merkel, and A. Salas, in Proceedings of the Third Australian-Japan Workshop on Plasma Theory and Computation Robertson, New South Wales, Australia, November 1995.
"The PIES 3D MHD Equilibrium Code: Theory and Applications,'' J. L. Johnson, D. A. Monticello, A. H. Reiman, and S. Deshpande, Proceedings of the Fourth Symposium on Plasma Dynamics: Theory and Applications, (Trieste, Italy, Nov. 1995).
"Finite-beta Equilibria for Wendelstein 7-X Configurations Using the Princeton Iterative Equilibrium Solver Code,'' S. Arndt, P. Merkel, D. A. Monticello and A. H. Reiman, Phys. Plasmas 6, 1246-1252 (1999).
"Reduction of Islands in Full-Pressure Stellarator Equilibria,'' S. Hudson, D. Monticello, A. Reiman, Phys. Plasmas 8, 3377 (2001).
"Eliminating Islands in High-Pressure Free-Boundary Stellarator Equilibrium Solutions," S.R. Hudson, D. A. Monticello, A.H. Reiman, A.H. Boozer, D.J. Strickler, S.P. Hiershman, and M.C. Zarnstorff, Physical Review Letters, Dec. 30, 2002.
"Constructing integrable full-pressure full-current free-boundary stellarator magnetohydrodynamic equilibrium solutions," S.R.Hudson, D.A.Monticello, A.H.Reiman, D.J.Strickler, S.P.Hirshman, L-P. Ku, E.Lazarus, A.Brooks, M.C.Zarnstorff, A.H.Boozer, G-Y. Fu and G.H.Neilson, Nuclear Fusion 43(10):1040,2003.
"Solving the 3D MHD Equilibrium Equations in Toroidal Geometry by Newton's Method," H. J. Oliver, A. H. Reiman and D. A. Monticello, Journal of Computational Physics, Volume 211, Issue 1, 1 January 2006, Pages 99-128. PIES Free Boundary Stellarator Equilibria with Improved Initial Conditions, Michael Drevlak, D. Monticello, A. Reiman, Nucl. Fusion 45 No 7 (July 2005) 731-740.
"Localized Breaking of Flux Surfaces and the Equilibrium beta Limit in the W7AS Stellarator," A. Reiman, M.C. Zarnstorff, D. Monticello, A. Weller, J. Geiger, and the W7-AS Team, Proceedings of the 21st IAEA Fusion Energy Conference, Chengdu, China, 2008.
"Equilibrium and Flux Surface Issues in the Design of the National Compact Stellarator Experiment," A. Reiman, S. Hirshman, S. Hudson, D. Monticello, P. Rutherford, A. Boozer, A. Brooks, R. Hatcher, L. Ku, E. A. Lazarus, H. Neilson, D. Strickler, R. White, and M. Zarnstorff, Fusion Science and Technology 51, pp. 145-165 (2007).
"Pressure-induced breaking of equilibrium flux surfaces in the W7AS stellarator," A. Reiman, M.C. Zarnstorff, D. Monticello, A. Weller, J. Geiger and the W7-AS Team, Nucl. Fusion 47, 572-578 (2007).
"Plasma Equilibrium in a Magnetic Field with Stochastic Regions," J. A. Krommes and A. H. Reiman, Phys. Plasmas 16, 072308 (2009).
"Efficient Numerical Calculation of MHD Equilibria with Magnetic Islands, with Particular Application to Saturated Neoclassical Tearing Modes," D. Raburn, Ph.D. thesis, April, 2011.
16. Codes it is similar to and differences (public version):
VMEC code (ORNL): Calculates 3D equilibria using representation for magnetic field which assumes nested flux surfaces, so that it cannot handle islands and stochastic regions.
HINT code (NIFS, Japan): 1) PIES is true free-boundary code; 2) PIES allows direct specification of current and pressure profiles, including in islands (needed e.g. for neoclassical tearing modes (NTMs)); 3) PIES has physical model for handling nonzero pressure gradients in stochastic regions; 4) PIES does not require differentiation of B to calculate localized currents with |j_\|| >> |j_\perp| near rational surfaces; 5) Very different algorithms make benchmarking between codes valuable for understanding limitations of each code.
SIESTA (ORNL): Fixed boundary version is under development. Does not yet have a refereed publication describing the algorithm.
17. Results of code verification and convergence studies (with references):
The PIES code has had extensive verification and convergence studies since it was first developed in the mid 1980s. As a first step in the verification process, each piece of code was tested as it was built. One method used for this was the method of manufactured solutions. For example, to test the Ampere solver, we calculated curl(B) for several different fields, B, to get a corresponding j. For each such j, we verified that the Ampere solver recovered the corresponding B.
Initial convergence studies: ``Convergence Properties of a Nonvariational 3D MHD Equilibrium Code,'' H. S. Greenside, A. H. Reiman, and A. Salas, J. Comput. Phys. 81, 102 (1989).
Since the initial implementation of the code, we have developed alternative algorithms for almost every piece of the code. Calculations with the new pieces of the code have been benchmarked against the previous versions. For example, we now have two different ways to calculate magnetic coordinates, one using a Newton algorithm, and the other requiring field-line following. As another example, we have developed an Ampere's Law solver that works in terms of the vector potential, while the original version solves for the field directly. To verify the overall code we directly check whether the equilibrium equations are satisfied. With our formulation of the equilibrium equations, the residual takes the form B - (curl)^{-1} j(B) , where (curl)^{-1} denotes the inverse of the curl operator. We monitor this residual as the code converges. For additional verification, we also check that the solution obeys the equilibrium equation in the more conventional form (curl(B)) \times B - \nabla p = 0 to within discretization error. We verify that this goes to zero as the resolution is increased. We have also compared with analytic test cases, and these comparisons have been published or have been presented at meetings:
  • Soloveev equilibria as well as the large aspect ratio stellarator expansion: (Greenside, Reiman, and Salas, J. Comput. Phys. 81 (1989) 102);
  • Helical force-free Bessel function equilibria with islands (Reiman and Greenside, J. Comput. Phys. 87 (1990) 349);
  • The narrow-island, saturated tearing mode solutions of White et al and Zakharov (Pomphrey and Reiman, Sherwood paper 3D29 (1989))
We have benchmarked the solutions of the PIES code against those from other codes, including cases that have been validated against experimental data:
  • 2D j-solver equilibria for TFTR and DIII-D;
  • Biot-Savart vacuum field solvers;
  • The linearized resistive time-dependent code developed by Mike Hughes (comparisons with marginal stability for tearing modes on TFTR);
  • A Helical Grad-Shafranov solver that can handle magnetic islands, including neoclassical effects on island widths, developed by our graduate student, Dan Raburn;
  • VMEC. ("Comparison of ATF and TJ-II Stellarator Equilibria as Computed by the 3D VMEC and PIES Codes", Johnson, Monticello, Reiman, Salas, Fraguas and Hirshman, Comp. Phys. Commun. 77, 1 ( 1993)).
The Johnson et al paper contains benchmarking between convergence studies with PIES and VMEC. One such set of studies is shown in the figure. We now routinely use VMEC solutions for an initial guess for PIES code calculations. Each such calculation provides a comparison of a VMEC solution with a PIES solution. In cases with small islands where the solutions are not approximately equal, we verify that the solutions converge to each other as the resolution is increased. A helical Grad-Shafranov solver that can handle magnetic islands was constructed by a graduate student, D. Raburn, and was benchmarked against the PIES code, including cases with current profiles in the islands to model the effects of neoclassical tearing modes (NTMs). A three-way collaboration has been established with NIFS (Japan) and IPP Greifswald (Germany) to benchmark the PIES and HINT codes.
18. Present and recent applications and validation exercises (with references as available):
Benchmarking against other codes over the last 20 years has included cases that were validated against experimental data.
  • 2D j-solver equilibria for TFTR and DIII-D;
  • The linearized resistive time-dependent code developed by Mike Hughes (comparisons with marginal stability for tearing modes on TFTR);
  • ("Comparison of ATF and TJ-II Stellarator Equilibria as Computed by the 3D VMEC and PIES Codes", Johnson, Monticello, Reiman, Salas, Fraguas and Hirshman, Comp. Phys. Commun. 77, 1 ( 1993)).
PIES was extensively used in the design of the NCSX stellarator, including guiding modifications in the shapes of the modular coils to improve flux surfaces in the calculated equilibria.
"Eliminating Islands in High-Pressure Free-Boundary Stellarator Equilibrium Solutions," S.R. Hudson, D. A. Monticello, A.H. Reiman, A.H. Boozer, D.J. Strickler, S.P. Hiershman, and M.C. Zarnstorff, Physical Review Letters, Dec. 30, 2002.
"Constructing integrable full-pressure full-current free-boundary stellarator magnetohydrodynamic equilibrium solutions," S.R.Hudson, D.A.Monticello, A.H.Reiman, D.J.Strickler, S.P.Hirshman, L-P. Ku, E.Lazarus, A.Brooks, M.C.Zarnstorff, A.H.Boozer, G-Y. Fu and G.H.Neilson, Nuclear Fusion 43(10):1040,2003.
A. Reiman, S. Hirshman, S. Hudson, D. Monticello, P. Rutherford, A. Boozer, A. Brooks, R. Hatcher, L. Ku, E. A. Lazarus, H. Neilson, D. Strickler, R. White, and M. Zarnstorff, Fusion Science and Technology 51, pp. 145-165 (2007).
Equilibrium reconstructions for the W7AS stellarator, including islands and stochastic regions. This work has allowed us to study the issue of stochastization of flux surfaces with increasing beta, and to validate the code against W7AS experimental data. The PIES calculations have been found to be consistent with the available experimental data, including the effect of the divertor control coils on the field line stochasticity, and the observed differences in the pressure profiles in the region calculated to be stochastic.
"Equilibrium and Stability of High-Beta Plasmas in Wendelstein 7-AS," M.C. Zarnstorff, A. Weller, J. Geiger, E. Fredrickson, S. Hudson, J.P. Knauer, A. Reiman, A. Dinklage, G.-Y. Fu, L.P. Ku D. Monticello, A. Werner, the W7-AS Team and NBI-Group, in Fusion Energy 2004 (Proc. 20th Int. Conf., Vilamoura, 2004) (Vienna: IAEA) CD-ROM file (EX/3-4)
"Localized Breaking of Flux Surfaces and the Equilibrium beta Limit in the W7AS Stellarator," A. Reiman, M.C. Zarnstorff, D. Monticello, A. Weller, J. Geiger, and the W7-AS Team, Proceedings of the 21st IAEA Fusion Energy Conference, Chengdu, China, 2008.
"Pressure-induced breaking of equilibrium flux surfaces in the W7AS stellarator," A. Reiman, M.C. Zarnstorff, D. Monticello, A. Weller, J. Geiger and the W7-AS Team, Nucl. Fusion 47, 572-578 (2007).
"Plasma Equilibrium in a Magnetic Field with Stochastic Regions," J. A. Krommes and A. H. Reiman, Phys. Plasmas 16, 072308 (2009).
19. Limitations of code parameter regime (dimensionless parameters accessible): Equilibrium equations are valid when the time evolution of the plasma is sufficiently slow that the ? dv/dt term can be neglected in the momentum equation. The PIES code presently assumes a scalar pressure, and it neglects the effects of flow. Upgrades to the code are planned to include kinetic and flow effects.
20. What third party software is used? (eg. Meshing software, PETSc, ...): At present, the code requires a few NAG routines. These could easily be replaced by equivalent routines in the public domain, but we have not done so because NAG is available on all of the platforms where PIES has been used.
21. Description of scalability: Our parallelization is still at an early stage, and we have not yet looked at the issue of scalability.
22. Major serial and parallel bottlenecks: The major serial bottleneck is generally the construction of magnetic coordinates, for which we generally use field line following. We have also implemented a Newton algorithm for this purpose, which is considerably faster, but the field line following algorithm tends to be more robust. We are presently implementing a more efficient field line following algorithm. Because our parallelization is still at an early stage, we do not yet know what the parallel bottlenecks will be.
23. Are there smaller codes contained in the larger code? Describe: The key components of the code are those which: (1) Calculate magnetic coordinates for a given magnetic field; (2) Solve a magnetic differential equation along the magnetic field lines; (3) Express the solution to Ampere's Law in three dimensions in terms of the solution of a 3D Poisson equation; (4) Solve a 3D Poisson in arbitrary specified coordinates.
24. Supported platforms and portability: PIES has been compiled and run on a number of platforms, including LINUX, SUN, some NERSC platforms, some platforms in Greifswald, some platforms in NIFS. It is usually an easy process to port the code to any platform as PIES is mostly platform independent.
25. Illustrations of time-to-solution on different platforms and for different complexity of physics, if applicable: We are beginning to implement and test our new JFNK scheme for stellarator calculations. For test cases, it has been found to be about an order of magnitude faster than the Picard algorithm that we are currently using. A typical Picard calculation for a W7X stellarator equilibrium takes about 10 hours on a single AMD Opteron 2500 mHz core. We expect the time advantage of the JFNK scheme to be particularly manifest in time-dependent transport calculations, where the changes in the equilibrium from one time step to the next are small. We are beginning to work on parallelizing the code.