NUBEAM NTCC Module Documentation
*** NUBEAM was developed at the Princeton Plasma Physics Laboratory ***
*** utilizing funding from the U.S. Department of Energy, 1979--today. ***
The note is illustrated with examples from the test program "nubeam_driver"
which is included in the NUBEAM NTCC module source distribution.
Written by D. McCune -- March 2002. Last update: September 2005.
Email comments & questions to: dmccune@pppl.gov
Contents
========
0. Introduction; Update Notices
1. Dependencies & Prerequisites.
2. Use of XPLASMA.
3. Code-generator & Fortran-90 Modules.
4. Inputs.
5. Executing the NUBEAM calculation.
6. Outputs.
7. 2d Outputs.
8. Odds & Ends: Sawteeth; State file save & restore; PPPL help.
9. Known limitations / problems.
10. Operation of Test Programs
A. Appendices.
The NUBEAM module is a Monte Carlo package for time dependent modeling
of fast ion species in an axisymmetric tokamak. Multiple fast ion
species can be present, due either to beam injection of energetic
neutral particles deposited as ions in the target plasma, or as a
product of nuclear fusion reactions.
The model self consistently handles classical guiding center drift
orbiting, collisional and atomic physics effects during the slowing
down of the fast specie population (represented by an ensemble of Monte
Carlo model particles), with options for sawtooth-induced or anomolous
radial diffusion, in a time evolving plasma with nested magnetic
flux surfaces represented by a numerical MHD equilibrium.
Timestepping is explicit, with the Monte Carlo deposition of new particles
and slowing down of the Monte Carlo ion ensemble against a fixed plasma
for a prescribed time interval, after which the Monte Carlo model state
is recorded, and control is returned to the plasma model, allowing the
plasma model to then evolve for the same time interval, before resuming
Monte Carlo deposition and slowing down in the next time interval, and so
on until the simulation is complete.
A numerical "goosing" technique is used to resolve the disparity in
timescale between fast guiding center drift orbit bounces and the next
next fastest physical timescale, be this collisions, atomic physics, or
anomolous diffusion. This is the key to the computational efficiency
of the model. The method is basically as described in the original
balance, plasma rotation, time-evolving numerical MHD equilibrium, and
fusion product ions) were added since the original paper was written.
These new features are described in a recently published CPC article [B].
The NUBEAM software package originally came about as an integral part of
TRANSP[C,D], but has now been extracted as an independent module for the
National Transport Code Collaboration (NTCC-- http://w3.pppl.gov/ntcc).
This document describes how to use NUBEAM, from the point of view of
a fortran programmer with a driver code capable of utilizing NUBEAM's
capabilities. Physics details are described elsewhere. It is presumed
that the user has access to a unix workstation, fortran-90, C and C++
compilers, and is familiar with the NTCC Modules Library website
( http://w3.pppl.gov/ntcc ), where this and related software can be
found.
references
----------
[A] R.J. Goldston et al., "New Techniques for Calculating Heat and
Particle Source Rates due to Neutral Beam Injection in Axisymmetric
Tokamaks", J. Comp. Phys. 43 (1981) 61.
[B] A. Pankin, D. McCune, R. Andre et al., "The Tokamak Monte Carlo Fast
Ion Module NUBEAM in the National Transport Code Collaboration Library",
Computer Physics Communications Vol. 159, No. 3 (2004) 157-184.
[C] R.J. Hawryluk, "An Empirical Approach to Tokamak Transport", in Physics
of Plasmas Close to Thermonuclear Conditions, ed. by B. Coppi, et al.,
(CEC, Brussels, 1980), Vol. 1, p. 19.
[D] J. Ongena, M. Evrard, D. McCune, "Numerical Transport Codes", in the
Proceedings of the Third Carolus Magnus Summer School on Plasma Physics,
(Spa, Belgium, Sept 1997), as published in Transactions of Fusion Technology,
March, 1998, Vol. 33, No. 2T, pp. 181-191.
NUBEAM uses an intermediate interpolation mechanism (XPLASMA) for
input/output communication of *profile* data. All such data is
in MKS units -- except "temperature" and "average energy" profiles
which in accordance with widely used plasma physics tradition are
given in KeV.
However, scalar quantities, communicated directly as members of
compound data types both for input and output, are NOT guaranteed
to be MKS. It should be born in mind that the original NUBEAM
code was written in cgs, and that some of this is exposed in the
interface, and cannot be easily changed now because of legacy
issues (compatibility with existing usage). The units of all
scalar quantities are given with the documentation of each
compound data type, which is present in the appendices of this
docuement.
There is an active program for development of NUBEAM capabilities.
Some of the long term goals of development are (a) MPI-parallelization,
(b) Handling of interaction with RF wavefields, (c) generalization of
geometrical framework to support stellarator applications. In addition,
minor updates and bugfixes are ongoing.
NUBEAM is developed within the TRANSP production and development system
with source code shared amongst authors via a cvs repository. Improvements
are mapped into the NTCC version of NUBEAM on a weekly basis.
Occasionally, specific upgrades will be summarized here.
(DMC Sept 2005)
Timesteps on the NUBEAM calculation can now be run in a separate process
on a client-server basis. This will ease the use of the (soon to be
released) MPI-parallel version of NUBEAM from a serial transport code
driver. The process for doing this can be described as follows:
(a) client writes NUBEAM inputs into files.
(b) client writes a namelist, readable by the "nubeam_driver"
program, which names the input and output files.
(c) client runs a script to invoke "nubeam_driver" as a sub-
process, which will read input files, calculate a NUBEAM
timestep, and write output files; nubeam_driver will also
maintain its own "state file" to track the state of the
fast ion model and assure accurate restartability.
(d) client reads NUBEAM outputs from files.
(e) at this point client can proceed as if it had used NUBEAM via
its subroutine library API.
The calls NBI_INTERP_PROFILES, NBSTART, DEPALL, ORBALL and NBFINISH
are replaced by the invokation of the server process. A more detailed
description is given in the main body of the help file.
For serial NUBEAM calculations, best performance is to be had by using
the API (subroutine calls). In a test case with a small number of test
particles requiring only 15 seconds per NUBEAM timestep, the use of
client-server communications imposed an extra wall clock cost of
about 6 seconds per step (this as measured in Sept. 2005 on an old
Dell Inspiron 8100 Linux laptop; such results are likely to be highly
system-dependent).
The NTCC module NUBEAM depends on the following NTCC modules which
must also be installed:
PREACT -- atomic/nuclear cross sections and reaction rate data.
XPLASMA --module for numerical representation of axisymmetric plasma
equilibrium geometry and plasma parameter profiles.
NSCRUNCH--tool for regularization of Fourier moments representation
of MHD equilibrium (used by XPLASMA).
PSPLINE --interpolation library (used by XPLASMA).
KDSAW -- a Kadomtsev-style MHD sawtooth "mixing model".
R8SLATEC--a library of mathematical subroutines.
RANDOM -- a portable, parallelizable, high quality random number
generator.
EZCDF -- easy interface to NetCDF -- portable binary files used for
save and restore of the state of the NUBEAM calculation.
PORTLIB --library of subroutines for portable low level utilities
(e.g. for naming files or for fetching current machine
precision).
NUBEAM comes with two test programs, which use additional modules not
needed by NUBEAM itself:
"nubeam_driver" -- non-interactive, self-contained namelist-driven
test program and example NUBEAM driver code (which can also be used
as a standalone NUBEAM timestep server).
"nubeam_test" -- NUBEAM debugging program, with command-line control
and interactive graphics. NUBEAM crashes can be diagnosed with this
program, using state data written from XPLASMA and NUBEAM.
The nubeam_test program uses the additional NTCC modules:
TRREAD -- read TRANSP database (used by nubeam_driver to support an
option to take an MHD equilibrium and profiles from an
existing TRANSP run).
TRGRAF -- command line interactive control of SGLIB plots (used by
nubeam_test).
UREAD -- command line user interface / scripting package (used by
nubeam_test).
SGLIB -- AG4010-compatible scientific graphics package (used by
nubeam_test).
All NTCC modules can be found in the catalog at
http://w3.pppl.gov/NTCC
Note: the NUBEAM download package may already include some of these
modules; others will have to be downloaded separately. Details are
given in the installation instructions (README file) provided with
the module code.
In addition to the NTCC module dependencies, the following freeware math
libraries should also be installed on systems where NUBEAM is to be run:
math libraries:
http://www.netlib.org -- blas & lapack (linear algebra)
http://www.fftw.org -- fftw (fast fourier transforms)
**version 2.1.5**
**cautions -- versions 3.x.x not backward compatible!**
portable binary file i/o library:
http://www.unidata.ucar.edu/packages/netcdf -- NetCDF
(Be sure to check the module README file; dependencies are subject to
change as technical requirements evolve).
NUBEAM is written in fortran-90, and, because portions of the interface
use fortran-90 compound types, some parts at least must be called from
fortran-90 user code. The author has not experimented with calling
NUBEAM from a C or C++ driver code.
Note that the atomic physics database software -- PREACT -- must be
properly configured. It uses a directory of files to store pre-
computed tables e.g. of beam-target maxwellian averaged reaction rate
coefficients. (If the PREACT NTCC module's test programs run correctly,
then, PREACT has been properly configured).
PREACT includes C++ components, and the system and compilers must be
chosen and properly configured to support mixed language fortran/C++
applications. If the PREACT fortran driver test program, "fpreact_test",
runs correctly, then this software prerequisite has been met.
If the NUBEAM test programs can be built and run correctly, then, all
the software prerequisites have been met.
In addition to the above listed components, upon which NUBEAM is
directly dependent, there is also an implicit dependency on an
external thermal neutral gas model (such as the FRANTIC or NUT NTCC
modules). A thermal neutral gas model is needed to compute transport
of neutral atoms in a plasma-- with concomittant ion and electron
source terms and transport of momentum and energy. The thermal
neutral gas model also computes profiles of average energy and
density of neutral atoms, which determine charge exchange loss of
fast ions; conversely, fast neutral deposition and charge exchange
recapture drive thermal neutral sources which should be tracked by
the thermal neutral gas model, and at the same time, the presence of
circulating fast ion species strongly depletes thermal neutrals. Thus
the NUBEAM and thermal neutral gas calculations are inter-dependent.
Author's note: this is an area where more work is needed. The FRANTIC
NTCC module, which, like NUBEAM uses PREACT for atomic physics data, is
fast and efficient; it supports simultaneous presence of multiple
isotopes of H and He, can handle NUBEAM's neutral sink rate profiles,
and computes all necessary transport terms-- but, all in an over-
simplified 1d geometry (flux surfaces remapped to concentric cylinders);
this is what TRANSP uses at present.
The current version of NUT supports a realistic 2d geometry, but the
currently released version only allows one plasma species to be present,
which is not sufficient for simulation of present experiments. Future
versions are expected to remove this restriction.
XPLASMA is a general tool for representing and sharing plasma geometry
and profiles, and for interpolation of such profiles between disparate
grids. Although coded to allow eventual generalization to 3d, XPLASMA
at present is restricted to axisymmetric configurations (i.e. tokamaks).
XPLASMA documentation can be found at:
http://w3.pppl.gov/~pshare/help/xplasma.htm
When setting up to call NUBEAM in a new application, the plasma MHD
equilibrium and various plasma parameter profiles must be loaded into
XPLASMA. The equilibrium can either be loaded directly via the
appropriate XPLASMA calls (see XPLASMA documentation), or it can be
loaded from a standard data source such as EFIT or TRANSP MDSplus trees
available from various tokamak experimental databases, using appropriate
software such as NTCC's I2MEX module. Equilibrium and plasma parameter
profiles are defined against the standard flux surface label coordinate
rho = sqrt([enclosed-toroidal-flux]/[toroidal-flux-within-boundary])
which is a dimensionless quantity ranging from 0 (the magnetic axis)
to 1 (the plasma boundary).
The complete specification of the equilibrium defines the plasma flux
surface geometry and magnetic fields. XPLASMA provides calls for
interpolating field and metric information, and for coordinate
transformations between magnetic, cylindric, and cartesian coordinate
systems used by various components of the plasma and neutral beam
model.
Over the space defined by the MHD equilibrium, XPLASMA allows an
arbitrary number of additional profiles to be defined, most of
which are 1d profiles of the form f(rho). Each profile is actually
an interpolating function, such as a spline (the PSPLINE NTCC module
is used to implement profile "objects"). In XPLASMA, each profile
is assigned a name and an integer "id" number on creation. When a
driver sets up a NUBEAM call, it provides NUBEAM with data structures
containing lists of XPLASMA "ids" pointing to the necessary data, such
as electron and ion temperatures, densities, etc. Similarly, NUBEAM
creates XPLASMA profile objects, and returns lists of XPLASMA "ids" to
the caller, containing the NUBEAM results-- electron heating, ion
heating, fusion rates, particle and momentum source profiles, etc.
By using the XPLASMA (PSPLINE) profile interpolation methods, neither
NUBEAM nor its driver code needs explicit knowledge of the other's
grids (although such knowledge is numerically available using
appropriate XPLASMA calls, if desired).
In addition to setting up the definition of the plasma equilibrium
geometry, fields, and profiles, XPLASMA also allows equilibria to be
"numerically extrapolated" to define a vacuum region and provide a
nested set of nominal surfaces beyond rho=1. The definition of the
vacuum region-- required by NUBEAM-- includes the specification of an
axisymmetric "wall" or "limiter", a closed path encompassing a region
that includes the entire core plasma as a sub-region. The simplest way
to define this region is to specify the wall location at a fixed distance
from the plasma boundary-- but for realistic simulations, a piece-wise-
linear limiter or wall description based on engineering data (such as
the one imbedded in EFIT G files used by many operating tokamaks) would
typically be more suitable. The limiter/wall location is important,
because in normal operation, NUBEAM defines a lost orbit as one coming
within a finite Larmour radius (FLR) of a limiter or wall. XPLASMA
will define (and NUBEAM use) an extrapolated mapped region big enough
to fill the rectangle [Rmin,Rmax]x[Zmin,Zmax] defined by the minimum
and maximum major radius and elevation of the given limiter/wall
configuration. Details of the setup are shown in the driver program
(nubeam_test) and are thoroughly described separately in the XPLASMA
NTCC module documentation, as well.
In its full generality, XPLASMA allows definition of profiles in any
of the forms f(rho), g(rho,theta), h(R,Z), each profile being assigned
its own unique name and XPLASMA id upon creation. XPLASMA profiles
are the mechanism for profile data input and output to/from NUBEAM.
Subsequent sections describing NUBEAM input and output will refer
frequently to XPLASMA profiles and profile ids.
Several components of the NUBEAM interface are written by a Python code
generator. The generator input specification file is meant to be
readable by humans and is a useful document in its own right, as is
described below.
XPLASMA supports an irregular 2d grid (rho,theta) specially tailored
to NUBEAM, and especially suitable for NUBEAM's Monte Carlo algorithm.
This grid bins the core plasma into 2d zones of roughly equivalent
cross-sectional area, organized as a nested sequence of radial "zone
rows", each with a different number of poloidal zones-- proportional
to the zone row's distance from the magnetic axis. These grids are
known by the acronym MCgrid -- Monte Carlo grid-- and are described
in detail in the XPLASMA documentation:
http://w3.pppl.gov/~pshare/help/xplasma.htm
(look for "MCgrid_Routines").
The set of spatially 2d data defined on the XPLASMA MCgrid includes:
2d neutral source profiles (of each neutral species from each
fast species)
2d neutral sink rate profiles (of each neutral species from each
fast species)
2d beam density, <Eperp>, <Epll> <toroidal_velocity> of each fast
species.
2d fast neutral density, average energy and toroidal velocity--
from charge exchange of partially slowed down fast ion species.
and for each fast specie:
4d fast ion distribution functions f(E,vpll/v,rho,theta)
A separate section, below, describes all 2d outputs in detail.
--------------------------
A note on units: all profile quantities passed through XPLASMA are
in MKS, except: specie temperatures and average energies are in KeV.
Almost all NUBEAM input and output quantities are passed through as
elements of fortran-90 compound data types. The compound data type
definitions are contained in the fortran-90 module "nbi_types.f90"
which is imported into user software by the statement
use nbi_types
This module defines data types only-- it does not contain any actual
data. To communicate with NUBEAM, the user code will declare an instance
of each appropriate type. When setting up input for NUBEAM, for example,
an instance of type "nbitype_sys" will be created, a NUBEAM routine
will be called to set default values for the type instance's data
elements, the user routine will modify individual elements appropriately,
and then a NUBEAM routine will be called, passing the instance of the
compound data type, to input the requested settings. The following
code example illustrates more specifically:
use nbi_types
type(nbitype_sys) :: zsys ! basic NUBEAM system inputs
call nbi_init_sys(zsys) ! load defaults in zsys
! user need only set elements for which non-default values are desired.
zsys%runid = my_runid ! set RUNID string
zsys%nonlin = my_msg_lun ! Fortran unit number for messages
zsys%nseed = my_seed ! random number seed
call nbi_set_sys(zsys,ierr) ! system inputs passed to NUBEAM.
if(ierr.ne.0) [...handle error...]
A python code generator, "nbigen.py", is used to maintain the NUBEAM
interface and state file i/o facility. The use of a code generator
makes it much easier to add new input and output options to NUBEAM.
By using compound types with defaults, newly added input options can
be safely defaulted without breaking existing user code-- i.e. a new
option can be added to NUBEAM, and the upated NUBEAM code distributed,
without breaking user implementations which might not choose to make
use of the new option right away.
The NUBEAM source is distributed with nbigen.py and also the file
"nbspec.dat" which is the NUBEAM module data specification file.
When nbigen.py is run, it reads nbspec.dat and from these specifications
writes interface routines such as "nbi_init_sys" and "nbi_set_sys", as
well as NUBEAM's internal data modules "nbi_dimensions" and "nbi_com".
The nbspec.dat file is important because it defines the actual input
and output data for the entire NUBEAM module, as well as the module's
internal data structures, as well as defining elements of the state
of the time dependent NUBEAM calculation. The nbspec.dat file is
heavily commented, and the user may refer to it for many details of
NUBEAM input/output. Comments in the nbspec.dat file that start with
"!" pertain to NUBEAM data elements; comments that start with "#"
pertain to nbspec.dat format and details of the python code generator.
All the input and output fortran-90 compound data types and associated
calls are shown in the appendix of this document (these appendices are
also written by nbigen.py from the definitions in nbspec.dat).
In summary, the code generator reads
nbspec.dat -- NUBEAM i/o and module data specification file
basic data types are: R -- REAL, I -- INTEGER,
L -- LOGICAL, D -- REAL*8, C*n -- CHARACTER*n.
An instance of each specified item is generated
in NUBEAM's internal fortran-90 module; in addition,
input/output items are made members of publicly
declared compound data types used for transfer of
information in and out of the NUBEAM internal module.
and the code generator writes the following NUBEAM components:
nbi_dimensions_mod.f90 --
nbi_com_mod.f90 -- NUBEAM's internal fortran-90 modules. These can
not be referenced directly by user code, but they
contain copies of comments from nbspec.dat, and may be
useful as documentation.
nbi_types.f90 -- public fortran-90 module defining input/output
data structures.
nbi_alloc.f90 -- allocate most NUBEAM data structures based on user
specified grid sizes-- called by user after a
successful "nbi_set_dims" call...
nbi_alloc2.f90 --allocate remainder of NUBEAM data structures based
on additional user input (e.g. desired number of
Monte Carlo particles in simulation). The user
does not call this directly.
nbi_dalloc.f90 --deallocate NUBEAM data structures. Since the NUBEAM
simulation has state, this should not be done unless
all NUBEAM fast ions have thermalized, and NUBEAM is
never again to be called.
nbi_init.f90 -- set defaults for input data structures
nbi_set.f90 -- pass inputs to NUBEAM. Since NUBEAM assumes state,
this is generally done ONLY ONCE per run per input
data structure.
nbi_get.f90 -- fetch current settings of input data structures.
nbi_update.f90 --update input. Most input data elements are not
updatable. For example, array dimensions can only
be set once at the beginning of a run. Attempts to
update non-updatable members are ignored. See the
Inputs section for more details.
nbi_state_io.f90 -- subroutines for NUBEAM state save/restore operations.
User callable routine: see nbi_states.f90.
nbi_inputs_io.f90 -- subroutines for reading/writing NUBEAM inputs files
(inputs not communicated via XPLASMA); used to enable
NUBEAM client-server implementation.
nbi_outputs_io.f90 -- subroutines for reading/writing NUBEAM outputs files
(outputs not communicated via XPLASMA); used to enable
NUBEAM client-server implementation.
nbi_ascii.f90 -- ascii dump of NUBEAM input and output quantities, for
debugging.
nbo_get.f90 -- retrieve NUBEAM output (into output data structures).
Portions of subroutines of the nubeam_driver test program are also
written by nbigen.py.
User settable controls for the NUBEAM calculation are defined in
nbspec.dat in the sections "Array_Dimensions" and "Inputs". The
definitions are replicated in the appendices of this document.
User settable array dimensions are given in the "Array_Dimensions"
section in lines of the form
U <name> = <default-value>[:<max-value>] [! comments...]
If no <max-value> is specified, then there is no limit, although in
no case would it be reasonable for the actual value to differ from
the default by many orders of magnitude.
The nbspec.dat file also contains a number of definitions of additional
array dimensions that are not directly user settable.
An example of a user settable array dimension is:
U mib = 32:200
! max number of neutral beams "nbeam" in simulation.
! beams may be added during a simulation but not removed.
which specifies the maximum number neutral beams in a given simulation.
All array dimensions are set once at the beginning of a run. Once set,
they cannot be modified. All array dimensions are constituents of the
state of the NUBEAM calculation.
As with all inputs, the user sets dimensions by initializing the
elements of a compound data structure and passing the structure to
NUBEAM. The subroutine "nbi_init_dims" sets the default values for
all elements; user code modifies individual elements as needed, and
then a call to "nbi_set_dims" actually defines the array dimensions
and dynamic memory allocation for NUBEAM in the context of the
current run.
However, array dimension inputs must be provided first, before any
other inputs can be accepted by NUBEAM.
-> a special note on the relationship of radial grids in NUBEAM.
There is a 1d radial grid which will have an active size "nzones"
as set in the "grid" input data type. The corresponding array
dimension "mj" should be set to nzones + 3. Also, the 2d
"MCgrid" XPLASMA structure which supports NUBEAM 2d f(rho,theta)
outputs will support a number of zone rows of which nzones must
be an integer multiple. Thus, if a MCgrid with 10 radial zone
rows is planned, the value of "nzones" must be chosen from the set
{10, 20, 30, 40, 50, ... }. This is illustrated in the working
example code nubeam_driver; the driver subroutine "nbdrive_checknaml"
enforces this by requiring nubeam_driver namelist controls "nzones"
to be an integer multiple of "nzones_fb".
Inputs (other than array dimensions) are defined and documented in the
"Inputs" section of nbspec.dat. Inputs are grouped in "blocks" each of
which corresponds to a single compound data type to be defined in
nbi_types.f90. The specification of an individual input element has
the form
<data-type>[A][updatability] <name>[(dims,...)] = <default-value>
[! comments...]
Individual element data types are defined by a code (R/I/L/D/C*n for
REAL, INTEGER, LOGICAL, REAL*8, and CHARACTER*n, respectively). The
presence of "A" denotes "array", in which case array dimensions must
also be given.
By default, input elements can ONLY BE SET ONCE and are not updatable.
Updatable inputs are indicated by an "M" in the [updatability] field.
Items with an "X" in the [updatability] field are members of "lists".
Lists can be extended with new list members added, but, old list members
once added can never be removed.
The nbspec.dat file contains detailed comments describing each and every
input element. Here are a few sample input elements from "nbitype_sys",
as given in nbspec.dat:
Block: sys ! basic system information for Monte Carlo code
C*32 runid ="NUBEAM"
! runid string -- alphanumeric string used as prefix for
! output files (e.g. state files, etc.). Should generally
! be set to a runid string associated with the driver code.
IM nonlin=6 ! unit number for messages (updatable).
I nseed=1 ! random number seed; default not recommended!
The use of the NUBEAM input control system is best illustrated by
example: the nubeam_driver test program, described below.
All inputs (except XPLASMA profiles and MHD equilibrium) are considered
part of the internal state of the NUBEAM calculation.
The following are the input compound data types:
nbitype_sys -- basic system information.
nbitype_times -- timestep.
nbitype_grid -- grid information.
nbitype_beams -- list of neutral beams with geometry information
also includes list of thermal plasma species; for
each beam specie there must be a corresponding
thermal specie into which the beam ions thermalize.
nbitype_impurity-- impurity species.
nbitype_minority-- RF minority species.
nbitype_powers -- neutral beam powers.
nbitype_fusion -- fusion product model controls.
nbitype_fpp -- controls for using NUBEAM to calculate a deposition
and charge exchange source for a Fokker Planck fast
ion slowing down calculation.
nbitype_num -- general numerical control options.
nbitype_atomic -- atomic physics model controls and options.
nbitype_collid -- collision operator controls and options.
nbitype_flr -- finite larmor radius correction model options.
nbitype_saw -- sawtooth model controls.
nbitype_adif -- anomolous diffusion operator controls and options;
model is activated by providing anomolous diffusivity
as an input (XPLASMA) profile.
nbitype_ripple -- TF field ripple loss (crude models).
nbitype_fishbone-- MHD fishbone induced loss (crude model).
nbitype_outcon -- output control options.
nbitype_box -- 3d neutral density "cartesian box" output options.
nbitype_tube -- neutral densities along chords-- used for charge
exchange eflux diagnostic simulations.
nbitype_misc -- miscellaneous options.
nbitype_profiles-- XPLASMA ids for input profiles (ion and electron
and neutral gas temperatures, densities, angular
velocities, loop voltage, anomolous diffusity, ...
See nbspec.dat or the appendices of this document for the member data
element descriptions.
--------------------------
A note on units: all profile quantities passed through XPLASMA are
in MKS, except specie temperatures and average energies are in KeV.
Scalar quantities have units as documented in "nbspec.dat" and in the
appendix.
This is best illustrated by example.
The program "nubeam_driver", included in the NUBEAM NTCC module, is
a complete working example of the sequence of calls needed to drive
NUBEAM.
This program is structured so that the main routine does very little.
It accepts a character string "runid" as a command line argument, forms
filenames (namelist input, ascii output files) from this, and passes
these names with fortran i/o unit numbers to subroutine "nbdrive_main"
in which the actual work is done. The routines which make XPLASMA
and NUBEAM calls to set up the simulation are all in the library
"nbdrive" subroutine library. The code is built in this manner so
that a call out to nbdrive_main can be executed from inside a larger
application, e.g. to test successful installation of NUBEAM. The
total "nbdrive" library is about 2500 lines of code.
The "main subroutine" nbdrive_main carries out the following operations.
The routines making XPLASMA and NUBEAM calls are indicated.
-- set namelist defaults (call nbdrive_defaults).
-- read namelist (call nbdrive_readnaml).
-- basic checks of namelist values (call nbdrive_checknaml).
-- make optional NUBEAM calls performing additional checks on
input data (call nbdrive_nubeam_checknaml).
-- also get energy grids for fast ion distribution functions.
-- open output files (call nbdrive_openfiles).
-- allocate nbdrive memory and sets up timebase (call nbdrive_alloc).
[loop over "nstep" NUBEAM timesteps]
|
| (if namelist indicates):
| -- restore XPLASMA state (XPLASMA eq_restore call).
| -- restore NUBEAM state (NUBEAM nbi_states call).
| -- check state (call nbdrive_ckstate).
| (previously interrupted simulation now ready to be restarted).
|
| (first timestep only):
| -- NUBEAM initialization: array dimensions; input controls that
| are set only once (call nbdrive_nubeam_init).
| -- load MHD equilibrium, fields and vacuum region definition into
| XPLASMA -- (these are not varying in time for some operational
| modes of nubeam_driver-- but in a real transport simulation, the
| equilibrium would evolve in time). (call nbdrive_mhdeq).
|
| ** all timesteps **
|
| -- update MHD equilibrium (when time dependent TRANSP input data is
| used-- a nubeam_driver option) (call nbdrive_mhdeq).
|
| -- update time varying data other than profiles, e.g. beam powers
| and voltages (these do not actually vary in time, in nubeam_driver).
| (call nbdrive_nubeam_update).
|
| -- update ion density profiles (subtract previous timestep fast
| ion density from electron density; use Zeff and quasineutrality
| equations). (If time dependent TRANSP input data is available,
| update other profiles as well). (call nbdrive_prof).
|
| -- load time varying profile data into XPLASMA; transfer XPLASMA
| ids to NUBEAM. (call nbdrive_xpload).
|
| -- *** execute NUBEAM timestep *** (call NBDRIVE_EXECUTE).
|
| -- append results to ascii output files (call nbdrive_summary;
| call nbdrive_details).
|
| (if namelist requests interrupt on current step):
|
| -- save XPLASMA state (eq_save call).
| -- save NUBEAM state (nbi_states call).
| -- interrupt exit
|
[end of loop]
--cleanup and normal exit.
The NBDRIVE_EXECUTE call contains some file i/o calls to support
debugging in case of a crash, and:
CALL NBSTART(ierr) ! NUBEAM call-- initialize timestep.
! (set up optimized atomic physics tables, etc.)
if(ierr.ne.0) then [...handle error...]
CALL DEPALL(ierr) ! NUBEAM call-- beam/fusion product deposition
if(ierr.ne.0) then [...handle error...]
CALL ORBALL(iorbtot,ierr) ! NUBEAM call-- orbiting and slowing down
! integer, intent(out) :: iorbtot returns
! zero if there were no orbits to follow.
if(ierr.ne.0) then [...handle error...]
CALL NBFINISH(ierr) ! NUBEAM call-- finish timestep.
if(ierr.ne.0) then [...handle error...]
The reason the main calculation is broken up into four calls is to
allow development of future module options. For example, it is
planned to allow the call to ORBALL, which uses most of the cpu
time in a NUBEAM call sequence, to be replaced by a Fokker-Planck
solver.
Author's note: An earlier version of NUBEAM did not have the output
arguments "ierr" in all four of these calls. These have been added
now, to support eventual removal of all "bad_exit" calls inside NUBEAM,
as has been requested by prospective customers.
NUBEAM calculations can also be run in a separate process (which will
be available soon as an MPI-parallelized process. In the context of
the example program "nubeam_driver", this involves replacing the NUBEAM
calls
nbi_interp_profiles (as in nbdrive_main.f90 or nbdrive_xpload.f90)
nbstart (as in nbdrive_execute.f90)
depall (as in nbdrive_execute.f90)
orball (as in nbdrive_execute.f90)
nbfinish (as in nbdrive_execute.f90)
with a call to a script that will actually launch nubeam_driver with
a namelist that puts the program into "proxy_exec" mode. An example
of the necesarry namelist follows:
&nbdrive_naml
! proxy namelist for TRANSP run: 112989Y03
mode = "proxy_exec"
sawflag = .FALSE.
proxy_inputs = "112989Y03_proxy_inputs.dat"
proxy_xplasma = "112989Y03_proxy_xplasma.dat"
proxy_state = "112989Y03_proxy_state.dat"
proxy_outputs = "112989Y03_proxy_outputs.dat"
/
which specifies mode = "proxy_exec_new" for the first timestep of a
NUBEAM simulation, or "proxy_exec" for a subsequent timestep, which
then will read state file information based on preceding steps.
The namelist flag "sawflag" specifies whether a Kadomtsev-style mixing
step should be applied to the Monte Carlo fast ions, after the state
file is read.
The namelist specifies the files containing (non-XPLASMA) NUBEAM inputs,
XPLASMA data (on input and modified on output), (non-XPLASMA) NUBEAM
outputs, and, a state file, which is read (unless mode="proxy_exec_new")
and then updated at the end of the step.
The code (nubeam_driver itself) manages the state file-- but see the
subtopic on state synchronization.
The "inputs file" can be written by calling a python-generated NUBEAM
subroutine:
call nubeam_inputs(<filename>,"W",ierr)
And the XPLASMA profile inputs are written by the XPLASMA call:
call eq_save(<filename>,ierr).
Then, a script is executed to launch the subprocess. This could involve
many steps, e.g. connecting to a remote parallel supercomputer, and
sending/receiving the input/output files.
After the server completes the NUBEAM timestep, the data is read back
into memory of the current process by using the following calls. The
NUBEAM "outputs file" is read by calling a python-generated subroutine:
call nubeam_outputs(<filename>,"R",ierr)
And the XPLASMA data, now augmented with NUBEAM output profiles, including
fast ion distribution functions, is read in with:
call eq_restore(<filename>,ierr).
Here is an example of code (from TRANSP) used to set up the inputs for
a NUBEAM call which is to be run as a service. This writes the namelist
for nubeam_driver in "proxy_exec" mode, writes the inputs file, and
writes the xplasma file containing MHD equilibrium and profile inputs.
! NUBEAM runs outside TRANSP as a parallel server process...
! (a) write proxy namelist
call nbi_chk_inputs(ierr)
if(ierr.ne.0) then
call errmsg_exit(' ?set_nbprofs: nbi_chk_inputs error.')
endif
zprefix1 = trim(fi_workpath)//trim(runid)//'_proxy_' ! for file i/o
zprefix2 = trim(runid)//'_proxy_' ! for namelist values
open(unit=lunphd, &
file=trim(zprefix1)//'namelist.dat', &
status='unknown',iostat=ierr)
if(ierr.ne.0) then
call errmsg_exit(' ?set_nbprofs: proxy_namelist file open failure.')
endif
write(lunphd,*) '&nbdrive_naml'
write(lunphd,*) ' '
write(lunphd,*) ' ! proxy namelist for TRANSP run: '//trim(runid)
write(lunphd,*) ' '
if(ifirst) then
write(lunphd,*) ' mode = "proxy_exec_new"'
else
write(lunphd,*) ' mode = "proxy_exec"'
endif
write(lunphd,*) ' proxy_inputs = "'//trim(zprefix2)//'inputs.dat"'
write(lunphd,*) ' proxy_xplasma = "'//trim(zprefix2)//'xplasma.dat"'
write(lunphd,*) ' proxy_state = "'//trim(zprefix2)//'state.dat"'
write(lunphd,*) ' proxy_outputs = "'//trim(zprefix2)//'outputs.dat"'
write(lunphd,*) '/'
close(unit=lunphd)
! --> record NUBEAM state filename in trcom -- needed by wrstrt/resume
pserve_state = trim(zprefix1)//'state.dat'
! (b) save xplasma state for proxy nubeam exec
call eq_save(trim(zprefix1)//'xplasma.dat',ierr)
if(ierr.ne.0) then
call errmsg_exit( &
' ?set_nbprofs: proxy NUBEAM XPLASMA eq_save failure.')
endif
! (c) save NUBEAM control inputs for proxy nubeam exec
call nbi_inputs(trim(zprefix1)//'inputs.dat','W',ierr)
if(ierr.ne.0) then
call errmsg_exit(' ?set_nbprofs: proxy NBI_INPUTS write failure.')
endif
The nubeam_driver execution is then controlled from a script which moves
to the directory containing the input/output data files and invokes the
code with an argument which identifies the namelist. (Much of the script
deals with error handling and is somewhat system dependent).
Then, when control is returned to TRANSP, it executes:
! fetch results of parallel NUBEAM step
zprefix = trim(fi_workpath)//trim(runid)//'_proxy_'
! (a) read xplasma, updated by PLL NUBEAM
call eq_restore(trim(zprefix)//'xplasma.dat',ierr)
if(ierr.ne.0) then
call errmsg_exit( &
' ?fimain: proxy NUBEAM XPLASMA eq_restore failure.')
endif
! (b) read scalar-reduced outputs of PLL NUBEAM
call nbo_outputs(trim(zprefix)//'outputs.dat','R',ierr)
From this point on, all the code for interpolating NUBEAM outputs (e.g.
from the NUBEAM grid back to the TRANSP grid) is identical with that
which would be used if NUBEAM had been driven in the traditional way
by subroutine calls.
In the context of TRANSP this architecture is a prototype for allowing
a serial driver code to invoke a parallel subprocess running on
separate computational hardware. Ultimately we expect to use this not
just for NUBEAM but for other computational intensive modules as well
(such as RF wavefield and ray tracing calculations).
If a time-dependent transport code maintains its own mechanism for
check-point restart, it is important to synchronize the state of the
NUBEAM simulation with the driver code's own state, so that, in case
of restart, the NUBEAM simulation will resume execution on the correct
timestep.
How this is done? I will describe what TRANSP does. I use an NTCC
portlib integer function "jsystem" to execute a system command. Then
character variables
character*120 nubeam_statefile ! NUBEAM's own state
character*120 nubeam_statesave ! NUBEAM state synchronized to TRANSP
When writing a TRANSP checkpoint state:
character*500 zcmd
zcmd = 'cp '//trim(nubeam_statefile)//' '//trim(nubeam_statesave)
istat = jsystem(zcmd) ! save a synchronized NUBEAM state.
When restarting TRANSP from a previously saved checkpoint state:
zcmd = 'cp '//trim(nubeam_statesave)//' '//trim(nubeam_statefile)
istat = jsystem(zcmd) ! restore a synchronized NUBEAM state.
NUBEAM output can be divided into three categories: scalars (and
arrays of scalars), 1d profiles f(rho), and 2d profiles f(theta,rho)
on the 2d irregular grid. The 2d profiles are considered in the
next section "2d_Outputs".
Scalars contains 0d information, such as terms in the global energy,
momentum, or particle balance for individual fast ion species--
typically an array of numbers dimensioned by the fast ion species
index array dimension (mibs). Examples will be given below.
1d profiles are of the form f(rho), are defined on the NUBEAM internal
grid but stored as XPLASMA objects, allowing interpolation to the caller's
own rho grid. For each 1d output profile, a unique integer XPLASMA id
code is returned.
All scalar and 1d profile outputs are defined in nbspec.dat; Outputs
are organized into several blocks, each of which corresponds to a
compound data type definable in the caller's code by the statement
use nbi_types
The caller's code creates an instance of each output data type, such as
type (nbotype_power_balance) :: zpbal
type (nbotype_powers) :: zpowers
and then, after the beam code timestep has been completed, issues a call
call nbo_get_power_balance(zpbal)
after which individual elements can be extracted into the caller's data
structures, as in:
do isi=1,nsfast
call my_imap(isi,isb) ! mapping from NUBEAM to user fast ion index
! (user defined).
pinjs(isb) = zpbal%pinjs(isi)
pshine(isb) = zpbal%bpshins(isi)
...
enddo
Each element of each compound data type is defined in nbspec.dat in
lines of the form:
<data-type>[A][S][code] <name>[(dims,...)]
[! comments...]
Where <data-type> usually specifies D for REAL*8, [A] is present and
dimensions given if the output type is an array, and [S] indicates an
item that is reused by NUBEAM in its next timestep and must therefore
be saved in the NUBEAM code's state file (the python generator handles
this). There is also a symbolic [code] specification which indicates
a profile subtype {blank, "^", "|", "~", or "@"} which will be described
below.
As an example of a scalar data definitions, the block "power_balance"
contains the specifications:
RA pinjs(mibs) ! watts
! injected power (or fusion product source power), by species.
RA bpshins(mibs) ! watts
! BEAM SHINE THRU POWER, BY BEAM SPECIES
which are used in the preceding "zpbal" pseudo-code example.
1d f(rho) profile outputs are similarly declared, but are identifiable
by the presence of the array dimension "mj", which is the rho-grid
dimension in NUBEAM. Thus for example, the block "powers" contains
(as seen in nbspec.dat):
RA pbe(mj) ! watts
! POWER TO ELECTRONS FROM BEAM HEATING -- for electron power balance
RA pbes(mj, mibs) ! watts
! POWER TO ELECTRONS FROM fast ion HEATING, BY FAST SPECIES
RA pbi(mj) ! watts
! POWER TO IONS FROM BEAM SLOWING DOWN -- for ion power balance
RA pbis(mj, mibs) ! watts
! POWER TO IONS FROM fast ion HEATING, BY SPECIES
RA pfe(mj) ! watts
! POWER TO ELECTRONS FROM FUSION ION HEATING -- for electron power bal.
RA pfi(mj) ! watts
! POWER TO IONS FROM FUSION ION SLOWING DOWN -- for ion power bal.
which define a set of XPLASMA output profiles. The corresponding compound
data type contains for each of these items an integer scalar or array which
is formed by deleting the "mj" dimension from the above specifications, i.e.
in this case:
integer :: id_pbe, id_pbes(mibs), id_pbi, id_pbis(mibs), id_pfe, id_pfi
which are members of a compound fortran-90 structure of type "nbotype_powers".
Thus, to retrieve the heating profiles, the user code employs:
use nbi_types
type (nbotype_powers) :: zpowers
...
call nbo_get_powers(zpowers) ! get profile ids
id = zpowers%io_pbi ! interpolate PBI to user grid
call my_xplasma_extractor(id, my_pbi_array, ...)
id = zpowers%io_pbe ! interpolate PBE to user grid
call my_xplasma_extractor(id, my_pbe_array, ...)
...etc...
Where the last calls are to a user written routine that uses XPLASMA to
interpolate the NUBEAM profile data to the caller's own rho-grid, with
possible normalizations or units transformations applied.
The majority of output profiles are integrated particle, power, torque,
source or energy density profiles in MKS units of (# of particles), watts,
Nt-m, #/sec, or Joules, respectively. To reconstruct local density on the
user's rho-grid, the typical method would be to use XPLASMA to interpolate
the integrated profile to the boundaries of the user's rho-grid, and then
employ a finite-difference formula such as
n(j) = (N(rho(j+1/2))-N(rho(j-1/2)))/dVol(j)
where n(j) = the density in the user's zone j, N(rho) is the integrated
profile interpolated to the boundaries of zone j, and dVol(j) is the volume
of zone j. XPLASMA interpolation defines N(rho).
Some output profiles are different, and these are indicated in nbspec.dat
by profile-type symbol codes in nbspec.dat. Examples of such profiles are
the fast ion driven currents:
RA@ curjbeam(mj) ! beam ion driven current (shielded), amps
RA@ curjfusi(mj) ! fusion ion driven current (shielded), amps
The "@" code indicates that these profiles, which specify enclosed
toroidal current as a function of flux surface label, are integrated
using the cross-sectional area rather than the volume of plasma flux
zones.
The complete set of profile type codes are as described in nbspec.dat
comments:
# for profiles there are some output category codes:
#
# <default> -- item is an integrated density, source density, torque
# density, or power density (n, n/sec, Nt-m, watts), starting
# with zero at the magnetic axis. The corresponding number in
# a user's zone [rho(j),rho(j+1)] is the difference in profile
# value between these two end points. Internally these are
# Monte Carlo sums done on a zone by zone basis; quantity will
# be smoothed (cf dbxsmoo input parameter).
#
# "@" -- item is an integrated profile of enclosed current density (A)
# starting with zero at the magnetic axis; quantity is still
# smoothed.
#
# "~" -- item is a characteristic parameter such as "average energy"
# not integrated from the axis; quantity is still smoothed.
#
# "^" -- item is a characteristic parameter such as "slowing down time"
# but does not need to be smoothed; piecewise linear output
# profile is provided.
#
# "|" -- item is a boundary oriented "flow" quantity; evaluated at
# rho(j) is the nominal flow across the specified surface; units
# vary. Internally these are Monte Carlo sums of flows across
# boundaries, generally using orbit averaging to reduce variance.
# quantity is fit using piecewise linear interpolation, NO
# smoothing.
The scalar and profile outputs are organized into several blocks, each
of which has a corresponding fortran-90 compound data type containing
scalars and XPLASMA integer id's for profiles:
nbotype_fusion -- computed fusion reaction rates.
nbotype_density -- fast ion particle and energy density data.
nbotype_deposition-- details on fast ion deposition.
nbotype_nclass -- fast ion statistics wanted by NCLASS
(W. Houlberg's neoclassical NTCC module).
nbotype_n0_fast -- fast neutral densities (flux surface averaged).
nbotype_trap_fraction -- trapping fraction information.
nbotype_erngfi_output -- densities and trapping fractions split
over a coarse user defined energy grid.
nbotype_average_energies -- average perpendicular and parallel
energies of fast ion species.
nbotype_excited_states -- excited states effect on deposition
and charge exchange recapture.
nbotype_rotation -- flux surface averaged angular velocities of
fast ion populations.
nbotype_compression -- effects of change in MHD equilibrium on fast species.
nbotype_powers -- plasma heating profiles.
nbotype_currents -- driven currents.
nbotype_sources -- fast ion driven thermal ion and electron sources.
nbotype_torques -- torques (for the thermal ion angular momentum
balance equation).
nbotype_neutral_sources -- thermal neutral atom sources.
nbotype_neutral_sinks -- thermal neutral atom sinks.
nbotype_recapture -- details of fast ion charge exchange loss and recapture.
nbotype_fokker_planck -- classical slowing-down and pitch-angle-scattering
times of newly deposited fast ions.
nbotype_mc_statistics -- some numerical data on the Monte Carlo
algorithm.
nbotype_radial_current -- radial currents (flows) (amps).
nbotype_power_balance -- scalar power balance data on fast species.
nbotype_momentum_balance -- scalar angular momentum balance data on
fast species.
nbotype_torque_work -- work done on rotating plasma by fast ion torques.
nbotype_particle_balance -- scalar particle balance data on fast species.
It should be noted that the indexing of fast ion species between NUBEAM
and the user's code might not be identical. NUBEAM provides an integer
mapping function
inubeam = index_nbfi(iZ,iA,iType)
which returns the NUBEAM index of the fast specie of indicated integer
Z (e.g. 1 for Hydrogen), A (e.g. 2 for Deuterium), and:
iType = 1 ... if seeking a Monte Carlo beam specie
iType = 2 ... if seeking a Fokker-Planck beam specie (expected
for the future: NUBEAM deposition combined with
Fokker-Planck with ICRF heating in a separate module
calculation... not yet available).
itype = 3 ... if seeking a Monte Carlo fusion product specie
This is probably best illustrated by an example: see nbdrive_mapfast.f90
in the nubeam_driver test program.
--------------------------
A note on units: all profile quantities passed through XPLASMA are
in MKS, except specie temperatures and average energies are in KeV.
Scalar quantities have units as documented in "nbspec.dat" and in the
appendix.
NUBEAM creates several outputs on the "irregular" 2d (rho,theta) MCgrid
structure described in the XPLASMA documentation. These are stored in
XPLASMA, associated with a MCgrid "id" and identified by character
strings:
SRBN02_<neutral-gas>_<fast-specie> -- thermal neutral source profiles
#/sec in each zone.
example: SRBN02_H0_D_MCBEAM
-- 2d H0 neutral source due to D Neutral Beam specie.
SNBN02_<neutral-gas>_<fast-specie> -- neutral sink <nb*sigma*v>, 1/sec
due to charge exchange
example: SNBN02_D0_HE4_MCBEAM
-- 2d D0 neutral sink due to charge exchange with
HE4 Neutral Beam specie.
SNB0I2_<neutral-gas>_<fast-specie> -- neutral sink <nb*sigma*v>, 1/sec
due to impact ionization
example: SNBN02_D0_HE4_MCFUSN
-- 2d D0 neutral sink due to impact ionization by
HE4 Fusion Product specie (fusion alphas).
BDENS2_<fast_specie> -- 2d fast ion density profile, #/zone
EBA2PP_<fast_specie> -- 2d profile of <Eperp>, KeV
EBA2PL_<fast_specie> -- 2d profile of <Epll>, KeV
FBM_<fast_specie> -- ** fast ion distribution function **
-- a function of pitch and energy at each 2d zone
-- pitch and energy grids as per user input.
XPLASMA also defines for any MCgrid structure:
BMVOL -- volume of each MCgrid zone (m**3).
which is frequently needed, e.g. for normalization to density.
A note on the definition of "pitch". Although often denoted as vpll/v,
NUBEAM assumes a convention where vpll/v > 0 means moving "with the
plasma current" and vpll/v < 0 means moving "against the plasma current".
This means:
pitch = (v.B)/(|v|*|B|)*sigma, where
sigma = 1 if Btoroidal and Itoroidal are parallel, and
sigma = -1 if Btoroidal and Itoroidal are anti-parallel.
---------------
It should be noted that the indexing of fast ion species between NUBEAM
and the user's code might not be identical. NUBEAM provides an integer
mapping function
inubeam = index_nbfi(iZ,iA,iType)
which returns the NUBEAM index of the fast specie of indicated integer
Z (e.g. 1 for Hydrogen), A (e.g. 2 for Deuterium), and:
iType = 1 ... if seeking a Monte Carlo beam specie
iType = 2 ... if seeking a Fokker-Planck beam specie (expected
for the future: NUBEAM deposition combined with
Fokker-Planck with ICRF heating in a separate module
calculation... not yet available).
itype = 3 ... if seeking a Monte Carlo fusion product specie
This is probably best illustrated by an example: see nbdrive_mapfast.f90
in the nubeam_driver test program.
---------------
In the following discussion, inubeam refers to NUBEAM's fast ion
specie index...
A set of routines are provided for indexing into the distribution
function's pitch and energy grids. The following code sample gives
the flavor. Working examples are in the nubeam_driver test program,
subroutine nbdrive_summary_2d.f90 ...
integer iZ,iA,iType ! fast specie information (set by user)
integer inubeam ! NUBEAM fast specie index
integer inumE ! number of energy zones
integer inumA ! number of vpll/v zones
real zE,zEmax ! energy (user supplied) & max energy
real zvpllv ! vpll/v (user supplied)
integer iEzone ! energy zone index
integer iAzone ! vpll/v zone index
inubeam = index_nbfi(iZ,iA,iType) ! get NUBEAM specie index (see [9])
call get_fbm_index_info(inubeam,zEmax,inumE,inumA)
...
iEzone = 1 + (zE/zEmax)*inumE ! init guess
iEzone = ndxfbe(zE,inubeam,iEzone) ! actual energy zone index
iAzone = ndxfba(zvpllv,inubeam) ! vpll/v zone index
! ... can now refer to individual elements of the fast ion distribution
! ... data array FBM(iEzone,iAzone,<spatial-zone-index>,<specie-index>)
subroutine get_fbm_index_info and integer functions index_nbfi, ndxfbe, and
ndxfbe are defined by the NUBEAM library.
-------------
To fetch 2d outputs, the appropriate character name strings have to
be generated. The names consist of a root concatenated with specie
identification label substrings and separated by underscore ("_")
characters.
The <fast_specie> label substrings can be generated by NUBEAM subroutine
calls:
integer inubeam ! the NUBEAM index of the fast specie wanted...
integer inum_fast,igot
character*20, dimension(:), allocatable :: zlbls
character*40 zlbl_2dobj
call nbfi_numfast(inum_fast)
allocate(zlbls(inum_fast))
call nbfi_lbl(zlbls,inum_fast,igot)
zlbl_2dobj = 'BDENS2_'//zlbls(inubeam)
call mcgrid_getobj(..., zlbl_2dobj, ...)
The above code uses the NUBEAM auxilliary routines "nbfi_numfast" and
"nbfi_lbl":
subroutine nbfi_numfast(inum_fast)
! return number of NUBEAM fast species...
integer, intent(out) :: inumfast
subroutine nbfi_lbl(lbls,idim,igot)
! return a label for each known fast specie, labels have this form...
! <id> will be H or D or T or He3 or He4 or (rarely) something from
! higher on the periodic table of elements
! <id>_MCBEAM -- Monte Carlo beam specie
! <id>_MCFUSN -- Monte Carlo fusion product specie
! <id>_FPBEAM -- Fokker Planck beam specie (not computed by MC code)
! <id>_FPRFI -- Fokker Planck RF minority ion specie (not by MC code)
! igot gives the number of labels set; igot=-1 if idim is too small
use nbi_com
implicit NONE
integer, intent(in) :: idim ! dimension of lbls character string array
character*(*), dimension(idim), intent(out) :: lbls ! label strings
integer, intent(out) :: igot ! number of label strings returned
........
---------------
The <neutral-gas> label can be generated from the table
Z A label
------------
1 1 "H0"
1 2 "D0"
1 3 "T0"
2 3 "He30"
2 4 "He40"
For which an NTCC library routine "namels" can also be used. Code
which generates and concatenates these label substrings is illustrated
by the example nubeam_driver subroutine nbdrive_summary_2d.f90; other
methods can of course be used.
C-----------------------------------------------------------------------
C NDXFBE -- FIND INDEX TO FAST ION DIST. FCN ENERGY ZONE
C
INTEGER FUNCTION NDXFBE(ZENERGY,ISBEAM,ISTART)
C
C INPUT
C ZENERGY = ENERGY, EV
C ISBEAM = NUBEAM FAST ION SPECIES INDEX
C ISTART = WHERE TO START SEARCH
C
C OUTPUT
C NDXFBE = INDEX TO FAST ION DIST FCN ZONE CONTAINING THE ENERGY
C ZENERGY, OR 0 IF ZENERGY IS NEGATIVE OR .GT. MAXIMUM
C
real, intent(in) :: zEnergy ! Energy, eV
integer, intent(in) :: isbeam ! NUBEAM specie index
integer, intent(in) :: istart ! search index (initial guess)
---> integer function NDXFBE_R8(ZENERGY_R8,ISBEAM,ISTART) also available.
real*8, intent(in) :: zenergy_r8
C--------------------------------------
integer function NDXFBA(ZCOSPA,ISBEAM)
C
C INPUT
C ZCOSPA = cos(pitch angle) = vpll/v relative to direction of Ip.
C ISBEAM = NUBEAM FAST ION SPECIES INDEX (not currently used because
C at present the vpll/v grid is the same for all species).
C
C OUTPUT
C NDXFBA = INDEX TO FAST ION DIST FCN ZONE CONTAINING THE ENERGY
C ZENERGY, OR 0 IF abs(ZCOSPA).gt.1
C
use nbi_com
implicit NONE
C
real, intent(in) :: zcospa ! vpll/v
integer, intent(in) :: isbeam ! specie index (currently not used).
---> integer function NDXFBA_R8(ZCOSPA_R8,ISBEAM,ISTART) also available.
real*8, intent(in) :: zcospa_r8
C-----------------------------------------------------------------------
subroutine GET_FBM_INDEX_INFO(ISBEAM,zEmax,inumE,inumA)
C
C return summary information on energy and pitch grids for
C fast specie distribution function.
C
use nbi_com
implicit NONE
C
integer, intent(in) :: isbeam ! NUBEAM specie index
C
real, intent(out) :: zEmax ! max energy in grid
integer, intent(out) :: inumE ! energy grid size (#zones)
integer, intent(out) :: inumA ! vpll/v grid size (#zones)
C
---> subroutine GET_FBM_INDEX_INFO_R8(ISBEAM,zEmax_r8,inumE,inumA)
also available...
C-----------------------------------------------------------------------
subroutine GET_FBM_EGRID(ISBEAM,zEgrid,imax,igotE)
C
C return energy grid associated with fast ion specie (ISBEAM)
C note-- these are the energy zone BOUNDARIES, and the number of
C boundaries is the number of zones + 1.
C
use nbi_com
implicit NONE
C
integer, intent(in) :: isbeam ! NUBEAM specie index
integer, intent(in) :: imax ! dimension of zEgrid array
real, intent(out) :: zEgrid(imax) ! energy grid
integer, intent(out) :: igotE ! energy grid size (#bdys)
C
C-------------------------------
C
---> subroutine GET_FBM_EGRID_R8(ISBEAM,zEgrid_R8,imax,igotE)
also available...
C-----------------------------------------------------------------------
subroutine GET_FBM_AGRID(ISBEAM,zAgrid,imax,igotA)
C
C return vpll/v grid associated with fast ion specie (ISBEAM)
C note-- these are the vpll/v zone BOUNDARIES, and the number of
C boundaries is the number of zones + 1.
C
use nbi_com
implicit NONE
C
integer, intent(in) :: isbeam ! NUBEAM specie index
integer, intent(in) :: imax ! dimension of zAgrid array
real, intent(out) :: zAgrid(imax) ! pitch angle grid
integer, intent(out) :: igotA ! pitch angle grid size (#bdys)
---> subroutine GET_FBM_AGRID_R8(ISBEAM,zAgrid_R8,imax,igotA)
also available...
NUBEAM supports a method, based on Kadomtsev mixing, to model the
effects of sawteeth. The mixing model is based on the NTCC module
KDSAW. The user requests that sawtooth effects be computed (but
not necessarily applied) by means the following means:
(before NUBEAM timestep is execute):
use nbi_types
type (nbitype_saw) :: zsaw
call nbi_get_saw(zsaw) ! get current sawtooth model settings
zsaw%nlsawb = < TRUE or FALSE > ! flag to compute mixing for beam ions
zsaw%nlsawf = < TRUE or FALSE > ! flag to compute mixing for fusion
! product ions
call nbi_update_saw(zsaw) ! pass information to NUBEAM
If the sawtooth flags are set, then the results of Kadomtsev mixing,
applied at the end of the NUBEAM timestep, are precomputed. However,
they are not APPLIED unless a subsequent call is made:
(after NUBEAM timestep executed, when the sawtooth actually occurs):
CALL SAWNBI
The effect of this call is to apply the precomputed mixing to actually
change the state of the NUBEAM Monte Carlo particle lists, and, change
the NUBEAM internally stored fast ion particle, momentum and energy
density profiles. NOTE-- a direct call to SAWNBI only works if NUBEAM
is being controlled from the subroutine library API. If it is being
run in client-server mode with "nubeam_driver", there is a namelist
switch "sawflag" which is set to .TRUE. to request a sawtooth mixing
operation on fast ions, on the next server call.
The modeling of a sawtooth event typically involves a discontinuous
breaking off and re-initialization of PDEs representing the evolution of the
plasma-- handling the timestepping issues and maintaining constraints such
as quasi-neutrality in the presence of discontinuously changed fast ion
density profiles-- is basically the responsibility of the user's code.
However, if the nlsawb and nlsawf options are set, NUBEAM does output
(via XPLASMA) post-sawtooth fast ion particle density, perpendicular and
parallel energy density, and angular velocity profiles, which can be used
by the caller IF a sawtooth DOES occur.
It is preferable if the beginning and end of NUBEAM timesteps can be
synchronized with sawtooth events, so that the timestep control flow
of the user model would look like this:
(a) use NUBEAM with sawtooth options set, compute fast ion evolution
from TB1 to TB2. **Save** post-sawtooth profiles, mapped to the
user's grid.
(b) evolve plasma from time TB1 to TB2
(c) at time TB2: does sawtooth occur? If not, do nothing; If yes:
1. CALL SAWNBI to update state of NUBEAM model (or set a flag
to cause SAWFLAG=.TRUE. to be set in the namelist of the
next nubeam_driver "proxy_exec" server mode call).
2. change fast ion densities, etc. using data saved from previous
call (a) to NUBEAM.
3. re-initialize PDE's of transport model, in conformance with
data constraints (e.g. Zeff, quasi-neutrality) taking changed
fast ion densities into account.
(d) TB1 <-- TB2, TB2 <-- TB2 + dtBEAM; go back to step (a).
If the user's code has a check-point/restart capability, i.e. if it
has the ability to restart itself from saved files the middle of a
run, e.g. in case of a system crash or for debugging purposes, then,
the NUBEAM state save / state restore feature must also be used. The
NUBEAM procedure for saving state must be synchronized with general
state-save operations of the user's code-- so, the user must choose
when a NUBEAM state save operation is to be done; saving state for
later restart cannot be part of the normal NUBEAM calculational cycle.
The NUBEAM library provides a single simple routine to do state save
and state restore operations:
call nbi_states(<filename>,'S',ierr) ! save state to file
and possibly, only as part of the procedure for restarting an
interrupted user calculation:
call nbi_states(<filename>,'R',ierr) ! restore state from file
In either case, an integer error code "ierr" is returned; ierr=0
means there was no error.
The NUBEAM state file (path) is named by the character string argument
<filename>. The format of the file is portable binary (NetCDF). Its
contents can be examined with `ncdump'.
CAUTION: the method described in this section is only correct if the
traditional subroutine-based call method is used for NUBEAM.
The NUBEAM module is distributed with a standalone test driver,
`nubeam_test'. This driver requires two (2) NetCDF files for
input:
XPLASMA state (just before execution of NUBEAM timestep).
-- MHD equilibrium and input profiles.
-- standard name: <workpath>/<runid>_xplasma_nubeam.cdf
NUBEAM state (just before execution of NUBEAM timestep).
-- NUBEAM state: e.g. accumulated fast ions from prior
timesteps.
-- standard name: <workpath>/<runid>_nbi_start.cdf
These files are generated automatically when NBSTART(...) is called
to initiate a NUBEAM timestep execution. The variables
<workpath>
<runid>
are character string data members of the "nbitype_sys" input block for
NUBEAM.
If NUBEAM crashes unexpectedly in NBSTART, DEPALL, ORBALL or NBFINISH,
the above state files contain all the information necessary to reproduce
the crash in the `nubeam_test' standalone code.
--caution-- some crashes may depend on random number sequence. In this
case, reproducibility in `nubeam_test' is guaranteed only if the exact
same version of code, compilers, etc., are used to build `nubeam_test'
as are used to build the driver program from which the crash occurred.
Even so, the NetCDF files are portable, and it may well be possible with
this data to reproduce a crash elsewhere, e.g. on the workstation of
a PPPL NUBEAM expert.
**IF PPPL debugging assistance is desired**
(a) put the 2 files <runid>_nbi_start.cdf and <runid>_xplasma_nubeam.cdf
on an anonymous ftp server (or email a compressed tar file) so that
a PPPL expert can access the data.
(b) send email to dmccune@pppl.gov, and include the following information:
-- name of tar file and anonymous ftp server.
-- nature of problem requiring debugging (try to be specific).
-- "tail of log file" or other output showing symptoms of problem,
including specific messages from NUBEAM; traceback information
in case of arithmetic error; whatever is available.
NUBEAM runs only in REAL*8 (64 bit) precision. No other precision
choices are available.
NUBEAM is being parallelized-- an MPI parallel "server" version of
NUBEAM will be available soon (written Sept 2005 dmc).
The legacy "dtn" orbit timestep control problem is now fixed; the
code will now adjust itself to find a correct value automatically.
The definition of a modular interface to a Fokker Planck (instead of
Monte Carlo) fast ion slowing down calculation (using NUBEAM for
deposition calculations) has not been completed. In this context,
a physics limitation of NUBEAM should be noted. NUBEAM has no
Monte Carlo RF operator-- therefore, effects of RF resonance on beam
injected ions cannot be calculated except by using a Fokker Planck
package. This has disadvantages, because the Fokker Planck package
makes other approximations (zero orbit width, zero FLR). Development
of an appropriate RF Monte Carlo operator is a difficult computational
research problem.
Statistical Variance of Results --
NUBEAM is a Monte Carlo code. This means (a) its results are statistical
in nature -- Monte Carlos sums, (b) the accuracy of these sums is subject
to statistical variance of magnitude dependent on the number of particles
used in the simulation (variance ~ nptcls**0.5), and (c) small changes to
code inputs produce generally small changes in code outputs, but only
within the accuracy of the inherent statistical variance. These accuracy
limitations have serious implications. For example, when the code is
ported to another machine, the reproducibility of results will generally
be good only within the statistical variance, i.e. not nearly as good as
the arithmetic accuracy of the machines. This is because even small
differences will cause at some point a change of sequence of Monte Carlo
events, after which the association of the random number sequence with
subsequent decision points in the code is shifted-- putting the entire
simulation on a different track. Such variance can be mitigated by
precomputing random number sequences separately for each particle track,
but, this would require major changes to the code and is not planned in
the short term.
The NUBEAM NTCC module comes with two test programs:
nubeam_driver -- NUBEAM demonstration program. This namelist driven
program can run NUBEAM through several timesteps
against a fixed target plasma (only the target ion
density is reduced to make room for the fast ions
as their density increases due to beam deposition).
This program gives a complete, self-contained
example of the set-up of a NUBEAM calculation and
use of its outputs. This program produces ascii
output files that can be compared to reference
versions, e.g. for regression testing.
(nubeam_driver now includes a "proxy_exec" mode
which allows it to perform NUBEAM timesteps as
a service for other codes which provide the
necessary namelist and input files).
nubeam_test -- NUBEAM standalone debugging program. This program
uses "state file" information, and so is "parasitic"
on some other successful NUBEAM installation within
a transport code (or a code like nubeam_driver). Its
purposes is to give a small, standalone tool that
can accurately reproduce problems that might arise
in the use of NUBEAM, for debugging purposes. This
program has a command line interactive interface and
produces NTCC "sglib" graphical output.
These programs can and should be used to test installation of NUBEAM NTCC
distributions. Summary instructions follow.
Recipe:
(a) have NTCC program nubeam_driver built and in your $PATH.
(b) copy the test data (input data and reference output data) from
the NTCC source directory containing nubeam_driver.f90, the
nubeam_driver main program routine, to a work directory.
(c) to exercise the code with an MHD equilibrium from a TRANSP run:
> nubeam_test nubeam
in addition to messages this produces "nubeam_summary.dat" and
"nubeam_details.dat", which can be examined, and, a NetCDF file
"nubeam_xplasma_state.cdf" which can be plotted and compared to
reference data:
> plot_xplasma nubeam_xplasma_state.cdf ! results of this run
> plot_xplasma nubeam_xplasma_state.ref ! reference results
append or type at prompt "@quick_nbplots" to get a standard set of
plots.
(d) to exercise the code with an MHD equilibrium from an EFIT run:
> nubeam_test nubeam2
in addition to messages this produces "nubeam2_summary.dat" and
"nubeam2_details.dat" which can be examined, and, a NetCDF file
"nubeam2_xplasma_state.cdf" which can be plotted and compared to
reference data:
> plot_xplasma nubeam2_xplasma_state.cdf ! results of this run
> plot_xplasma nubeam2_xplasma_state.ref ! reference results
append or type at prompt "@quick_nbplots" to get a standard set of
plots.
(e) for a somewhat longer running test with more statistics:
> nubeam_test nubeam3
in addition to messages this produces "nubeam3_summary.dat" and
"nubeam3_details.dat" which can be examined, and, a NetCDF file
"nubeam3_xplasma_state.cdf" which can be plotted and compared to
reference data:
> plot_xplasma nubeam3_xplasma_state.cdf ! results of this run
> plot_xplasma nubeam3_xplasma_state.ref ! reference results
append or type at prompt "@quick_nbplots" to get a standard set of
plots.
NOTE:
"plot_xplasma" is a test program, an interactive fortran program that
produces xterm/tektronix-compatible plots, which is available with the
xplasma NTCC module.
CAUTION:
Comparison of results to reference files from non-identical platforms
may show significant statistical variation.
nubeam_driver is provided mainly as a coding example for use of NUBEAM.
Recipe:
(a) have NTCC programs nubeam_test and tek2ps built and in your $PATH.
(tek2ps is an sglib utility-- see the NTCC sglib documentation).
(b) copy the test data, scripts, and references results sample_*.* to
a work directory.
(c) run the "sample_run.csh" script.
(d) compare the post-script file "foo.ps" produced with the reference
post-script file "sample_run_output.ps". The results should be
very similar, but, small differences may be visible due to
Monte Carlo statistical variance.
Alternative:
The program can be run interactively in an `xterm' window with
AG4010 tektronics graphics emulation capability. In this case, it
will send plots to the screen, and `tek2ps' is not needed. To
invoke the program, enter
> nubeam_test
at the first program prompt, enter
@sample_run
this will launch a NUBEAM timestep using state file data specified
inside the script. The plots produced can be compared with the
plots contained in "sample_run_output.ps".
Explanation:
The program is distributed with test data: two NetCDF files,
sample_nbi_start.cdf
sample_xplasma_nubeam.cdf
containing the NUBEAM and XPLASMA "state file" information needed to
drive nubeam_test. A pair of files of this type is produced anytime
NUBEAM's NBSTART subroutine is called-- the corresponding filenames
will be
<runid>_nbi_start.cdf
<runid>_xplasma_nubeam.cdf
where <runid> is a NUBEAM control input, a character string element
in the "nbiytpe_sys" data structure.
If built with the same compilers and libraries and run on the same
platform, nubeam_test will precisely reproduce results from NUBEAM
calls in production applications, for which the appropriate state
data has been saved.
The "nubeam_driver" program -- distributed with the NTCC module --
contains a short, well-commented driver program
nubeam_driver/nubeam_driver.f90
which calls the subroutine
nbdrive/nbdrive_main.f90
which is the gateway to a well-commented library of routines
nbdrive/*.f90
which implement the test program. The test program is meant to serve
as an example driver for prospective users of the NUBEAM module.
The subtopics here give summary explanation of the subroutines called
from member routines of the nbdrive library. All calls are resolved
by the NUBEAM module or its dependencies.
The routines nbi_init_*, nbi_set_*, nbi_update*, nbi_get_*, and
nbo_get* are not listed here; these generated interfaces and associated
compound data types are described in Appendix_B and Appendix_C.
The calls for fetching the default values and/or current values, and
for setting or updating the values of inputs to the NUBEAM code are
described in Appendix B. All such routines take an instance of a
fortran-90 compound data type "nbitype_<tname>" as an argument, and
have names of the form:
nbi_init_<tname> -- get <tname> default values.
nbi_set_<tname> -- set <tname> values in NUBEAM internal storage
(can only be done once for many types).
nbi_get_<tname> -- get <tname> currently set values from NUBEAM
internal storage.
and for updatable types,
nbi_update_<tname> -- update (modify) <tname> values inside NUBEAM.
or
nbi_add_<tname> -- add new members to list of values inside NUBEAM.
(the set of neutral beams, impurity species, and RF minority species
are treated as lists, to which members can be added, but once added,
cannot be removed).
The calls for fetching the results of a NUBEAM calculation, including
scalar outputs and arrays of scalars, and XPLASMA ids for radial profiles,
are all described in Appendix C. All such routines take an instance
of a fortran-90 compound data type "nbotype_<tname>" as an argument,
and have names of the form:
nbo_get_<tname> -- retrieve <tname> output data.
These routines are executed only after all inputs have been specified
and all input profiles loaded and interpolated.
for all of these calls, integer ierr contains the status or
completion code; 0 indicates "OK".
call NBSTART(ierr) ! set up timestep; prepare atomic physics
! tables for current plasma parameters, etc.
call DEPALL(ierr) ! beam & fusion product deposition calculation
call ORBALL(iorbtot,ierr)
! orbit calculation; integer iorbtot returns the
! total number of orbits followed.
call NBFINISH(ierr) ! finish up: normalize output profiles and store
! results in XPLASMA; set up output data structures.
After this sequence is complete, "nbo_get_*" calls can be used to fetch
specific output information.
The routines
integer function index_nbfi(iZ,iA,iType)
subroutine nbfi_numfast(inum_fast)
subroutine nbfi_lbl(zlbls,idim,igot)
can be used to find indices and labels for fast species used by NUBEAM;
these are explained fully in the main text under the sections on output.
There is also a set of auxilliary routines for getting information about
the pitch and energy grids of the Monte Carlo fast ion distribution
function outputs -- ndxfbe, etc. -- also described in the section on 2d
outputs.
A note of caution is in order. The pitch and energy grids for the
fast ion distribution functions will not be defined until the first
NUBEAM NBSTART call has successfully completed. The pitch grid is
evenly spaced in (vpll/v) over a known range [-1,1], but the energy
grid is generally not evenly spaced, and its maximum energy can vary
for different fast species.
A method for getting the energy grid in advance, before the first
NUBEAM call, is illustrated in nbdrive/nbdrive_nubeam_checknaml.f90.
Although the maximum energy for distribution functions for beam ions
is controlled by the namelist, the maximum energy for fusion product
ions is determined by a NUBEAM subroutine call (based on knowledge of
the fusion product birth energies):
subroutine nbfusn_emax(iZ,iA,emax)
implicit NONE
integer, intent(in) :: iZ ! Z of fusion product
integer, intent(in) :: iA ! A of fusion product
real, intent(out) :: emax ! a number safely above the birth energy
--or--
subroutine r8_nbfusn_emax(iZ,iA,emax)
implicit NONE
integer, intent(in) :: iZ ! Z of fusion product
integer, intent(in) :: iA ! A of fusion product
real*8, intent(out) :: emax ! a number safely above the birth energy
(These calls can be made any time, e.g. before the main NUBEAM
calculation has been initialized).
when the energy grid size and maximum energy have been set, and arrays
to hold the energy grids allocated, the following NUBEAM call can fill
in the energy grids, which will then match the grids NUBEAM uses
internally to calculate the distribution functions:
subroutine nbi_efbm(ilfprod,inum,zemax,ilin,zelin,efbmb,efbm)
!
! compute the energy grid for a MC fast ion specie dist. fcn.
! (code taken from old trcore/datckb.for -- dmc 10/04/01)
implicit NONE
logical,intent(in) :: ilfprod ! .TRUE. if this is a fusion product specie
integer,intent(in) :: inum ! total no. of zones
real,intent(in) :: zemax ! max energy
integer,intent(out) :: ilin ! no. of linearly spaced zones; ilin.ge.1
real,intent(out) :: zelin ! max. linearly spaced energy (min is zero)
real,intent(out) :: efbmb(inum+1) ! energy grid bdys
real,intent(out) :: efbm(inum) ! energy grid zone ctrs
(and there is also the real*8 precision version, subroutine r8_nbi_efbm).
Take note that the "energy grid bdys" array has one more element than
the "energy grid zone ctrs" array.
These calls are described in detail in the XPLASMA NTCC module
documentation,
http://w3.pppl.gov/~pshare/help/xplasma.htm
and so their usage is merely summarized here.
call eqm_msgs(...) ! redirect XPLASMA messages
call eq_save(filename,ierr) ! save XPLASMA state to a NetCDF file.
call eq_restore(filename, ierr) ! restore XPLASMA state from a
! NetCDF file written by a prior eq_save
! call; any existing state is over-written.
call eqm_fromgeqdsk(...) ! build complete XPLASMA equilibrium,
! fields, and vacuum region description
! using G-eqdsk file or MDSplus path;
! any existing state is over-written.
(if eqm_fromgeqdsk is not used, the following calls will be used
to feed in data from a TRANSP run; this also serves as an example
usable by any code with an axisymmetric MHD equilibrium with
geometry of the form R(rho,theta),Z(rho,theta)).
call eqm_select(...) ! re-initialize XPLASMA, in preparation
! for new MHD equilibrium data
call eqm_rho(...) ! define basic rho grid
call eqm_chi(...) ! define poloidal angle grid
! with apologies:
! sometimes "chi" and sometimes
! "theta" is used to name this
! coordinate.
call eqm_uaxis(...) ! define a user grid -- typically with
! an association to the basic rho grid.
call eqm_rhofun(...) ! define a named XPLASMA profile function
! f(rho). A unique XPLASMA id is returned.
call r4_eqm_rhofun(...) ! as eqm_rhofun, with REAL instead of
! REAL*8 precision arguments.
call r4_eqm_rzmag(...) ! specify MHD equilibrium geometry
call eqm_bset(...) ! specify orientation (sign) of toroidal
! and poloidal magnetic fields; the fields
! are calculated if sufficient data has
! been provided.
call eqm_dbdy_grid(...) ! specify vacuum vessel wall at fixed
! distance from plasma boundary, and
! create vacuum region [R,Z] grid
external eqm_brz_adhoc
call eqm_brz(eqm_brz_adhoc,...) ! create adhoc extrapolation for
! B fields and pseudo-flux-coordinates
! for the vacuum region.
call eq_rholim(...) ! return limits of XPLASMA rho grid --
! in practice, rho_min=0, rho_max=1 always.
call eq_glimrz(...) ! find Rmin,Rmax,Zmin,Zmax of a flux surface.
call eq_rzget(...) ! [mag. coords] -> [cyl. coords] map
call eq_inv(...) ! [cyl. coords] -> [mag. coords] inverse map
call eq_volume(...) ! volume (m**3) enclosed by flux surfaces.
call eq_area(...) ! cross-sectional area (m**2) enclosed by
! flux surfaces.
call eq_rgetf(...) ! return interpolated values of XPLASMA
! profile f(rho); XPLASMA id is input.
call eq_ganum(axis_name,id) ! get XPLASMA axis id from name
Note: call eq_ganum('__RHO',id_rho) fetches the id for the
XPLASMA rho or sqrt-normalized-toroidal-flux grid. This
is useful for allowing the user to define his/her own
rho grid (with an eqm_uaxis call) and to have it associated
with the XPLASMA rho grid. User rho grids need not be
evenly spaced.
call eq_gfnum(fcn_name,id) ! get XPLASMA profile id from name
Note: usually for f(rho) profiles in NUBEAM context.
call mcgrid_define(...) ! create a MCgrid; return a MCgrid id.
call mcgrid_getsym(...) ! return updown symmetry of MCgrid.
call mcgrid_getnumr(...) ! return no. of zone rows of a MCgrid.
call mcgrid_getnumz(...) ! return total no. of zones in a MCgrid.
call mcgrid_getnumzns(...) ! return array of no. of poloidal zones
! per radial zone row
call mcgrid_getobj(...) ! get data defined over a MCgrid.
! (it was put there by NUBEAM calls to
! mcgrid_putobj(...)
call mcgrid_getarr2(...) ! get array data over a MCgrid: the
! fast ion distribution function.
set random number seed from system clock
(used by nbdrive if the seed is defaulted in the namelist).
integer function irns(idummy)
C******************** START FILE IRNS.FOR ; GROUP RNG ******************
FUNCTION IRNS(IDUMMY)
C 11/29/78
C THIS FUNCTION WAS WRITTEN BY HARRY H. TOWNER OF THE PRINCETON PLASMA
C PHYSICS LAB. IRNS WILL RETURN A UNIQUE INTEGER SUITABLE FOR SETTING
C THE RANDOM NUMBER SEED.
c 5/01 RGA -- replace with f90 standard call
integer iv(8)
call date_and_time(values=iv)
secnds = 60*(60*iv(5)+iv(6))+iv(7)+iv(8)/1000.
IRNS=1.E4*SECNDS
IRNS = 2*IRNS + 1 !MAKE SURE IT'S ODD, RTM 26 FEB 86
RETURN
END
C******************** END FILE IRNS.FOR ; GROUP RNG ******************
subroutine namels(iZ,iA,cname)
integer, intent(in) :: iZ,iA
character*(*), intent(out) :: cname
return short string to identify an isotope:
in out
iZ iA
1 1 "H"
1 2 "D"
1 3 "T"
2 3 "He3"
2 4 "He4"
higher values of iZ and iA are also allowed, but typically not needed
in the NUBEAM context.
A small set of routines are used to categorize and label fusion reactions
that can be computed by the NUBEAM model. The fusion rate data going into
these estimates comes from the PREACT NTCC module.
Experimentally, it is of considerable interest to estimate not just
the total neutrons, or, the total number of fusion reactions on a
reaction-by-reaction basis, but also, the break-down of rates based
on the reagent types (fast ions vs. thermal ions). That was the
motivation for the "catfus" categorization scheme-- a TRANSP legacy.
The catfus* routines are useful for indexing into the fusion related
output data of NUBEAM (or in the case of 2d outputs, generating the
necessary character label for the desired MCgrid data object), as
shown in the example code, below.
The catfus routines use a category index to distinguish among reactions
and partially among reagent types. Then, within each category index
there are two "reagent type" codes specifying:
code meaning
0 thermal ion
1 beam ion
2 fusion product ion
3 RF tail ion
For each category, reagent#1 can be either thermal or one of {1,2,3};
reagent#2 can be either thermal or {1: beam ion}.
This allows a separate category index / reagent code combinations for
the following fusion reactions
DD:P DD:N D-T D-He3 TT
(T-He3 category also exists but cross-section data is suspect at
this time).
with the following reagent combinations ("target" means thermal ion
target):
thermonuclear (indices 1-6)
beam-target (indices 1-8)
beam-beam (indices 1-6)
(fusion_product)-target (indices 9-10: T & He3 fusion products)
(fusion_product)-beam (indices 9-10: with D beams)
(minority_tail_ion)-target (indices 11-14)
(minority_tail_ion)-beam (indices 11-14)
(erratum: indices 13 & 14 list "D" as an RF minority which is
unlikely to occur, experimentally).
Certain fast ion -- fast ion combinations are not supported. For
example, the fusion_product - fusion_product T-He3 fusion in a D beam
D thermal target T/He3 fusion product burnup simulation -- known to
be small -- cannot be calculated by NUBEAM at present.
To avoid the situation where a physical reaction would be associated
with more than one index / reagent-type combination, some categories
require that the first reagent be non-thermal, and/or that the
second reagent be thermal. For example, the fusion rates associated
with D beams interacting with T thermal plasma, i.e. beam-target D->T
(index 1) is calculated separately from the rates associated with
T beams interacting with D thermal plasma, i.e. beam-target T->D
(index 7). But the thermonuclear and beam-beam DT fusion reactions
are symmetric with respect to reagent type, and so are reserved to
category index 1.
The code
use nbi_types
type(nbotype_fusion) :: zfusion
integer i,inum_react,i1,i2
character*20 label
call nbo_get_fusion(zfusion) ! get NUBEAM fusion rate results
call catfus_num(inum_react) ! get number of known reactions
do i=1,inum_react
call catfus_mtype(i,i1) ! fetch available fast ion type, this index
do i2=0,1
! 0: reagent 2 is thermal; 1: reagent 2 is beam ion
call catfus_name(i,i1,i2,label)
if(label.ne.'error') then
! read fusion data vs (rho,theta) for this reaction/reagent
! type combination ... use -id_mcgrid, so that if the item
! is not present e.g. because not all reagents are present,
! the result is silently set to zero.
label = trim(label)//'_2d' ! append this suffix
call mcgrid_getobj(-id_mcgrid,label,...)
if (i2.eq.0) then
write(6,*) label,' fusion reaction with thermal ions', &
zfusion%btnta(i),' per second.'
else
write(6,*) label,' fusion reaction with beam ions', &
zfusion%bbnta(i),' per second.'
endif
else
! the answer is "zero"
endif
enddo
enddo
The index "i" is also to be used to specify the "reaction index" in
array data members of the nbotype_fusion data type (described in the
appendix showing NUBEAM output data structures). Data members with
the substring "bt" refer to "beam-target", or more generally, "fast
ion -- thermal target" reactions; data members with the substring "bb"
refer to "beam-beam", or, more generally, "fast ion -- fast ion" fusion
reactions. The above example code illustrates the point.
subroutine catfus_num(inum_react)
integer, intent(out) :: inum_react ! no. of fusion reaction categories
Get total number of fusion reaction categories known by NUBEAM.
subroutine returns a fusion reaction category index, based on the
input description.
subroutine catfus_btindx(c3beam,c3target,itypbeam,indx)
character*(*), intent(in) :: c3beam ! "D" "T" or "HE3"
character*(*), intent(in) :: c3target ! "D:P" "D:N" "T" or "HE3"
integer, intent(in) :: itypbeam ! 1 for beam, 2 for fusion prod, 3 for
! RF tail; 0 for thermal
integer, intent(out) :: indx ! reaction index (or 0 if none found)
(note the target type has not been specified -- there are two
possibilities: beam fast ions, or thermal ions).
given the reaction index and reagent types, generate a label string
for the fusion reaction.
In the NUBEAM context, this call generates the label needed to name
and fetch a profile of 2d (rho,theta) resolved fusion rate data--
example in nbdrive/nbdrive_summary_2d.f90 ...
subroutine catfus_name(ireact,irtyp1,irtyp2,name)
integer,intent(in) :: ireact ! reaction index
integer, intent(in) :: irtyp1 ! type of 1st reagent
integer, intent(in) :: irtyp2 ! type of 2nd reagent
! rule: irtyp1.ge.irtyp2 .and. irtyp2.le.1 .and. irtyp2.ge.0
! additional rules: see the module data
! generally, reaction indices .gt.8 are for FP and RFI ions only.
!
! irtyp1|2 = 0 -- thermal irtyp1=0, irtyp2=0: thermonuclear reaction
!
! in what follows, "target" means thermal target population.
!
! irtyp1|2 = 1 -- beam ion irtyp1=1, irtyp2=0: beam-target reaction
! irtyp1=1, irtyp2=1: beam-beam reaction
!
! irtyp1|2 = 2 -- fusion product (FP): T or He3
! irtyp1=2, irtyp2=0: FP-target reaction
! irtyp1=2, irtyp2=1: FP-beam reaction
!
! (D is not a fusion product; P and He4 are non-reacting fusion products)
!
! irtyp1|2 = 3 -- RF tail ion (RFI)
! irtyp1=3, irtyp2=0: RFI-target reaction
! irtyp1=3, irtyp2=1: RFI-beam reaction
! (generally, T or He3 tail ions are of interest).
character*(*), intent(out) :: name ! name (label string) for reaction
The name string "error" is returned if any of the input arguments
are out of the expected range, or if the reagent types specified
are not supported for the given reaction index.
given a reaction index, return the type of fast species available
for the 1st reagent; the 2nd reagent can be a beam ion or thermal
target plasma species.
subroutine catfus_mtype(ireact,itypfast)
integer,intent(in) :: ireact ! reaction index
integer,intent(out) :: itypfast ! fast-ion type
! 1=beam, 2=fusion product, 3=RF tail
! generally, the categorization allows fast-ion -- fast-ion reactions
! to be considered, but the "second" fast-ion specie must be a beam
! ion specie; the "first" specie can be beam-ion (B),
! fusion-product-ion (FP), or RF-tail-ion (RFI).
!
! categories for FP-FP, RFI-RFI, and FP-RFI are not supported;
! categories B-B, FP-B, and RFI-B *are* supported.
!
return in character format reagent isotope information
subroutine catfus_reagent(ireact,reagent1,reagent2)
integer,intent(in) :: ireact ! reaction index
character*(*),intent(out) :: reagent1,reagent2 ! reagents
Typical values returned: "D", "T", "HE3". For the DD reaction, the
second reagent can be "D:P" or "D:N" indicating the proton or neutron
branch of that fusion reaction.
This is the f90 module data (set of constants) used by the subroutines
catfus*...
module catfus_mod
! Fusion reaction indexing data
! the categorization comes originally from TRANSP at PPPL
! ( http://w3.pppl.gov/transp)
integer, parameter :: nreact=14 ! no. of indexed reactions
! reagents isotope species:
character*3, parameter, dimension(nreact) :: r1 = &
(/ 'D ','D ','D ','D ','T ','T ','T ','HE3', &
'T ','HE3','T ','HE3','D ','D ' /)
character*3, parameter, dimension(nreact) :: r2 = &
(/ 'T ','HE3','D:P','D:N','T ','HE3','D ','D ', &
'D ','D ','D ','D ','T ','HE3' /)
! reagent categories:
! 0->thermal, 1->beam, 2->fusion product, 3->RF tail
! valid reagent1 categories, by reaction index (lower indicies
! allow beam or thermal;
! higher indices are special purpose, generally speaking).
integer, parameter, dimension(nreact) :: typmin = &
(/ 0, 0, 0, 0, 0, 0, 1, 1, &
2, 2, 3, 3, 3, 3 /)
integer, parameter, dimension(nreact) :: typmax = &
(/ 1, 1, 1, 1, 1, 1, 1, 1, &
2, 2, 3, 3, 3, 3 /)
! maximum valid reagent2 categories (1 means beam or thermal;
! 0 means thermal only).
integer, parameter, dimension(nreact) :: typ2max = &
(/ 1, 1, 1, 1, 1, 1, 0, 0, &
1, 1, 1, 1, 1, 1 /)
end module catfus_mod
These calls are described in detail in the TRREAD NTCC module
documentation,
http://w3.pppl.gov/~pshare/help/trread_hlp.html
and so their usage is merely summarized here.
These routines are used in nbdrive_mhdeq.f90, only if the nbdrive
namelist specifies a TRANSP run as the source for MHD equilibrium
data to be used in the upcoming NUBEAM simulation.
call plc_msgs(lun,' ') ! set FORTRAN i/o unit number for TRREAD messages
call rconnect(runpath, ierr) ! connect to TRANSP output database
call rplabel(...) ! get information on a TRANSP function
call rpdims(...) ! get information on a radial TRANSP grid
call t1mhdeq(...) ! read the TRANSP MHD equilibrium at a
! selected time...
call t1scalar(...) ! read a TRANSP scalar function value at
! a selected time...
human readers need not be concerned...
maintainers: do not change this line unless you know what your are doing:
!-->PYTHON_START
! (nbigen.py: python generated .hlp code, do not edit)
Array dimensions must be set before any other
input structure.
Methods and structure contents are shown, with
comments copied from nubeam/nbspec.dat ...
To define the referenced fortran-90 structures,
user code should include the statement:
use nbi_types
use nbi_types
type(nbitype_dims) :: zdims ! instantiate object
integer :: ierr ! status code, 0=OK
call nbi_init_dims(zdims) ! set defaults
[...user code...] ! modify elements...
! this can only be done ONCE:
call nbi_set_dims(zdims,ierr) ! send to NUBEAM
if(ierr.ne.0) [...handle error...]
... ... ... ...
call nbi_get_dims(zdims,ierr) ! retrieve from NUBEAM
if(ierr.ne.0) [...handle error...]
! note: most nbi_get methods do not return ierr.
! nbi_get_dims can be used to check if
! NUBEAM dimensions have been set.
members of nbitype_dims data type:
integer :: mib
! -> default value: 32
! -> maximum user value allowed: 200
! = nbeam
! max number of neutral beams "nbeam" in simulation.
! beams may be added during a simulation but not removed.
integer :: mj
! -> default value: 55
! = nzones + 2
! mj is related to the max number of radial zones in the
! simulation: mj must exceed <desired-number-of-zones> + 2
integer :: mig
! -> default value: 7
! -> maximum user value allowed: 10
! = ngmax
! max number of non-impurity thermal plasma species.
! the number of knownn non-impurity species is "ngmax";
! ngmax can increase up to the max number, mig, during a run.
! (inside the Monte Carlo code, the number of species with
! non-zero density is ng, ng.le.ngmax at all times).
integer :: mis
! -> default value: 21
! -> maximum user value allowed: 50
! = nsmax
! number of externally supplied thermal neutral sources:
! leading to neutral densities & temperatures for CX evaluation
integer :: mibs
! -> default value: 7
! -> maximum user value allowed: 15
! = nsfast
! max total number of fast ion species including beam ions,
! fusion product ions, and RF minority ions
integer :: mimpt
! -> default value: 100
! -> maximum user value allowed: 100
! = nrhix
! max number of impurity thermal plasma species "nrhix".
! nrhix can increase up to the max number during a run.
integer :: mmini
! -> default value: 3
! -> maximum user value allowed: 10
! = nmini
! max number of RF minority species; the actual number "nmini"
! can increase up to the max number during a run.
integer :: mtheta
! -> default value: 64
! = nthets
! number of poloidal grid points in bicubic spline
! representations of [1/B] and [J] used in orbiting
integer :: miefi
! -> default value: 5
! -> maximum user value allowed: 10
! = nerngfi
! coarse energy grid (cf erngfi(1:nerngfi)) for select
! output profiles
integer :: mimxba
! -> default value: 50
! = nznbma
! number of pitch angle zones (evenly spaced in vpll/v) in
! lab frame fast ion distribution functions.
! note: pitch angle grid the same for all species
integer :: mimxbe
! -> default value: 100
! = nznbme
! maximum number of energy zones in fast ion distribution
! functions; actual number and zone spacing can vary from
! fast specie to fast specie.
! definition of spatial zones: see xplasma MC grid.
integer :: miboxm
! -> default value: 12
! -> maximum user value allowed: 100
! = nbbox
! max number of "beam in box" 3d neutral density profiles to
! compute-- actual number "nbbox" can vary from timestep
! to timestep
integer :: mitubem
! -> default value: 3
! -> maximum user value allowed: 200
! = nbtube
! max number of "beam in tube" 3d neutral density profiles to
! compute
integer :: mituben
! -> default value: 100
! -> maximum user value allowed: 500
! = maxval(NSEGtube)
! max number of segments per tube for the "beam in tube"
use nbi_types
type(nbitype_sys) :: zsys ! instantiate object
integer :: ierr ! status code, 0=OK
call nbi_init_sys(zsys) ! set defaults
[...user code...] ! modify elements...
! this can only be done ONCE:
call nbi_set_sys(zsys,ierr) ! send to NUBEAM
if(ierr.ne.0) [...handle error...]
... ... ... ...
call nbi_get_sys(zsys) ! retrieve from NUBEAM
members of nbitype_sys data type:
! basic system information for Monte Carlo code
CHARACTER*32 :: runid
! --> default value: "NUBEAM"
! runid string -- alphanumeric string used as prefix for
! output files (e.g. state files, etc.). Should generally
! be set to a runid string associated with the driver code.
CHARACTER*120 :: workpath
! --> default value: " "
! work directory path -- for debug statefile...
! blank means current working directory
INTEGER :: nonlin
! --> default value: 6
! unit number for messages.
INTEGER :: nseed
! --> default value: 1
! random number seed; default not recommended!
LOGICAL :: nlbout
! --> default value: .FALSE.
! ***DISABLED*** DMC Apr 2007 -- incompatible with
! mpi parallelization
! (NUBEAM maintainer will evaluate cost of
! restoration of this feature, if requested).
! .TRUE. to turn on "lost particles" ascii output file
! filename will be RUNID(1:LRUNID)//'NB.OUT' where LRUNID
! is the nonblank length of RUNID.
INTEGER :: lunnbx
! --> default value: 98
! unit number for "lost particles" ascii output file
! (not used unless switch NLBOUT is set .TRUE.; seldom used).
! (will be used if ref_namelist=T is set, see below).
INTEGER :: lunres
! --> default value: 99
! auxilliary unit number for ascii file i/o -- supports
! debug restarts and restarts of "lunnbx" file.
! **** this LUN will be used **** even if .NOT.NLBOUT.
INTEGER :: mrstrt
! --> default value: 1
! set to zero if there are *no* check-point / restart files
! (a property of the driver code).
! this changes the way the NB.OUT file is handled, and may
! have other effects in the future.
INTEGER :: nclass
! --> default value: 0
! set =1 to activate beam-by-beam full:half:third
! energy species resolved profiles of n, Pbe, Pbi,
! <v.B>, <(vcrit**3/v**3)v.B> (lab frame v.B) for NCLASS
! fast ion interface -- see nclass output block.
! *** warning *** dmc Dec 2003 -- definition and
! construction of nclass=1 outputs is not yet complete ***
LOGICAL :: only_io
! --> default value: .FALSE.
! set .TRUE. when calling NUBEAM from a serial driver that
! will use a separate MPI-parallel process to compute the
! fast ion evolution. In this mode, only the input and
! output structures are allocated-- NOT the large internal
! computational structures. The following routines are
! NOT available in this mode: NBI_INTERP_PROFILES,
! NBSTART, DEPALL, ORBALL, FI_FINISH; instead the timestep
! advance sequence is:
! CALL NUBEAM_WRITE_INPUT(...)
! (invoke local procedure for calling parallel NUBEAM
! as separate process)
! CALL NUBEAM_READ_OUTPUT(...)
! these calls are more fully explained elsewhere.
LOGICAL :: ref_namelist
! --> default value: .FALSE.
! set .TRUE. to write a reference namelist to enable use of
! "nubeam_driver" (and nbdrive subroutine library) to produce
! a standalone calculation close to the one actually being
! performed, provided the profiles used to drive this step
! can be acquired independently. The namelist is written as
! if the profiles are available in a TRANSP "NetCDF" output
! dataset, but, nubeam_driver supports other options. The
! ascii namelist file <runid>_nbdrive_namelist.dat, written
! on fortran LUN LUNNBX, is commented automatically and can
! be edited with a text editor.
use nbi_types
type(nbitype_times) :: ztimes ! instantiate object
integer :: ierr ! status code, 0=OK
call nbi_init_times(ztimes) ! set defaults
[...user code...] ! modify elements...
! (nbi_set method not recommended)
... ... ... ...
call nbi_get_times(ztimes) ! retrieve from NUBEAM
! usually, all elements are reset every timestep:
[...user code...] ! modify elements...
call nbi_update_times(ztimes) ! send to NUBEAM
members of nbitype_times data type:
REAL*8 :: tbm1 ! *updatable*
! --> default value: 0.0d0
! start time of current timestep -- seconds
REAL*8 :: tbm2 ! *updatable*
! --> default value: 0.0d0
! stop time of current timestep -- seconds
! the start time of timestep j+1 should match
! the stop time of timestep j (else a warning will be issued).
use nbi_types
type(nbitype_grid) :: zgrid ! instantiate object
integer :: ierr ! status code, 0=OK
call nbi_init_grid(zgrid) ! set defaults
[...user code...] ! modify elements...
! this can only be done ONCE:
call nbi_set_grid(zgrid,ierr) ! send to NUBEAM
if(ierr.ne.0) [...handle error...]
... ... ... ...
call nbi_get_grid(zgrid) ! retrieve from NUBEAM
! the following should be preceeded by a nbi_get call.
! only items marked as updatable are affected
! inside NUBEAM by this call (see element description)
[...user code...] ! modify elements...
call nbi_update_grid(zgrid) ! send to NUBEAM
members of nbitype_grid data type:
! basic grid information
INTEGER :: id_mcgrid ! *updatable*
! --> default value: 0
! xplasma id of 2d irregular "Monte Carlo" spatial grid.
! **must be supplied each time and shape must not change**
INTEGER :: nzones
! --> default value: 0
! number of radial zones (equally spaced in sqrt(tor.flux))
! this determines the grid for calculation of 1d outputs
! such as Pbe and Pbi; must be set to a number .le. mj-3
! where "mj" is a user supplied array dimension.
INTEGER :: nthets
! --> default value: mtheta
! number of poloidal angle grid points in [J] and [1/B] 2d
! splines (cf nubeam/magcf.for)
INTEGER :: nznbma
! --> default value: mimxba
! number of pitch angle zones for distribution functions
! this is the same for all fusion product species
INTEGER :: nznbme
! --> default value: mimxbe
! number of energy zones for beam ion distribution functions
! note for fusion product species the number of energy zones
! is always equal to the maximum, mimxbe, regardless of the
! nznbme specification for beam species.
REAL*8 :: ebdmax
! --> default value: 0.0
! maximum energy in distribution function for beam ion species
! the maximum for fusion product species is set by the fusion
! product source energy (e.g. 3.5e6 eV for DT alphas-- adjusted
! upward to 5.0e6 allow for possible energy diffusion).
! **energy in eV**
use nbi_types
type(nbitype_beams) :: zbeams ! instantiate object
integer :: ierr ! status code, 0=OK
call nbi_init_beams(zbeams) ! set defaults
[...user code...] ! modify elements...
! set integer values to match the number
! corresponding list elements to *add* on this call
call nbi_add_beams(zbeams,ierr) ! append to NUBEAM list(s)
if(ierr.ne.0) [...handle error...]
... ... ... ...
call nbi_get_beams(zbeams) ! retrieve from NUBEAM
members of nbitype_beams data type:
! here *both* the beams and the thermal plasma species
! are defined together, because, for each beam specie
! there must be a thermal plasma specie for it to
! thermalize into.
!
! except where noted, actual values must be given...
! in the course of a run, beams and species can be added
! but once added can never taken away.
INTEGER :: ngmax ! *extensible list*
! --> default value: 0
! the number of thermal plasma species
REAL*8 :: backz(10) ! (1:ngmax) *extensible list*
! --> default value: 0.0
! the Z of each thermal specie (atomic no.)
REAL*8 :: aplasm(10) ! (1:ngmax) *extensible list*
! --> default value: 0.0
! the A of each thermal specie (AMU)
INTEGER :: nbeam ! *extensible list*
! --> default value: 0
! the number of distinct neutral beam sources
REAL*8 :: xzbeama(200) ! (1:nbeam) *extensible list*
! --> default value: 0.0
! Z of injected ion specie (atomic no.)
REAL*8 :: abeama(200) ! (1:nbeam) *extensible list*
! --> default value: 0.0
! A of injected ion specie (AMU)
LOGICAL :: nlco(200) ! (1:nbeam) *extensible list*
! --> default value: .TRUE.
! .TRUE. if injected "co", i.e. in direction
! of the toroidal plasma current.
INTEGER :: ntrace(200) ! (1:nbeam) *extensible list*
! --> default value: 0
! non-zero value means this beam (j) is a
! "trace element" of the beam indicated by
! the value of ntrace(j); it inherits the
! geometry of that beam and the interpretation
! of energy fractions is affected.
! (of course, usually, the default ntrace=0 is OK)
REAL*8 :: rtcena(200) ! (1:nbeam) *extensible list*
! --> default value: 0.0