1 XPLASMA Xplasma -- D. McCune Apr. 2000 -- PPPL -- dmccune@pppl.gov Xplasma2 -- D. McCune Sept. 2006 -- PPPL -- dmccune@pppl.gov This document describes xplasma -- a set of routines for numerical representation of fusion plasma MHD equilibrium and plasma parameter profiles. Abstract (Sept 2006 updated version) -- The xplasma software offers a representation standard for tokamak plasma MHD equilibria and plasma parameter profiles defined over a variety of coordinates. Each xplasma dataset defines a "dictionary", a single namespace, in which each name refers to a specific data item-- typically, a list, a grid discretizing a coordinate, or a profile defined over one, two, or three gridded coordinates. Xplasma provides efficient methods for accessing the data associated with these items, or for adding any number of new items to the dictionary and collective dataset. An xplasma dataset (dictionary and all item contents) can be written to or restored from a NetCDF file by means of simple subroutine calls. The original implementation (xplasma 1.0) provided a single global xplasma dataset that was accessible by a set of fortran-77 style calls. An executable program could only contains one xplasma dataset at a time. The full fortran-95 implementation (xplasma 2.0) allows multiple instantiation of xplasma datasets in a single process, while preserving the fortran-77 interface with an implied reference to a "golden" xplasma pointer set aside for this purpose. Layered over efficient spline interpolation software (pspline [1]), xplasma has been used for several years to communicate data between integrated simulation codes (e.g. TRANSP, ONETWO) and NTCC modules (e.g. NUBEAM). There are also widely used applications for accessing experimental data (TRANSP & EFIT MDSplus trees) that use xplasma for internal storage of data items. It is planned to use xplasma for implementation of a plasma state component in the Fusion Simulation Project (SciDAC-2 SWIM project). [1] http://w3.pppl.gov/NTCC/PSPLINE/pspline.html 2 Organization_of_help_document The main topics are: Introductory Information F95 Upgrade (changes from xplasma vsn 1 to xplasma vsn 2). F95 (vsn 2) xplasma basic concepts Using F95 xplasma data Here, methods for accessing existing xplasma data are described. Initializing F95 xplasma data Here, first steps for creating a new xplasma object from scratch. Modifying F95 xplasma data Here, methods for writing data into xplasma are introduced. Error Handling Time evolution of xplasma Here, a discussion of what can and what cannot be easily changed. Creation and update of named xplasma data items More information on writing data into xplasma Global Items A few bits of modifiable globally accessible xplasma data Saving to File -- the remainder of the document: the original xplasma 1 (f77 interface) documentation. 2 Xplasma_Introductory_Information Xplasma is currently being used to communicate numerical representations of axisymmetric tokamak plasma timeslices between experimental databases and simulation codes, and between components or processes of multi- component or multi-process (distributed) simulation codes. Xplasma defines a sorted dictionary of named data items. It provides efficient mechanisms for looking up and accessing data. And it allows new data to be named and added in a straight-forward manner. The details are given in the subtopics. The contents of each xplasma dataset is thought of as a snapshot of the plasma state at a particular time in an experiment or simulation. 3 Xplasma_Data_Item_Identification Each xplasma data item has a dictionary entry with: a) a unique NAME[1] b) a unique integer ID code c) an AUTHOR name [1] d) a LABEL [2] e) a UNITS label [2] f) an integer data TYPE code The xplasma public module defines: integer, parameter, public :: xplasma_listType = 1 ! LIST data type code integer, parameter, public :: xplasma_coordType= 2 ! COORDinate data type integer, parameter, public :: xplasma_gridType = 3 ! GRID data type code integer, parameter, public :: xplasma_profType = 4 ! PROFILE data type code integer, parameter, public :: xplasma_blkbxType= 5 ! BLACK BOX data type code Notes: [1] up to 32 characters long, uppercase converted alphanumeric characters and underscore ("_"); first character must not be a numeric digit. [2] any ascii character string, or blank. Although xplasma does not restrict the length, these labels may be used by graphical display codes which may impose their own length restrictions e.g. by truncation. Recommendations: labels <~ 50 characters, units <~ 16 characters. 3 Xplasma_Data_Types COORDINATES-- an abstract coordinate, discretized by zero or more grids. There are two types: periodic, and non-periodic. Periodic coordinates are discretized by grids spanning [-pi:pi] or [0:2*pi]. For purposes of interpolation, there are no bounds: shifts of 2*pi*N are applied to bring interpolation coordinate data into bounds. For very large values of |N|, accuracy of interpolation will be lost. Non-periodic coordinates have fixed minimum and maximum values. Xplasma enforces prescribed limits for some coordinates; otherwise the minimum and maximum values are defined when the first grid discretizing the coordinate is provided. The xplasma public module defines (this is a partial list): integer, parameter, public :: xplasma_rho_coord = 1 ! radial flux coord. integer, parameter, public :: xplasma_theta_coord = 2 ! poloidal angle coord. integer, parameter, public :: xplasma_phi_coord = 3 ! toroidal angle coord. integer, parameter, public :: xplasma_R_coord = 4 ! R coord integer, parameter, public :: xplasma_Z_coord = 5 ! Z coord Xplasma methods exist for finding the periodicity attribute, min and max value (when defined), and number of GRIDs discretizing the coordinate. GRIDS-- specific discretization of a coordinate. A grid is a strictly ascending finite sequence of numbers, the first of which precisely matches the minimum value of the corresponding coordinate (or -pi or 0 in the case of a periodic coordinate), and the last of which precisely matches the maximum value of the corresponding coordinate (or pi or 2*pi in the case of a periodic coordinate). In gridded calculations, the GRID in the xplasma sense of the word is usually the set of boundaries of the finite zones of the calculation. The set of zone centers only constitutes a grid in the xplasma sense of the word, when it is augmented with the inner and outer boundaries of the gridded region. Xplasma methods exist for finding the size of any grid and retrieving its values. PROFILES-- 1d, 2d, or 3d interpolating functions based on array data defined over xplasma GRIDs. Interpolation method is defined when the profile is defined; there are at present 4 options: type: storage requirement: 1d 2d 3d NC step function N-1 (N1-1)* (N1-1)*(N2-1)*(N3-1) (N2-1) C0 piecewise linear N N1*N2 N1*N2*N3 f continuous C1 cubic Hermite 2*N 4*N1*N2 8*N1*N2*N3 df/dx continuous C2 cubic Spline 2*N 4*N1*N2 8*N1*N2*N3 d2f/dx2 continuous (N, N1, N2, N3 refer to the sizes of GRIDs discretizing known coordinates). When profile data is accessed by interpolation, it is only necessary to know the coordinates over which it is defined, not the specific grids. On the other hand, if the original data is needed, it is straightforward to acquire the specific grid information and the original data upon which the profile interpolating function is based. Cubic interpolation is affected by boundary conditions. For data defined over periodic coordinates, periodic boundary conditions are employed. For data defined over non-periodic coordinates, the boundary conditions are set when the data is defined. Xplasma methods exist for interpolating profiles as well as finding the underlying coordinates and grids. For multidimensional profiles, input coordinate data are tagged by COORDINATE or GRID ids; it is not necessary for the caller to know the order of storage of the data. Xplasma methods also allow straightforward means of acquiring information on the interpolation method and underlying grids. LISTS-- a list, identified by a single name, can contain 1 or more elements; each element consists of: -- a 32 character name (not part of the xplasma dictionary namespace). -- an integer value -- a (real*8) floating point value -- character data, any length, trailing blanks trimmed away. List lengths, element names, and contents are readily acquired by available xplasma methods. BLACK BOXES-- a named BLACK BOX consists of: -- an integer type code. The meaning is user defined meaning, but, it should be noted that xplasma itself uses black box type codes K, 0 <= K <= 101, to aid internal calculations and methods. User defined type codes should be outside this range. -- an integer 1d array of length specified by the user. -- a floating point (real*8) 1d array of length specified by the user. Black box datasets are used by xplasma for various purposes. Examples are: data caches for flux surface integration, representation of quantities defined over non-rectilinear grids. 3 Conventions Xplasma was developed for communication of data on axisymmetric tokamaks. Physics units: MKS; KeV for temperatures and particle energies. Flux coordinates: Radial flux coordinate "rho" -- rho = sqrt([enclosed-toroidal-flux]/[toroidal-flux-enclosed-by-boundary]) Poloidal angle coordinate "theta" -- theta=0 on large major radius side of plasma; theta increases along flux surfaces in counter-clockwise direction when plasma cross-section is viewed to the right of the axis of symmetry. But: data access methods provide optional controls to enable user data to be interpreted as theta_rev = 2*pi - theta grad(theta_rev) = -grad(theta), i.e. a clockwise poloidal angle coordinate. Toroidal angle coordinate "phi" -- increasing phi draws a circle in a counter-clockwise direction about the axis of symmetry, when the machine is viewed from above. ------------------- Note: adaptation of xplasma for non-tokamak applications may well be possible, but has not been attempted as of September 2006. A major effort would be required, e.g. to support stellarators or other non-axisymmetric plasma configurations. 4 Setting_theta_orientation Setting the orientation of the poloidal angle coordinate "theta" (aka "chi" in the fortran-77 interface) can now ONLY be done on a call-by-call basis. In the xplasma (vsn 2) fortran-95 interface, use the ccwflag[n] optional logical arguments. I.e. if the 2nd coordinate of a 2d profile evaluation call is a poloidal angle which is to be interpreted as having reversed orientation, set the optional argument ccwflag2=.FALSE.; if these arguments are omitted the normal orientation (.TRUE.) is assumed to apply. In the fortran-77 interface, most routines affected by the poloidal angle orientation include a vector size argument (ivec). To signal that the poloidal angle is to be considered reversed, set ivec to [the vector size] x (-1). I.e. a negative value of ivec signals poloidal coordinate interpretation is to be reversed. This will be verified in the documentation of the specific routines. 3 Interpolation_Methods xplasma supports several interpolation methods for data. The interpolation method is defined when an interpolating function is set up. In the xplasma software documentation, the desired interpolation method is indicated by an argument "iorder" in the setup routine. (Sometimes also referred to as the "spline type"). Generally the following methods are available: iorder name comment -1 zonal data step function 0 piecewise linear 1st derivatives defined, not continuous 1 Akima Hermite 1st derivatives continuous 2 Cubic Spline 1st & 2nd derivatives continuous Hermite is a piecewise cubic method that is C1 given 1st derivatives at the grid points. Akima Hermite uses a numerical method [1] for computing these derivatives which minimizes "ringing" in the case of noisy data. Cubic splines are C2 and preferable for smooth data, but can yield "ringing" artifacts in the case of noisy data. [1] Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1 3 Axisymmetric-core-plasma The core plasma is represented on a user defined flux coordinate grid (rho,theta), where "rho" is a flux surface label-- by convention, rho=sqrt(normalized-toroidal-flux) in xplasma itself, although physics applications can use any coordinate as long as the mapping to rho is known. After completion of the setup sequence, xplasma's core plasma model will contain * the flux surface geometry R(rho,theta), Z(rho,theta), in m. * equilibrium related magnetic profiles g(rho) (Tesla*m), and the poloidal flux profile psi(rho) (Wb/rad). Other magnetic profiles such as q(rho) are optional but often present. * the field components BR(rho,theta),BZ(rho,theta), and [mod(B)](rho,theta), all in Tesla. * any number of user defined 1d plasma parameter profiles f(rho) * any number of user defined 2d plasma parameter profiles g(rho,theta). The flux surface geometry and field components are represented as bicubic splines. User defined plasma parameter profiles are represented according to the user's specification, with splines, Hermite, piecewise linear, and zonal step function representations being available. For each profile there is an associated id number and name string. For reasons of efficiency, the id numbers are used for profile identification in the interpolation routines. Routines for fetching the id number from the name are provided and are efficient to log(N) of the number of names, but still should not be used in applications' innermost loops. 3 Limiter-description There are two options for the "limiter" or "vessel wall" specification: * a closed sequence of (R,Z) points outside the core plasma, and, * a list of bounding circles and infinite lines enclosing a space that contains the core plasma. The "scrape off plasma" is considered to be in the space inside the limiters but outside the core plasma. Generally, after the limiter is specified the user will want to set up an (R,Z) rectilinear grid which covers a rectangular space enclosing both the core plasma and the scrape-off-layer region between the core plasma and the limiters or vessel wall as given. Automatic generation of evenly spaced (R,Z) meshes is available as an option. Once this is done, xplasma will construct a piecewise bilinear function d(R,Z), where d gives (in R,Z units) the distance of (R,Z) from the nearest limiter, and * a negative number denotes a position inside the limiter * a zero denotes a position on the limiter * a positive number denotes a position outside the limiter. 3 Scrape-off-plasma-coverage After limiters are defined, a rectilinear (R,Z) grid can be constructed which covers both the core plasma and the scrape-off plasma. An (R,Z) representation of equilibrium, fields, and plasma profiles can then be built up as needed. 3 Intended_applications It is envisioned that xplasma might be used as follows: (1) a transport simulation code evolves a plasma MHD equilibrium and a set of plasma parameters. (2) "physics modules" are programmed to acquire data from the xplasma module and write data back into it. (3) the transport simulation code uses routines that add data to xplasma, to specify the plasma equilibrium and necessary plasma parameter profiles, and then it calls the physics module. (4) the physics modules reads its input data using the xplasma read routines, and then does its modeling calculations, developing a set of output profiles. (5) the physics module writes its results back into xplasma using more data definition routines, then exits. (6) the transport simulation code reads back the data created by the physics module. This reading and writing through xplasma will handle interpolation so that the main transport code and the physics module need not share COMMON blocks nor use the same numerical grids. 3 Physical_units_conventions Since the purpose of xplasma is to facilitate communication of data between separately developed physics models, it is necessary to recommend a convention regarding physical units. The convention adopted is based on emerging consensus in the tokamak physics modeling community: MKS for all quantities, except for temperatures and energies where KeV are used. Thus: MKS/KeV NOT cgs/eV (sorry, TRANSP...) "g" Tesla*m NOT gauss*cm "psi" Webers/rad NOT gauss*cm**2/rad "R,Z" m NOT cm "Te,Ti" KeV NOT eV "ne" m**-3 NOT cm**-3 [note it is customary to give the poloidal flux function normalized per unit toroidal angle; then Bpol = (1/R)*grad(psi), note (1/R) not (1/(2*pi*R))]. 3 Performance_considerations Considerable effort has gone into optimization of the performance of xplasma subroutines. In particular, xplasma is layered over the pspline library, which has been shown to work well on modern vector and super- scalar computer architectures. But, general purpose routines (e.g. for interpolation) necessarily have a fair amount of control logic and error checking to support general purpose use. This will inevitably result in performance somewhat less than what is optimally achievable for a specific code. Codes for which e.g. field or metric interpolation are critical performance considerations likely already have optimized code for these purposes. In many cases it will be best simply to use xplasma as a way to load the input data structures for these optimized routines at the start of the calculation, and use the optimized routines (not xplasma calls) in the innermost loops of the calculation. 3 Non-axisymmetry Although the xplasma interface is designed to allow eventual support of non-axisymmetric field equilibria, the September 2006 version of the software is limited to axisymmetric configurations only. 2 F95_Upgrade A project has been completed to upgrade xplasma for the Fusion Simulation Project. The upgrade includes: -- simplified fortran-95 interface; simplified underlying data structures. -- object-oriented implementation, allowing multiple xplasma "instances" to co-exist independently in memory. (This is fortran-95; the term "object-oriented" is used loosely). -- ability to update equilibrium without deleting entire object contents and rebuilding from the ground up. The fortran-77 style original xplasma interface is maintained for the most part, but a few routines have been dropped-- see subtopic. The f77 routines all refer implicitly to a reserved xplasma object pointer set aside for this purpose. To use f77 routines to modify any other xplasma object, it is necessary to reset the f77 object pointer: see the subtopic 3 Reset_F77_pointer To reset the f77 xplasma pointer-- which defines the xplasma object to which all legacy f77 xplasma calls refer-- use the following module and its pair of contained routines grab_f77_xplasma(swant,ierr) release_f77_xplasma(ierr) -------------------------------------- details: use xplasma_obj_instance ... ! the following points to an xplasma object, to which future ! f77 xplasma calls should refer. This pointer has been initialized ! elsewhere...! type (xplasma), pointer :: swant integer :: ierr ! status code, 0=normal on exit ! reset the f77 xplasma pointer with this call: call grab_f77_xplasma(swant,ierr) if(ierr.ne.0) [handle-unlikely-error] To later restore the f77 xplasma pointer to its prior state, use: call release_f77_xplasma(ierr) if(ierr.ne.0) [handle-unlikely-error] Note that these calls should occur in pairs! If (e.g. due to error handling logic) the pairing is broken, there are two likely serious consequences: (a) F77 xplasma calls interact unexpectedly with a different xplasma object, with unpredictable results; (b) eventual f77 xplasma pointer stack underflow or overflow error. 3 Deprecated_F77_Routines The "C" interface to fortran-77 xplasma will be frozen. The "R4" interface to fortran-77 xplasma will be frozen. All the fortran-95 interfaces are in REAL*8 precision. The following fortran-77 style routines will be dropped in the upgraded version: eqm_flxbdy / eqm_rbdy -- "move limiter to plasma boundary" / "restore original limiter". dropped because: (a) no evidence of use; (b) error-prone; if "eqm_flxbdy" call is not followed by "eqm_rbdy", xplasma suffers a side effect; distance to limiter calculations, used e.g. for fast ion orbit loss calculations, are modified. advice: codes wishing to treat the plasma boundary as the limiter can use routines like "eq_dfast" to compute distance-to-boundary. Such routines return a positive number for locations beyond the plasma boundary, and a negative number for locations within the plasma boundary. eqm_chi_cwdir_set, eqm_chi_cwdir_restore, eqm_chi_cwdir, eq_get_chi_cwdir control clockwise/counter-clockwise orientation of poloidal angle variable. Dropped due to being error prone-- side effect of a subroutine's reset of chi_cwdir seen by all subsequent calls to xplasma routines. advice: In the xplasma fortran-95 interface, optional arguments are provided to specify angle coordinate orientations. In affected fortran-77 interface routines, the poloidal angle coordinate orientation can be reversed, separately on each call, by setting the argument vector length to a negative number. Details are provided in the updated subroutine documentation, below. eqm_phi_cwdir_set, eqm_phi_cwdir_restore, eqm_phi_cwdir, eq_get_phi_cwdir dropped: no evidence of use. Optional arguments in fortran-95 interface only. eqm_spordr_set, eqm_spordr_get dropped: xplasma feature to control spline storage mode (compact or non-compact) was found to be unnecessary and has been removed. eqm_nfast dropped: control of resolution of a 2nd (R,Z)->(rho,theta) inverse map, no longer necessary. 3 Nomenclature The original fortran-77 xplasma API used "chi" to designate the poloidal angle coordinate in a 2d flux coordinate representation of a tokamak plasma equilibrium. The new fortran-95 interface uses "theta" to designate this coordinate. However, since the fortran-77 and fortran-95 interfaces are both maintained, both names persist in the xplasma API. Hopefully this will not cause much confusion. 2 F95_Xplasma_Basics In order to use the Fortran-95 Xplasma Interface in code, it is necessary to be able to refer to xplasma objects. A consequence of the design goal of allowing multiple xplasma objects to co-exist, is that each subroutine that accesses or manipulates xplasma data must include an argument which specifies which instance of xplasma is to be affected or used. This argument is a "pointer to an xplasma object". A pointer is used, because each xplasma object can be quite large, such that it is sometimes advantageous to reuse existing xplasma code by means of resetting a pointer, rather than having to explicitly copy the entire contents from one xplasma object to another. Advanced users may want for various reasons to invent additional fortran-95 modules which declare xplasma objects and pointers. However, for starting out it is recommended to use the public module provided with the software: use xplasma_obj_instance This causes the following pointers to be defined, with suggested semantics: s -- pointer to currently active xplasma object F77 xplasma calls implicitly refer to this one. s1-- pointer to previous xplasma object (e.g. from a prior timestep in a time-dependent simulation). s2-- spare. The module xplasma_obj_instance also imports all public xplasma modules and definitions. INITIALIZATION-- xplasma_obj_instance defines SAVEd global data. Because fortran-95 does not yet allow pointers to be initialized at compile time to point to any particular object, the pointers are initially NULL. To be activated, the pointers in this module must be initialized with by calling a public subroutine defined in the module: call xoi_init(ierr) This call simply associates the pointers with specific xplasma objects inside the module; no data initialization is performed. It also sets a flag inside the module indicating that initialization has been performed, so that if the call is repeated, no action is taken. The integer argument "ierr" returns a status code, 0=normal. Almost all xplasma subroutine calls return such a status code. (NOTE: spelling variant: some xplasma routines spell this "ier" in their argument lists). Nearly all f95 xplasma calls take an xplasma object pointer as their first argument. In the documentation, the symbol "s" is used for this pointer. In some applications of xplasma, the symbol can be taken as short hand for "plasma state". Almost all xplasma calls also have a status code "ierr" as their final required argument. Ordering of calling arguments in f95 xplasma calls: Example: call xplasma_eval_prof(s,id,x,result,ierr,ideriv1=1,n_out_of_bounds=iout) 1st argument -- s -- xplasma object pointer next set of arguments -- id,x -- required inputs (here, the profile id and the x values at which to evaluate it) next set of arguments -- result,ierr -- required outputs (here, the result of the evaluation, and an error status code) ("ierr" is almost always the last required output argument). next set of arguments -- ideriv1 -- optional inputs (here, requesting the 1st derivative of the profile w.r.t. x). (the documentation will indicate that the default is to return the profile value unless otherwise directed by optional arguments). final set of arguments -- n_out_of_bounds -- optional outputs (here, the number of x values that were out of bounds, with respect to the coordinate over which the profile identified by "id" is defined). Of course, not all xplasma subroutines contain arguments in every one of these categories. Nevertheless, the ordering of arguments is consistent with this ordering for all f95 xplasma public interfaces. 3 Xplasma_public_modules xplasma_obj_instance uses "xplasma_obj_instance_mod" (where xplasma objects, pointers and xoi_init are defined). uses "xplasma_definitions", which... uses "xplasma_profs" -- methods for creating certain profiles. uses "xplasma_mcgrid"-- methods for creating and accessing a type of fast ion distribution function "black box" and related items. uses "xplasma_rzgeo" -- methods for defining MHD equilibrium geometry and related profiles. uses "xplasma_flxint"-- methods for computing flux surface integrals and averages. uses "xplasma_sol" -- methods for defining limiter and scrape-off region. uses "xplasma_ctran" -- coordinate transformation methods and evaluation of metric jacobian. uses "xplasma_obj" -- the kernel -- dictionary and supported data types; error handling support; NetCDF i/o. 3 Xplasma_public_constants The following are public parameters that are defined in the public xplasma modules, which may be helpful in user applications: Data types known to the xplasma dictionary: integer, parameter, public :: xplasma_listType = 1 ! LIST data type code integer, parameter, public :: xplasma_coordType= 2 ! COORDinate data type integer, parameter, public :: xplasma_gridType = 3 ! GRID data type code integer, parameter, public :: xplasma_profType = 4 ! PROFILE data type code integer, parameter, public :: xplasma_blkbxType= 5 ! BLACK BOX data type code Commonly used coordinates: integer, parameter, public :: xplasma_rho_coord = 1 ! radial flux coord. integer, parameter, public :: xplasma_theta_coord = 2 ! poloidal angle coord. integer, parameter, public :: xplasma_phi_coord = 3 ! toroidal angle coord. integer, parameter, public :: xplasma_R_coord = 4 ! R coord integer, parameter, public :: xplasma_Z_coord = 5 ! Z coord Designation for coordinate subsets (usable as argument in a handful of xplasma routines): integer, parameter, public :: xplasma_flux_coords = 1 ! generic flux coord. integer, parameter, public :: xplasma_cyl_coords = 2 ! generic cyl. coord. Author name usable to designate quantities that should be deleted when the plasma MHD equilibrium changes: character*32, parameter, public :: xplasma_xmhd = '__XPLASMA_MHDEQ_DERIVED' The following arguments define algorithm choices for the [R,Z] -> [rho,theta] inverse map, in order from the slowest but most accurate to the fastest but most approximate... integer, parameter, public :: xp_imap_newton = 1 integer, parameter, public :: xp_imap_polar = 2 integer, parameter, public :: xp_imap_rzlinear = 3 3 Reserved_names For reasons of compatibility with the original xplasma implementation and existing uses of the software, certain item names related to the plasma MHD equilibrium are reserved. These are: name -- definition (all are profiles unless otherwise specified) G -- g(rho) = R*B_phi PSI -- psi(rho); (1/R)*grad(Psi) = B_pol P -- P(rho), plasma pressure for MHD equilibrium Q -- q(rho), rotational transform d[toroidal_flux]/d[poloidal_flux] R -- R(rho,theta), (with Z) position of flux surfaces Z -- Z(rho,theta), (with R) position of flux surfaces PSI_RZ -- Psi(R,Z) BR -- BR(rho,theta), field component in R direction BZ -- BZ(rho,theta), field component in Z direction BPHI -- Bphi(rho,theta), field component in phi direction BMOD -- mod(B)(rho,theta), total field magnitude BR_RZ -- BR(R,Z), field component in R direction BZ_RZ -- BZ(R,Z), field component in Z direction BPHI_RZ -- Bphi(R,Z), field component in phi direction BMOD_RZ -- mod(B)(R,Z), total field magnitude MAG_AXIS -- list, [R,Z] at magnetic axis BDY_RZ_MINMAX -- list,{Rmin,Rmax,Zmin,Zmax} of plasma boundary (nominal last closed flux surface). User applications should avoid using these names for storage of items in xplasma. Several additional named items created by XPLASMA start with "__", two consecutive underscores. Users should avoid defining names that start with this prefix. 3 User_defined_modules Although the public module xplasma_obj_instance defines sevaral xplasma instances and associated pointers, the user may want to define his own xplasma instances at some point. The existing modules old_xplasma/xplasma_obj_instance.f90 old_xplasma/xplasma_obj_instance_mod.f90 can be used as examples for how to do this. The user's code will need to see use xplasma_definitions in order to make use of xplasma calls. 2 Using_F95_Xplasma_Data In this section, the new Fortran-95 Xplasma interface is described, with emphasis on methods for using data in an existing dataset. The xplasma pointers (s,s1) defined in xplasma_obj_instance are presumed to be available in the coding examples. The argument "ierr" is always a return status code, 0=normal. For simplicity of examples, error checking / error handling code is usually omitted -- but real codes should check the error codes! To access any xplasma call, the calling code must declare use xplasma_obj_instance or else use a user-defined module that otherwise provides access to xplasma object and interface definitions. In all f95 code examples it is assumed that this has been done. 3 Example_code The xplasma module comes with a debug library, which supports SGlib/ElVis graphical examination of xplasma datasets. The library is called xplasma_debug.a; the main source code for this library is xplasma_debug_module.f90-- this code has examples of pretty much every kind of access to xplasma data. Given an xplasma file path, or a URL pointing to a publicly accessible xplasma file, the data can be plotted running "plot_xplasma", a command line interactive fortran program, which uses the debug library. The f77 interface to xplasma has been re-implemented by using the f95 xplasma routines internally. This interface resides in the library old_xplasma.a; the source routines for this library, e.g. eqm_rhofun.f90, contain examples of use of xplasma f95 routines. This may be useful especially for anyone already familiar with the f77 routines. 3 Copy_and_read ! copy contents of "s" to "s1". Prior contents of "s1" are removed; all ! associated memory is freed; "s" is not changed. call xplasma_copy(s,s1,ierr) ! read in new contents of "s" from a file. Prior contents of "s", if any, ! are removed, and associated memory is freed. "s" and "s1" are not sharing ! any memory, so, the removal of the prior contents of "s" has no effect on ! "s1". call xplasma_read(s,filename,ierr) ! "filename" is a character string variable containing a file path to ! a read accessible xplasma NetCDF file. !-------------------------- ! example: use xplasma_obj_instance use my_module ! defines "file_path" and "runid" character string variables. implicit NONE character*100 filename filename = trim(file_path)//trim(runid)//'_my_xplasma.cdf' call xplasma_read(s,filename,ierr) if(ierr.ne.0) then write(6,*) ' ? xplasma open or read failure: '//trim(filename) call xplasma_error(s,ierr,6) call bad_exit endif Note: if the filename argument has the form of a URL (i.e. starts with "http://"), xplasma_read will attempt to spawn a /usr/bin/curl process to download the indicated URL to a local disk file (in the process current working directory), and will then attempt a read of the downloaded file. 3 Global_information Certain global information, such as a global label, the current simulation or experiment time in seconds, sign or orientation of the toroidal field and toroidal plasma current, are stored in each xplasma. The accuracy of this information depends on the software that created the original xplasma object or file. To retrieve this information use the following routine. Other than the status code, all arguments are optional, allowing the user to select which piece or pieces of information are desired. subroutine xplasma_global_info(s,ier, & initLabel,time,axisymm,scrapeoff,bphi_ccw,jphi_ccw,nitems, & prof_counter,eq_counter,bdytol,ajac_maxVaR,rzOrder,kmom) ! this routine returns global xplasma object information, according ! to which optional arguments are provided... type (xplasma), pointer :: s integer, intent(out) :: ier ! error set iff s is uninitialized. ! this is the label passed in an xplasma_(re)Init call... character*(*), intent(out), optional :: initLabel ! can also use xplasma_time_get for this: real*8, intent(out), optional :: time logical, intent(out), optional :: axisymm ! axisymmetry flag logical, intent(out), optional :: scrapeoff ! SOL flag integer, intent(out), optional :: bphi_ccw ! Bphi counter-clockwise flag ! 0 means: never specified; 1: counter-clockwise (CCW); -1: CW integer, intent(out), optional :: jphi_ccw ! Jphi counter-clockwise flag ! 0 means: never specified; 1: counter-clockwise (CCW); -1: CW integer, intent(out), optional :: nitems ! number of items in xplasma ! dictionary at time of call. integer, intent(out), optional :: prof_counter ! profile counter ! global counter incremented for each creation or update of each profile integer, intent(out), optional :: eq_counter ! profile counter ! global counter incremented for each update of MHD equilibrium real*8, intent(out), optional :: bdytol ! out-of-range error tolerance ! relative tolerance applied for interpolations over non-periodic ! coordinates; for periodic coordinates shifts of 2*n*pi are applied ! as necessary to bring interpolations in range real*8, intent(out), optional :: ajac_maxVar ! det[J] variation tol. ! relative tolerance for variations of det[J] within a flux surface; ! used to check for singular or near-singular metric tensor in ! R(rho,theta),Z(rho,theta) flux surface definitions. integer, intent(out), optional :: rzOrder ! spline fit order for MHD ! equilibrium geometry R(rho,theta) z(rho,theta) data ! 1 for C1 Hermite; 2 (the default) for C2 Spline. integer, intent(out), optional :: kmom ! #moments for Fourier ! representation of equilibrium Also available: subroutine xplasma_time_get(s,ztime,ier) ! fetch the time (seconds) affiliated with xplasma object "s". type (xplasma), pointer :: s real*8, intent(out) :: ztime integer, intent(out) :: ier 3 Equilibrium_derived_information The plasma MHD equilibrium defines the transformation between magnetic flux coordinates (rho,phi,theta) and cylindrical coordinates (R,phi,Z). It also defines the magnetic field and flux surface geometry. Xplasma provides a routine for transforming coordinate sets between cartesian, cylindrical, and flux coordinate spaces. It also provides for retrieval or calculation of extrema-- e.g. minimum and maximum R, Z, and mod(B) on any flux surface, and for finding the distance to the nearest point on any flux surface. These are described in the subtopics... 4 Magnetic_axis To get the magnetic axis postion (m) and field strength(T): real*8 :: R_axis, Z_axis, B_axis Use this call: call xplasma_mag_axis(s,ierr, raxis=R_axis, zaxis=Z_axis, baxis=B_axis) Note that the last 3 arguments are keyworded and optional-- any user selected subset of the information may also be had using the same subroutine. The B field information is mod(B); to get the sign use xplasma_global_info and fetch "bphi_ccw" (see preceding topic). Note that this is the actual field strength, not the "vacuum" field-- i.e. g(0)/R_axis, not g(1)/R_axis, where g = R*B_phi, constant on each flux surface-- xplasma profile name "G". 4 Flux_surface_extrema To find the minimum and maximum R and/or Z and/or |B| occuring on a plasma flux surface: real*8 :: R_min, R_max, Z_min, Z_max ! spatial extrema, m real*8 :: B_min, B_max ! extrema in mod(B), T For the plasma boundary surface: subroutine xplasma_RZminmax_plasma(s, rmin,rmax, zmin,zmax, ier, & ccw_theta,i2pi, thRmin,thRmax, thZmin,thZmax) ! return the plasma boundary R & Z min & max type(xplasma), pointer :: s real*8, intent(out) :: Rmin,Rmax real*8, intent(out) :: Zmin,Zmax integer, intent(out) :: ier ! optionally return the poloidal angle locations of the extrema logical, intent(in), optional :: ccw_theta ! Theta CCW in poloidal plane, default=T integer, intent(in), optional :: i2pi ! angle coord normalization options ! i2pi=1 (default) returned value in range [0,2pi] ! i2pi=2 returned value in range [-pi,pi] real*8, intent(out), optional :: thRmin,thRmax real*8, intent(out), optional :: thZmin,thZmax The optional controls ccw_theta and i2pi only have an effect if theta location of extrema are being returned. For any plasma surface identified by "rho": subroutine xplasma_RZminmax(s,rho,ier, phi, & ccw_theta,i2pi, & rmin,rmax, zmin,zmax, & thRmin,thRmax, thZmin,thZmax) ! very accurately determine Rmin/max, Zmin/max of a flux surface. type (xplasma), pointer :: s real*8, intent(in) :: rho ! required input: flux surface. integer, intent(out) :: ier ! completion code real*8, intent(in), optional :: phi ! toroidal angle location ! ignored for axisymmetric cases real*8, intent(out), optional :: rmin,rmax ! Rmin & Rmax of surface real*8, intent(out), optional :: zmin,zmax ! Zmin & Zmax of surface ! optionally return the poloidal angle locations of the extrema logical, intent(in), optional :: ccw_theta ! Theta CCW in poloidal plane, default=T integer, intent(in), optional :: i2pi ! angle coord normalization options ! i2pi=1 (default) returned value in range [0,2pi] ! i2pi=2 returned value in range [-pi,pi] real*8, intent(out), optional :: thRmin,thRmax real*8, intent(out), optional :: thZmin,thZmax The optional controls ccw_theta and i2pi only have an effect if theta location of extrema are being returned. For the minimum and maximum mod(B) on any surface "rho": subroutine xplasma_Bminmax(s,rho,ier, phi, & ccw_theta,i2pi, bmin,bmax, thbmin,thbmax) use eqi_rzbox_module ! very accurately determine Bmin/max of a flux surface type (xplasma), pointer :: s real*8, intent(in) :: rho ! required input: flux surface. integer, intent(out) :: ier ! completion code real*8, intent(in), optional :: phi ! toroidal angle location ! ignored for axisymmetric cases logical, intent(in), optional :: ccw_theta ! Theta CCW in poloidal plane, default=T integer, intent(in), optional :: i2pi ! angle coord normalization options ! i2pi=1 (default) returned value in range [0,2pi] ! i2pi=2 returned value in range [-pi,pi] real*8, intent(out), optional :: Bmin,Bmax ! Bmin & Bmax of surface real*8, intent(out), optional :: thBmin,thBmax The optional controls ccw_theta and i2pi only have an effect of thBmin and/or thBmax are being returned. 4 Computational_domain_extrema For the computational domain (this generally encompasses a rectangle in [R,Z] space that covers and extends beyond the plasma boundary; if there is no [R,Z] rectangle defined, the plasma boundary extrema are returned: subroutine xplasma_RZminmax_extended(s, rmin,rmax, zmin,zmax, ier, & sol, id_Rgrid, id_Zgrid) ! return the R & Z min & max of the extended mapped region; if ! there is no such region return the R & Z min & max of the ! plasma type(xplasma), pointer :: s real*8, intent(out) :: Rmin,Rmax real*8, intent(out) :: Zmin,Zmax integer, intent(out) :: ier 4 Coordinate_transform A single interface "xplasma_ctrans" enables computation of any coordinate transformation between cartesian, cylindric, and flux coordinates. The interface xplasma_ctrans aliases two routines xplasma_ctran1 and xplasma_ctrann which accept scalar and vector arguments, respectively. The transformations are defined as follows: (R,phi,Z) --> (x,y,z) x = R*cos(phi) y = R*sin(phi) z = Z (x,y,z) --> (R,phi,Z) R = sqrt(x**2+y**2) if R > 0, phi = atan2(y,x), otherwise phi = 0 Z = z (rho,theta) --> (R,Z) bicubic splines R(rho,theta), Z(rho,theta) are evaluated. An extrapolation of the system defines (rho,theta) --> (R,Z) for rho > 1 -- it is continuous at the boundary rho=1 with the interior splines, but not continuously differentiable. (R,Z) --> (rho,theta) various options for estimating the inverse map are provided. The options control an accuracy/speed trade off. If many evaluations are needed (e.g. for straight line tracking in flux coordinate space), the fastest option is to compute a bilinear profile interpolant pair rho(R,Z) and theta(R,Z) on some grid. The map is calculated accurately at the grid points; the interpolation of theta(R,Z) is handled carefully in the vicinity of the theta branch cut (2*pi discontinuity). This method is supported by the "maptype=xp_imap_rzlinear" option and is sufficient for many purposes. For greater cost in cpu time, the (R,Z) splines can be inverted using a Newton method with accuracy tolerances set down near machine precision if necessary: "maptype=xp_imap_newton", the default, combined with "tol="; code comments: real*8, intent(in), optional :: tol ! map tolerance ! tol applies to Newton map (maptype = xp_imap_newton = default) ONLY ! (R_in,Z_in)->(rho_out,theta_out): ! |R_in-R(rho_out,theta_out)| < tol*max(R_in,|Z_in|) ! |Z_in-Z(rho_out,theta_out)| < tol*max(R_in,|Z_in|) ! default: s%bdytol ! s%bdytol's own default value is 1.0d-10 ! cf xplasma_bdytol_set; xplasma_global_info. An option of intermediate speed that is generally accurate to ~1.0d-4 and requires no pre-tabulation of map data, can be had by setting "maptype=xp_imap_polar". Almost all arguments are optional-- the presence of arguments defines what transformation is desired. Examples: INTEGER, PARAMETER :: R8=SELECTED_REAL_KIND(12,100) real(R8) :: r,z,rho,theta ! real*8 also works... real(R8) :: r2(2),z2(2),rho2(2),theta2(2) ! scalar (rho,theta) -> (R,Z) call xplasma_ctrans(s,ierr,rho_in=1.0_R8,theta_in=0.0_R8, & r_out=r,z_out=z) ! scalar (R,Z) -> (rho,theta) call xplasma_ctrans(s,ierr,r_in=r,z_in=z,rho_out=rho,theta_out=theta) r2(1)=r r2(2)=r-0.01_R8 z2(1:2)=z ! vector (R,Z) -> (rho,theta) call xplasma_ctrans(s,xp_use_coord_vecsize,ier, & r_in=r2, z_in=z2, rho_out=rho2, theta_out=theta2) "xp_use_coord_vecsize" is a logical parameter (value is .TRUE.) defined by the xplasma public module; it means the vector size is inferred from the arguments (all sizes must match). The alternative is "xp_set_coord_vecsize" (.FALSE.) which allows a size less than or equal to the smallest passed vector size to be used. Full interface: interface xplasma_ctrans module procedure xplasma_ctran1 ! scalar module procedure xplasma_ctrann ! vector end interface Scalar interface: subroutine xplasma_ctran1(s,ier, & x_in,y_in,z_in, r_in,phi_in, rho_in,theta_in, & tol, maptype, ccw_theta,ccw_phi, i2pi, & x_out,y_out,z_out, r_out,phi_out, rho_out,theta_out, & nregion,cosphi,sinphi) Vector interface: subroutine xplasma_ctrann(s,vsize_flag,ier, & isize, & x_in,y_in,z_in, r_in,phi_in, rho_in,theta_in, & tol, maptype, ccw_theta,ccw_phi, i2pi, & x_out,y_out,z_out, r_out,phi_out, rho_out,theta_out, & nregion,cosphi,sinphi) ! VECTOR coordinate mapping interface type (xplasma), pointer :: s logical, intent(in) :: vsize_flag ! vector size flag ! .TRUE. means size is to be inferred from passed vector sizes, ! all of which must match; ! .FALSE. means size is to be set by optional argument "isize"; ! all vectors must size.ge.isize then. integer, intent(out) :: ier ! completion code 0=OK !------------------------ ! optional inputs... integer, intent(in), optional :: isize ! vector size if(vsize_flag) real*8, intent(in), dimension(:), optional :: x_in,y_in,z_in ! Cartesian Coords in real*8, intent(in), dimension(:), optional :: r_in,phi_in ! Cylindric Coords (z also) real*8, intent(in), dimension(:), optional :: rho_in,theta_in ! Flux Coords (phi also) !------------------------------------------------------------ ! the following applies to (R,Z) -> (rho,theta) "inverse" map: real*8, intent(in), optional :: tol ! map tolerance ! tol applies to Newton map (maptype = xp_imap_newton = default) ONLY ! (R_in,Z_in)->(rho_out,theta_out): ! |R_in-R(rho_out,theta_out)| < tol*max(R_in,|Z_in|) ! |Z_in-Z(rho_out,theta_out)| < tol*max(R_in,|Z_in|) ! default: s%bdytol; s%bdytol's own default value is 1.0d-10 integer, intent(in), optional :: maptype ! = 1 = xp_imap_newton = default -- iterative Newton root finder ! with initial guess generator using polar map (slowest but ! most accurate, accuracy can be controlled with tol. ! = 2 = xp_imap_polar -- polar map inversion method ! much faster, not as accurate as Newton map ! = 3 = xp_imap_rzlinear -- bilinear inverse map -- after ! initialization (using polar map), the fastest, but less ! accurate still (depends on R,Z grid resolution). This ! method requires a user defined (R,Z) rectangle grid; if ! none is available the polar map is used instead. !------------------------------------------------------------ logical, intent(in), optional :: ccw_phi ! Phi CCW view from above, default=T logical, intent(in), optional :: ccw_theta ! Theta CCW in poloidal plane, default=T integer, intent(in), optional :: i2pi ! angle coord normalization options ! i2pi=1 (default) returned value in range [0,2pi] ! i2pi=2 returned value in range [-pi,pi] !------------------------- ! optional outputs... real*8, intent(out), dimension(:), optional :: x_out,y_out,z_out ! Cartesian Coords in real*8, intent(out), dimension(:), optional :: r_out,phi_out ! Cyl. Coords (z also) real*8, intent(out), dimension(:), optional :: rho_out,theta_out ! Flux Coords (phi also) real*8, intent(out), dimension(:), optional :: cosphi,sinphi ! phi sin & cos integer, intent(out), dimension(:), optional :: nregion ! region code ! 0=not checked; 1=inside core plasma; 2=outside plasma but inside ! mapped region; 3=beyond mapped region. 4 Metric_jacobian The 2x2 form of the metric jacobian [J] [ dR/drho dR/dtheta ] [ dZ/drho dZ/dtheta ] This is useful for various calculations. For example, 2*pi*R*det[J]*drho*dtheta is the differential volume element; grad(rho) and grad(theta) are readily derived from [J]; for f=f(rho,theta), grad(f) = (df/drho)*grad(rho) + (df/dtheta)*grad(theta). There is a generic interface xplasma_rzjac that maps to two specific routines: xplasma_rzjac1 for all scalar arguments; xplasma_rzjacn for all vector arguments. The metric jacobian array elements, or R, or Z, or det[J], or components of grad(rho) or grad(theta), or any combination thereof, can be acquired with xplasma_rzjac. The full interface follows... Full interface: interface xplasma_rzjac module procedure xplasma_rzjac1 ! scalar module procedure xplasma_rzjacn ! vector end interface Scalar interface: subroutine xplasma_rzjac1(s,rho_in,theta_in,ier, ccwflag1, & r1,z1,drdrho1,dzdrho1,drdtheta1,dzdtheta1, & rzdetj1,drhodr1,drhodz1,dthdr1,dthdz1) Vector interface: subroutine xplasma_rzjacn(s,rho_in,theta_in,ier, & ccwflag, r,z,drdrho,dzdrho,drdtheta,dzdtheta, & rzdetj, drhodr,drhodz, dthdr,dthdz) ! evaluate terms associated with 2d Jacobean -- vector version ! [dR/drho dR/dtheta] ! [dZ/drho dZ/dtheta] ! *** all vectors must be of the same length *** type (xplasma), pointer :: s real*8, dimension(:) :: rho_in ! vector of radial coordinate in real*8, dimension(:) :: theta_in ! vector of poloidal angle coordinate in ! (see ccwflag below) integer, intent(out) :: ier ! completion code 0=OK ! Poloidal angle orientation: TRUE (default) if dZ/dtheta is ! positive on the large major radius size and negative on the ! small major radius side of each flux surface; FALSE for the ! reverse. logical, intent(in), optional :: ccwflag real*8, intent(out), dimension(:), optional :: r,z ! (R,Z) values if desired real*8, intent(out), dimension(:), optional :: drdrho,dzdrho ! rho derivatives if desired real*8, intent(out), dimension(:), optional :: drdtheta,dzdtheta ! theta derivatives if desired real*8, intent(out), dimension(:), optional :: rzdetj ! determinant of 2x2 jacobian matrix real*8, intent(out), dimension(:), optional :: drhodr,drhodz ! grad(rho) components if desired real*8, intent(out), dimension(:), optional :: dthdr,dthdz ! grad(theta) components if desired 4 Distance_to_flux_surface For any point [R,Z] or a vector of points {R(1:n),Z(1:n)}, it is possible to very accurately find the closest point of approach on a chosen plasma flux surface. The public interface xplasma_bdfind aliases two routines, xplasma_bdfind1 with scalar input and output arguments, and xplasma_bdfindn with vector input and output arguments, to compute these distances. Three methods are available, which trade accuracy for speed-- as described in the interface comments. Distances returned are in meters; positive numbers denote points outside the flux surface; negative numbers denote points inside. The user can optionally receive information on the location of the point on the flux surface which is closest to the input point(s). Both these routines are aliased to "xplasma_bdfind" in the module interface. Full interface: interface xplasma_bdfind module procedure xplasma_bdfind1 ! scalar module procedure xplasma_bdfindn ! vector end interface Scalar arguments: subroutine xplasma_bdfind1(s,zR,zZ,ier, & maptype, phi_in, rho_in, ccw_theta, ccw_phi, i2pi, outside_only, & theta_out, phi_out, dist) Vector arguments: subroutine xplasma_bdfindn(s,isize,zR,zZ,ier, & maptype, phi_in, rho_in, ccw_theta, ccw_phi, i2pi, outside_only, & theta_out, phi_out, dist) ! accurately find the distance of points (R,Z) from a flux surface ! and the closest point on the flux surface ! although phi arguments are supplied to support a possible future ! upgrade to 3d, these have no effect on axisymmetric calculations. use eqi_rzbox_module type (xplasma), pointer :: s integer, intent(in) :: isize ! vector size ! all array arguments' sizes must be .ge.isize ! the target points whose distances from a surface are to be ! calculated: real*8, dimension(:), intent(in) :: zR ! R values real*8, dimension(:), intent(in) :: zZ ! Z values integer, intent(out) :: ier ! completion code 0=OK !------------ ! optional arguments integer, intent(in), optional :: maptype ! = 1 = xp_imap_newton = default -- iterative Newton root finder ! with initial guess by searching method, accurate w.r.t. ! spline fit of R & Z, to machine precision, but costly... ! = 2 = xp_imap_polar -- polar "circle-fit" distance estimate: ! The boundary is treated as a piecewise linear interpolation ! between circles fit to nearby points on the boundary ! much faster, not as accurate as Newton map ! = 3 = xp_imap_rzlinear -- bilinear inverse map -- after ! initialization (using polar map), the fastest, but less ! accurate still (depends on R,Z grid resolution). This ! method requires a user defined (R,Z) rectangle grid; if ! none is available the polar map is used instead. real*8, dimension(:), intent(in), optional :: phi_in ! toroidal angle ! (not needed for axisymmetric equilibria) real*8, intent(in), optional :: rho_in ! **SCALAR** -- the surface from ! which to calculate distances. DEFAULT is the boundary surface rho=1.0 logical, intent(in), optional :: ccw_phi ! Phi CCW view from above, default=T logical, intent(in), optional :: ccw_theta ! Theta CCW in poloidal plane, default=T integer, intent(in), optional :: i2pi ! angle coord normalization options ! i2pi=1 (default) returned value in range [0,2pi] ! i2pi=2 returned value in range [-pi,pi] logical, intent(in), optional :: outside_only ! flag that caller only ! expects to need data from beyond plasma boundary (default .FALSE.) ! (this switch only useful for optimization of maptype=3 calls). !------------- ! results-- optional but at least one should be present real*8, dimension(:), intent(out), optional :: theta_out ! theta of ! nearest point on surface real*8, dimension(:), intent(out), optional :: phi_out ! phi of ! nearest point on surface real*8, dimension(:), intent(out), optional :: dist ! distance to ! surface. If dist(j) > 0.0 then (zR(j),zZ(j)) is outside the test ! surface; if dist(j) < 0.0 then the point is inside the test surface; ! if dist(j) = 0.0 then the point is on the test flux surface. A caveat: when the fastest (bilinear) map is used: the distance information is good to within the resolution of the map, but... the location-on-surface (theta_out) data can be inaccurate in the vicinity of discontinuities in the location of the nearest point. These discontinuities occur when passing through a plane where two widely separate nearest points are equidistant. The bilinear interpolation of the fast map data then yields a theta value somewhere between the two nearly equidistant points, which is not correct. The solution, when the accuracy of the angle coordinate of the nearest point on the surface is of concern, is to use the slower but more accurate mapping options. 3 Volume_and_Area When xplasma is loaded with an equilibrium, two spline functions Volume(rho) -- volume, m**3, enclosed by rho surface, 0 <= rho <= 1. Area(rho) -- area, m**2, of poloidal cross section enclosed by rho surface, 0 <= rho <= 1. These are splines defined on the rho grid of the equilibrium. Two convenience routines are for evaluating these functions: xplasma_volume and xplasma_area. These are interface names that map to two actual routines, one which just returns the total plasma enclosed volume (or area), the other accepting vector of rho values for which either the volume (area) or its derivative w.r.t. rho may be returned. This may best be illustrated by example: real*8, parameter :: rho_bdy = 1.0d0 real*8 :: plasma_volume, plasma_area real*8 :: my_rho_grid(5) = (/ 0.0d0, 0.25d0, 0.5d0, 0.75d0, 1.0d0 /) real*8 :: my_vols(5),my_areas(5) call xplasma_volume(s,plasma_volume, ierr) call xplasma_volume(s,my_rho_grid,my_vols, ierr) call xplasma_area(s,plasma_area, ierr) call xplasma_area(s,my_rho_grid,my_areas, ierr) The variant with vector arguments accepts an optional argument with which the derivative of the volume are area may be retrieved: interface xplasma_volume module procedure xplasma_vol1 module procedure xplasma_voln end interface subroutine xplasma_voln(s,xs,ans,ier, ideriv) ! ** public ** ! return Volume(x) or dVol/dx -- calculate if needed. type (xplasma), pointer :: s real*8, dimension(:), intent(in) :: xs ! evaluation vector real*8, dimension(:), intent(out) :: ans ! result of evaluation integer, intent(out) :: ier ! completion code 0=OK integer, intent(in), optional :: ideriv ! derivative option 0,1,2 (Here "x" is used for "rho" -- radial flux coordinate). The accuracies of volumes and areas are to 64 bit machine precision for the given equilibrium representation (see next subtopic). 3 Flux_surface_integration There are many flux surface averaged metric and/or field quantities that may be needed in a plasma simulation, beyond area and volume. So, xplasma includes a highly efficient and accurate facility for computing flux surface integrated or averaged metrics on any user specified grid. To do this, a program using the xplasma library needs write access (see subtopic). The steps involved are: a) enable write access to xplasma (xplasma_author_set). b) define integration grid & have xplasma initialize data structures to support integration. c) have xplasma compute each desired integration, one at a time. If integrand data caching is enabled, the efficiency is greatly improved, when multiple integrals or averages need to be computed. d) clean up temporary data structures e) release write access (xplasma_author_clear). 4 Algorithm_and_accuracy The integration method used is a 10 point gauss-kronrod-patterson non- adaptive method, segmented to evaluate a separate 10 point sum between each pair of grid points on the equilibrium grid. This means that that the 10 point summation is applied to a sequence of C-infinity cubic spline segments, or, algebraic sums, products, and quotients of cubic spline segments for the geometry and magnetic fields as interpolated over the equilibrium grid. For non-singular integrands, this yields accuracy to 64 bit machine precision, as has been verified by comparing to higher-N gauss-kronrod-patterson summations. The code to carry this out is derived from www.netlib.org quadpack/dqng.for. Of the standard integrations, only the "neoclassical" flux surface averages "S" and "S" involve integrand singularities. These are of the form sqrt(1-B/Bmax). At present, for the sake of convenience and speed, these are evaluated using the same non-adaptive method with a simple patch to the singularity, resulting in accuracy of about 1 part in 10**4 compared to a slow, precise method. This is thought to be sufficient for the purposes of these quantities. User defined spline profiles f(rho,theta) can also be integrated or averaged, but the machine precision accuracy will only be obtained if the user splines are defined over the equilibrium grid-- i.e. the same grid as used by "R", R(rho,theta), "Z", Z(rho,theta), and the magnetic field splines BMOD, BZ, etc. The accuracies claimed are with respect to these spline representations of the equilibrium geometry and fields. 4 Write_access_to_xplasma Using the flux surface integration facility requires writing working data to be written into xplasma-- so, it is necessary to introduce routines to enable this. By default, xplasma is opened "write locked"-- dictionary items may not be added. To open xplasma to allow new items to be written, it is necessary to call the following pair of routines: character*32 :: author_name = call xplasma_author_set(s,author_name,ierr) [... execute code requiring write access to xplasma ...] call xplasma_author_clear(s,author_name,ierr) The code actually maintains a stack of author names-- i.e. the code executed could include subroutine calls to another library which writes xplasma data under a different authorname, with its own pair of xplasma_author_set and xplasma_author_clear calls. (If a routine fails to issue the trailing xplasma_author_clear call, the next routine that does issue such a call will see an error, because the author_name passed will not match the author_name at the top of the author stack). 4 Define_integration_grid The user's integration grid spans a subspace of rho [0,1] containing at least one grid zone. The entire [0,1] space need not be covered. The subdivided space is defined by passing the rho zone boundaries to the setup routine: real*8 :: my_rho_grid(5) = (/ 0.0d0, 0.25d0, 0.5d0, 0.75d0, 1.0d0 /) integer :: id_integ ! integrator object ID call xplasma_create_integ(s,'my_integrator',my_rho_grid,id_integ,ierr, & cache_enable = .TRUE.) Subsequent integrations using this grid will refer to the integer ID "id_integ" to reference the integrator. The integrator defined is now ready to compute surface averages at each of the grid points (my_rho_grid(1:5)), or, zonal averages covering the grid zones my_rho_grid(1:2), my_rho_grid(2:3), etc. For many applications, the 1d integrator is all that is needed. However, it is also possible to set up a 2d grid, with an additional subdivision of space by the poloidal angle coordinate: integer ::id_integ_2d ! 2d integrator object ID real*8 :: my_theta_grida(n) my_theta_grid(1) = 0.0d0 ... my_theta_grid(n) = 2*pi ! can also use [-pi,pi] ! or perhaps not span the entire theta coordinate range. call xplasma_augment_integ(s,'my_2d_integrator',id_integ,my_theta_grid, & id_integ_2d,ierr) Note that now two integrators are defined, with IDs: id_integ -- 1d integrator id_integ_2d -- 2d integrator. 4 Execute integrations The typical call to evaluate a flux surface average looks like this: real*8 :: answer(nn) ! nn the grid size specified earlier to real*8 :: ans_zon(nn-1) ! "xplasma_creat_integ"... real*8 :: ans_2d(nth,nn) ! nth the theta grid size specified earlier real*8 :: ans_2dzon(nth,nn-1) ! to "xplasma_augment_integ"... ! get the differential volume averaged <1/R^2> on my flux surfaces... call xplasma_rho_zonint(s,id_integ,'<1/R^2>S',answer,ierr) ! one can also evaluate call xplasma_rho_zonint(s,id_integ,'<1/R^2>',ans_zon,ierr) ! BUT it is much slower and usually the differential average is enough ! for modeling applications... ! <1/R^2> = zonal integral[dVol*(1/R**2)]/integral[dVol] over entire zone ! <1/R^2>S = surface integral[dVol*(1/R**2)]/integral[dVol] ! = surface integral[dVol*(1/R**2)]/(dVol/drho at surface) ! dVol = 2*pi*R*[(dl/dtheta)/grad(rho)]*dtheta*drho ! (available from the equilibrium metric tensor) ! integrations can also be done separately for each theta zone. The ! 1st dimension of ans_2d or ans_2dzon must match the number of theta ! zones: ! get volume of each zone in a 2d partitioning of flux space: call xplasma_2d_zonint(s,id_integ_2d,'dVOL',ans_2dzon,ierr) Note that the 3rd argument (character string) identifies the integration to be computed. The following integrals and averages are available. The intregral name passed argument must match one of these exactly, after uppercase conversion: surface-oriented zone-oriented remark "VOL" enclosed volume, m**3 "dVOL" volume in each zone "AREA" enclosed cross sectional area, m**2 "dAREA" area in each zone "DVDRHO" dVol/drho at each surface, m**3 "LPOL" poloidal path length, m "SURF" area of flux surface, m**2 "ITOR" enclosed toroidal current, A "<1/R^2>S" "<1/R^2>" m**-2 "<1/R>S" "<1/R>" m**-1 "<1/R^3>S" "<1/R^3>" m**-2 "S" "" m**2 "S" "" m "<1/B^2>S" "<1/B^2>" T**-2 "<1/B>S" "<1/B>" T**-1 "S" "" T**2 "S" "" T "S" "" T**2 "S" "" m**-1 "S" "" m**-2 "S" "" m**-4 "S" "" m**-5 "S" "" dimensionless "<1/(R*grad(RHO))>S" "<1/(R*grad(RHO))>" dimensionless "S" "" m**-2*T**-2 the following have sqrt singularities; the result is less accurate than with continuously differentiable integrands: "S" <(1-sqrt(1-h))(1+h/2)/(h*h)>; h=B/Bmax "S" <(1/B**2)*sqrt(1-h)*(2+h)/3>; h=B/Bmax if "USER_FUNC" is the name of a profile vs. (rho,theta) existing in the current xplasma: "S" "" surface or zone averaged USER_FUNC "I(USER_FUNC)S" "I(USER_FUNC)" surface or zone integral of USER_FUNC Any profile of (rho,theta) can be integrated or averaged in this way. 5 xplasma_rho_zonint subroutine xplasma_rho_zonint(s,idi,int_name,result,ier, & dvol_out,iwarn_dvol,dvdrho_out,iwarn_dvdrho) ! integrator -- over rho zones result of form f(rho) indexed by zone ! 1 result per zone (#zones) ! or rho surfaces, also result of form f(rho) ! 1 result persurface (#zones + 1) type (xplasma), pointer :: s integer, intent(in) :: idi ! integrator data structure ("black box") id character*(*),intent(in) :: int_name ! name of desired integration real*8, dimension(:), intent(out) :: result integer, intent(out) :: ier ! status code 0=OK ! ** optional outputs ** real*8, dimension(:), optional :: dvol_out ! zone volumes integer, intent(out), optional :: iwarn_dvol ! =0 if zone volumes OK ! iwarn_dvol=1: zone volumes not available; =2: array dimension error. real*8, dimension(:), optional :: dvdrho_out ! dV/drho @ surfaces integer, intent(out), optional :: iwarn_dvdrho ! =0 if dV/drho OK ! iwarn_dvdrho=1: dV/drho not available; =2: array dimension error. 5 xplasma_2d_zonint subroutine xplasma_2d_zonint(s,idi,int_name,result,ier, & dvol_out,iwarn_dvol,dvdrho_out,iwarn_dvdrho) ! integrator -- 2d domain; return integrations or averages ! over theta^rho zones or across theta zones at rho surfaces. type (xplasma), pointer :: s integer, intent(in) :: idi ! integrator data structure ("black box") id character*(*),intent(in) :: int_name ! name of desired integration real*8, dimension(:,:), intent(out) :: result integer, intent(out) :: ier ! status code 0=OK ! result is dimension [#theta-zones,#rho-zones] or ! [#theta-zones,#rho-surfaces] according as the integrand ! specified is zone oriented or surface oriented. ! The #zones are defined when the integrator dataset (idi) ! is defined. ! ** optional outputs ** real*8, dimension(:,:), optional :: dvol_out ! zone volumes integer, intent(out), optional :: iwarn_dvol ! =0 if zone volumes OK ! iwarn_dvol=1: zone volumes not available; =2: array dimension error. real*8, dimension(:,:), optional :: dvdrho_out ! dV/drho @ surfaces integer, intent(out), optional :: iwarn_dvdrho ! =0 if dV/drho OK ! iwarn_dvdrho=1: dV/drho not available; =2: array dimension error. ! Options for integrands is the same as for xplasma_rho_zonint; ! Surface oriented results e.g. LPOL are on rho surfaces at each theta ! zone; ! Zone orented results e.g. dVol are on rho^theta zones. 4 Cleanup After completion of use of the integrator or integrators, it is a good idea to clean up, i.e. free the memory and temporary data associated with the integrators. call xplasma_remove_item(s,id_integ,ierr) ! remove integrator (rho) ...and if integrations segmented in theta space were done: call xplasma_remove_item(s,id_integ_2d,ierr) ! remove integrator (rho,theta) (In general, "xplasma_remove_item" will report an error in case of an attempt to remove a non-existent item, or an item not owned by the current author). Finally, as previously mentioned, it is important to close the "write" connection to xplasma: call xplasma_author_clear(s,author_name,ierr) (author_name the same character string as used in a prior xplasma_author_set call). 3 Poloidal_flux_grid For codes that use a grid based on poloidal flux, it can be useful to get a mapping to rho=sqrt(toroidal_flux/toroidal_flux_at_bdy). The following xplasma subroutine provides a precise mapping: subroutine xplasma_rhopsi_find(s,psivals,rhovals,ierr, tol, iwarn) ! find rho values corresponding to a specified set of Psi values ! Psi -- Poloidal flux, Wb/rad ! rho -- sqrt(Tor_flux/Tor_flux_at_bdy) ! all input Psi values should be in the range [Psimin,Psimax] where ! Psimin corresponds to the magnetic axis and Psimax-Psimin ! corresponds to the (unsigned) poloidal flux, Wb/rad, enclosed ! within the core plasma. Usually Psimin=0 is set. type (xplasma), pointer :: s real*8, dimension(:), intent(in) :: Psivals ! Psi values in any order real*8, dimension(:), intent(out) :: rhovals ! rho values output ! sizes of Psivals and rhovals must match integer, intent(out) :: ierr ! status code, =0 on normal exit ! error occurs if xplasma is unitialized or contains no MHD equilibrium; ! Psi-out-of-range is handled (see iwarn, below). real*8, intent(in), optional :: tol ! accuracy tolerance ! on output, rho values satisfy ! abs(Psi(rhovals(i))-Psivals(i)) <= tol*[Psi at bdy] ! (1:npsi) -- sqrt(toroidal_flux/toroidal_flux_at_bdy) ! 0 on axis, 1 at the edge. integer, intent(out), optional :: iwarn ! #of Psi values out of range ! if Psi <= Psimin, rho=0 is returned; if Psi >= Psimax rho=1 is ! returned. Caution-- note that the input Psi values are in physical units, Wb/rad, not normalized to 1.000 at the boundary. But the following conventions are used: Psi=0 at the magnetic axis Psi is "unsigned", always non-negative; rho > 0 ==> Psi'(rho) > 0. The range of currently known Psi values within the plasma can be had with the following call: real*8 psimin,psimax call xplasma_psi_range(s,psimin,psimax) The sign of the poloidal flux can be had by fetching jphi_ccw, as in: integer :: jphi_ccw call xplasma_global_info(s,ierr, jphi_ccw=jphi_ccw) The convention is: jphi_ccw=1 means that toroidal current flows counter- clockwise in the tokamak as viewed from above. 3 Find_data_items Acquiring xplasma data generally involves two steps: a) translate name to integer ID (look up in dictionary) b) access the data, using the ID. This section describes step (a). If data of a known type is expected to exist under a known name, the following calls are the easiest to use: character*32 :: name integer :: id name = call xplasma_profId(s,name,id) ! look up a PROFILE --or-- call xplasma_gridId(s,name,id) ! look up a GRID --or-- call xplasma_coordId(s,name,id) ! look up a PROFILE --or-- call xplasma_listId(s,name,id) ! look up a LIST --or-- call xplasma_blkbxId(s,name,id) ! look up a BLACK BOX Note there is no "ierr" code. If the item does not exist, or if it is of the wrong type, or if there is some other error (e.g. "s" is a NULL pointer), id=0 is returned. 4 Alternate_method: ! this returns the ID of an item, no matter what its type is. ! setting nf_noerr=.TRUE. causes non-existence of the item not to ! be considered an error; this is an optional argument. call xplasma_find_item(s,name,id,ierr,nf_noerr=.TRUE.) ! itype is an integer... call xplasma_get_item_info(s,id,ierr, itype=itype) if(itype.eq.xplasma_profType) then else if(itype.eq.xplasma_gridType) then ...etc... 4 xplasma_common_ids The following call is used mainly within xplasma to retrieve the IDs of commonly used items; it is publicly accessible: subroutine xplasma_common_ids(s,ier, & id_g,id_psi,id_P,id_R,id_Z,id_Bmod,id_BR,id_BZ, & id_axis,id_RZminmax) ! return commonly requested profile IDs ! g(rho), psi(rho), R(rho,theta), Z(rho,theta) type (xplasma), pointer :: s integer, intent(out) :: ier integer, intent(out), optional :: id_g integer, intent(out), optional :: id_psi integer, intent(out), optional :: id_P integer, intent(out), optional :: id_R integer, intent(out), optional :: id_Z integer, intent(out), optional :: id_BR integer, intent(out), optional :: id_BZ integer, intent(out), optional :: id_BMOD integer, intent(out), optional :: id_axis integer, intent(out), optional :: id_RZminmax 3 Information_on_items Use "xplasma_get_item_info" to retrieve one or more of the following: character*32 :: name character*32 :: author character*(*) :: label character*(*) :: units integer :: itype ! type of item (list, profile, grid, etc.) All the return arguments are OPTIONAL and keyworded outputs. An optional input argument nf_noerr can be set .TRUE. to ignore "invalid ID" errors, i.e. not to create a message or set an error code in this case. Example: character*32 :: my_name,my_author call xplasma_get_item_info(s,id,ierr, name=my_name, author=my_author) Full interface: subroutine xplasma_get_item_info(s,id,ier, & nf_noerr, & name,label,units,itype,author) ! return info on item-- what info determined by presence of ! optional arguments type (xplasma), pointer :: s integer, intent(in) :: id ! id of item for which info is requested integer, intent(out) :: ier ! return status code (0=OK) ! optional control: logical, intent(in), optional :: nf_noerr ! .TRUE. to allow ivalid id ! here is what can be returned... character*(*), intent(out), optional :: name ! name of item character*(*), intent(out), optional :: label ! item label character*(*), intent(out), optional :: units ! item units label integer, intent(out), optional :: itype ! item type character*(*), intent(out), optional :: author ! name of author 3 Tables_of_contents It is possible to retrieve an entire table of contents of any category of xplasma data type. Since the interface uses optional arguments, the user's code decides what type of table of contents list is to be fetched. In addition, there is special treatment of lists containing elements that identify profiles (list element name and integer value match name and id of profile). The following routine fetches the lengths of the ID lists for tables of contents: subroutine xplasma_num_items(s,ier, & author_only, & num_lists, num_plists, num_coords, num_grids, num_profs, num_blkbxs) ! return the total numbers of items of various types in xplasma ! container object. type (xplasma), pointer :: s integer, intent(out) :: ier ! error code -- only set if s not initialized character*(*), intent(in), optional :: author_only ! set this to restrict the list to items owned by the specified author integer, intent(out), optional :: num_lists ! no. of lists integer, intent(out), optional :: num_plists ! no. lists refering ! to profiles integer, intent(out), optional :: num_coords ! no. of coordinates integer, intent(out), optional :: num_grids ! no. of grids integer, intent(out), optional :: num_profs ! no. of profiles integer, intent(out), optional :: num_blkbxs ! no. of black boxes The following routine fetches the actual lists of ids of the indicated type. Note that ierr=60 is returned if a passed integer array meant to receive the contents list is too small... subroutine xplasma_contents(s,ier, & author_only, & id_lists, id_plists, id_coords, id_grids, id_profs, id_blkbxs) ! return sorted list(s) of ids of objects of indicated type(s). ! lists of lists, lists of profiles, etc. type (xplasma), pointer :: s integer, intent(out) :: ier ! error code, 0=OK character*(*), intent(in), optional :: author_only ! set this to restrict the list to items owned by the specified author integer, dimension(:), intent(out), optional :: id_lists ! lists integer, dimension(:), intent(out), optional :: id_plists ! lists that ! refer to profiles integer, dimension(:), intent(out), optional :: id_coords ! coordinates integer, dimension(:), intent(out), optional :: id_grids ! grids integer, dimension(:), intent(out), optional :: id_profs ! profiles integer, dimension(:), intent(out), optional :: id_blkbxs ! black boxes The returned contents lists are sorted in alphabetic order of the names corresponding to each item ID. 3 Lists The LIST is one of the data types supported within xplasma. Each list has a name, and contains a named set of elements. Each element in addition to its name can have an integer, a floating point value, and a character string associated with it. Element names must be unique within the list but are not part of and do not conflict with the xplasma namespace-- indeed, one of the uses of lists is to gather the names and IDs of related xplasma data items. Some example applications of lists: -- storage of commonly referenced data, such as the plasma magnetic axis location; -- storage of scalar data, such as plasma ion species lists with the Z and A of each ion species specified. -- storage of lists of names and IDs of useful sets of xplasma profiles. If a list name is known, use the following call to get its integer ID: call xplasma_listId(s,name,id) ! look up a LIST Methods for accessing list data are described in the subtopics. 4 xplasma_list_info This routine optionally returns the size of a list and labeling information for the list (the list itself, not the member elements): subroutine xplasma_list_info(s,id,ier, & nelems,name,label,units,author) ! ! get the size of an existing list ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item integer, intent(out) :: ier ! completion code integer, intent(out), optional :: nelems ! no. of list elements in list character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author 4 xplasma_getList_size To just get the size of a list: subroutine xplasma_getList_size(s,id,nelems,ier) ! ! get the size of an existing list ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item integer, intent(out) :: nelems ! no. of list elements in list integer, intent(out) :: ier ! completion code 4 xplasma_getList_names To get the names of list elements, and, optionally, element data and/or global labels for lists, this routine can be used: subroutine xplasma_getList_names(s,id,enames,ier, & name,label,units,author,ivals,r8vals,chvals) ! ! get the element names from an existing list ! optional arguments allow additional data to be fetched as well. ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item character*(*), dimension(:), intent(out) :: enames ! the names... integer, intent(out) :: ier ! completion code character*(*), intent(out), optional :: name ! name (of whole list) character*(*), intent(out), optional :: label ! label (for whole list) character*(*), intent(out), optional :: units ! units label (whole list) character*(*), intent(out), optional :: author ! author integer, intent(out), dimension(:), optional :: ivals ! integer values real*8, intent(out), dimension(:), optional :: r8vals ! real*8 values character*(*), intent(out), dimension(:), optional :: chvals ! str vals 4 xplasma_getList_ivals subroutine xplasma_getList_ivals(s,id,ivals,ier) ! ! get integer values in list (r4 precision) ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item integer, dimension(:), intent(out) :: ivals ! the integer values integer, intent(out) :: ier ! completion code 4 xplasma_getList_r8vals subroutine xplasma_getList_r8vals(s,id,r8vals,ier) ! ! get real*8 values in list ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item real*8, dimension(:), intent(out) :: r8vals ! the real*8 values integer, intent(out) :: ier ! completion code 4 xplasma_getList_r4vals subroutine xplasma_getList_r4vals(s,id,r4vals,ier) ! ! get floating point values in list (r4 precision) ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item real, dimension(:), intent(out) :: r4vals ! the floating point values integer, intent(out) :: ier ! completion code 4 xplasma_getList_chvals subroutine xplasma_getList_chvals(s,id,chvals,ier) ! ! get character string values in list ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of list item character*(*), dimension(:), intent(out) :: chvals ! the character string values integer, intent(out) :: ier ! completion code 3 Coordinates The COORDINATE is one of the data types supported within xplasma. It contains little actual data-- it exists mainly as a mechanism to define the semantics of grids, with the following associations: (a) profiles are defined over grids; (b) grids are associated with coordinates. Each specific grid is associated with a single coordinate. Multiple grids can be associated with the same coordinate. Interpolation of profiles is organized by means of coordinates. Like all xplasma data items, coordinates have names, labels, physical units, and authors. Additional attributes of coordinates are: (a) whether or not a coordinate is periodic; (b) the minimum and maximum value (range) of a coordinate. When an xplasma object is initialized, several coordinates are predefined; their IDs are public parameters supplied by xplasma f90 modules: integer, parameter, public :: xplasma_rho_coord = 1 ! radial flux coord. integer, parameter, public :: xplasma_theta_coord = 2 ! poloidal angle coord. integer, parameter, public :: xplasma_phi_coord = 3 ! toroidal angle coord. integer, parameter, public :: xplasma_R_coord = 4 ! R coord integer, parameter, public :: xplasma_Z_coord = 5 ! Z coord integer, parameter, public :: xplasma_rhox_coord =6 ! "rho" extrapolation integer, parameter, public :: xplasma_thx_coord = 7 ! "theta" extrapolation integer, parameter, public :: xplasma_B_coord = 8 ! mod(B) grid integer, parameter, public :: xplasma_vpv_coord = 9 ! vpll/v coord. These parameters are frequently used, e.g. to identify the coordinates of input arguments in multivariate profile interpolation calls. The above list is sufficient for some applications. However, codes may define additional coordinates. For example, in the NUBEAM fast ion Monte Carlo model, a separate energy coordinate is defined for each fast ion species distribution function, because, the energy range desired for distribution functions varies from species to species, e.g. fusion product ion species vs. beam ion species. Several routines exist for making inquiries about coordinates, as shown in the subtopics. 4 xplasma_coord_info subroutine xplasma_coord_info(s,id,ier, & ngrids,periodic,xmin,xmax,name,label,units,author) type (xplasma), pointer :: s integer, intent(in) :: id ! coordinate object id integer, intent(out) :: ier ! completion code, 0=OK integer, intent(out), optional :: ngrids ! # of grid discretizations logical, intent(out), optional :: periodic ! periodicity flag real*8, intent(out), optional :: xmin ! min value real*8, intent(out), optional :: xmax ! max value character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author 4 xplasma_coord_gridId It is possible to retrieve the IDs of specific grids discretizing a coordinate. This is probably not how it will be done in applications, which are more likely to find grids by starting from a profile ID. Nevertheless, this public interface is supplied. The value of the input argument "indx" should be between 1 and "ngrids", the upper limit being the number of grids known to discretize the coordinate (see the xplasma_coord_info call). subroutine xplasma_coord_gridId(s,id,indx,idgrid,ier) ! get grid id from coordinate grid list index type (xplasma), pointer :: s integer, intent(in) :: id ! coordinate object id integer, intent(in) :: indx ! which grid is desired integer, intent(out) :: idgrid ! grid id of desired grid integer, intent(out) :: ier ! completion code, 0=OK 4 xplasma_coord_isPeriodic This public interface is provided; the periodicity attribute can also be obtained through an optional argument in the xplasma_coord_info call. subroutine xplasma_coord_isPeriodic(s,id,periodic,ier) ! get information whether coordinate object is a periodic angle coord. type (xplasma), pointer :: s integer, intent(in) :: id ! coordinate object id logical, intent(out) :: periodic ! T if this is a periodic coordinate integer, intent(out) :: ier ! completion code, 0=OK 3 Grids The GRID is one of the data types supported within xplasma. Each grid is a strict ascending sequence of numbers covering a specific range, and is associated with and discretizes a specific COORDINATE. Some coordinates have a prescribed range-- e.g. the normalized radial flux coordinate xplasma_rho_coord always covers the range 0 (the magnetic axis) to 1 (the plasma boundary). Grids over periodic coordinates cover a range of either [0,2pi] or [-pi,pi]. Some coordinates have their ranges defined by the first discretizing grid-- e.g. xplasma_R_coord and xplasma_Z_coord-- the dimensions of the [R,Z] box defining a computational region that will vary from tokamak to tokamak. With the exception of periodic coordinates, all grids discretizing a given coordinate must cover exactly the same range-- i.e. the first and last elements of the sequence of numbers defining the grid are constrained by the coordinate discretized by the grid. If the name of a grid is known and its ID is needed, the following call can be used: call xplasma_gridId(s,name,id) ! look up a GRID 4 xplasma_grid_size Get the size of a grid... subroutine xplasma_grid_size(s,id,nx,ier) ! get grid size type (xplasma), pointer :: s integer, intent(in) :: id ! grid object id integer, intent(out) :: nx ! grid size integer, intent(out) :: ier ! completion code, 0=OK 4 xplasma_grid Get the grid itself. The array "x" must precisely match the size of the grid. subroutine xplasma_grid(s,id,x,ier, ccwflag) ! get grid type (xplasma), pointer :: s integer, intent(in) :: id ! grid id object real*8, intent(out) :: x(:) ! the grid returned integer, intent(out) :: ier ! completion code, 0=OK logical, intent(in), optional :: ccwflag ! orientation flag ! default value is .TRUE., signifying counter-clockwise orientation If the grid is periodic, ccwflag=.FALSE. causes the grid to be returned in reverse order, corresponding e.g. to a clockwise drawn angle coordinate. The internal grid is always stored with clockwise orientation. The ccwflag=.FALSE. transform of a periodic grid theta(...) with ntheta points covering range [-pi:pi] is: theta_reverse(j) = -pi + (pi - theta(ntheta+1-j)) for j=1:ntheta 4 xplasma_grid_info This routine returns information on a grid and its associated coordinate. Other than the xplasma object pointer, the ID, and return status code, all arguments are optional, so, the caller can select the information desired. subroutine xplasma_grid_info(s,id,ier, & size,xmin,xmax,perio,coord,nrefs,name,label,units,author) ! get grid reference count (# of profile items using this grid) type (xplasma), pointer :: s integer, intent(in) :: id ! grid id object integer, intent(out) :: ier ! completion code, 0=OK integer, intent(out), optional :: size ! size of grid real*8, intent(out), optional :: xmin ! minimum grid value real*8, intent(out), optional :: xmax ! maximum grid value logical, intent(out), optional :: perio ! T for periodic grids integer, intent(out), optional :: coord ! coordinate associated w/grid integer, intent(out), optional :: nrefs ! no. of references to grid character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author 3 Profiles The PROFILE is one of the data types supported within xplasma. Indeed, xplasma was largely invented to support the storage, retrieval, and interpolation of profiles used in plasma simulations. Profiles are defined over grids (and hence over coordinates as well). For purposes of interpolation, it is generally enough to know the coordinates over which a profile is defined, but, it is of course possible to retrieve the underlying grids, and the original un- interpolated data, if desired. Profiles can be 1d, 2d, or 3d-- i.e. defined over 1, 2, or 3 grids (coordinates). The interpolation method for a profile is defined by the creator of the profile. The possibilities are the methods supported by the PSPLINE NTCC module-- piecewise linear (continuous), cubic Hermite (continuous and continuously once differentiable), and cubic spline (continuous and continuously twice differentiable along any coordinate). In addition, XPLASMA allows zonal step function data to be defined. For more information on interpolation, see the PSPLINE NTCC web page, http://w3.pppl.gov/NTCC/PSPLINE. There is a generic interface for profile interpolation: xplasma_eval_prof. Using this interface, retrieval of profile information basically comes down to making an appropriate xplasma_eval_prof call. This method will be described in detail. If the name of a profile is known, the following xplasma call retrieves its ID: call xplasma_profId(s,name,id) ! look up a PROFILE 4 xplasma_prof_info The following call allows access to detailed information on a given profile, as identified by its ID. Many of the arguments are self- explanatory, but further information on those that might not be is given below. subroutine xplasma_prof_info(s,id,ier, & rank,splineType,gridId1,gridId2,gridId3,profId1,gridType,counter, & name,label,units,author) ! base argauments in/out type (xplasma), pointer :: s integer, intent(in) :: id ! profile object id integer, intent(out) :: ier ! completion code, 0=OK !--------------------------------------------------- ! select desired information with optional arguments... integer, intent(out), optional :: rank ! rank integer, intent(out), optional :: splineType ! -1 for step function; 0 for piecewise linear, ! 1 for C1 Hermite, 2 for C2 cubic spline integer, intent(out), optional :: gridId1 ! id of 1st grid integer, intent(out), optional :: gridId2 ! id of 2nd grid (or zero) integer, intent(out), optional :: gridId3 ! id of 3rd grid (or zero) integer, intent(out), optional :: profId1 ! id of associated profile ! (sometimes f(R,Z) has an associated f(rho,theta) profile available, ! and vice versa). integer, intent(out), optional :: gridType ! generic grid type code ! 0=unknown or mixed ! 1=flux coords (rho) or (rho,theta) [phi someday] ! 2=cylindric coordinates (R,Z) [phi someday] integer, intent(out), optional :: counter ! profile counter ! each new version of this profile gets a new counter value. character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author Notes on argumentes: rank -- the dimensionality of the profile: 1 for 1d, 2 for 2d, 3 for 3d. profId1 -- ID for quantity defined over "alternate" grids. For some profiles there are two versions-- one defined over flux coordinates and one defined over [R,Z] coordinates. If for the input profile ID an alternate version exists, the alternate version ID can be returned with this argument. When interpolations are performed, if the coordinates supplied are not appropriate for a given profile ID, the profile's alternate ID, if available, will be tried before giving up. gridType -- this can be used to determine if the profile is defined over flux coordinates or cylindrical coordinates. The following public parameters are available for reference: integer, parameter, public :: xplasma_flux_coords = 1 ! generic flux coord. integer, parameter, public :: xplasma_cyl_coords = 2 ! generic cyl. coord. 4 xplasma_prof_gridInfo An alternate method of acquiring a profile's grid information is shown here; xplasma_prof_info can also be used. subroutine xplasma_prof_gridInfo(s,id,icoord,idgrid,ierr, & id_actual) type (xplasma), pointer :: s integer, intent(in) :: id ! id of profile about which info is sought integer, intent(in) :: icoord ! coordinate sought integer, intent(out) :: idgrid ! grid over this coordinate, or zero ! zero means profile identified by "id" is not a function of the ! specified coordinate; this is not considered an error. integer, intent(out) :: ierr ! error status code returned, 0=OK ! errors: profile id is invalid; coordinate id is invalid... integer, intent(out), optional :: id_actual ! "actual" profile id ! explanation: for some quantities there are two representations, ! e.g. Bphi(R,Z) and Bphi(rho,theta). The passed id can refer to ! either of these; the coordinate sought is found for the appropriate ! representation. Thus, if id points to Bphi(R,Z) and info on a rho ! grid is sought, id_actual (if present) will be set to the id of ! Bphi(rho,theta). 4 xplasma_eval_prof xplasma_eval_prof is a generic interface that can be used in any of the following modes: (1) interpolate 1d, 2d, or 3d profiles to given target coordinates. (2) interpolate 1d, 2d, or 3d profiles using previously calculated lookup information. (3) retrieve the original profile data (no interpolation). Each type of use is considered in the subtopics. 5 standard_interpolation Interpolations are supported on 1d, 2d, and 3d profiles, and, in a single call it is possible to perform interpolations on 1 profile or multiple profiles, at 1 target location or at a vector of target locations. Function values or derivatives can be returned. Angle coordinates may optionally be reversed with transformations theta -> 2pi - theta. When the rank of the profiles being evaluated exceeds one, the coordinate arguments can be given in any order, with grid or coordinate ID information given to indicate which is which. Generally for non-periodic coordinates it is an error to provide target coordinate data that is out of bounds; this behavior can be overridden with an optional argument "force_bounds". The number of out of bounds target points can be returned in the optional output argument "n_out_of_bounds". The relative tolerance for out-of-bounds tests is "bdytol", available with the call: call xplasma_global_info(s,ierr, bdytol=bdytol); interpolation target coordinates that exceed the coordinate minimum xmin or maximum xmax by more than bdytol*[xmax-xmin] are considered out of bounds. The result "ans" returned is a scalar if there is a single profile id and a scalar target point. It is a vector if there is either a list of profile ids or a vector of target points. If there is both a list of profile ids and a vector of target points, ans is a 2d array with size(ans,1) matching the target vector dimension, and size(ans,2) matching the number of profile ids; ans(:,1) will be the results of the evaluation of the 1st profile; ans(2,:) will be the results for all profiles at the 2nd target coordinate location. Here are some example calls: integer :: id1,id2(2) real*8 :: x,xvec(3),yvec(3) real*8 :: ans00,ans2(2),ans3(3),ans2d(3,2) ! 1 id, 1 target, eval 1st derivative call xplasma_eval_prof(s,id1,x,ans00,ier, ideriv1=1) ! 2 ids, 1 target, eval profile value call xplasma_eval_prof(s,id2,x,ans2,ier) ! 1 id, vector target, eval profile values call xplasma_eval_prof(s,id,xvec,ans3,ier) ! 2 ids, vector target, eval value of 1st profile, derivative of 2nd call xplasma_eval_prof(s,id2,xvec,ans32,ier, ideriv1s = (/0,1/)) ! 2d profile evaluation (no derivatives; xvec is "rho"; yvec is "theta"): call xplasma_eval_prof(s,id2, & xplasma_rho_coord,xvec, xplasma_theta_coord,yvec, ans32, ier) When the 2d or higher evaluation routine is called, it is possible for some of the profiles evaluated to be 1d, provided that the coordinate scalar or vector argument needed for evaluation of the 1d profile is present. What follows is the full interface for the calls mapped to when multiple ids and a vector target is provided to xplasma_eval_prof. For 1d profile evaluation: subroutine xplasma_eval_1dprofsxs(s,ids,xs,ans,ier, & ideriv1,ideriv1s,ccwflag1,force_bounds,n_out_of_bounds) ! evaluate multiple profiles f(x) 1d at a vector of points type (xplasma), pointer :: s integer, dimension(:), intent(in) :: ids ! profile ids real*8, dimension(:), intent(in) :: xs ! evaluation vector real*8, dimension(:,:), intent(out) :: ans ! result of evaluations integer, intent(out) :: ier ! completion code 0=OK ! ans(1:size(xs),1) -- result of eval of profile ids(1) ! ans(1:size(xs),2) -- result of eval of profile ids(2) ...etc... ! default: 0= fcn value integer, intent(in), optional :: ideriv1 ! derivative control ! 1 for df/dx; 2 for d2f/dx2; 3 for d3f/dx3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv1s ! as ideriv1 but specified separately for each profile ! default: T -- CCW orientation, no CW->CCW transform logical, intent(in), optional :: ccwflag1 ! default: F -- T to force all points in bounds with min & max logical, intent(in), optional :: force_bounds ! optional output: return no. of target points out of range integer, intent(out), optional :: n_out_of_bounds .............. For 2d profile evaluation: subroutine xplasma_eval_2dprofsxs(s,ids,idx1,x1s,idx2,x2s,ans,ier, & ideriv1,ideriv1s,ideriv2,ideriv2s, & ccwflag1,ccwflag2,force_bounds,n_out_of_bounds) ! evaluate multiple profiles f(x1,x2) 2d at a vector of points type (xplasma), pointer :: s integer, dimension(:), intent(in) :: ids ! profile ids integer :: idx1 ! id of x1 grid or coord. real*8, dimension(:), intent(in) :: x1s ! evaluation vector x1 integer :: idx2 ! id of x2 grid or coord. real*8, dimension(:), intent(in) :: x2s ! evaluation vector x2 real*8, dimension(:,:), intent(out) :: ans ! result of evaluations integer, intent(out) :: ier ! completion code 0=OK ! size(x1s)=size(x2s)=size(ans,1) expected... ! ans(1:size(x1s),1) -- result of eval of profile ids(1) ! ans(1:size(x1s),2) -- result of eval of profile ids(2) ...etc... ! default: 0= fcn value integer, intent(in), optional :: ideriv1 ! derivative control ! 1 for df/d[x1]; 2 for d2f/d[x1]2; 3 for d3f/d[x1]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv1s ! as ideriv1, 1 per evaluation profile id ! default: 0= fcn value integer, intent(in), optional :: ideriv2 ! derivative control ! 1 for df/d[x2]; 2 for d2f/d[x2]2; 3 for d3f/d[x2]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv2s ! as ideriv2, 1 per evaluation profile id ! default: T -- CCW orientation, no CW->CCW transform logical, intent(in), optional :: ccwflag1 ! for x1 logical, intent(in), optional :: ccwflag2 ! for x2 ! (these are ignored except for theta and phi angle coordinates ! default: F -- T to force all points in bounds with min & max logical, intent(in), optional :: force_bounds ! optional output: return no. of target points out of range integer, intent(out), optional :: n_out_of_bounds For 3d profile evaluation: subroutine xplasma_eval_3dprofsxs(s,ids,idx1,x1s,idx2,x2s,idx3,x3s, & ans,ier, & ideriv1,ideriv1s,ideriv2,ideriv2s,ideriv3,ideriv3s, & ccwflag1,ccwflag2,ccwflag3,force_bounds,n_out_of_bounds) ! evaluate multiple profiles f(x1,x2) 2d at a vector of points type (xplasma), pointer :: s integer, dimension(:), intent(in) :: ids ! profile ids integer :: idx1 ! id of x1 grid or coord. real*8, dimension(:), intent(in) :: x1s ! evaluation vector x1 integer :: idx2 ! id of x2 grid or coord. real*8, dimension(:), intent(in) :: x2s ! evaluation vector x2 integer :: idx3 ! id of x3 grid or coord. real*8, dimension(:), intent(in) :: x3s ! evaluation vector x3 real*8, dimension(:,:), intent(out) :: ans ! result of evaluations integer, intent(out) :: ier ! completion code 0=OK ! size(x1s)=size(x2s)=size(x3s)=size(ans,1) expected... ! ans(1:size(x1s),1) -- result of eval of profile ids(1) ! ans(1:size(x1s),2) -- result of eval of profile ids(2) ...etc... ! default: 0= fcn value integer, intent(in), optional :: ideriv1 ! derivative control ! 1 for df/d[x1]; 2 for d2f/d[x1]2; 3 for d3f/d[x1]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv1s ! as ideriv1, but separately specified for each evaluation ! default: 0= fcn value integer, intent(in), optional :: ideriv2 ! derivative control ! 1 for df/d[x2]; 2 for d2f/d[x2]2; 3 for d3f/d[x2]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv2s ! as ideriv2, but separately specified for each evaluation ! default: 0= fcn value integer, intent(in), optional :: ideriv3 ! derivative control ! 1 for df/d[x3]; 2 for d2f/d[x3]2; 3 for d3f/d[x3]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv3s ! as ideriv3, but separately specified for each evaluation ! default: T -- CCW orientation, no CW->CCW transform logical, intent(in), optional :: ccwflag1 ! for x1 logical, intent(in), optional :: ccwflag2 ! for x2 logical, intent(in), optional :: ccwflag3 ! for x3 ! (these are ignored except for theta and phi angle coordinates ! default: F -- T to force all points in bounds with min & max logical, intent(in), optional :: force_bounds ! optional output: return no. of target points out of range integer, intent(out), optional :: n_out_of_bounds 5 interpolation_using_prior_lookup Every interpolation involves two steps: zone lookup, and evaluation. Especially in the case of profiles defined over unevenly spaced grids, the zone lookup is a significant portion of the computational cost of the interpolation. Sometimes it is possible, and worth the effort, to pre-evaluate zone lookup and reuse this information through several interpolation calls. Xplasma provides a way to do this. In order for the user application to take advantage of this capability, it needs to declare special data structures, given public definitions in the xplasma modules, to store the results of a zone lookup calculation in such a way that it can be passed on to an xplasma interpolation routine built for this purpose. Such a declaration will look something like this: type (xpeval) :: x1_intrp,x2_intrp,x3_intrp Then, a generic interface xplasma_x_lookup can be used which will accept evaluation of a lookup either for a single scalar value or for a vector of values-- see the subtopic for information on the lookup pre-evaluation routines "xplasma_x_lookup". Once the lookups have been evaluated, xplasma_eval_prof can be called, supplying this information, to perform actual interpolations. Each xpeval item contains the lookup information either for a single point on a single coordinate, or a vector of points on a single coordinate. The evaluation routines expect to do a vector evaluation; at least one of the xpeval items must contain a vector of n lookup results. Those with only scalar lookup data will be replicated n times on evaluation. xplasma_eval_prof will make the following call for interpolation of a 1d profile: subroutine xplasma_xpeval_1dprof(s,id,xinfo,ans,ier, ideriv) ! evaluate a single 1d profile function f(x) at a vector of ! target points for which lookup has already been executed ! and stored in "xinfo". ! the x axis grid id of the profile and of the xinfo must match ! the size of the output and the no. of points in the xinfo must ! match. ! xinfo was set up by a prior xplasma_x_lookup call. type(xplasma), pointer :: s integer, intent(in) :: id ! 1d function to evaluate type (xpeval) :: xinfo ! prepared interpolation information real*8, dimension(:), intent(out) :: ans ! interpolation result integer, intent(out) :: ier ! status code, 0=OK integer, intent(in), optional :: ideriv ! derivative selection The following call is generated for interpolation of a 2d profile: subroutine xplasma_xpeval_2dprof(s,id,xinfo1,xinfo2,ans,ier, & ideriv1, ideriv2) type(xplasma), pointer :: s integer, intent(in) :: id ! 2d function to evaluate type (xpeval) :: xinfo1,xinfo2 ! prepared interpolation information sets real*8, dimension(:), intent(out) :: ans ! interpolation result integer, intent(out) :: ier ! status code, 0=OK integer, intent(in), optional :: ideriv1 ! derivative selection, x1 integer, intent(in), optional :: ideriv2 ! derivative selection, x2 And the following call is generated for interpolation of a 3d profile: subroutine xplasma_xpeval_3dprof(s,id,xinfo1,xinfo2,xinfo3,ans,ier, & ideriv1, ideriv2, ideriv3) type(xplasma), pointer :: s integer, intent(in) :: id ! 3d function to evaluate type (xpeval) :: xinfo1,xinfo2,xinfo3 ! prepared interpolation information sets real*8, dimension(:), intent(out) :: ans ! interpolation result integer, intent(out) :: ier ! status code, 0=OK integer, intent(in), optional :: ideriv1 ! derivative selection, x1 integer, intent(in), optional :: ideriv2 ! derivative selection, x2 integer, intent(in), optional :: ideriv3 ! derivative selection, x3 In all of these calls, the optional arguments default to zero-- meaning, interpolate the function value, not a derivative. ideriv*=1 means to evaluate the 1st derivative along the indicated coordinate. ideriv*=2 means to evaluate the 2nd derivative along the indicated coordinate. With these methods one interpolates one profile at a time. 6 lookup_routines These routines must be called to load the xpeval structures BEFORE calling one of the xplasma_xpeval routines (via xplasma_eval_prof): When done it is important to FREE MEMORY associated with the xpeval structure, called "xinfo" in the interfaces below. Failure to free memory after use will result in memory leaks (pointers are involved). !------------------------------------------------------ subroutine xpeval_free(xinfo) type(xpeval) :: xinfo ! deallocate & clean up xinfo !------------------------------------------------------ interface xplasma_x_lookup module procedure xplasma_x_lookup1 module procedure xplasma_x_lookupn end interface subroutine xplasma_x_lookup1(s,idx,x,xinfo,ier, & ccwflag,force_bounds,n_out_of_bounds) ! do x lookup -- scalar x value type (xplasma), pointer :: s integer, intent(in) :: idx ! x grid (grid) id real*8, intent(in) :: x ! single scalar x value type (xpeval) :: xinfo ! *** lookup results integer, intent(out) :: ier ! status code, 0=OK ! default: T -- CCW orientation, no CW->CCW transform logical, intent(in), optional :: ccwflag ! default: F -- T to force all points in bounds with min & max logical, intent(in), optional :: force_bounds ! optional output: return no. of target points out of range integer, intent(out), optional :: n_out_of_bounds ............... subroutine xplasma_x_lookupn(s,idx,x,xinfo,ier, & ccwflag,force_bounds,n_out_of_bounds) ! do x lookup -- scalar x value type (xplasma), pointer :: s integer, intent(in) :: idx ! x grid (grid) id real*8, intent(in), dimension(:) :: x ! vector of x values type (xpeval) :: xinfo ! *** lookup results integer, intent(out) :: ier ! status code, 0=OK ! default: T -- CCW orientation, no CW->CCW transform logical, intent(in), optional :: ccwflag ! default: F -- T to force all points in bounds with min & max logical, intent(in), optional :: force_bounds ! optional output: return no. of target points out of range integer, intent(out), optional :: n_out_of_bounds ............... Some further notes on arguments: idx -- this must be a grid ID, not a coordinate ID. x -- if scalar, the generic interface maps to xplasma_x_lookup1; if vector, the generic interface maps to xplasma_x_lookupn. xinfo -- dummy argument name for xpeval object, where the results of the interpolation lookup are stored. optional arguments: ccwflag -- if idx points to a grid over an angle coordinate, e.g. the poloidal angle flux coordinate, setting this flag to .FALSE. will cause the input data to be transformed by the formula theta -> 2pi - theta which would be appropriate if the caller's code wants to treat the poloidal angle coordinate as drawn clockwise around the flux surface (inside xplasma, the internal representation is always drawn counterclockwise). force_bounds -- force all input arguments to be in range using min and max function calls. This has no effect on periodic coordinate arguments, which always brought in range with 2*n*pi shifts as needed. n_out_of_bounds (out) -- number of x values not in bounds. Normally it is an error for out of bounds x values to be provided. 5 data_retrieval The original profile data can be retrieved using xplasma_eval_prof with the arguments shown for the following non-generic interfaces. In all cases, the arrays receiving the data must EXACTLY match the dimensioning of the profile being retrieved. The dimensioning can be reconstructed (i.e. use xplasma_prof_info to retrieve rank and grid IDs; use xplasma_grid_size to get sizes of individual grids, which are also the correct profile array dimensions). Error codes will be set if there is a mismatch in array dimension sizes or array rank; if the return code "ier" is set to a non-zero value, use: call xplasma_error(s,ier,6) to cause the error report to be printed on fortran unit 6 (stdout). Optional arguments "ccwflag*" are used to specify orientation of periodic coordinates. These have the same meaning as the corresponding arguments used when the profile is created. Setting ccwflag to .FALSE. will reverse the indexing order of the data along the th dimension, if that dimension corresponds to a periodic coordinate. If the th coordinate is not periodic, ccwflag is ignored. --------------- Retrieve 1d profile (this method in the xplasma_eval_prof generic interface): subroutine xplasma_getprof_1data(s,id,zdata,ier, & ccwflag1) ! return the original data associated with this profile type (xplasma), pointer :: s integer, intent(in) :: id ! profile id real*8, dimension(:), intent(out) :: zdata ! data returned integer, intent(out) :: ier ! completion code 0=OK logical, intent(in), optional :: ccwflag1 ! CCW flag (default T) Retrieve 2d profile (this method in the xplasma_eval_prof generic interface): subroutine xplasma_getprof_2data(s,id,zdata,ier, & ccwflag1,ccwflag2) ! return the original data associated with this profile type (xplasma), pointer :: s integer, intent(in) :: id ! profile id real*8, dimension(:,:), intent(out) :: zdata ! data returned integer, intent(out) :: ier ! completion code 0=OK logical, intent(in), optional :: ccwflag1 ! dim. 1 CCW flag (default T) logical, intent(in), optional :: ccwflag2 ! dim. 2 CCW flag (default T) Retrieve 3d profile (this method in the xplasma_eval_prof generic interface): subroutine xplasma_getprof_3data(s,id,zdata,ier, & ccwflag1,ccwflag2,ccwflag3) ! return the original data associated with this profile type (xplasma), pointer :: s integer, intent(in) :: id ! profile id real*8, dimension(:,:,:), intent(out) :: zdata ! data returned integer, intent(out) :: ier ! completion code 0=OK logical, intent(in), optional :: ccwflag1 ! dim. 1 CCW flag (default T) logical, intent(in), optional :: ccwflag2 ! dim. 2 CCW flag (default T) logical, intent(in), optional :: ccwflag3 ! dim. 3 CCW flag (default T) 4 xplasma_RZeval_2d This is a specialized routine for evaluating R(rho,theta) and/or Z(rho,theta), with automatic failover to extrapolated profiles of R and Z for rho > 1, when these are available. In general, if xplasma has scrap off layer information, i.e. limiter locations, it will, when updating the MHD equilibrium, create extrapolated (rho,theta) -> (R,Z) map which however is not differentiable across the rho=1 boundary. This routine gives access to both the original and extrapolated R and Z in the context of a single call. subroutine xplasma_RZeval_2d(s,ids,idx1,x1s,idx2,x2s,ans,ier, & ideriv1,ideriv1s,ideriv2,ideriv2s, & ccwflag1,ccwflag2,force_bounds,n_out_of_bounds) ! use xplasma_eval_2dprofsxs to ! evaluate multiple profiles f(x1,x2) 2d at a vector of points ! IF all the profiles are "R" and "Z", ! and IF there is a scrape-off region, ! split the evaluation so that points inside the plasma use ! the standard R & Z, and points outside use the extrapolated ! profiles. type (xplasma), pointer :: s integer, dimension(:), intent(in) :: ids ! profile ids integer :: idx1 ! id of x1 grid or coord. real*8, dimension(:), intent(in) :: x1s ! evaluation vector x1 integer :: idx2 ! id of x2 grid or coord. real*8, dimension(:), intent(in) :: x2s ! evaluation vector x2 real*8, dimension(:,:), intent(out) :: ans ! result of evaluations integer, intent(out) :: ier ! completion code 0=OK ! size(x1s)=size(x2s)=size(ans,1) expected... ! ans(1:size(x1s),1) -- result of eval of profile ids(1) ! ans(1:size(x1s),2) -- result of eval of profile ids(2) ...etc... ! default: 0= fcn value integer, intent(in), optional :: ideriv1 ! derivative control ! 1 for df/d[x1]; 2 for d2f/d[x1]2; 3 for d3f/d[x1]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv1s ! as ideriv1, 1 per evaluation profile id ! default: 0= fcn value integer, intent(in), optional :: ideriv2 ! derivative control ! 1 for df/d[x2]; 2 for d2f/d[x2]2; 3 for d3f/d[x2]3 ! 2 & 3 available only for C2 spline fits integer, intent(in), dimension(:), optional :: ideriv2s ! as ideriv2, 1 per evaluation profile id ! default: T -- CCW orientation, no CW->CCW transform logical, intent(in), optional :: ccwflag1 ! for x1 logical, intent(in), optional :: ccwflag2 ! for x2 ! (these are ignored except for theta and phi angle coordinates ! default: F -- T to force all points in bounds with min & max logical, intent(in), optional :: force_bounds ! optional output: return no. of target points out of range integer, intent(out), optional :: n_out_of_bounds 3 Generic_blackBoxes The BLACK BOX is one of the data types supported within xplasma. A black box is a storage mechanism for data not fitting into other types. It consists of the following: a) a user defined integer "type" code; b) a user specified 1d array or buffer of any size, of integers; c) a user specified 1d array or buffer of any size, of 64 bit (real*8) floating point numbers. The black box is used within xplasma for a number of purposes-- caches for flux surface integration, data to support rapid methods for evaluation of coordinate transformations, "MCgrid" data, ... Integer black box type codes in the range 0 to 101 are used by xplasma library software; it is recommended that user defined black boxes stay away from using type codes in this range. In order to avoid unnecesary and costly copying of large arrays, the black box data access allows pointers to the black box integer and floating point array data to be returned. This means however that the contents of black box datasets can be modified directly by assignment of (pointer) array elements, without going through an interface call. User caution is advised. 4 xplasma_blackBox_info The black box type code, buffer sizes, and generic xplasma item labeling can be retrieved with this call. Except for the status code, output arguments are optional, allowing the caller to control what information is to be returned. subroutine xplasma_blackBox_info(s,id,ier, & type,isize,r8size,name,label,units,author) ! ! get the size of an existing list ! type (xplasma), pointer :: s integer, intent(in) :: id ! id of black box item integer, intent(out) :: ier ! completion code integer, intent(out), optional :: type ! "type" of black box ! generally, meaning of iitype is user defined... integer, intent(out), optional :: isize ! number of integer words integer, intent(out), optional :: r8size ! number of floating pt words character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author 4 xplasma_blackBox_retrieve Black box data, or pointers to data, can be retrieved with this call. As usual, most output arguments are optional: subroutine xplasma_blackBox_retrieve(s,id,ier, & itype,iarray,r8array,ia_ptr,r8a_ptr,name,label,units,author) ! retrieve blackBox data ! optional arguments; user controls what data is wanted. type (xplasma), pointer :: s integer, intent(in) :: id ! id of black box item integer, intent(out) :: ier ! completion code, 0=OK integer, intent(out), optional :: itype ! "type" of black box ! meaning is user defined. ! if copies of arrays are wanted: integer, dimension(:), intent(out), optional :: iarray real*8, dimension(:), intent(out), optional :: r8array ! if pointers to arrays are wanted: integer, dimension(:), pointer, optional :: ia_ptr real*8, dimension(:), pointer, optional :: r8a_ptr ! optional labeling info... character*(*), intent(out), optional :: name ! item name character*(*), intent(out), optional :: label ! item label character*(*), intent(out), optional :: units ! item units label character*(*), intent(out), optional :: author ! author 3 Limiter_information Xplasma supports the notion of an axisymmetric "limiter" or vessel wall. The region between the limiter and the core plasma contains plasma and neutral gas on open field lines and is referred to here as the scrape off layer. For reasons related to legacy code, xplasma supports two methods for defining an axisymmetric limiter. Therefore, the routine for retrieving limiter information supports both representations. However, the routines "xplasma_bdycon" and "xplasma_limcon" support methods for returning the plasma boundary and limiter represented as a closed contour (to good approximation), regardless of the internal representation. The interface xplasma_lim_distance allows quick access to the distance from any point in space to the nearest limiter point. Specific call interfaces are in the subtopics. 4 xplasma_lim_info subroutine xplasma_lim_info(s,ier, & itype, npts,rpts,zpts, & dist, & nlines, Rl, Zl, thl, & ncircs, Rc, Zc, rad) ! retrieve limiter information type(xplasma), pointer :: s integer, intent(out) :: ier integer, intent(out), optional :: itype ! itype=100 means piecewise linear closed contour (axisymmetric limiter) ! itype=101 means list of circles & infinite lines integer, intent(out), optional :: npts ! #pts in contour (itype=100 only) real*8, dimension(:), intent(out), optional :: rpts,zpts ! the closed contour; rpts(npts)=rpts(1) and zpts(npts)=zpts(1) (m). real*8, intent(out), optional :: dist ! itype=101 limiters only... ! if .ge.0.0d0, the plasma boundary + dist is considered a limiter ! if .lt.0.0d0-- no plasma boundary based limiter. integer, intent(out), optional :: nlines ! #lines in itype=101 limiter real*8, dimension(:), intent(out), optional :: Rl,Zl,thl ! point (m) through which line passes, and line orientation in DEGREES ! e.g. 45.0 means up and to the right... integer, intent(out), optional :: ncircs ! #circles in itype=101 limiter real*8, dimension(:), intent(out), optional :: Rc,Zc,rad ! center location and radius of each circle (m) 4 xplasma_lim_RZminmax subroutine xplasma_lim_rzminmax(s,zRmin,zRmax,zZmin,zZmax,ier, & itype) ! get the Rmin,Rmax,Zmin,Zmax of the vacuum region ! enclosed by the mechanical limiter ! optional: get a limiter type code also. type (xplasma), pointer :: s real*8, intent(out) :: zRmin,zRmax,zZmin,zZmax ! R & Z min & max (m) integer, intent(out) :: ier integer, intent(out), optional :: itype ! limiter representation ! within xplasma: 100 for piecewise contour; 101 for "circles and ! lines" possibly augmented by a fixed distance from the plasma ! boundary. 4 xplasma_limcon subroutine xplasma_limcon(s,rlim,zlim,inum_got,ier, & tol,maptype) ! return a closed contour of length size(rlim)=size(zlim) or less ! describing the limiter. ! if tol.le.0.0d0 -- the length of the description is size(rlim) exactly ! if tol.gt.0.0d0 -- length may be shortened to remove colinear points. ! also, if the limiter was originally given in contour form, and, ! size(rlim) is sufficient, the original specification will be ! returned. type (xplasma), pointer :: s real*8, dimension(:), intent(out) :: rlim,zlim ! limiter contour ! the sizes of these vectors must be the same. ! unless the original specification is being returned, the code ! requires a minimum size of 30 points. integer, intent(out) :: inum_got ! contour length returned ! if inum_got.lt.size(rlim), rlim(inum_got+1:size(rlim))=0 on exit. ! similarly for zlim. integer, intent(out) :: ier ! status code, 0=OK real*8, intent(in), optional :: tol ! tolerance specification ! if zero or negative, no shortening of the limiter description is ! attempted. If .gt. 0.0, adjacent segments that are colinear to within ! tol*[Rmax] are replaced by a single segment. DEFAULT: s%bdytol integer, intent(in), optional :: maptype ! mapping option ! 1=slowest but most exact; 2=intermediate; 3=bilinear interpolation ! 2 or 3 recommended. DEFAULT: 2 4 xplasma_bdycon Return a closed contour describing the plasma boundary; length of contour taken from size of vector arguments provided: subroutine xplasma_bdycon(s,rbdy,zbdy,ier) ! return a closed contour around the plasma boundary, similar to ! xplasma_limcon... type (xplasma), pointer :: s real*8, dimension(:), intent(out) :: rbdy,zbdy ! bdy contour ! the sizes of these vectors must be the same. integer, intent(out) :: ier 4 xplasma_lim_distance interface xplasma_lim_distance module procedure xplasma_lim_dist1 module procedure xplasma_lim_distn end interface subroutine xplasma_lim_dist1(s,r,z,d,ier, maptype, rlim1,zlim1) ! scalar version of xplasma_lim_distn -- axisymmetry assumed. type (xplasma), pointer :: s real*8, intent(in) :: R,Z ! input location (m) real*8, intent(out) :: d ! output distance from nearest point ! on limiter or wall (m). integer, intent(out) :: ier integer, intent(in), optional :: maptype ! see xplasma_lim_distn real*8, intent(out), optional :: rlim1,zlim1 ! nearest limiter point end subroutine xplasma_lim_dist1 subroutine xplasma_lim_distn(s,rvec,zvec,dist,ier, & maptype,rlim,zlim) ! find distance from each element of a vector of (R,Z) pairs ! to the nearest point on a limiter. ! optionally return the location of that point on the limiter ! in (R,Z). ! the routine assumes axisymmetry. Optional phi arguments may ! be added someday. type (xplasma), pointer :: s real*8, intent(in), dimension(:) :: rvec,zvec ! (R,Z) input vector real*8, intent(out), dimension(:) :: dist ! distance values returned integer, intent(out) :: ier integer, intent(in), optional :: maptype ! distance map option ! if it is required to compute the distance from the (R,Z) point to ! the plasma boundary, this specifies the option (for xplasma_bdfind). ! Default is the circle-fit option (2). (1) is slower but more ! precise; (3) is less accurate but would be faster if a bilinear ! distance map has already been computed. real*8, intent(out), dimension(:), optional :: rlim,zlim ! locations ! of nearest contact points on limiter, for each ! input location. 3 MCgrid_Data The MCgrid, or "Monte Carlo grid", is an irregular 2d spatial grid that was developed originally for capture of 2d binned data in a Monte Carlo calculation (NUBEAM). It is used for fast ion distribution functions and for certain spatially 2d profiles output by NUBEAM-- such as, beam halo thermal neutral sources, beam-target and beam-beam fusion rates. The grid is aligned with flux coordinates. It is constructed of "zone rows" that are equally spaced in rho=sqrt(Psi_tor/Psi_tor_at_bdy); each zone row is subdivided into a different number of zones equispaced in the equilibrium poloidal angle coordinate theta, with fewer subdivisions for zone rows near the axis, and more for zone rows in the edge. The result is a set of zones all with roughly equal volume and cross sectional area, as is desirable for consistency of Monte Carlo summation statistics. The number of poloidal zones per zone row is linear in the zone row index: typically, 4 in the 1st zone row adjacent to the magnetic axis, 8 in the next row out, 12 in the next row, etc. There are two variants, according as the underlying MHD equilibrium is updown symmetric, or updown asymmetric. As internally stored, the first row layout is this for the updown symmetric variant: _____ / | \ / 2 | 1 \ |___|___| with theta=0 on the large major radius side; increasing to theta=pi on the large major radius side covering the upper half of the plasma cross section above the midplane. The next row out would have 4 zones, the next 6 zones, etc. And this for the updown asymmetric variant: _____ / | \ / 4 | 3 \ |___|___| | | | \ 1 | 2 / \__|__/ with theta=-pi on the lower branch on the small major radius side, theta=0 on the large major radius side, and theta=pi on the upper branch on the small major radius side. The next zone row contains 8 zones, the next 12 zones, and so on. Poloidal zones are stored contiguously, with the poloidal zone index increasing with increasing theta coordinate, theta being oriented counter-clockwise in the plasma cross section drawn to the right of the machine axis of symmetry. This is how the data is stored internally. Note, however, that on retrieval of MCgrid data, there are options to reverse the storage order to be consistent with a clockwise oriented poloidal angle coordinate (as preferred by some codes), and, for both updown symmetric and updown asymmetric data, the start point for theta zone indexing can be specified as either 0, -pi, or the default (which is 0 for updown symmetric data and -pi for updown asymmetric data). These choices are controlled by optional arguments in the MCgrid data access routines. The MCgrid itself is a named "black box" item in the xplasma dataset; in principle there could be multiple MCgrids defined, but in practice there is usually only one. Each profiles defined over an MCgrid is itself a named "black box" item. When a MCgrid is set up, the volume of each MCgrid zone is computed, and this is itself stored as a profile over the MCgrid; the name is the MCgrid name with the suffix "_DVOL" appended. These integer parameters define the black box type codes used for MCgrids and their associated profiles: integer, parameter, public :: xplasma_bbgtype=17 ! MC grid type integer, parameter, public :: xplasma_bbftype=18 ! MC profile type Although information on these grids and profile could be obtained by using generic black box data access calls, it is better to use the calls specifically provided to support MCgrid items, which are described in the subtopics. Calls for access to MCgrid data are shown here; calls to create MCgrid data are given under the section heading Modifying_F95_Xplasma_Data. 4 Finding_MCgrid_ID Usually there is only one MCgrid per dataset, so, the optional argument id_mcgrid1 can be used in the public interface shown below. If the MCgrid is known to have been created by a named author, the search can be constrained by that name: e.g "author_only='NUBEAM'". subroutine xplasma_mcgrid_find(s,ierr, author_only, & id_mcgrid1,num_mcgrids,id_mcgrids) ! find MCgrid grid definitions that are available in the current ! xplasma-- usually there will be only one we think. type(xplasma), pointer :: s integer, intent(out) :: ierr ! status code returned (0=OK) character*(*), intent(in), optional :: author_only ! restrict search to MCgrid objects written by this author or code integer, intent(out), optional :: id_mcgrid1 ! MCgrid id returned, if... ! if exactly one exists. If none exist, id_mcgrid1=0 is returned; if ! >1 exists,id_mcgrid1=- is returned. integer, intent(out), optional :: num_mcgrids ! number of MCgrids found. integer, dimension(:), optional :: id_mcgrids ! MCgrid ids returned. 4 MCgrid_info The following interface enables access to details on the MCgrid. It should be noted that for historical reasons the MCgrid supports an extension beyond the plasma boundary, but, known applications of the MCgrid are restricted to the plasma core at present. Other than the error code, the return arguments are optional, so, the caller can select just the information that is needed. subroutine xplasma_mcgrid_info(s,id_mcgrid,ierr, & nzrow,nzrow_ext, nzons,nzons_ext, nth0,nphi, nths, udsym, & name,label) ! fetch information on xplasma_mcgrid grid type(xplasma), pointer :: s integer, intent(in) :: id_mcgrid ! id of MCgrid object to be queried integer, intent(out):: ierr ! status code 0=OK integer, intent(out), optional :: nzrow ! # of zone rows in plasma integer, intent(out), optional :: nzrow_ext ! # inside and outside integer, intent(out), optional :: nzons ! # of zones in plasma integer, intent(out), optional :: nzons_ext ! # inside and outside integer, intent(out), optional :: nth0 ! # of zones @axis cover [0:pi] integer, intent(out), optional :: nphi ! # of phi zones (1:axisymmetry) integer, dimension(:), intent(out), optional :: nths ! number of theta zones in each zone row logical, intent(out), optional :: udsym character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: author ----------- It is important to notice that nth0 gives only the number of zones in the upper (or lower) halfplane of the first zone row, regardless of whether the underlying MCgrid is updown symmetric or updown asymmetric. Thus, an updown symmetric grid with 2 zones in the first row, and an updown asymmtric grid with 4 zones in the first row, would both return a value of 2 for nth0. 5 MCgrid_zone_indexing_example [A working example is in xplasma_debug_module.f90] integer :: id_mcgrid real*8 :: rho = 0.375d0 ! rho where index is wanted real*8 :: th = 0.2d0 ! theta angle where index is wanted real*8, parameter :: rhomax = 1.0d0 ! only looking at core plasma region... integer :: nrow,irow ! zone rows integer, dimension(:), allocatable :: nths,nthsum logical :: udsym ! updown symmetry flag integer :: indx call xplasma_find_mcgrid(s,ierr, id_mcgrid1=id_mcgrid) if(ierr.ne.0) [handle xplasma error] if(id_mcgrid.eq.0) [no MCgrid exists] if(id_mcgrid.lt.0) [multiple MCgrids exist] call xplasma_mcgrid_info(s,id_mcgrid,ierr, nzrow=nrow, udsym=udsym) allocate(nths(nrow),nthsum(0:nrow)) call xplasma_mcgrid_info(s,id_mcgrid,ierr, nths=nths) nthsum(0)=0 do irow=1,nrow nthsum(irow)=nthsum(irow-1)+nths(irow) enddo ! find row index using subroutine. The indexing is counting from -pi, ! regardless of the symmetry flag-- there is an implicit assumption that ! MCgrid profile data is fetched with this option!!! call mcindx(rho,th,udsym,nrows,rhomax,nths,nthsum,indx) !------------------------- subroutine mcindx(rho,th,udsym,nrow,rhomax,nths,nthsum,indx) ! private ! this routine assumes that the theta indexing is oriented ! counterclockwise and that it starts at -pi for each row, ! *even* for updown symmetric data. real*8, intent(in) :: rho,th ! flux coords in logical, intent(in) :: udsym ! updown symmetry flag integer, intent(in) :: nrow ! #rows real*8, intent(in) :: rhomax ! max rho in range integer, intent(in) :: nths(nrow) ! # theta zones / row integer, intent(in) :: nthsum(0:nrow) ! # in prior rows integer, intent(out) :: indx !---------------------------------------- integer :: irho,ith real*8 :: zth,zthlim real*8, parameter :: CPI = 3.1415926535897931D+00 real*8, parameter :: ZERO = 0.0d0 !---------------------------------------- indx = 0 if(udsym) then zth = -abs(th) zthlim = ZERO else zth = th zthlim = CPI endif irho = 1 + rho*nrow/rhomax if((irho.ge.1).and.(irho.le.nrow)) then ith = 1 + (zth+CPI)*nths(irho)/(zthlim+CPI) if((ith.ge.1).and.(ith.le.nths(irho))) then indx = nthsum(irho-1) + ith endif endif end subroutine mcindx 4 MCgrid_zone_volumes The following interface allows retrieval of the volumes of all the MCgrid individual zones: subroutine xplasma_mcgrid_volumes(s,id_mcgrid,vols,ierr, & lstart0_th,ccwflag_th) ! return the zone volumes associated with a MCgrid. Compute them ! if necessary. type(xplasma), pointer :: s integer, intent(in) :: id_mcgrid ! MC grid id real*8, dimension(:) :: vols ! the array of volume elements integer, intent(out) :: ierr ! status code (0=OK) ! theta indexing options logical, intent(in), optional :: lstart0_th ! default depends on symmetry logical, intent(in), optional :: ccwflag_th ! default: T ! lstart0_th=T means user's first theta zone lower bdy is 0; F means it ! is -pi. If defaulted: updown symmetric data assumed to start at 0; ! updown asymmetric data assumed to start at -pi. ! ccwflag_th=T means zone index increases as one goes along a row of ! zones at fixed radial index in a counter-clockwise direction; F means ! the opposite. T is the default. 4 List_of_MCgrid_profiles The following routine allows retrieval of a list of IDs of available profiles defined over a MCgrid. In typical usage, a first call might get the number of such profiles, followed by allocation of an integer array to hold the list of profile IDs, followed by a second call to actually fetch the IDs. subroutine xplasma_mcgrid_findProfs(s,id_mcgrid,ierr, & author_only, num_profs, id_profs) ! find profiles defined over given MCgrid type(xplasma), pointer :: s integer, intent(in) :: id_mcgrid ! MC grid id integer, intent(out) :: ierr ! status code returned (0=OK) character*(*), intent(in), optional :: author_only ! restrict search to MCgrid profiles written by this author or code integer, intent(out), optional :: num_profs ! number of profiles found. integer, dimension(:), optional :: id_profs ! profile ids returned. 4 Information_on_MCgrid_profile This routine returns information on an individual MCgrid profile. Except for the status code, all output arguments are optional, allowing the caller to select the desired information: subroutine xplasma_mcgrid_profInfo(s,id_mcprof,ierr, & id_mcgrid,nzons,irank,idim1,idim2, & gridId1,gridId2,normcode, & name,label,units,author) ! fetch information on a single profile defined over an MCgrid... type (xplasma), pointer :: s integer, intent(in) :: id_mcprof ! profile ID integer, intent(out) :: ierr ! status code returned: 0=OK integer, intent(out), optional :: id_mcgrid ! MCgrid ID integer, intent(out), optional :: nzons ! #spatial zones integer, intent(out), optional :: irank ! #non-spatial dimensions integer, intent(out), optional :: idim1,idim2 ! non-spatial dim. sizes ! grid IDs for non-spatial dimensions (if available) integer, intent(out), optional :: gridId1 integer, intent(out), optional :: gridId2 ! normalization code integer, intent(out), optional :: normcode !------------- character*(*), intent(out), optional :: name character*(*), intent(out), optional :: label character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author For multidimensional profiles with velocity space coordinates (i.e. distribution functions), the arguments gridId1 and gridId2 can be used to retrieve the IDs of the velocity space coordinates used. 4 Fetch_MCgrid_profile_data This routine retrieves the actual profile data into a user provided array. An array with correct dimensioning must be supplied by the user. The last dimension of the user supplied array is the spatial dimension, matching the MCgrid size-- the value of the "nzons" argument returned by xplasma_mcgrid_info(...), for profiles defined just over the core plasma region. subroutine xplasma_mcgrid_getobj(s,id_p,ierr, & lstart0_th,ccwflag_th, & data_1d,data_2d,data_3d,label,units) ! retrieve previously stored MCgrid profile ! (see xplasma_mcgrid_putobj) type (xplasma), pointer :: s integer, intent(in) :: id_p ! MCgrid profile ID integer, intent(out) :: ierr ! exit status code (0=OK) ! theta indexing options logical, intent(in), optional :: lstart0_th ! default depends on symmetry logical, intent(in), optional :: ccwflag_th ! default: T ! lstart0_th=T means user's first theta zone lower bdy is 0; F means it ! is -pi. If defaulted: updown symmetric data assumed to start at 0; ! updown asymmetric data assumed to start at -pi. ! ccwflag_th=T means zone index increases as one goes along a row of ! zones at fixed radial index in a counter-clockwise direction; F means ! the opposite. !----------------- ! one of the optional arguments data_1d, data_2d, data_3d must ! be provided to receive the data. ! the provided data array sizes must match the stored data. real*8, dimension(:), intent(out), optional :: data_1d real*8, dimension(:,:), intent(out), optional :: data_2d real*8, dimension(:,:,:), intent(out), optional :: data_3d !----------------- character*(*), intent(out), optional :: name character*(*), intent(out), optional :: labels character*(*), intent(out), optional :: units character*(*), intent(out), optional :: author 2 Initializing_F95_Xplasma_Data The previous sections discussed methods for accessing data in an xplasma object that was already created. This and the following sections describe how to build up an xplasma object from scratch, and how to evolve it in a time dependent simulation. The following pieces of an xplasma dataset normally do not evolve in time; rather, they are initialized once at the start of the lifetime of an xplasma object: (a) the flux coordinate grids "rho" and "theta" over which equilibrium profiles are defined: flux surface label "rho" = sqrt(toroidal_flux/toroidal_flux_at_bdy) poloidal angle "theta" -- in the internal xplasma representation, moving in the direction of +theta at fixed rho traces out a flux surface in the counter-clockwise direction, when the flux surface cross section is drawn to the right of the tokamak axis of symmetry (R=0). (b) the location of an axisymmetric limiter or vacuum vessel wall, usually specified as a closed contour sequence of [R,Z] pairs. (c) the R and Z grids covering a rectangle in space that encloses the limiters and will enclose the core plasma at all times. (d) additional coordinates and grids as may be defined by the user, e.g. in velocity space; multiple grids may discretize the same coordinates. The items that do evolve over the lifetime of an xplasma object are: (e) the equilibrium over flux coordinates R(theta,rho), Z(theta,rho) -- the geometry of flux surfaces (m). g(rho) = R*|B_phi|, a flux surface constant, m*T; B_phi = bphi_ccw*(g/R); bphi_ccw=+1 if toroidal field points counter-clockwise -1 if clockwise, viewed from above Psi(rho) = poloidal flux function, Wb/rad, zero on axis, dPsi/drho > 0, with relation to the poloidal field Bpol = -jphi_ccw*(1/R)*grad(Psi); jphi_ccw=+1 if toroidal current flows counter-clockwise -1 if clockwise, viewed from above P(rho) = a profile representing the pressure used in the Grad- Shafranov equation to produce the equilibrium The derived profile BR(theta,rho), BZ(theta,rho), BPHI(theta,rho), BMOD(theta,rho) give spline approximations to the field components. (f) the equilibrium over (R,Z) coordinates Psi(R,Z), Wb/rad, Bpol = -jphi_ccw*(1/R)*grad(Psi) BR_RZ(R,Z), BZ_RZ(R,Z), BPHI_RZ(R,Z), BMOD_RZ(R,Z) are provided as spline approximations to the field components over (R,Z) coordinates. (g) any collection of named data items-- * profiles-- zonal data, piecewise linear data, or splines over up to 3 independent coordinates * lists -- each containing its own collection of named, labeled elements with associated scalar integer and floating point values. * black boxes -- generic storage of "unstructured" integer and floating point data. The INITIALIZATION of an xplasma object consists of the definition of the elements that do not (normally) vary in time. The most direct method of initialization is to start from an EFIT G-eqdsk file containing a free boundary equilibrium reconstruction: use xplasma_obj_instance character*32 :: author_name = 'My_code' character*60 :: filename = call xoi_init(ier) ! init xplasma pointers in module call xplasma_author_set(s,author_name,ier) call xplasma_fromgeqdsk(s,filename,ier, & nrho = , & ntheta = ) ! also a good idea: real*8 :: curtime curtime = call xplasma_time_set(s,curtime,ier) call xplasma_author_clear(s,author_name,ier) The optional arguments specify the desired equilibrium flux coordinate grid resolutions; if defaulted, the code uses nrho=51, ntheta=101. Loading xplasma from a G-eqdsk file for the first time, items (a)-(c), are defined, and initial values for items (e)-(f) are also provided. The equilibrium can subsequently be modified with later calls to xplasma_fromgeqdsk, reading modified G-eqdsk files, but then the limiter and [R,Z] grids found must be the same as were given in the first G-eqdsk file read; this is checked and an error will be reported if there is a mismatch. Further initializations may involve the definition of additional coordinates and grids, e.g. in velocity space, which are application dependent. The xplasma public interfaces include the following routines likely to be used in initialization: xplasma_create_coordinate xplasma_create_grid which will be described below. 3 Alternate_initialization_method The alternative to using a G-eqdsk file is to use a series of calls to spell out each of the components of the equilibrium. In this context, initialization consists in the definition of the flux coordinate grids, the limiters, and the [R,Z] grids. Code to do this is sketched below. Details on the interfaces used will be given below. !-------------------------- ! on the FIRST call, the rho and theta grids need to be defined... character*32 :: author_name = 'My_code' call xplasma_author_set(s,author_name,ier) allocate(rho(inrho),th(inth)) call xplasma_create_grid(s,'my_rho',xplasma_rho_coordinate,rho,id_rho,ier, & label='my rho grid for the equilibrium') call xplasma_create_grid(s,'my_theta',xplasma_theta_coordinate,th,id_th, & ier, label='my theta grid for the equilibrium') ! the field SIGNS also have to be set, once, at the beginning of the ! lifetime of the xplasma object: ! bphi_ccw=1 means the toroidal field points counter clockwise; ! jphi_ccw=-1 means the toroidal current flows clockwise ! these orientations as viewed from above the tokamak ! values of +/-1 are allowed for each of these arguments call xplasma_field_ccw_set(s,ier, bphi_ccw=1, jphi_ccw=-1) ! example ! also on the first call, a LIMITER should be defined... allocate(rpts(nlim),zpts(nlim)) call xplasma_mklim_contour(s,rpts(1:nlim),zpts(1:nlim),iwarn,ier) ! and an [R,Z] rectangular space covering the plasma and limiter should be ! defined: allocate(rgrid(nrgrid),zgrid(nzgrid)) < i.e. rgrid(1) <= minval(rpts); rgrid(nrgrid) >= maxval(rpts) > < i.e. zgrid(1) <= minval(zpts); zgrid(nzgrid) >= maxval(zpts) > call xplasma_create_rzgrid(s,rgrid(1:nrgrid),zgrid(1:nzgrid), & id_rgrid,id_zgrid,ier) ! these arguments returned... ! done for now... call xplasma_author_clear(s,author_name,ier) 2 Modifying_F95_Xplasma_Data In this section, the new Fortran-95 Xplasma interface is described, with emphasis on methods for adding data to an existing dataset, or, creating one from scratch. The xplasma pointers (s,s1) defined in xplasma_obj_instance are presumed to be available in the coding examples. The argument "ierr" is always a return status code, 0=normal. For simplicity of examples, error checking / error handling code is usually omitted -- don't try this at home! 3 Code_authorship Since xplasma is envisioned as a means of communication between disparate code modules or even separately executing processes, it is useful to be able to trace the origin of xplasma data items to specific codes or modules. In practice, this means a code must "identify itself" to access xplasma. A code that writes into xplasma must provide an "author id", which is any alphanumeric name, not starting with a digit, up to 32 characters long. Failure to do this will usually result in an "xplasma locked" error report. Codes which write to xplasma should use the following pair of calls: character*32 author_name author_name = "my_code" ! or whatever name you choose... call xplasma_author_set(s,author_name,ier) ! push my name onto author stack call xplasma_author_clear(s,author_name,ier) ! pop my name off the stack An error could conceivably occur on the latter call, if an intervening call (after xplasma_author_set) was made to another module, which issued another xplasma_author_set but failed to issue a matching call to xplasma_author_clear; the software would then detect that the passed name does not match the name at the top of the author stack. To prevent such errors, it is important to always pair these calls. 3 Update_the_MHD_equilibrium FREE BOUNDARY MHD equilibria are best read into xplasma from a G-eqdsk file: call xplasma_author_set(s,author_name,ier) call xplasma_fromgeqdsk(s,filename,ier, [...optional arguments...]) ! also a good idea: real*8 :: curtime curtime = call xplasma_time_set(s,curtime,ier) call xplasma_author_clear(s,author_name,ier) The first time this is done, the resolution of the flux coordinate "rho" and "theta" grids can be specified-- note that only evenly spaced grids are allowed for equilibria built from G-eqdsk at present. If the grid size optional arguments "nrho" and "nth" are omitted, the default values nrho=51 and nth=101 are used-- sufficient for many purposes. On subsequent calls, the equilibrium flux coordinate grids are reused from the prior state. PRESCRIBED BOUNDARY MHD equilibria are updated through a sequence of calls which redefine the equilibrium profiles-- as described in the subtopic. 4 xplasma_fromgeqdsk The full xplasma_fromgeqdsk interface is copied below. Except possibly for specifying "nrho" and "nth" on the first call, most users should leave the optional arguments unspecified-- the default values are well tested and reliable. The argument "filename" can be either a unix file path to a read-accessible G-eqdsk file in the standard GA-documented EFIT format, or, it can be a path to an MDS+ tree (access to experimental database). MDS+ path strings take the following form: MDS+[/]:[:]:(;t=