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
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.
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.
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.
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.
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.
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.
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
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.
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.
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.
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.
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))].
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.
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.
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
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.
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.
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.
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.
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.
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
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.
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.
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.
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.
! 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.
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
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...
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".
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.
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
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=<desired relative accuracy tolerance>"; 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.
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
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.
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).
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).
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
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 "<F(H)NCLASS>S" and "<F(H)NC_TSC>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.
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 = <name-of-my-program-in-quotes>
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).
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.
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
"<R^2>S" "<R^2>" m**2
"<R>S" "<R>" m
"<1/B^2>S" "<1/B^2>" T**-2
"<1/B>S" "<1/B>" T**-1
"<B^2>S" "<B^2>" T**2
"<B>S" "<B>" T
"<BZ^2>S" "<BZ^2>" T**2
"<grad(RHO)>S" "<grad(RHO)>" m**-1
"<grad(RHO)^2>S" "<grad(RHO)^2>" m**-2
"<grad(RHO)^2/R^2>S" "<grad(RHO)^2/R^2>" m**-4
"<grad(RHO)^2/R^3>S" "<grad(RHO)^2/R^3>" m**-5
"<R^2*grad(RHO)^2>S" "<R^2*grad(RHO)^2>" dimensionless
"<1/(R*grad(RHO))>S" "<1/(R*grad(RHO))>" dimensionless
"<grad(RHO)^2/B^2>S" "<grad(RHO)^2/B^2>" m**-2*T**-2
the following have sqrt singularities; the result is less accurate than
with continuously differentiable integrands:
"<F(H)NCLASS>S" <(1-sqrt(1-h))(1+h/2)/(h*h)>; h=B/Bmax
"<F(H)NC_TSC>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:
"<USER_FUNC>S" "<USER_FUNC>" 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.
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.
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.
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).
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.
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 = <the-name-you-want>
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.
! 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
<it is a profile>
else if(itype.eq.xplasma_gridType) then
<it is a grid>
...etc...
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
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
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.
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.
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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.
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).
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.
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
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
i