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 /Rext [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 -lesc ! on Alpha-Linux: ! f90 -assume no2underscore ex1.f90 -L -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 -lesc ! on Alpha-Linux: ! f90 -assume no2underscore ex1.f90 -L -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 ..//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.