orbitnotes

Stored for public access on machine a.nersc.gov in /tmp/wk1/u454/ORBIT

Start by taking all files from ORBIT:

Files for the code orbits.f (spline version):

Begin with the PEST equilibrium file, the example given is pesteq, an ITER equilibrium. For information how to produce this file see Manickam.

First type gmake, which will construct the mapmc files. The mapmc files take as input fort.17.

1. copy pesteq into fort.17, the input file.

2. Type mapmc, which will produce fort.18 and fort.19. This is the PEST mapping code, set to give the equilibrium in Boozer coordinates, jacobian 1/B^2.

3. Rename these files for your own identification, copy fort.18 as m1iter1, and fort.19 as m0iter1 m0 and m1 prefixes are universal, the rest is a name of your choice.

4. Inside spline.f, subroutine defolt, find the mp0 and mp1 names and set them to be m0iter1 and m1iter1. Note that according to the setting of map, several different equilibria can be accessed. Set map at the beginning of spline.f to give the equilibrium you want. You can also produce analytic second order Shafranov equilibria by selecting numeric = 0 at the beginning of spline.f. See subroutine shaf for analytic equilibria input.

5. Set the major radius (magnetic axis!) in centimeters at the beginning of spline.f, and the ripple choice krip. I have ripple models for ITER, Tore Supra, TFTR, and NSTX. Sorry, you will have to construct one for DIIID. If no ripple is desired krip=0 skips it. The magnetic axis location is essential since the code units are defined by the major axis and the gyro radius. The gyro radius is calculated later in orbit from the value of z, energy, B, and proton mass. 5. Type shs. shs is a shell, stored in ORBIT, which will compile and run spline.f, producing the spline data set spdata. It also produces an output file spout and, in the case of an analytic Shafranov equilibrium, a plot file gmeta, giving equilibria characteristics. Copy and store spdata as eq1001, or whatever, because otherwise it will be written over the next time you run spline.f

6. If the outside surface of the equilibrium is not well represented, the spline dimension is too small. At the beginning of spline.f, the parameters lsp and lst are respectively the poloidal flux and poloidal angle spline dimensions. Increase lst. Limits on these dimensions are governed by the common blocks spline, used by both orbits.f and spline.f. Likewise, if the magnetic axis location or plasma position are poor, increase lsp. The accuracy of the representation depends on the plasma shape. The most detailed plots of the equilibrium is provided by orbits.f

7. Orbits.f, (the spline version) takes as input spdata. In the main of orbits.f there are many options, for single particle run, the addition of perturbations, different particle distributions, collisions, drag, subsidary options such as stochastic loss calculations, ripple contours, the trapped passing boundary, etc. The guiding center equations used are given in White, Phys Fluids 2, 845 1990.

8. The shell sho compiles and runs orbits.f. The file kam0 is a universal stochastic threshold data file, see White et. al. Physics of Plasmas 3, 3043, 1996. Output is a plot file, gmeta, a data file orbout, and other data files for lost particle distributions, etc., which can be constructed and set as desired, usually in the routine pdist (particle distribution) or plost (lost particle distribution).

9. There are all kinds of working switches in the code, for example in the alpha particle distribution routine alphdep, ntrap = 1 can be used to produce only a trapped distribution. These switches have not been moved into an input file because I modify them all the time, and produce new ones constantly. I try to put write statements with them, so always check orbout to be sure you are doing what you want to do, i.e. that the switches are set the way you want. A good expedient is to always do a short run with few particles to check what you are doing before starting a long run.

10. Single particle orbits nplot = 1 runs single particle orbits, giving plots of time history. The first plot shows the accuracy of energy conservation. The time step is controlled by the energy conservation, using the limiter dele, normally set to about 5.e-8. If time dependent MHD modes are used (set npert = 1, and set the amplitudes and frequencies in subroutine amp1) energy is not conserved, and a fixed time step must be used. dele = 10 accomplishes this, the time step is dt0, which can be adjusted.

11. Particle loss analysis nplot = 5 is for doing particle loss analysis. There are several subroutines for loading different types of Monte-Carlo distributions. The subroutine pdist will give plots of particle distributions. It can be called at any time during a run, but normally only at the beginning and end. Inside the subroutine eject, near the beginning, there is a write statement (commented out) which writes out individual particle loss data, including the particle number. To observe a particular loss orbit leave the initial nplot =5 load parameters the way the were, and activate runone using this particle number, using the call in the main, just before the time step loop. The routine runone loads the particle you have selected into position 1, and then switches to nplot =1, giving a plot of that particular orbit. Inside plot5 there are several optional calls call sigma - gives a statistical error analysis of particle loss call pdist - plots the particle distribution functions call plost - plots the lost particle distribution functions call mupzeta(1) initial particle positions, pzeta, mu plane call mupzeta(0) final particle positions, pzeta, mu plane call dump0 - writes out a file for lost particle analysis, to be used with lost.f