LAST REVISION: # $Id: readme.inf,v 1.12 2003-07-18 18:25:58 xshare Exp $

Toroidal magnetohydrodynamic equilibrium code ESC
=================================================

ESC is an easy to use, light-weight and numerically efficient fixed 
boundary magnetohydrodynamic (MHD) equilibrium code, which solves
the toroidal Grad-Shafranov equation for the poloidal flux. ESC computes
the cylindrical grid coordinates R(s, tau) and Z(s, tau) on flux
surfaces s=const and along poloidal rays tau=const, as well as
a number of radial profiles.

ESC's key features are:

* Newton scheme to iterate the right hand side of the Grad Shafranov
equation. (Most equilibrium codes are based on the slower Piccard scheme.)

* compact Fourier representation in the poloidal direction (only two modes
for Z are needed).

* radial cubic spline expansion.

The geometry is prescribed by two lists of points representing the plasma 
boundary coordinates. In addition, two input profile must be specified. 
These can be either:

* the pressure [mu0 Pa] and the parallel current <J. B>/Rext <B. grad phi> [mu0 A/m^2]

or

* the pressure [mu0 Pa] and the safety factor profile q

All profiles are defined on a uniform radial mesh in sqrt(toroidal flux).


o KEYWORDS

Equilibrium, MHD, toroidal geometry, flux, pressure, fixed bondary, 
parallel current, safety factor


o INSTALLATION

See INSTALL file provided with the distribution. 

o LANGUAGE

ESC is entirely written in C (C99). However, all interface routines are 
also callable from Fortran.

o PRECISION

8-byte reals (C-double) and default (4 or 8 byte) integers (C-int).

o PREREQUISITE

None, ESC is self contained.

o PLATFORM/OS

alpha/OSF1 (cc), alpha/Linux (gcc), Solaris 2 (cc), Linux (gcc), IRIX 6.x (cc), Cray/Unicos (cc)


o AUTHOR

L. E. Zakharov.

o CONTACT

A. Pletzer (pletzer@pppl.gov)

o USAGE

Because ESC comes as a list of procedures (subroutines) the code can easily
be embedded into another package. Thus, communication is at the application 
programming interface (API) level; there are no namelists, input or output 
files, except for dump files which allow ESC to be restarted from a previous
run. For an extensive documentation of the ESC API, refer to 

http://w3.pppl.gov/NTCC/ESC/html/globals.html

If you prefer the quick start approach, the following will get you up to
speed. In a few words, a typical ESC run involves (in Fortran): 

* initializing 

  call esc0(label) ! label is string: eg 'nstx.12345Z00'

* setting up the profiles

  ! j//mode
  call esc_p_jparallel(pressure, jparallel, rbound, zbound, nt1, &
     & rbtor, r0)

or

  ! q mode
  call esc_p_q(pressure, q, rbound, zbound, nt1, &
     & rbtor, r0)


* executing the code

  call esc(ifail, iflag, mpol, time)

