grideir.txt
===========
This document describes the program that generates the Eirene grid for TJ-II,
and also gives the list of subroutines, and intended action, that could be used
for any device different from TJ-II (program grideir_dev) including a simple
example.
INDEX
-----
1. Introduction.
1.1 Program grideir for TJ-II.
1.2 Namelist "namsup".
1.3 Contents of the Eirene geometry file.
2.1 Program grideir_dev for any device.
2.2 Subroutines needed for the user.
2.3 Listing and comments for the subroutines needed.
1. Introduction
------------
The description of the Eirene geometry model adopted at Ciemat for the TJ-II
helical axis stellarator is discussed in the Power Point document "TJ-II Eirene
geometry", where many graphical ilustrations appear, here we will present a
written explanation.
Only a single semiperiod of the device (a toroidal angle span of 45 degrees in
TJ-II) is required to model the geometry because periodicity and stellarator
simmetry can complete the full 4 periods of the device.
The Eirene cells are disposed in nphi toroidal layers (typically the value of
nphi is around 40). These layers are all of them vertical and lie exactly on
toroidal planes. For the first toroidal layer the "anterior" cell walls are in
the plane Phi=0, while the "posterior" walls of the last layer are in the plane
Phi = 45 degrees. This means that each layer (for nphi=40) extends over a Phi
interval dphi (1.125 degrees). The midle plane of the first layer is at Phi =
dphi/2 (0.5625 degrees), etc. Those planes, for every sucesive layer, play an
important role in the model, in all this document they will always be called
the "middle toroidal plane" of the corresponding layer.
Now in every "middle toroidal plane" an equally spaced cartesian grid is
established with nhor intervals (nax=nhor+1 "Nodes") along both the horizontal
and vertical directions (a "straight" position). The center of this uniform
grid is always placed just on the magnetic axis position of the corresponding
"middle toroidal plane". As TJ-II has an helical axis with a very wide
elongation (around 28 cm for a major radius of 150 cm and plasma average radius
between 10 and 20 cm, depending of the configuration) this means that the
center of the horizontal-vertical cartesian grid (The "Grid" from now on)
describes a helical path in space (it rotates poloidaly 180 degrees,
clock-wise, along the 45 degrees of one semiperiod).
The number of intervals of the grid is always an odd number (this is enforced
by the program), so that the magnetic axis (the "Grid" origin) is always at the
center of the middle "Grid" interval. The cells constructed from these middle
intervals as explained below will also play an important role in the model and
>from now on they will be named "Axial" cells. These "Axial" cells enclose the
magnetic axis inside them. In other words: the axis is always at the center of
the "middle toroidal plane" of the "Axial" cells.
The cell contruction starts at the "Grid" vertices (the "Nodes") of each
"middle toroidal plane". At each Node a straight horizontal line is drawn
perpendicular to that plane and prolonged until it crosses the next toroidal
layer wall (an angular distance of dphi/2), at which intersection point it
ends. These intersection points will be the vertices of the cells. With this
construction the cells are hexahedra (6 face polyhedra with 8 vertices and 12
edges) whose walls (the faces) are planar quadrilaterals. Four of these
quadrilaterals are vertical parallelograms: the two "anterior" and "posterior"
walls lie on vertical toroidal planes and are square, equal, but not parallel;
the two lateral walls are rectangular with horizontal and vertical sides,
parallel but of different sizes. The remaining two walls (top and bottom) are
equal horizontal, symmetrical trapezia (quadrilaterals with two parallel
sides).
The construction makes it obvious that the set of cells is contiguous, with
neither holes nor overlappings, and that the calculation of the external normal
vectors to the walls (the "Normals", needed to calculate trajectory
intersections in Eirene) is very easy.
Up to now nothing has been said of the "Grid" size. This size is determined by
the variable rxlim of the subroutine g3d_getlim_gr in the g3d TJ-II local
library, which is defined, for any given configuration, as the maximum poloidal
distance in cm from the plasma border to the magnetic axis in all the torus,
increased by 2. In this way a set of squares with half side equal to rxlim
centered on the magnetic axis at each "middle toroidal plane" and with
horizontal and vertical sides ("straight" position) enclose the plasma
completely. These squares will be called "The Big Squares" from now on.
More generally, a possible supplementary margin (dxmarg) can be added to rxlim,
so that the total half width of the "Grid" is dmaxcm = rxlim + dxmarg. In
consequence the distance between the "Nodes" (in other words the horizontal and
vertical width of the cells) is dxcm = 2*dmaxcm/nhor, the same for all layers
and cells.
In spite of this, since then total number of cells is at a premium, not all the
cells that could be contained inside those "Big Squares" are retained. For each
toroidal layer only those cells are retained that are inside a rectangle
defined as the minimum size "straight" rectangle such that ALL the cells that
are OUTSIDE have a average plasma radius normalised to 1 at border (from now on
the "normalised plasma radius") greater than some given limit "slim". Moreover
that rectangle is extended on every side by "kmarg" cells, although never
trespassing the "Big Square" limits.
In this way, although all the cells are of the same horizontal and vertical
width, the number of cells along the horizontal and vertical directions do not
need to be the same, and this is also true for the different layers. In summary
the region ocupied by the cells is always a straight rectangle, but with
varying width and heigth from layer to layer.
In TJ-II dxmarg is usually 0.0 (the results given by the g3d library outside
that distance from the axis are not accurate but only extrapolations) while
slim = 1.2 anf kmarg = 5. Nevertheless all these values could be changed in the
namelist "namsup" ($1.2).
For example: for the 100_44_64 configuration nhor is 35, the maximum half width
dmaxcm (= rxlim in TJ-II) is 34.75 cm (the average plasma radius is 19.25 cm
for this configuration) and the resulting cell width (the side of the two
square walls) is 1.186 cm. These are typical values for most of the
configurations.
In this way these cells, called from now on the "Central" cells, include, with
certainty, the totality of the plasma, but not necessarily the totality of the
TJ-II Vacuum Vessel (VV), although its most important parts from the point of
view of neutral physics and orbit behaviour (the Hard Core in particular)
remain inside. This is one of the reasons to include margins (the 2 cm added in
rxlim and also kmarg) farther from the plasma, and why the visualisation of the
geometry is convenient (cf. $1.3).
Nevertheless the VV is totally enclosed by a 2nd set of cells (the "Vacuum"
cells). This set, only 8 in each layer, is contiguous (without overlapping)
with the extreme borders of the "Central" cell region and it is arranged as a
ring wrapping that region. They are constructed in the same way as the first
set, but this time the size is calculated so as to include the VV totally (but
see later), also with some margin (half a cell width, i.e., dxcm/2). This
produces hexahedral cells also, similar to the "Central" ones, but generally
larger and usually without any square sides, although the sides are still plane
quadrilaterals.
It is evident that a wide enough plasma margin could include ALL the VV, but we
consider that in TJ-II very far away parts of the VV (i.e. the ports) are not
too relevant for the physics of neutrals and the price to pay in number of
cells for the inclusion of these relatively useless parts should be too high,
or, alternatively, for a given total number of cells --this number of cells is
limited inside Eirene to 999999-- the "resolution" inside the "Central" cell
region (which is the essential part for physics) would be deteriorated. For
this reason we have established the possibility to clip off the most distant
parts of the VV, using the variable vvrad (see $1.2). In our model every point
of the VV that is at a distance from the TJ-II TF coil centers greater then
vvrad centimeters (usualy vvrad=50) is ignored and the point shifted to the
said distance. This eliminates the most extreme parts of the TJ-II ports,
somewhat rounding those remote and wild territories.
Once the two sets of cells have been constructed, physical magnitudes can be
attached to them: plasma normalised radius, magnetic field, drifts, and also
other data such as volume, etc. For each physical variable the choice has been
to assign to the full cell the value of the variable at its center of mass. In
particular a very important role is played by the normalised plasma radius,
which is defined for every cell. If the cell is inside the plasma (center and
all vertices inside the plasma) that value is the true normalised radius of its
center, and is ² 1. For cells partially or totally outside the plasma the
associated cell value is always greater than 1 and with a codification which
indicates the relative situation of its center and vertices with respect to the
VV, for example, how many vertices are outside the VV, etc (see the
"radii.list" file). These codified ficticious "normalised plasma radius" are
listed below:
Codification for the ficticious normalised plasma radius (s)
============================================================
"Central" cells (0.0 ² s ² 22.0):
s < 1.0 Inside plasma. All vertices inside VV.
1.0 ² s ² 2.99 Inside plasma, but some corners outside VV.
The "real" radius is the fractional part of s.
s > 2.99 Outside plasma:
s = 6.0 8 corners and center inside VV.
s = 11.0 Some corners outside, center inside.
s = 12.0 Some corners outside, center outside.
s = 21.0 8 corners outside, center inside (improbable ??).
s = 22.0 8 corners outside, center cutside.
"Vacuum" cells (101.0 ² s ² 822.0):
For these cells the first digit of s indicates the number of order
of the cell (from 1 to 8, left to right first, then bottom to top).
For the 1st "Vacuum" cell the code is:
s = 101.0 Some corner outside, center inside.
s = 102.0 Some corner outside, center outside.
s = 121.0 8 corners outside, center inside.
s = 122.0 8 corners outside, center outside.
And similar codes for the rest of these cells.
Cells outside the geometry limits (s > 822.0):
s = 1.e6 Cells totally outside the geometry.
--------------------------
Points placed outside the "Vacuum" cell extremes correspond the "External
Darkness" and are illegal cells.
Finally the VV has been modeled (including the clipping) by a series of
nphvv=nphi+1 toroidal cross sections, in coincidence with the toroidal planes
of the walls. In each cut, using the TF coil center as the origin (different
for each toroidal angle), a poloidal angle has been defined, starting at -180
degrees on the inside horizontal side, going counter-clock wise, until closing
the loop at the same point for +180. Then the vertex position at each of nths
points (usualy 91) is calculated. These toroidal and poloidal grids form a kind
of "mosaic" that covers all the VV. The facets of this mosaic are
quadrilaterals, sometimes very nonplanar and even often convex with respect to
the poloidal angle origin. These quadrilaterals have been bisected, in a
consistent way, in two triangles (which, of course, are planar) and the
external normal vectors of these triangles have been calculated. These vertex
positions and normal vectors allow the localisation of the intersection of
trajectories with the VV mosaic in a fast and effective way.
1.1 Program grideir for TJ-II
-------------------------
This is the program to generate the Eirene grid for TJ-II. It is located (as
all the other Eirene files) in the directory $ei4n, that is:
"/r2/fusion/guasp/EIRENE-2004_new" in the jen50,
"/disco02/fusion/guaspx/EIRENE-2004_new" in the fenix.
This program must be linked with the local g3d library of TJ-II (see script
"grideir"), which provides the needed data for TJ-II: magnetic axis position,
plasma radius at any given point, magnetic field vector and drifts at arbitrary
points, cross sections for the TJ-II Vacuum Vessel (VV), checking if a point is
inside the VV, etc.
The execution and compilation script "grideir" has a single argument, which is
the name of the configuration (the default is "100_44_64"). There is a separate
Eirene geometry file for each configuration, a binary file, rather big
(20054296 bytes for the geom_opt50.1oct that corresponds to the TJ-II 100_44_64
configuration). The program takes about 5 minutes to generate the geometry on
both Ciemat computers.
With this information and the data given in the namelist "namsup" (see $1.2)
the program starts looking, in the TJ-II configuration data base, for the
coresponding g3d file (these files are placed at "/r4/temporal/guasp/griddir"
with names such as 100_44_64.grd, etc.).
If that configuration file does not exist the program aborts. Once that file is
read (subroutine g3d_geom_nbi) the program finds, by calling the subroutine
g3d_getlim_gr, the main characteristics of the configuration i.e.: rotational
transform, plasma radius, maximal rectangular size of the plasma cross sections
(rxlim), etc. Some of these data are written to the Eirene binary file (unit
30).
Then the program calculates the magnetic axis points (subroutine g3d_findax)
and writes them into the file.
The next step is to loop along the toroidal angles, calculating the position of
the center of the hexahedra and obtaining the plasma normalised radius at these
centers (subroutine g3d_rad_gr) and then the 8 vertices of each hexahedron are
calculated in the way described at $1 (subroutine genhex). As commented before
($1) the horizontal and vertical extension of the cells are calculated by
finding the minimal "straight" rectangle such that all the cells that are
outside have a normalised radius greater than slim and then allowing a margin
of a few cell widths (kmarg) but without trespassing the "Big Square" limits.
Then the resulting cell volumes are calculated by the volum subroutine and the
magnetic field and drifts at the center of the cell are calculated with the
g3d_B_Field_cyl and g3d_drift subroutines.
All these values are written in the file, as well as the vertices and center
coordinates. At the same time several checks are performed: possible illegal
values, correct toroidal angles for the cell vertices, too small volumes, etc.,
etc.
A particularly important check is made in order to ensure that the
intersections of the magnetic axis with the anterior and posterior axial cell
faces remain inside them (by construction the center of the axial cells is
always at the magnetic axis). If the axis falls outside one of these faces, or
too near the borders, the program aborts and reccomends either to increase the
number of toroidal cells (nphi) or, alternatively, to decrease the maximum
number of horizontal and vertical cells (nhor).
When all the "Central" cells have been constructed and written the program
starts to construct the "Vacuum" cells, 8 of them for each toroidal layer, with
width and height variable so as to totally occupy the space that exists between
the limit of the "Central" cell region and the VV, extended by a small margin
(half a cell width).
To do this the toroidal cross sections of the TJ-II VV for each toroidal cut
along nths poloidal angles (subroutines intvvgr and getvac) are obtained and
stored, and then hexahedric cells similar to but bigger than the "Central" ones
are produced and checked.
Finally the VV vertices are written and VV triangles are consructed by
bisecting the "facets" formed by the VV vertex "mosaic", the external normal
vectors of the resulting triangles are calculated (gentri subroutine), checked
for possible null triangles and written.
With this the program stops. Any error detected is commented and produces a
premature abort.
1.2 Namelist "namsup".
-----------------
This namelist contains only 7 variables (and some other more for tests that are
not explained here) and determines the maximum number of horizontal, vertical
and toroidal cells. These variables are the following (Fortran convention for
data types is used):
namelist /namsup/ nhor, nphi, nths, dxmarg, slim, kmarg, vvrad
Where:
nhor: Number of maximum horizontal and vertical "Central" cells in each
---- toroidal layer. It must be an odd number (this is enforced by the
program) so that the center of the middle cells (the "Axial" cells)
falls exactly on the magnetic axis. It cannot be larger than 47
and the default value is 35.
nphi: Number of cell layers along the toroidal angle in a single
---- semiperiod (45 degrees for TJ-II). The anterior walls of the first
toroidal layer fall exactly at toroidal angle 0 degrees. The posterior
walls of the last layer fall exactly at 45 degrees (one semiperiod).
It cannot be larger than 51 and the default value is 40.
nths: Number of poloidal points taken for each toroidal cross section
---- of the TJ-II VV. The points of the cross section begin at poloidal
angle -180. degrees and finish at +180. (completely closing the loop
in the counterclockwise direction).
An odd value is recomended (but not forced) in order that the poloidal
angle 0. degrees be included. It cannot be larger than 91 and this is
also the default value.
The toroidal cuts for the VV are always in coincidence with the
cell layers. Therefore there are nphi+1 VV toroidal cuts.
The point used as the radial origin to measure the poloidal angle is
always the the TF coil center corresponding to each toroidal
plane (this is calculated by subroutine intvvgr).
dxmarg: The extension of the plasma margin that added to rxlim (cf.$1.)
------ gives the total maximal "Grid" width. In cm, default value = 0.
slim: The limit for the maximum normalised plasma radius.
---- For each toroidal layer all retained cells are in a straight
rectangle such that all the cells that are outside have
normalised plasma radius that is larger than slim. No units.
Default value 1.2 .
kmarg: The number of cells that extend all sides of the rectangle
----- calculated with slim. Default value 5. Nevertheless the extended
rectangle never trespass the "Big Square", that is the sides
never are greater than 2*(rxlim + dxmarg).
vvrad: Limit for the VV minor radius: for TJ-II it seems convenient to
----- "clip" the more external parts of the VV (ports, etc), so that
all VV points that fall at a distance, from the local coil center,
greater than vvrad are ignored, and their distance limited to vvrad.
The default value is 50. cm (remember that the TF coil radius is
40. cm). This choice assures that all the "useful" parts of the
VV are retained and only the most remote parts of the ports are
eliminated.
1.3 Contents of the Eirene geometry file.
------------------------------------
This file (fort.30) is binary and has been written (as usual in Eirene) with
the real variables as real*8 type and the integer ones as integer*4.
a) Headings
1st. Record: Name of the configuration (character*80)
-----------
2nd. Record: nper, nsemiper
----------- Where nper in the number of periods of the device (4 in TJ-II)
and nsemiper is always 1, since only one semiperiod of the device
is used. Periodicity and stellarator simmetry allow the code to
complete the remaining parts.
3rd. Record: nhor, nphi, nax, nd, nphvv, nths1, ndph
----------- Where nhor, nphi are the same as in the namelist ($1.2)
nax is nphi+1 (the number of toroidal cuts)
nd is nhor+1 (the maximum number of nodes along the horizontal
or vertical directions)
nphvv is nphi+1 (the number of toroidal cuts for the VV)
nths1 is nths-1 (the number of poloidal intervals in the VV)
ndph is always 1
It is evident that most of these values are redundant, they
are a "residual" of former geometric models and have been
retained by reasons of compability and checking.
4rd. Record: dmaxcm,dxcm, dphi,dphih,dphvv,dths,R0cm,ravcm
----------- These correspond to the grid sizes and device properties
(all of them in cm or degrees)
dmaxcm: Maximum half width of the "Central" cell zone
dxcm: Horizontal and vertical cell width
dphi: Toroidal interval angle (45./nphi)
dphih = dphi/2
dphvv = dphi
dths = Poloidal interval angle (360./(nths-1))
R0cm: Major device radius (150. in TJ-II)
ravcm: Average plasma minor radius (Configuration dependent)
b) Magnetic axis data
next nphi+1 records: Pihr, Rax, Zax, Xax, Yax
------------------- Data for the magnetic axis at Axial cell walls
where:
Phir: toroidal angle (in radians)
Rax, Zax: cylindrical coordinates of the axis (cm)
Xax, Yax: remaining cartesian coordinates of the
axis (cm)
next nphi records: Pihr, Rax, Zax, Xax, Yax
----------------- Ibid. for the nphi centers of the Axial cells
c) Horizontal and vertical "Grid" coordinates
next nd (= nhor+1) records: zz
-------------------------- zz (in cm) are the horizontal (and vertical)
coordinates of the "Nodes" on each "middle
toroidal plane" (the same for all the layers),
taking as origin the magnetic axis position
(which is the center of the Axial cell)
d) Cells
Now come the cells in a loop over the nphi toroidal layers (whose
index is iph) containing:
First the "Central" cells
next num records: where "num" is the number resulting of doing two
---------------- loops. The first one (the more external) is "ih", the
horizontal coordinate along the "middle toroidal plane". The
2nd one (and most internal) is "iz", the vertical coordinate.
As not every point is retained and, in addition, there are more
data inserted, this number "num" can not be calculated "a
priori". Each record contains:
ih,iz,iph, Br,Bphi,Bz, s, Vol, Rvt,Xvt,Yvt,Zvt,
Vdrift,Vpol
where:
ih: number of order of the horizontal loop
iz: ibid. for the vertical one
iph: ibid. for the toroidal loop
Br,Bphi,Bz: cylindrical components of the magnetic
field vector at the center of the cell (in T)
s: normalised plasma radius corresponding to the center
of the cell, or the ficticious radius if it is outside
the plasma.
Vol: volume of the cell (in cm3)
Rvt(1:13): radius (in cylindrical coordinates and cm) of
the 13 "key" points of the hexahedra. Such 13 points
are: the 8 vertices, first the anterior faces, then
the posterior ones and last (9-12) the intersections
of the 4 toroidally directed edges with the midle
toroidal plane of the cell layer, always in the
counterclockwise direction, and finally (number 13)
the cell center of mass
Xvt(1:13): ibid. for the cartesian X coordinate
Yvt(1:13): ibid. for the cartesian Y coordinate
Zvt(1:13): ibid. for the cartesian Z coordinate
Vdrift(1:3): cylindrical components of the magnetic drift due
to gradients and curvatures without the "kinetic"
factor (see $2.3) (MKS units)
Vpol(1:3): ibid. for the poloidal rotation velocity (see $2.3)
The last record of this series has a negative value for "ih".
Then come 3 more records
1st Record: kplim(1:4,iph)
2nd Record: pclim(1:4,iph)
3rd Record: pvlim(1:4,iph)
where:
kplim: the 4 extreme limits of ih and iz for this toroidal layer
(ih_min, ih_max, iz_min, iz_max) (remember that not all the
cells inside the "Big Square" are retained).
pclim: the spatial limits of the retained "Grid" region, along the
"middle toroidal plane" taking as origin the center of the
Axial cell. In cm.
pvlim: ibid. for the extremes of the Vacuum cells. These are the
absolute limits of the Eirene geometry for this toroidal
layer (in cm). Everything that falls outside these limits
is outside all the legal cells.
Now come the 8 "Vacuum" cells in a series of records completely similar to
those described above for the "Central" cells, although this time the ih and
iz indices correspond to the "Vacuum" cell numeration (3x3: first left to
right, then bottom to top, with the [2,2] not considered).
As before, the last record of this series has a negative value for "ih".
Finally the toroidal loop ends with a record where both "ih" and "iz" are
negative.
e) Vacuum Vessel
The next series corresponds to the Vacuum Vessel mosaic.
First a single record containing the value of nphi (just a checking
for confirmation).
Then comes a loop for the first nphi toroidal cross sections (index ip).
Inside this loop there appears first a record containing
nths1, id1, phvv
where nths1 is just nths-1 (the number of poloidal intervals)
id1 is allways 1
phvv is the toroidal angle of the cross section (radians).
All these values are used for checking purposes.
The next series is a poloidal loop over the nths angles (index ith)
that contains, for each poloidal angle:
Rvvcm,Zvvcm,Xvvcm,Yvvcm, (triang(1:4,j),j=1,2)
where Rvvcm and Zvvcm are the cylindrical coordinates of the
VV vertex (in cm)
Xvvcm and Yvvcm the remaining cartesian coordiantes
while (triang(1:3,j),j=1,2) are the 3 cartesian components of the
VV triangle external normal vectors for
each of the two triangles of this VV facet
(these vectors are normalised to 1)
and (triang(4,j),j=1,2) is the scalar product of the normal vector
with the central triangle position vector.
All these variables are needed for the VV localisation procedure.
Once the toroidal VV loop is finished there is one more single record for
checking purposes containing:
nths1, id1, phvv
as before.
And then, finally, the poloidal loop is repeated once more for the last
toroidal cross section (that is, the cross section number nax = nphi+1,
coresponding to phi = 45. in TJ-II) but this time the records contain only:
Rvvcm,Zvvcm,Xvvcm,Yvvcm
since the triangles are not necessary for this isolated cross section.
And here the file terminates.
Along all the file reding procedure (subroutine leegeomv of deck leegeom.f of
the Eirene geometry library libgeom), all the quoted check and many more are
performed and at the least failure the lecture aborts.
There exists a program ("testgeo", executed with the sentence runlibb testgeo)
that makes exhaustive checks (around 10**7 in total) of all the geometry
functions: cell localisation, periodicity and stellarator symmetry
transformations, interception of lines with the cell walls, normal vectors
checking, VV triangle localisation, interception of lines with the VV etc.,
etc.
Anoher set of programs allows the graphical visualisation of the geometry
cells: "vertor" represents toroidal cross sections at a given angle of the
cell structure, including the plasma; "verlat" represents an ortographic
projection, either from a top point of view or from a lateral external one, of
a single semiperiod. Both programs are executed by means of the sentence
"runlibb" and use a Disspla-CA emulation running over the graphical library
NCARG.
There is also a program "tradasc" (executed also with the sentence "runlibb")
that translates the binary geometry file to an ascii file, that then could be
ported to any other computer. For this program the Eirene binary geometry is
expected to be (as usual) in fort.30 and the resulting ascii file appears in
fort.2. It must be noted that this last file can be very big: for example in
the case of the 100_44_64 configuration (geom_opt50.1oct file), the original
binary file has 20 MBy while the resulting ascii file (fort.2) has 40.4 MBy.
A second program "tradbin" (executed also with the sentence "runlibb") performs
the reverse operation. It expects to find the ascci translated file at fort.2
and generates the binary geometry in fort.31
In this way the porting of any geometry file between different computers can be
done without using the "grideir" program which needs the TJ-II g3d local
library that, perhaps, could not be available in the second computer.
2.1 Program grideir_dev for any device
----------------------------------
Although the Eirene geometry at Ciemat has been especifically devised for the
TJ-II stellarator, care has been taken that this Eirene implemetation were
"device covariant", so that, in principle, any other device, different from
TJ-II in dimensions, plasma or VV shape, number of periods, etc., could execute
propperly, once the correct geometry file has been generated.
As described above ($1) the grideir program is exclusively prepared for TJ-II,
linking with the g3d local especific library.
In order to allow the possible use of the same geometric model for other
devices (not TJ-II) we have prepared an alternative program: grideir_dev,
entirely similar to grideir but that does not connect to the local g3d library.
This means that the missing subroutines (10 in total) must be constructed by
the user to obtain the required data from the eventual specific alternative
libraries of the device in question.
2.2 Subroutines needed from the user
--------------------------------
As a guide to help construct those subroutines, and to explain what each one of
them is expected to do, a fortran deck with name subr_dev.f has been prepared
containing the 10 required subroutines, fully commented, and applied to a very
simple example: a 5-period tokamak-like stellarator with circular planar
magnetic axis, circular plasma cross section without shift and a circular
Vacuum Vessel. This is one of the simplest examples that can be conceived and,
of course, the use of so complex a geometry model for such a trivial case would
be ridiculous overkill (Eirene has resources enough to deal easily with a
device as simple as this), and the example is presented here only as an
illustrative guide.
The example subroutines are listed below ($2.3) and, as will be seen, are
expected to provide such things as the determination, for a given arbitrary
point, of the plasma normalised radius, the magnetic field vector, the drifts,
determine if the point is inside or outside the Vacuum Vessel, etc. Or also,
for a given toroidal angle, the localisation of the magnetic axis, the VV cross
section for several poloidal angles, etc., etc.
Special care should be taken with the definition of limiters, puffings, NBI
lines, pumps and diagnostic chords in devices different from TJ-II, in order to
avoid possible geometrical inconsistencies or unfortunate choices. As all these
definitions are given to Eirene either by data on namelist or, in the case of
the TJ-II limiters, reading an external file (limiter.dat file) these eventual
problems can be easily avoided.
The program has been checked, not only with the available tests, but also
running the Eirene code for the example case using 1,2,3,4,5,7 and 12 periods,
including the visualisation of the perticle trajectories (program "vertray") in
both computers without any aparent problems.
In all the subroutines cylindrical coordinates are used and, unless explicitely
stated, the units are MKS and radians.
The listing of these subroutines is given below.
2.3 Listing and comments for the subroutines needed
-----------------------------------------------
c subr_dev.f
ccccc
c Ficticious subroutines to replace the TJ-II g3d local library
c In this example the helical axis stellarator TJ-II is replaced
c by a 5 period Tokamak-like stellarator with circular plasma
c planar axis and circular Vacuum Vessel
c (the simplest case possible)
c Unless explicitelly stated all units are MKS and radians
ccccc
subroutine g3d_geom_nbi(ns,file)
character*(*) file
ccccc
c Initialises g3d geometry file in TJ-II
c Both variables are input and are NOT used here in the example.
ccccc
common/com_dev1/ npert
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
ccccc
ccc Sets data for this simple example, transmitted by commons
R0t = 2.0 ! Major radius (in m)
ravt = 0.22 ! Plasma minor radius (in m)
rvvt = 0.30 ! VV minor radius (in m)
aiplt = 100. ! Plasma current (in kA)
Btort = 2.0 ! Toroidal field at axis (in T)
npert = 5 ! Number of periods
ccccc
return
end
subroutine g3d_curr(aicc,aihc,aivf,aitf)
ccccc
c In TJ-II it provides the coil currents (circular and helical coils
c VF coils and TF coils, all in tenths of kA)
c All 4 arguments are output and have only an informative role
ccccc
common/com_dev1/ npert
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
ccccc
c Here, in this example, "aicc" is put to the Plasma current in kA
aicc = aiplt ! Plasma current (in kA) in this example
aitf = 0.
aihc = 0.
aivf = 0.
ccccc
return
end
subroutine g3d_getlim_gr(rxlim,R0,rav,aiot0,Rsw0,
$ nh,nph,nz,nper)
ccccc
c Dimensions and limits for the TJ-II geometry
c All variables are output and in MKS units
ccccc
common/com_dev1/ npert
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
R0 = R0t ! Major radius
rav = ravt ! Plasma radius
Rsw0 = 0. ! Helical axis swing (none here)
nper = npert ! Number of periods
rxlim = ravt + 0.02 ! Maximum plasma extension added 2 cm
ccccc
c The value of rxlim is used to calculate the "Grid" maximum size,
c it is defined as the maximal distance from the plasma border to the
c magnetic axis, all along the torus, and with 2 cm added.
c It is a very important parameter.
ccccc
ccccc
c The 4 values that folow have only informative value and are not used here
c in the example, they are only relevant for TJ-II.
c aiota is iota at plasma border
aiota = R0t/ravt * 0.2e-4*(aiplt*1.e3)/(ravt*100.) / Btort
c aiot0 is iota at magnetic axis. A parabolic J profile assumed here
c in fact any value could do because, after all, it is never used.
aiot0=2.*aiota ! iota at axis (parabolic J radial profile)
nh = 65 ! Horizontal points in g3d grid
nph = 48 ! Toroidal points by semiperiod
nz = 65 ! Vertical points in g3d grid
cccc
return
end
subroutine g3d_findax(Rx,Zx,Phr)
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
ccccc
c Magnetic axis position at a given toroidal angle
c Input: Phr = Toroidal angle in radians
c Output: Rx, Zx = cylindrical coordinates of axis in m
ccccc
Rx=R0t
Zx=0.
return
end
subroutine g3d_rad_gr(Rm,Phr,Zm,s)
ccccc
c Gives the local plasma minor radius, normalised to 1 at border,
c at a given point (cylindrical coordinates)
c Input: Rm,Phr,Zm = cilyndrical coordinates of the point (m and radians)
c Output: s = normalised local plasma radius
ccccc
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
dlm = sqrt((Rm-R0t)**2 + Zm**2) ! distance to axis (m)
s = dlm/ravt
return
end
subroutine g3d_B_Field_cyl(Rm,Phr,Zm,Br,Bphi,Bz)
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
ccccc
c Gives the magnetic field vector components at a given point
c Input: Rm,Phr,Zm = cylindrical coordinates (m and radians)
c Output: Br,Bphi,Bz = cylindrical components of the magnetic field (T)
ccccc
Bphi = Btort*R0t/Rm ! Toroidal Field at Rm (T)
dhm = Rm-R0t ! horizontal distance to axis (m)
dlm = sqrt(dhm**2 + Zm**2) ! distance to axis (m)
s = dlm/ravt ! Normalised plasma radius
if(s.ge.1.e-10) then
thr = atan2(Zm,dhm) ! Poloidal angle
c A parabolical density current profile is assumed here
aiamp = 2.*aiplt*1.e3*s**2*(1.-0.5*s**2) ! Current in Amp
dlcm=dlm*100.
Bp = 0.2e-4*aiamp/dlcm ! Poloidal field in T
Br = Bp*sin(thr)
Bz = -Bp*cos(thr) ! Plasma current goes along increasing Phr
else
Br=0. ! Point on axis: Bp = 0.
Bz=0.
end if
return
end
Subroutine g3d_drift (Rm,Phr,Zm, vdrift, vpol)
dimension vdrift(3),vpol(3)
ccccc
c Gives the drift vector components at a given point
c Input: Rm,Phr,Zm = point cylindrical coordinates (m and radians)
c Output: vdrift = drift vector originated by the gradient and curvature
c that is: [ B x grad(Bm) ]/Bm**3 where [... x ...] denotes
c external product and Bm is the B modulus.
c Note that the "kinetic" factor (1+gam**2)*Ekin/q is
c absent ("gam" is pitch, "Ekin" kinetic energy and
c "q" is charge (see Miyamoto book) (MKS units)
c vpol = Poloidal rotation velocity induced by radial electric
c field. That is : [ grad(s) x B ]/Bm**2
c Note that it is not multiplied by the radial electric
c field intensity -d[V(s)]/ds where V(s) is the electric
c potential radial profile and s the normalised minor radius.
c (MKS units)
c
c The components of both vectors are cylindrical.
ccccc
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
c By simplicity we take both drifts null in this example
c (in any case we do not follow the ions in Eirene)
vdrift(1:3) = 0.
vpol(1:3) = 0.
return
end
subroutine intvvgr(dummy,Phr,rrdum,Rc,Zc)
ccccc
c Any "reasonable" point (Rc, Zc) interior to the Vacuum Vessel
c provided that the choice does not produce "reentrancy" effects at the given
c toroidal angle "Phr" (if reentrancy efects happen the VV could be "clipped")
c
c For TJ-II gives the position of the TF coil center at a given toroidal angle
c
c Input: Phr = Toroidal angle in radians
c Output: Rc, Zc = Cylindrical coordinates of the TJ-II TF coil center
c (or of the chosen "suitable" point) in m.
c "dummy" and "rrdum" are never used.
c
ccccc
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
Rc=R0t ! Here we choose just the center
Zc=0.
return
end
subroutine getvac(nths,Rv,Zv,Phidg, Tht)
dimension Rv(*),Zv(*)
ccccc
c Obtains "nths" points (Rv,Zv) of the Vacuum Vessel toroidal cut
c at toroidal angle Phidg (cylindrical coordinates).
c Input: nths number of required points (greater than 1)
c Phidg = Toroidal angle of the cut (in degrees this time !!!)
c
c Output: arrays Rv and Zv with the "nths" resulting cylindrical
c coordinates (in m)
c
c In TJ-II "Tht" (usually = -180.) is used as the starting poloidal angle
c (measured from the coil centers), and the points are equally spaced
c in poloidal angle (this is not a requirement) and go counterclockwise
c (this is a requirement for normal calculation consistency)
ccccc
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
pi = 4.*atan(1.)
dth = 2.*pi/float(nths-1) ! Poloidal angle interval
Th0 = -pi ! Initial poloidal angle
do ith = 1,nths
Thr = Th0 + (ith-1)*dth ! Poloidal angle (equispaced)
Rv(ith) = R0t + rvvt*cos(Thr)
Zv(ith) = rvvt*sin(Thr)
end do
return
end
subroutine vessel(Rm,Phr,Zm, ives)
ccccc
c Checks if a given point is inside the Vacuum Vessel (VV) or not.
c Input: Rm,Phr,Zm = cylindrical coordinates of the point (m and radians)
c Output: ives. If ives = 1 the point is totaly inside the VV
c if ives = 0 the point is outside or touching the VV
ccccc
common/com_dev2/ R0t,ravt,rvvt,aiplt,Btort
dlm = sqrt( (Rm-R0t)**2 + Zm**2 ) ! Distance at center VV in m
if(dlm.lt.rvvt) then
ives=1
else
ives=0
end if
return
end
----------------------------
This text is in file grideir.txt