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