* extracting the data

  ! get the geometry data on grid size np1*na1 (user's choice)
  ! np1: no of poloidal rays
  ! na1: no of radial surfaces
  call escGetR(R, np1, na1)  
  call escGetZ(Z, np1, na1) 

  ! many other data are available

* freeing the memory

  call esc1 ! clean up

More detailed examples can be found below. 

A note on efficiency:

Since the Newton scheme converges particularily fast when starting from 
a "good" equilibrium approximation, it is thus strongly advised to use any 
information available about neighbouring equilibria to accelerate the 
search of the final solution. This is generally possible when performing
for instance parametric studies (i.e. a parameter is gradually modified)
or in time dependent runs when the time intervals are small and hence
previous, old equilibrium solution can be used as initial guess. From our
experience of ESC in the TRANSP code, we found that ESC's CPU execution
time routinely runs at a fraction of a second on Pentium III PCs (10x faster
than any other equilibrium code within TRANSP) under these circumstances.

A note on robustness:

A drawback of the Newton scheme is that it is in general less robust than 
the Picard method. The following can be used to improve the robustness:

* When prescribing q, first run in j// mode and then switch to q mode, as
illustrated in the exmaple runs below.

* In some high beta cases, it may help to gradually increase the 
pressure. 

Note that increasing the pressure in a loop, or switching from j// to
q mode is hardly a performance issue; ESC's CPU time can still be 
expected to be of the order of a few seconds.


o EXAMPLE 1: SIMPLE FINITE PRESSURE, J// RUN

program ex1

  ! simple esc demo code
  ! f90 ex1.f90 -L<where libesc.a resides> -lesc
  ! on Alpha-Linux: 
  ! f90 -assume no2underscore ex1.f90 -L<where libesc.a resides> -lesc

  implicit none
  integer, parameter :: r8 = selected_real_kind(12,100)
  integer, parameter :: clockwise = 1, anticlockwise=0

  ! input data 

  integer, parameter :: ns = 21 ! number of radial profile points 
  ! (fixed by ESC)
  integer, parameter :: nt1= 17 ! number of boundary points (arbitrary)
  character*32, parameter :: label='TFTR.12345A67'
  real(r8) :: pressure(ns), jparallel(ns), rbound(nt1), zbound(nt1)
  real(r8) :: p(ns)

  ! output data

  integer, parameter :: ns_out=10  ! number of radial points 
  integer, parameter :: nt1_out= 33 ! number poloidal points
  real(r8) :: q(ns_out), r(nt1_out, ns_out), z(nt1_out, ns_out)

  integer i, ifail, ijob, mpol
  real(r8) :: t(nt1), s(ns), pi, rbtor, r0, time

  pi = acos(-1.0_r8)
  t = 2.0_r8*pi*(/ (real(i-1,r8)/real(nt1-1,r8), i=1, nt1) /)
  r0 = 1.0_r8    ! estimate of the magnetic axis position
  rbtor = 1.0_r8 ! r B_toroidal at vacuum
  rbound = r0 + 0.3_r8 * cos(t + 0.4_r8*sin(t))
  zbound =      0.5_r8 * sin(t)
  s = (/ (real(i-1,r8)/real(ns-1,r8), i=1, ns) /)
  pressure = 0.05_r8*(1.0_r8 - s**2)
  jparallel= 2.0_r8*(1.0_r8 - s**2)

  ! initialize

  call esc0(label)

  do i = 1, 3
     p = pressure*real(i-1,r8)/real(2,r8)

     ! set the geometry and profiles

     call esc_p_jparallel(p, jparallel, rbound, zbound, nt1, &
          & rbtor, r0)
     ijob = 0
     mpol = 4 ! number of poloidal modes (<=12)

     ! execute the code

     call esc(ifail, ijob, mpol, time)
     if(ifail /= 0) then
        print*,'ESC-FAILURE ',ifail
     endif
  end do

  ! access the data

  call escGetQ(q, ns_out)
  write(*,*)'Safety factor profile'
  write(*,'(5f10.5)') q
  call escGetR(r, nt1_out, ns_out, anticlockwise)
  call escGetZ(z, nt1_out, ns_out, anticlockwise)

  ! clean-up

  call esc1

end program ex1

o EXAMPLE 2: SIMPLE Q-MODE RUN

program ex2

  ! simple esc demo code
  ! f90 ex2.f90 -L<where libesc.a resides> -lesc
  ! on Alpha-Linux: 
  ! f90 -assume no2underscore ex1.f90 -L<where libesc.a resides> -lesc

  implicit none
  integer, parameter :: r8 = selected_real_kind(12,100)
  integer, parameter :: clockwise = 1, anticlockwise=0

  ! input data 

  integer, parameter :: ns = 21 ! number of radial profile points 
  ! (fixed by ESC)
  integer, parameter :: nt1= 17 ! number of boundary points (arbitrary)
  character*32, parameter :: label='TFTR.12345A67'
  real(r8) :: pressure(ns), qin(ns), jparallel(ns), rbound(nt1), zbound(nt1)
  real(r8) :: p(ns)

  ! output data

  integer, parameter :: ns_out=10  ! number of radial points 
  integer, parameter :: nt1_out= 33 ! number poloidal points
  real(r8) :: q(ns_out), r(nt1_out, ns_out), z(nt1_out, ns_out)

  integer i, ifail, ijob, mpol
  real(r8) :: t(nt1), s(ns), pi, rbtor, r0, time

  pi = acos(-1.0_r8)
  t = 2.0_r8*pi*(/ (real(i-1,r8)/real(nt1-1,r8), i=1, nt1) /)
  r0 = 1.0_r8    ! estimate of the magnetic axis position
  rbtor = 1.0_r8 ! r B_toroidal at vacuum
  rbound = r0 + 0.3_r8 * cos(t + 0.4_r8*sin(t))
  zbound =      0.5_r8 * sin(t)
  s = (/ (real(i-1,r8)/real(ns-1,r8), i=1, ns) /)
  pressure = 0.05_r8*(1.0_r8 - s**2)
  qin = 1.0_r8 + 2.0_r8*s**4
  jparallel = 2.0_r8*rbtor*(1.0_r8 - s**2)/(r0**2 *qin(1))

  ! initialize

  call esc0(label)

  ! Initial run: zero pressure & j//

  write(*,*) ' j// mode...'

  p = 0.0_r8 
  call esc_p_jparallel(p, jparallel, rbound, zbound, nt1, &
       & rbtor, r0)
  ijob = 0
  mpol = 4 ! number of poloidal modes (<=12)
  call esc(ifail, ijob, mpol, time)
  if(ifail /= 0) then
     print*,'ESC-FAILURE ',ifail
  endif

  ! switch to q mode

  write(*,*) ' q mode...'

  do i = 1, 5
     p = pressure*real(i-1,r8)/real(2,r8)
     ! set the geometry and profiles
     call esc_p_q(p, qin, rbound, zbound, nt1, &
          & rbtor, r0)
     ijob = 0
     mpol = 4 ! number of poloidal modes (<=12)
     call esc(ifail, ijob, mpol, time)
     if(ifail /= 0) then
        print*,'ESC-FAILURE ',ifail
     endif
  end do

  ! access the data

  call escGetQ(q, ns_out)
  write(*,*)'Safety factor profile'
  write(*,'(5f10.5)') q
  call escGetR(r, nt1_out, ns_out, anticlockwise)
  call escGetZ(z, nt1_out, ns_out, anticlockwise)

  ! clean-up

  call esc1

end program ex2

o TEST RUN

A test driver, drive.f90, is provided with the distribution. By default,
"gmake all" will produce the test program ../<MACHINE>/test/esc. This 
program will, when executed, generate the files prof.mtv and mesh.mtv which
can be compared against prof.mtv_ref and mesh.mtv_ref, respectively.

o DOCUMENTATION

http://w3.pppl.gov/NTCC/ESC/html/globals.html

o ERROR AND WARNING MESSAGES

Error and warning flags are raised when ifail differs from zero 
shortly after calling esc:

0 - normal operation;
1 - some problems;
2 - problems near the boundary;
4 - problems near the axis;
8 - problems in the middle;
3, 11, or 15 are symptomatic of localized inaccuracies (eg at the centre 
or near the edge)

In general, ESC fails when the plasma boundary exhibits a sharp V point 
near a separatrix.

o KNOWN ISSUES, PROBLEMS, LIMITATIONS AND BUGS

Results may differ between platforms at the 10th significant digit.

We encourage users to check independently the accuracy of the Grad-Shafranov solution when an error flag is returned. Depending on whether the error is localized at the edge or spread across the plasma, the solution can be salvaged or should be discarded. 

There is at present a hard-wired limit for the maximum number of points in the profile representation. It is the experience of the author that too many spline nodes can give rise to ringing, i.e. spurious oscillations in the profiles. The author recommends 20 points.

Due to some convergence problems in the number of poloidal modes, Esc at present does not take more than 12 modes. The author recommends 4 to 6 modes depending on the plasma boundary geometry. Periodically increasing this number from 4 to 6 in Transport codes has proved to be both efficient and robust. 

There is no control of the accuracy other than by specifying the number of poloidal (Fourier) modes.

A user reported problems near the edge, which could be reproduced by running Esc in pressure and j// mode, taking the q profile and running in pressure and q mode. The resulting j// profile showed structures near the edge.


o REFERENCES

L E Zakharov and A Pletzer, Theory of pertubed equilibria for solving the 
Grad-Shafranov equation, Phys Plasmas 6, 4693.
