Figure
User's Guide for DEGAS 2
Release V. 4.6^{1}
Daren Stotler and Charles Karney
Princeton Plasma Physics
Laboratory
Randall Kanzleiter
Rensselaer Polytechnic Institute
Jaishankar S.
Institute for Plasma Research, India
Abstract
This is the user's manual for DEGAS 2  A
Monte Carlo code for the study of neutral atom and molecular transport
in confined plasmas. It is intended to provide an
introduction to DEGAS 2 from the
user's pointofview: obtaining
the files, compiling them, setting up the input files, executing the
code, and interpreting the output.
References will be provided for the underlying physics, but the
essential aspects will be highlighted.
More detailed documentation of the procedures used to
create input files is contained in the typeset
source code of the relevant preprocessors; links
to those files are provided.
Chapter 1
Introduction
1.1 Purpose of Monte Carlo Neutral Transport
The interaction of neutral atoms and molecules with the background
plasma is an important component in the physics of fusion research
plasmas. These
neutrals affect both the particle and energy balance of the plasma,
providing a source of new plasma and a channel for heat transport across
the magnetic field that confines the charged particles. In addition, hot
neutrals which leave the plasma volume can interact with the wall of the
plasma chamber, sputtering impurities into the plasma and (in reactors)
possibly damaging the wall. Finally, these hot neutrals can also be used
for plasma diagnostics. Therefore, it is worth developing an accurate
description of the neutral behaviour.
The theory of neutral particle kinetics[1]
treats the transport of mass, momentum, and energy in a plasma due to
neutral particles that are themselves unaffected by magnetic
fields. This transport affects the global power and particle balances in
fusion devices, as well as profile control and plasma confinement
quality, particle and energy fluxes onto device components, performance
of pumping systems, and the design of diagnostics and the interpretation
of their measurements. Though analytic models (solving the
Boltzmann equation for the neutrals or using the diffusion approximation),
both single and multispecies, have been used to study neutral transport,
a variety of approximations must be made. Some numerical methods for solving
the Boltzmann equation use simplified geometries and atomic physics. On
the other hand, the Monte Carlo approach to integrating the Boltzmann
equation can treat in great detail asymmetric geometries and complex
atomic physics and wall interactions. Because of the similarity in the
problems, Monte Carlo neutral transport codes can build on the techniques
developed for neutron transport in reactor design[2].
1.2 Historical Background
The predecessor to DEGAS 2 was DEGAS[3]. The origin of the
"DEGAS" name can best be explained by pointing out that the first
Monte Carlo neutral transport code written by D. B. Heifetz was
named SEURAT. Its name was a direct reference to the artist Georges
Seurat. This appellation was inspired by the analogy between the
artist's pointillist style and the nature of the Monte Carlo algorithm.
Up close, one sees that Seurat's paintings consist of a huge number of very
tiny dots of paint of a variety of colors. Only when viewed from a distance
is the image and its coloring clearly apparent. In much the same way,
the track of an individual neutral flight in the Monte Carlo code
gives relatively little insight into the problem being studied. Yet, when the
results of many flights are compiled together, an accurate solution can
be obtained.
Heifetz' fondness for Impressionist artists was also manifested
in succeeding codes like MONET and DEGAS. The latter name was further
motivated by two associated puns. One is that "DEGAS" is an
abbreviation of "DEuterium GAS"; DEGAS was written to satisfy a
request for a code that could treat deuterium as well as hydrogen.
The other is that "DEGAS" is just short for "degassing"; another
proposed application for the code was in assessing measurements of
wall degassing in fusion devices.
For those who insist that a code's name must be an acronym, we will
reply that "DEGAS" stands for Divertor and Edge Gas Analysis System.
By the time the authors began contemplating the concepts that led
to DEGAS 2, the original DEGAS code was in widespread use. It also
enjoyed a very good reputation as being easy to use and well documented.
For this reason, we were reluctant to abandon the DEGAS name, even
though we knew we would be rewriting the code largely from scratch.
Nonetheless, much of the atomic and surface physics from
the original code has been carried over. Likewise, many of the
numerical techniques of DEGAS served as inspiration for
those in DEGAS 2.
1.3 Need for DEGAS 2
Two events pointed out limitations in the original DEGAS code[3].
 The combined DEGASB2 code was too slow
(about 100 Cray hours)[4].
The stronger the coupling between neutral and plasma species, the more
difficult it is to arrive at a converged solution. In this case which
looked at a high recycling
divertor for a compact, high field reactor, hundreds or thousands of complete
Monte Carlo profiles must be carried out in order to couple effectively
to a fluid code.
 The problem physics was difficult to modify. In the radiative and
detached divertor
regimes being investigated for future reactors, it is believed that the list
of important atomic processes will grow to include elastic collisions and
molecular hydrogenic species other than ground state H_{2}.
Hence, it was decided that a new code should be written which draws
on the extensive experience at PPPL and that from Garching and
Jülich.
The philosophy behind DEGAS 2 was defined during 19921993:
 High speed,
 Flexibility,
 Easeofuse,
 Well documented.
1.4 DEGAS 2 Features
An effort was made to realise the above mentioned objectives of
DEGAS 2 by the following means:
 Improve speed relative to DEGAS by utilizing the
"track length" estimator.
 Ensure flexibility by designing the code using quasiobject oriented
programming techniques. The modifier "quasi" refers to the fact that DEGAS 2
is written in FORTRAN (compatible with FORTRAN 77 and FORTRAN 90). Even though
the basic "classes" in the problem have been abstracted, the degree of
encapsulation actually achieved (via preprocessor macros) in DEGAS 2 is
limited. Furthermore, there is no provision for inheritance within this
programming environment.
 Provide an easytouse and well documented code by writing it with the
FWEB[8] package. FWEB provides the powerful preprocessing and
macro facility which underlies the "quasiobject oriented" implementation
and internal T_{E}X documentation.
Other significant features of DEGAS 2 include:
 Dynamic memory allocation for optimal performance in large problems,
 Operation on a number of UNIX platforms, compilable with both
FORTRAN 77 and FORTRAN 90,
 Input and output data stored in selfdescribing, platformindependent
binary files.
Presently, the main code uses exclusively
netCDF[9] (Network Common Data
Format) files for input and output. However, a simple postprocessing
utility is included which generates graphics and data in
HDF[10] (Hierarchial Data Format) or Silo[11] files.
Hopefully, a single format will be settled upon eventually.
 Parallel execution on distributed workstations and massively
parallel computers. This is the most promising
technique for improving the speed of DEGAS 2. The reason is that most of
the future advancements in computer speed will result from increases in the
number of processors. Parallel processing is practical to Monte Carlo since
the amount of communication required per computation is small.
 High quality coding standards.
Good programming practices such as use of the implicit none statement
and an assertion facility help to ensure that the code works properly.
This document will eventually describe a basic test suite which actually
demonstrates that the code is correct. In the code design process, clarity
has been given a high priority as well, facilitating readability and
maintainability.
 Code versions numbered and documented with the "CVS" (concurrent
version system) package.
Chapter 2
Background
2.1 Generic Monte Carlo Algorithms
Monte Carlo methods generally involve sampling from some "parent population".
The sampling process is usually facilitated by pseudorandom numbers generated
via a numerical algorithm. The quantity or quantities sampled are then
used as input to a subsequent (deterministic) calculation which results in the
desired output. This process is repeated a number of times, so as to
thoroughly sample the parent population, and a distribution of the output
data is compiled. In some cases, only the mean of this distribution is of
interest (as well as some estimate of its error), while in other cases, the
distribution itself is the objective.
Note that the process being simulated need not be random. For example, Monte
Carlo techniques may be the only way of computing the volume of some complex
threedimensional object. This is just one example of the versatility of
Monte Carlo techniques. But, Monte Carlo is not always the most efficient
solution for a problem because of the significant computational
effort required. Two reasons for using Monte Carlo are:
 No other solution is available.
 The details of the problem are sufficiently important that no
approximations
can be tolerated.
In the case of particle transport, the latter reason is usually invoked.
Examples of Monte Carlo particle transport include:
 Neutron tranport in fission devices,
 Radiation transport,
 Electron transport in semiconductors,
 Neutral particle transport in plasmas.
2.1.1 Sampling a Distribution
The most fundamental process of a Monte Carlo simulation is the sampling
of probability distributions. Typical examples for a particle simulations
would be:
 Choosing a velocity vector from a Maxwellian distribution,
 Picking a particle from a spatially varying source,
 Determining the distance to the next collision.
2.1.2 Sampling Distance to Next Collision
To illustrate the third case, consider particles with a total,
macroscopic cross
section Σ_{t}. We write the probability density function
for having
the next collision between x and x + dx, p(x) dx, as:
p(x) dx = Σ_{t} exp(−Σ_{t} x) dx. 
 (2.1) 
Then, the cumulative probability distribution P(y),
is the probability of having a collision for any x ≤ y. Note that
P(0) = 0 and P(∞) = 1.
Figure
Figure 2.1: A uniformly sampled ξ value is translated into a sample of y
distributed according to the "Probability Density" curve using the integrated
"Probability" curve.
As suggested by the above figure, P(y) can be sampled by choosing a
uniformly distributed random number ξ, with 0 ≤ ξ < 1. For
this example, P(y) is easily computed and inverted so that
where we have surreptitiously replaced 1 − ξ with ξ since both
are uniform over the unit interval. Again the logic of this process,
as depicted in the figure, is that if ξ is
uniformly distributed, then numbers between ξ to ξ+ ∆ξ
are chosen with the same frequency as particles having collisions between
y and y + ∆y. Thus, y is the sampled
distance to next collision.
2.1.3 Sampling More Complicated Distributions
Because of the simplicity of the functional form of the collision probability
density function, the distance to the next collision can be easily sampled
using a single, uniform random number. If this were not the case, more
general methods, such as this rejection technique, would have
to be employed.
As an example, consider the sampling of a velocity v from a Maxwellian
distribution,
p(v) = 
1

exp  ⎡ ⎣

− 
(v −v_{0})^{2}
v_{th}^{2}
 ⎤ ⎦

. 
 (2.4) 
For simplicity, limit the v values to v − v_{0} < M v_{th} for some
appropriate M.
Figure
Figure 2.2: Values of x distributed according to p(x) can be sampled by
uniformly selecting (x,ξ) pairs and tossing out ones with ξ > p(x).
Now, we choose a pair of random numbers (x, ξ) with x uniformly
distributed over the interval
v_{0} − M v_{th} ≤ x < v_{0} + M v_{th} and ξ uniform over
0 ≤ ξ < p_{max}, where p_{max} is the maximum value of p(v).
If ξ < p(x), accept this x as the sampled velocity. Otherwise, if
ξ > p(x), reject these (x, ξ) and sample another pair. That the
resulting set of x values is distributed as required can be seen by studying
the figure above. This figure also makes the case that the rejection approach
may not be very efficient for sampling a Maxwellian (indeed, DEGAS 2
uses the wellknown BoxMueller method for this task).
In general, a rejection technique factors the probability density p(x) into
two functions, p(x) = g(x) h(x), where it is known how to sample x from
the function g(x). The acceptance condition is just ξ < h(x). The
Maxwellian example simply sets g(x) = 1.
2.2 Particle Transport
One of the main benefits of doing Monte Carlo particle transport is that
most of the time the algorithm can be understood directly from a physical
picture. Namely, we think of the code as tracking particles through a
medium and having collisions with probabilities determined from the cross
sections of the various collision processes. The outcome of those collisions
is the same as for the actual physical particles, and so on. This is
known as the analog description of Monte Carlo particle transport.
In this section, we want to give a rough idea of how this picture relates
to the integral equations which provide the actual mathematical basis
for Monte Carlo techniques. Only from these equations can we understand
the nonanalog techniques (the rejection technique described in
the previous subsection is a nonanalog technique) used to make the Monte
Carlo process more efficient. We hope to familiarize the reader with
the integral equations to the point where more detailed
references[2,12] can be readily digested.
From an analog pointofview, we imagine a Monte Carlo particle simulation
as consisting of four parts:
 A way to sample an initial particle distribution,
 A way to track particles from one point in space to another,
 Tools to treat "interesting things" that might happen along the way,
 Some useful quantities to compute and record at each stage of the
process.
The last of these, referring to "scores" and their "estimators" will
be addressed in Sec. 2.3.
2.2.1 Motivating the Integral Equation
The initial particle distribution, or source, is given by a probability
density function Q(x) > 0. We take x here to represent a point in
phase space so that this function provides both the velocity and spatial
distribution of the source.
The tracking of particles is described in the integral equation by a function
L(y → x) which gives the probability that a particle starting
at a point y
in phase space ends up at point x. This function includes the "interesting
things" (collisions)
that might happen along the way, as we will describe in the next
subsection.
Using these two functions, we can write the density of particles
arriving at x, χ(x) as the sum of two pieces:
 Arrivals directly from source: Q(x).
 Particles transported from elsewhere:
And we have
χ(x) = Q(x) +  ⌠ ⌡

χ(y) L(y → x) dy. 
 (2.6) 
2.2.2 Transport and Collision Kernels
In order to describe the kernel L(y → x) in more detail, substitute
for y and x the phase space variables (→r,→v) and a species
index i. The kernel L can then be broken into two pieces:
L( 
→
r

′

, 
→
v

′

, i^{′} → 
→
r

, 
→
v

, i) = T( 
→
r

′

→ 
→
r

; 
→
v

, i) C( 
→
v

′

, i^{′} → 
→
v

, i; 
→
r

), 
 (2.7) 
where T(→r^{′} → →r; →v, i)
is the (probability) density of particles leaving a
collision at →r^{′} and having a subsequent collision at →r.
The ";" is used in
the arguments of T and C to separate changing parameters (on the left) from
those held fixed (on the right). Basically, sampling from T is as
described in
Section 2.1.2.
Collisions are described by C(→v^{′}, i^{′} → →v, i; →r) which provides the probability
that a particle of velocity →v^{′}, species i^{′}
will have a collision resulting in
a particle of velocity →v, species i. In practice, C is a
sum over the physical
reactions in which species i^{′} participates.
2.2.3 Integral Equation Details
One of the problems with reading Monte Carlo particle references is that
there are several slightly different "quantities of interest" equivalent
to χ(x). Furthermore, there is no wellestablished notation for those
quantities. In this subsection, we attempt to list these various quantities
and to provide a detailed form of an integral equation for later reference.
The quantity χ(x) is called the postcollision density. The
origin of this designation should be clear from the previous discussion.
Equivalently, the integral equation can be written in terms of the
precollision density ψ(x). The two are related by
ψ( 
→
r

, 
→
v

,i) =  ⌠ ⌡

d^{3}r^{′} χ( 
→
r

′

, 
→
v

,i) T( 
→
r

′

→ 
→
r

; 
→
v

,i). 
 (2.8) 
It is this precollision density which is most closely related to the
more familiar particle flux Φ and particle
distribution function f by the definitions:
ψ( 
→
r

, 
→
v

,i) = Σ_{t}( 
→
r

, 
→
v

,i) Φ( 
→
r

, 
→
v

,i) = Σ_{t}( 
→
r

, 
→
v

,i)  
→
v

 f( 
→
r

, 
→
v

,i). 
 (2.9) 
The integral equation for ψ bears a strong resemblance to that for
χ. We write it out here in detail to illustrate the details as well as
to provide a reference point for the rest of the document and other
reading:
 

 
 

+  ⌠ ⌡

∞
0

dR  ⌠ ⌡

d^{3}v C( 
→
v^{′}

,i^{′} → 
→
v

,i; 
→
r

− 
→
Ω

R) 
 
 

×T( 
→
r

− 
→
Ω

R → 
→
r

; 
→
v

,i)ψ( 
→
r

− 
→
Ω

R, 
→
v^{′}

,i^{′}), 
  (2.10) 

where
S( 
→
r

, 
→
v

,i) =  ⌠ ⌡

∞
0

dR Q( 
→
r

− 
→
Ω

R, 
→
v

, i)T( 
→
r

− 
→
Ω

R → 
→
r

; 
→
v

,i), 
 (2.11) 
and
T( 
→
r

− 
→
Ω

R → 
→
r

; 
→
v

,i) = Σ_{t}( 
→
r

, 
→
v

,i) exp  ⎡ ⎣

−  ⌠ ⌡

R
0

dR^{′} Σ_{t}( 
→
r

− 
→
Ω

R^{′}, 
→
v

,i)  ⎤ ⎦

. 
 (2.12) 
The actual form of the collision kernel C need not be written down
here. Instead, in this document we'll only talk about its implementation
from an analog point of view. The only additional piece of information
about it that we could add to the above equations is that C consists
of a sum over the relevant collision processes involving species i^{′}.
The approach of writing the path integrals as above instead of as general
phasespace volume integrals with the path enforced via
deltafunctions[12]
was copied from Case and Zweifel[13]
We find that this leads to expressions which are easier to work
with and think about.
2.3 Estimators
The objective of the Monte Carlo calculation is to evaluate integrals over the
distribution function:
I(g) ≡  ⌠ ⌡

d^{3}r d^{3}v f( 
→
r

, 
→
v

)g( 
→
r

, 
→
v

). 
 (2.13) 
For example, if we wish to estimate a reaction rate
ρ ≡ Σ→v,
I(ρ) =  ⌠ ⌡

d^{3}r d^{3}v ΣΦ( 
→
r

, 
→
v

). 
 (2.14) 
Values of ρ are accumulated along the random walk of the neutral flight
and used to estimate I(ρ). The mathematical description of this
accumulation process is called the estimator of I.
2.3.1 Simple Estimators
Consider a random walk γ which consists of k steps,
(x_{1}, x_{2}, …, x_{k}). I.e., the
particle is absorbed or otherwise terminated at the kth step.
The simplest estimator involves writing down the value of ρ
at the point of absorption. Thus, the absorption estimator is
ξ_{a}(γ) = 
ρ
ρ_{a}

w χ_{V}, 
 (2.15) 
where ρ_{a} is the absorption rate in the volume of interest
V and w is the statistical weight of the particle (constant for
this random walk process). The function χ_{V}
just permits the score to be localized to the volume of interest.
The expectation value of the estimator is just the sum over all
possible random walks γ of the probability of each
walk P(γ) times the corresponding estimator,
E(ξ_{a}) = 
∑
γ

P(γ) ξ_{a}(γ). 
 (2.17) 
Obviously, only one data point is collected for each random walk.
We can instead accumulate a score at each step (scattering collision)
of the random walk.
ξ_{c1}(γ) = 
k ∑
m=1

 ⎡ ⎣

ρ_{a} (x_{m})
ρ_{t} (x_{m})
 ⎤ ⎦

wχ_{V} (x_{m}), 
 (2.18) 
where ρ_{t} = ρ_{a} + ρ_{s} is the total reaction
rate for the particle, with ρ_{s} being the scattering
rate. The subscript c1 denotes the first collision estimator
described here. The expectation value of the estimator is again
formally written as
E(ξ_{c1}) = 
∑
γ

P(γ) ξ_{c1}(γ). 
 (2.19) 
Because of the greater amount of data accumulated along each particle
track, the variance of this estimator will usually be smaller
than that of the absorption estimator.
Consider now a random walk process in which absorption is forbidden.
This is our first departure from the usual analog processes.
Instead of an absorption, the weight of the particle is reduced
at each scattering collision by its survival probability.
So, at the mth step of the walk the weight of the particle is
w_{m} = w 
m ∏
j=1


ρ_{s} (x_{j})
ρ_{t} (x_{j})

. 
 (2.20) 
The flight must be terminated in some manner; the usual practice
is to have it undergo Russian roulette[2] when
w_{m} < w_{min}, where w_{min} is some specified
minimum weight. The collision estimator for this nonanalog
random walk process is then
ξ_{c2}(γ) = 
k ∑
m=1

w_{m}  ⎡ ⎣

ρ_{a} (x_{m})
ρ_{t} (x_{m})
 ⎤ ⎦

χ_{V} (x_{m}). 
 (2.21) 
For formal proofs that these expectation values do indeed yield I(ρ),
see Spanier and Gelbard[2].
2.3.2 Generalized Collision and Tracklength Estimators
Two principal
estimators are used in DEGAS 2. One is a generalization of the
collision estimator described in the previous subsection; the other
is the tracklength estimator.
These estimators can be obtained by making a slight modification of the
derivation given by Macmillan[51] (who had improved upon an
earlier approach described by Spanier[52]). To make the connection
to DEGAS 2 clearer, we replace the Σ and d variables of Macmillan
with ρ and t.
Macmillan's paper, like the earlier one by Spanier, aims to derive the
tracklength estimator. Collisions are described as either
scatterings, rate ρ_{s}, or absorptions ρ_{a}. The
starting point is the nonanalog random walk process described in
Sec. 2.3.1. The probability of an event occurring in
a time interval dt is (ρ_{s} + ρ_{a}) dt. Again, the weight
reduction factor applied at each scattering
is ρ_{s} / (ρ_{s} + ρ_{a}).
The trick used by Spanier is to introduce a fictitious "delta scattering"
that has no impact on the particle's trajectory or weight. These
"pseudocollisions" are not unrelated to those of the original
DEGAS code[3]. Additional flexibility in the final results is
obtained by associating some fraction of the absorptions α with
the delta scatterings. To obtain the estimators used in DEGAS 2, we
work with a total fictitious scattering reaction rate
ρ_{t,δ} ≡ ρ_{s,δ} + αρ_{a}, 
 (2.22) 
where ρ_{s,δ} is equivalent to Macmillan's Σ_{s,δ}.
Then, the probability of a delta scattering in a time interval dt is
ρ_{t,δ} dt. At each delta scattering, the weight is
reduced by a factor (ρ_{t,δ} − αρ_{a}) / ρ_{t,δ}.
The probability of a real scattering in the same interval is
(ρ_{s} + (1 − α) ρ_{a}) dt; at each such event, the
weight is reduced by ρ_{s} / [ρ_{s} + (1 − α) ρ_{a}].
If a particle travels a time t in a region of constant parameters
before having a real collision, the conditional probability of having
n pseudocollisions is
P(n  d) = 
( ρ_{t,δ} t)^{n}
n!

exp(− ρ_{t,δ} d). 
 (2.23) 
At each pseudo or real collision, we use the standard collision estimator
for the nonanalog random walk, taken from Eq. (2.21):
where w is the current weight of the particle and
ρ_{t}^{∗} = ρ_{t,δ} + ρ_{s} + (1 − α) ρ_{a} 
 (2.25) 
is the overall total reaction rate.
The expected value of this estimator over the n collisions in the
region is then
E(ξ  t) = 
∞ ∑
n=0

P(n  t) 
ρ
ρ_{t}^{∗}


n ∑
i=0

 ⎛ ⎝

ρ_{t,δ} − αρ_{a}
ρ_{t,δ}
 ⎞ ⎠

i

, 
 (2.26) 
where we have assumed that the initial weight of the particle is 1. The
last factor represents the subsequent weight reductions with each
pseudocollision.
Following Macmillan's evaluation of E(ξ  t), we find that we
can write
E(ξ  t) = c^{′} 
ρ
ρ_{s} + (1 − α) ρ_{a}

+ (1 − c^{′}) ρ 
1 − exp(− αρ_{a} t)
αρ_{a}

, 
 (2.27) 
where
c^{′} ≡ 
ρ_{s} + (1 − α) ρ_{a}
ρ_{t,δ} + ρ_{s} + (1 − α) ρ_{a}


 (2.28) 
should be compared with Macmillan's c.
Macmillan's expression for the weight reduction after time t
is exactly the same in our case:
E(w  t) = exp( − αρ_{a} t) 
ρ_{s}
ρ_{s} + (1 − α) ρ_{a}

. 
 (2.29) 
The factor at the end represents the weight reduction by the real collision
at time t. If the particle travels a time t through the region without
having a real collision, its weight should be reduced by just
exp( − αρ_{a} t).
In this latter situation in which there is not a real collision at time
t, only the pseudocollisions contribute to the score. An interesting
difference from MacMillan's approach is that if we consider the limit
ρ_{t,δ} → 0, we should get no score at all. In
contrast, MacMillan's pure "collision estimator" still compiles
a nonzero score as Σ_{s,δ} → 0 through the
αΣ_{a} contribution to the pseudocollision cross section.
Thus, MacMillan's "collision estimator" behaves somewhat like
a tracklength estimator if there are no collisions in a region.
For simple scores, this behavior is desirable: lower variances will
result in regions of small cross section. However, a few scores are
sufficiently complex that averaging along the particle path cannot
be done, and scoring can only be done at collisions. A pure
collision estimator is required for handling such scores.
In the case of no real collision at time t, the expectation value
for our estimator becomes:
E(ξ  t) = 
ρ
ρ_{t}^{∗}

ρ_{t,δ} 
1 − exp(− αρ_{a} t)
αρ_{a}

. 
 (2.30) 
Clearly, E(ξ  t) → 0 as ρ_{t,δ} → 0.
The limiting expressions above for ρ_{t,δ} → 0
can be obtained from MacMillan's with Σ_{s,δ} →− αΣ_{a}.
We now consider the expressions used in DEGAS 2. We start by using
α = 1 so that absorptions occur only at the delta scatterings.
This is consistent with our use of the nonanalog random walk, i.e.,
suppressed absorption. The collision estimator is then given by
ρ_{t,δ} = 0,
ξ_{c} ≡ 
ρ
ρ_{s}

exp(− ρ_{a} t) w χ_{V}, 
 (2.31) 
where w is the weight of the particle at the beginning of timestep
t and we have inserted the characteristic function of V χ_{V}
for completeness.
At the end of the timestep (after the collision is scored) the weight
is reduced by
w^{′} = w exp( − ρ_{a} t). 
 (2.32) 
The particle track continues from there, and multiple scores may be made
inside the volume V. If the flight travels through the volume V without
a collision, no contribution to the score is made. However, the
weight is still reduced by the factor in Eq. (2.32).
The tracklength estimator is obtained with ρ_{t,δ} →∞,
ξ_{t} ≡ ρ 
1 − exp(− ρ_{a} t)
ρ_{a}

w χ_{V}. 
 (2.33) 
Note that Eq. (2.30) yields this same result in the case
of no collisions. The same weight reduction, Eq. (2.32)
is applied at the end of t. Conveniently, this allows us to use the
same weight factor for both estimators (true in MacMillan's case as well).
Molecules and some other species are not treated with suppressed absorption.
Instead, an ionizing collision terminates the flight (typically changing
the species). The estimators for this case are obtained by setting
α = 0 in the above expressions:
ξ_{c} ≡ 
ρ
ρ_{s} + ρ_{a}

w χ_{V}, 
 (2.34) 
and
There is no weight reduction in this case, and w retains its initial
value for the duration of the flight.
2.3.3 Specification of Tallies
This section briefly covers the practical implementation of the estimators
described in the previous subsection.
DEGAS 2 currently features three primary types of estimators. These estimators
are used in accumulating three types of tallies. The tally types are
 test,
 reaction,
 sector.
These designations can also be viewed as characterizing the nature of the
information contained in the function g of Eq. (2.13).
The following discussion of the three estimators applies to the first
two of these; there is only one estimator for the sector tallies.
The first estimator, associated with the accumulating
test particle information, is just the track length estimator (see also
Sec. 2.3.2),
ξ_{g}^{TLE} = w 
1 − exp(−σ_{i} t_{V})
σ_{i}

g, 
 (2.36) 
where t_{V} is the time the particle track spent in volume V,
w = weight at start of the time step. Compare this with Eq. (2.33).
The function g can be a test particle attribute, such as
mass, momentum, energy, or reactionrelated, e.g., ion
momentum source due to charge exchange. The latter g quantity
enter the tracklength estimator as averages over the background
distribution[17].
The second type of estimator is the collision estimator,
ξ_{g}^{CE} = 
∑
m_{s}

w_{ms} 
g
σ_{s}

, 
 (2.37) 
where the sum is over m_{s}, all scattering collisions within V.
The quantity w_{ms} is the weight at m_{s}th collision (by using the
instantaneous weight rather than the initial weight,
the exponential weight reduction factor that appears explicitly in
Eq. (2.31) is incorporated automatically).
Again, g is a test particle attribute.
There is an analogous expression for the collision estimator associated
with a particular reaction,
ξ_{g}^{CE} = 
∑
m_{j}

w_{mj} 
g
σ_{j}

, 
 (2.38) 
where the sum is now over m_{j} collisions of reaction j within V,
The function g is some quantity related to the reaction and is likely
something known only at collisions. One example of such a quantity would
be the energy source due to H_{2} dissociation (i.e., the energy
exchange is determined only once the product velocities have been
determined).
The third type of estimator arises from the observation that
many integrals of interest can be written as:
I(g) ≡  ⌠ ⌡

d^{3}r g( 
→
r

)  ⌠ ⌡

d^{3}v f( 
→
r

, 
→
v

) =  ⌠ ⌡

d^{3} r g( 
→
r

) n( 
→
r

). 
 (2.39) 
This leads to the postprocessing estimator. The requirement is that
g not depend upon the test or background velocities. An example of
this is the H_{α} emission rate.
Information accumulated when a flight crosses a diagnostic sector (surface)
is compiled in the sector tallies. The
stimator is just w g, where
w = is the weight when the track hits the sector.
Examples of the g = functions include mass, incident angle, and energy.
This is effectively a collision estimator.
There is also a "null" or "none"
estimator that can be used, say, to eliminate the
contribution of a particular reaction to a tally without having to
remove that reaction from the problem.
For each tally, we need to specify
 A descriptive label,
 The type of tally (test, reaction, or sector),
 Its geometry, [volume (per zone), surface (diagnostic sector),
or detector (e.g., chordintegrated)],
 The dependent variable (e.g., energy, background momentum,
H_{α} emission rate). The requested variables
are matched directly against the labels in the atomic physics data files.
If a new reaction and data file containing a new variable are added to
the code, that variable can be used as the basis for a corresponding tally
without having to modify the code.
 The rank (or dimensionality) of the tally,
 Its independent variables (1 → rank) [e.g.,
zone, test particle species, energy bin (diagnostic sector),
wavelength bin (detector).
 The estimator to be used (tracklength, collision, postprocessed,
or null). For reaction tallies, the estimator chosen can vary with
reaction. It may be advantageous to combine results
from than one estimator, but no effort to implement such a capability
has been made yet.
 Any conversions that must be performed on the tally prior to
output.
Usually, these conversion will do some sort of scaling, e.g., by
mass or volume. An important example is
converting →v from Cartesian
to cylindrical coordinates when scoring.
2.4 Tracking Procedure
The two previous subsections are tied together with an outline of the
DEGAS 2 subroutine for following flights in Fig. 2.3.
Figure
Figure 2.3: Simplified flow chart of DEGAS 2 subroutine
follow, showing the steps
involved in tracking a flight from its birth to its termination. Note
the weight reduction, and collision and tracklength estimators
(Sec.
2.3). The operators from the integral equations are
shown along with the pre and postcollision densities. These labels
tie the practical algorithmic operations to the equations in
Sec.
2.2.3. Some details have been omitted
for clarity.
The flight is initiated as it is sampled from the spatial and velocity
distribution of a source. Conceptually, this step amounts to
specifying χ_{0}(→r^{′},→v^{′},i) using the
source term Q(→r^{′},→v^{′}). The subscript
0 denotes the initial track for the flight.
The next stage of the algorithm involves computing the reaction
rates ρ_{s}, ρ_{a}, etc. for the flight. Continuing
flights which require a recomputation (due to crossing from one
plasma region to another or to a change in species or velocity)
repeat the main loop starting at this same point. The
postcollision density χ_{n}
in Fig. 2.3 is intended to refer to them.
The time to the next collision is randomly sampled
(Sec. 2.1.2) from the sum of the nonabsorbing
reaction rates using a random number ζ.
The weight w is compared with some minimum weight to guarantee
that the flight terminates at some point. The "End flight"
process indicated in Fig. 2.3 actually represents
a Russian Roulette procedure[2] which may allow the flight
to live a little longer.
The transport kernel T translates the flight from position
→r^{′} to →r. At the same time, this step
represents the conceptual transition from postcollision density
to precollision density ψ, Eq. (2.8). The tracking
algorithm stops when t = t_{max}, the flight
crosses a zone boundary (denoting a change in plasma parameters), or
a "sector" (a material or diagnostic surface; these will be
discussed in Sec. 2.6) is reached. If a zone boundary
is crossed, the flight continues as is, although the rates must
be recomputed. Processing a sector may result in the continuation or the
termination of the flight, depending on the type of sector. Again,
the details have been omitted for simplicity.
If the flight reaches t = t_{max}, the implication is that
it has had a collision. A particular collision is randomly chosen from
amongst those possible for this particular flight.
At this point, precollision scores are recorded
and the collision itself is processed. Conceptually, this represents the
transition via the operator C from ψ_{n+1} to χ_{n+1}.
If the reaction chosen results in no particles to track (e.g.,
in no neutral species; particles which DEGAS 2 has been instructed
to track are called "test" particles), the flight is terminated.
Test particles resulting from the collision are tracked in the same
manner, starting with the computation of the reaction rates.
2.5 Random Numbers
The concepts behind the random number generator used in DEGAS 2 are
described in the file random.web. The
documentation at the top of the file
provides a nice introduction.
2.6 Geometry
The hiearchy of geometryrelated objects in
DEGAS 2 is, from the lowest level to the highest level, is:
 surface,
 cell,
 polygon,
 zone.
Individual cells are composed of surfaces, as discussed in the documentation
for the internal geometry (e.g., see geometry.web; follow
this link to the PDF version of this file).
The next level up
in two (and three, to a limited extent) dimensions
is a polygon.
Finally, a zone may consist of one or more polygons; properties
are constant across a zone. For example, a single zone might be used to
represent the vacuum region around the plasma which is comprised of several
(possibly disconnected) polygons. Or, a plasma flux surface on which density
and temperature are constant might be a single zone.
Note that polygons are used only for external interfaces to the
geometry such as that provided by definegeometry2d or
readgeometry and are not retained in the classes of the
main code. Cells and surfaces
are essentially used only by the code itself; the user would have to deal
with them only in debugging. Zones are used primarily for the
specification of input and output data. For more detailed
background on the DEGAS 2 geometry and some guidance on how to
set one up for your particular problem, see the introductions
to the external geomtry interface codes definegeometry
(PDF file here),
readgeometry
(PDF file here), or
boxgen (less well documented).
Another important component of the geometry are the "sectors";
they are defined at interfaces between adjacent zones of
different types to facilitate tracking or between zones of the same
type for diagnostic purposes.
Namely, flights are halted at sectors during tracking. If the next
zone along the track represents, say, a solid material, the
main code hands that flight off to the routines that
effect the interaction of the current particle with that material.
While stopped at that sector,
scores to corresponding "diagnostics" (groups of sectors)
are also tabulated. The essential identifying information
of a sector are its surface and zone numbers.
See the
sector class documentation for a more
detailed description. Default diagnostic sector groups defined
by all DEGAS 2 geometry setup programs are discussed in
Sec. 3.11.
2.7 Symmetry and Coordinate Systems
The coordinate system in the main DEGAS 2 code is purely 3D
Cartesian; that's essential to an efficient tracking algorithm. There
are, however, two coordinate system related concepts.
 The geometry's symmetry. In most problems, the user's objectives
can be met with an approximate 2D solution possessing an ignorable
coordinate. The two situations considered in DEGAS 2 are invariance
to translation along the y axis ("plane" symmetry) and
invariance to rotation about the z axis ("cylindrical" symmetry).
In both cases, the tracking is still done in 3D, but the output
results are averaged over the third dimension. As will
be discussed in more detail below, 1D and 3D "symmetries"
are also permitted. Most geometry related variables refer to
Cartesian coordinates. However, in cylindrically symmetric or nearly symmetric
(3D) cases, some geometry related variables are in cylindrical
coordinates.
 Coordinate system of input and output velocity fields.
The main need for this is in cylindrically symmetric or
nearly symmetric
problems, e.g., when DEGAS 2 is coupled to a 2D
plasma code like UEDGE. In these cases, the plasma flow velocities
are effectively specified at a toroidal angle ϕ = 0, corresponding
to y = 0.
To use these velocities in a collision at some particular location
for which ϕ ≠ 0 (equivalently, y ≠ 0),
the velocities
need to be rotated through that angle ϕ.
Likewise, a code like UEDGE expects to get from DEGAS 2
toroidally symmetric ion momentum
sources specified at ϕ = 0. So, the scores for those sources
need to be rotated from the local toroidal angle ϕ ≠ 0
to ϕ = 0 as they are
being compiled. These two transformations are handled by macros
in the background class.
In problems with
plane or no symmetry, these macros have no effect.
The currently permissible symmetry types (see
also geometry.hweb)
are:
 None
 Currently indicates that the symmetry has not been defined,
not that the geometry has no symmetry.
 Plane
 A 2D geometry in which the y coordinate is ignored.
 Cylindrical
 A 2D geometry in which the toroidal angle
ϕ is ignored.
 1D
 A 1D (plane) geometry in which the y and z
coordinates are ignored.
 3D Plane
 This geometry is close to be planar symmetric;
all of the geometrical objects are defined using 2D shapes
defined in a single plane. The problem is bounded in the y
direction, and the y variation of the results can be
computed and output. The nature (solid vs. plasma or vacuum)
of the objects can also vary in
the y direction. The problem is divided in the y direction
by y = constant surfaces. This might be described as a
"bread slice" approach to a 3D geometry.
 3D Cylinder
 This is the cylindrical equivalent of the
3D Plane symmetry. The full range of toroidal angles (360 degrees
or 2 π radians) is included in this case. Obviously, the
problem is divided in the ϕ direction by ϕ = constant
surfaces. The analogous simplified characterization of this approach
is the "pie slice" technique.
 3D Cylindrical Section
 This is also a cylindrical equivalent
of the 3D plane symmetry. The difference is that in this case
only a finite range of toroidal angles is considered. The problem
is bounded by ϕ = constant surfaces at the extrema of this
range. Extending the "pie slice" analogy, this case represents
some number of slices, but less than the whole pie!
2.8 Time Dependence
The default and most frequently used mode of running DEGAS 2 simulates a
steady state or time independent system. The reason for this is physically
obvious: unconstrained by magnetic field lines, neutral particles can
travel through a plasma much more quickly than charged species, and
in many cases, the neutral profile equilibrates on a time scale
short compared with global transport times.
For a rapidly evolving system, however, a
time dependent calculation may be needed. The implementation
of time dependence in DEGAS 2 follows
that used in EIRENE [5].
Only the defineback input preprocessor is presently set up
to handle direct user specification of a time dependent run;
see the corresponding source file
for details on how this is done.
However, the user can in principle also manually
set the controlling
parameters in the source class.
The user specified sources in a time dependent run are
either fixed in time during a time
interval (uniform) or start at the beginning of the interval
(δ function). More complex time variation of the sources
can be created by breaking up the source time dependence into a
set of discrete time intervals with the sources being uniform
over each.
To provide continuity between consecutive time intervals, an
additional "snapshot" source group is set up to sample
from the neutral particle distribution function (PDF)
that was in the volume at
the end of the end of a previous interval.
Again following [5], we represent this PDF via the
particle positions,
velocities, and weights; this is the purpose of the
snapshot class. These
data are written to the file indicated by the
snapshotfile in degas2.in (Sec. 3.5.1).
With this
approach, the amount of kinetic detail in the PDF is
maximized.
2.9 Atomic Physics
The interaction of electrons with atoms and ions in a divertor
plasma is complicated by having comparable collision and radiative
decay times. The techniques for treating such systems were first
described by Bates and McWhirter[6,7]. The physical
processes included in such a collisionalradiative model are:
 Electron collisional excitation K_{e,mn} and deexcitation
K_{d,mn},
 Spontaneous radiative transitions A_{mn},
 Electron collisional ionization K_{i,mn},
 Threebody recombination K_{r,mn},
 Radiative recombination β_{rad,m},
 Dielectronic recombination (for multielectron atoms)
β_{dia,m},
where m and n are principal quantum numbers. For the two state
terms, the initial state is the second number; the final is the
first. This apparently backward notation allegedly has a mathematical
origin.
A set of m equations describes the balance between these processes that
determines the density of the mth excited state of the atom,
N_{m},
 

− N_{m}  ⎡ ⎣

n_{e}  ⎛ ⎝

K_{i,m} + 
∑
n > m

K_{e,nm} + 
∑
n < m

K_{d,nm}  ⎞ ⎠


 
 

 
 

+ n_{e}  ⎛ ⎝

∑
n > m

N_{n} K_{d,mn} + 
∑
n < m

N_{n} K_{e,mn}  ⎞ ⎠


 
 

+ 
∑
n > m

N_{n} A_{mn} + n_{e} n_{i} ( β_{rad,m} + β_{dia,m} + n_{e}K_{r,m} ), 
  (2.40) 

where n_{e} is the electron density and n_{i} is the density of ions
corresponding to the atom under consideration (m → ∞,
effectively).
The basic assumption of the collisionalradiative model[6,7]
is that states above ground state will decay much more rapidly than
the ground state will change.
The evolution of the ground state is considered to be comparable to
the transport time scales. This physical process is the one considered by
the neutral transport code. In other words, the m = 1 equation will
essentially be solved by DEGAS 2 and N_{1} can be treated as an
externally known quantity. The relatively rapid equilibration of the
excited states permits the time derivatives of their densities
to be set to zero. The set of equations for m > 1 can then be
solved,
N_{m} =  ⎛ ⎝

N^{(i)}_{m}
N_{1}
 ⎞ ⎠

N_{1} +  ⎛ ⎝

N^{(ii)}_{m}
n_{i}
 ⎞ ⎠

n_{i}. 
 (2.41) 
In practice, the solution of these equations involves:
 Obtaining the various reaction rates, K_{e,mn},
K_{d,mn}, A_{mn}, K_{i,mn}, K_{r,mn}, β_{rad,m}, β_{dia,m},
as a function of n_{e} and T_{e},
 Solving the m > 1 equations for each n_{e}, T_{e}
pair with N_{1} = 1 and
n_{i} = 0, providing the "coupling to the ground state",
 Solving the equations again with N_{1} = 0 and
n_{i} = 1, yielding the "coupling to the continuum".
Tables of the population coefficients N^{(i)}_{m}/N_{1} and
N^{(ii)}_{m}/n_{i} are thus obtained as a function of n_{e} and
T_{e}.
With Eq. (2.41), the m = 1 equation can becomes

d N_{1}
d t

= − N_{1} n_{e} S_{eff} + n_{e} n_{i} R_{eff} + transport, etc., 
 (2.42) 
with
S_{eff} = 
∑
m ≥ 1

 ⎛ ⎝

N^{(i)}_{m}
N_{1}
 ⎞ ⎠

K_{i,m}, 
 (2.43) 
and
R_{eff} = − 
∑
m ≥ 2

 ⎛ ⎝

N^{(ii)}_{m}
n_{i}
 ⎞ ⎠

K_{i,m}+ 
∑
m ≥ 1

( β_{rad,m} + β_{dia,m} + n_{e} K_{r,m} ). 
 (2.44) 
S_{eff} is thus the effective ionization rate of the ground state
neutral, incorporating the multistep processes associated with all
of the excited states. Likewise, R_{eff} is the effective recombination
rate. Both are functions of n_{e} and T_{e}. Selected
values of
the population coefficients N^{(i)}_{m}/N_{1} and
N^{(ii)}_{m}/n_{i} may also be tabulated or used to compute
photon emission rates for specific transitions of interest
(e.g., Balmerα).
The above assumes that just the ground state evolves on the transport
time scales relevant to DEGAS 2. If metastable excited states
are present in the system (e.g., helium's 2^{1}S and 2^{3}S
states or hydrogen's n = 2 in an optically thick region),
they may need to be treated on this same footing. The need
for such a treatment and its formulation go beyond the
scope of this presentation. See the discussion below on
helium for some additional references.
Computing the electron energy sources and sinks associated with
Eq. (2.40) involves accounting for the energy exchanges
in each of the processes appearing in Eq. (2.40):
 

E_{1} ( − N_{1} n_{e} S_{eff} + n_{e} n_{i} R_{eff} ) 
 
 

− 
∑
n < m^{m,n}

N_{m} A_{nm}( E_{n} − E_{m} ) 
 
 

− n_{e} n_{i} 
∑
m ≥ 1

β_{rad,m} ( 〈E_{rad,m} 〉+ E_{m}) 
 
 

− n_{e} n_{i} 
∑
m ≥ 1

β_{dia,m} ( 〈E_{dia,m} 〉+ E_{m}), 
  (2.45) 

where
 E_{1} is the ionization energy,
 the first term is the loss due to ionization,
 the second is the gain due to recombination.
 The sums represent lost photons from radiative transitions and
recombination,
 〈E_{rad,m} 〉 and
〈E_{dia,m} 〉 are the average energies of electrons in
radiative and dielectronic recombination, respectively.
The two average energies are computed using
the expression given by Reiter[17],
〈E_{x,m} 〉 = T_{e}  ⎡ ⎣

3
2

+ 
T_{e}
β_{x,m}


d β_{x,m}
d T_{e}
 ⎤ ⎦

, 
 (2.46) 
where x denotes either rad or dia.
Using Eq. (2.41), Eq. (2.45) becomes

d E_{e}
d t

= − E_{1} ( N_{1} n_{e} S_{eff} − n_{e} n_{i} R_{eff} ) − N_{1} E_{loss}^{(i)} − n_{i} E_{loss}^{(ii)}, 
 (2.47) 
where the terms E_{loss}^{(i)}, E_{loss}^{(ii)} represent the loss
rates due to coupling to the ground state and continuum, respectively.
A different arrangement of these terms which can be
generalized to the case involving multiple "transported
states" has also been used:

d E_{e}
d t

= − N_{1} n_{e} E_{loss}^{(i)′} +n_{i} n_{e} E_{source}^{(ii)′} − n_{i} n_{e} E_{loss}^{(ii)′}. 
 (2.48) 
The relationships between the terms are nearly obvious, but
potentially prone to confusion. For clarity, they are:
 

E_{1} S_{eff} + E_{loss}^{(i)}/ n_{e}, 
  (2.49) 
 

  (2.50) 
 

  (2.51) 

Note that the "unprimed" energies have dimensions of
power, and the "primed" versions have dimensions of
power times volume. All are stored in tabular form
as functions of n_{e} and T_{e} for
use in DEGAS 2.
2.9.1 Hydrogen Collisional Radiative Model
The code we use to solve these equations for
the hydrogen atom, collrad,
is based on the one described
by Weisheit[20]. The original
code assumed that ionization and recombination were in balance
so that the neutral density could be expressed
in terms of the ion density. While this might be valid in
an astrophysical context, it is horribly wrong in fusion devices.
The first major change to the code is to have
it use Eq. (2.41) instead. The other
change is to replace the cross sections
with those provided in the JanevSmith database[18].
The resulting collisionalradiative data are believed to be
in good agreement with ADAS. Although the collrad code is not
itself distributed with DEGAS 2, it can be obtained from the DEGAS 2
authors (see Sec. 5.1).
The output data from collrad are contained in the file ehr2.dat.
Earlier versions of this same file, eh.dat
(as generated by Weisheit's original code) and ehr1.dat
[generated with Weisheit's code modified to use Eq. (2.41)] are
also included in the DEGAS 2 distribution since they were both used
extensively with the original DEGAS code.
The ehr2.dat file contains
 Ionization rate S_{eff} in cm^{3} s^{−1},
 Recombination rate R_{eff} in cm^{3} s^{−1},
 Neutral electron losses E_{loss}^{(i)} in erg s^{−1},
 Continuum electron losses E_{loss}^{(ii)} in erg s^{−1},
 Neutral "n=3 / n=1", N_{3}^{(i)} / N_{1},
 Continuum "n=3 / n=1", N_{3}^{(ii)} / n_{i},
 Neutral "n=2 / n=1", N_{2}^{(i)} / N_{1},
 Continuum "n=2 / n=1", N_{2}^{(ii)} / n_{i}.
Figure 2.4 shows the
variation of the effective ionization and recombination rates
with electron density and temperature.
Figure
Figure 2.4: Effective hydrogen ionization and recombination rates as a function
of electron temperature for selected electron densities.
The portion of the file for each quantity consists of 15 separate sections
(labeled by jn), corresponding to the various electron density values,
n_{e}(jn) = 10^{10 + (jn − 1)/2} cm^{−3}. 
 (2.52) 
Each section consists of 60 data values, one for each T_{e} (jt),
T_{e}(jt) = 10^{−1.2 + (jt − 1)/10} eV. 
 (2.53) 
The index jt starts at 1 in the first column and row and proceeds
rowbyrow. There are 10 rows of 6 entries each.
The density of the n = 3 excited state can be found with
N_{3} =  ⎛ ⎝

N^{(i)}_{3}
N_{1}
 ⎞ ⎠

N_{1} +  ⎛ ⎝

N^{(ii)}_{3}
n_{i}
 ⎞ ⎠

n_{i}. 
 (2.54) 
The corresponding emission rate for
the Balmerα line (a.k.a. "H_{α}", n = 3 → 2)
is given by
N_{3} A_{23} = 4.41 ×10^{7} N_{3} cm^{−3} s^{−1}.
Note that the energy of this transition is
E_{23} = 13.595 (2^{−2} − 3^{−2}) eV.
Similar expressions give the
n=2 density and the Lymanα (n = 2 → 1) emission
rate.
The data in ehr2.dat are transformed into netCDF files
containing the particular variables DEGAS 2 expects to find
(see Sec. 3.9) by the reactionwrite
code. Separate files for ionization and recombination
are written. In these files, S_{eff} and R_{eff} become
the "reaction rates". The "emission rate" variable
is obtained by rewriting the above expression
for the Balmerα emission rate as a power,
per atom (for ionization) or ion (for
recombination) and per electron:
 

E_{23} A_{23}  ⎛ ⎝

N_{3}^{(i)}
N_{1}
 ⎞ ⎠

/ n_{e}, 
  (2.55) 
 

E_{23} A_{23}  ⎛ ⎝

N_{3}^{(ii)}
n_{i}
 ⎞ ⎠

/ n_{e}. 
  (2.56) 

These quantities have dimensions of power times volume.
The wavelengths of this transition for each of the
hydrogen isotopes are also inserted into the
file, λ_{Hα} = 6562.80 Å,
λ_{Dα} = 6561.04 Å, and λ_{Tα} = 6560.45 Å.
The electron energy exchange terms in ehr2.dat correspond
to those in Eq. (2.47); reactionwrite transforms
them into those of Eq. (2.48).
2.9.2 Helium Collisional Radiative Model
The code used to solve the collisional radiative equations for
the helium atom was developed by Goto[21],
based on Fujimoto's original[22]. This code actually
generates data for two collisional radiative models. In the
first, only the ground state (1^{1}S) atom is considered slowly
varying (i.e., must be treated explicitly as a transported
species in DEGAS 2); this model was designated as
"formulation II" by Fujimoto. In the second, the metastable
2^{1}S and 2^{3}S states are also considered
slowly varying; Fujimoto called this "formulation I".
The appropriate species, reactions, and
data needed to do simulations based on this model have been
incorporated into DEGAS 2. However, investigations using it
have not demonstrated that the added detail yields any significant
differences in the final results. Moreover, a separate analysis
of the helium collisional radiative equations [Eq. (2.40)]
suggests that this particular
model may not even be valid for typical fusion edge conditions.
For these reasons, we will not discuss this model any further and
will focus only on the simpler "one state" model.
The primary modification made to Goto's code in preparation
for using its output with DEGAS 2 were to add the
computation of the electron energy exchange terms in
Eq. (2.48).
The resulting formulation II model and data are, not surprisingly, very
similar in structure to those of the hydrogen collisional
radiative model.
The file he_crII_alad.dat (in directory Aladdin/data) is
an Aladdin formatted version of these data. These same data
written in a format more like that of ehr2.dat are
contained in the file he_crII.txt
that can be downloaded
from the "Data" section of the DEGAS 2 web site. This file has
been provided as a convenience for those who would just like to
use the data, but not download the entire DEGAS 2 code.
The contents of this file are:
 Ionization rate S_{eff} in cm^{3} s^{−1},
 Neutral electron losses E_{loss}^{(i)′}
in eV cm^{3} s^{−1},
 Neutral 5877 Å (3^{3}D → 2^{3}P) emission
rate P_{5877}^{(i)} in eV cm^{3} s^{−1},
 Neutral 6680 Å (3^{1}D → 2^{1}P) emission
rate P_{6680}^{(i)} in eV cm^{3} s^{−1},
 Recombination rate R_{eff} in cm^{3} s^{−1},
 Continuum electron losses E_{loss}^{(ii)′}
in eV cm^{3} s^{−1},
 Continuum electron sources E_{source}^{(ii)′}
in eV cm^{3} s^{−1},
 Continuum 5877 Å (3^{3}D → 2^{3}P) emission
rate P_{5877}^{(ii)} in eV cm^{3} s^{−1},
 Continuum 6680 Å (3^{1}D → 2^{1}P) emission
rate P_{5877}^{(i)} in eV cm^{3} s^{−1},
The portion of the file for each quantity consists of 30 separate sections
(labeled by jn), corresponding to the various electron density values,
n_{e}(jn) = 10^{8 + 8*(jn − 1)/29} cm^{−3}. 
 (2.57) 
Each section consists of 65 data values, one for each T_{e} (jt),
T_{e}(jt) = 10^{(jt − 1)/16} eV. 
 (2.58) 
The index jt starts at 1 in the first column and row and proceeds
rowbyrow. There are 10 rows of 6 entries each, followed by a single
row of 5 entries.
The "emission rates" here correspond to the power in those lines
and are analogous to Eqs. (2.55) and (2.56) above.
The total emitted power (per unit volume) is obtained at run
time by combining them with the ground state and ion densities:
P_{5877} = P_{5877}^{(i)} N_{11S} n_{e} + P_{5877}^{(ii)} n_{i} n_{e}. 
 (2.59) 
Note that E_{33D → 23P} = 2.10955 eV,
A_{33D → 23P} = 7.069 ×10^{7} s^{−1}, and
E_{31D → 21P} = 1.85606 eV,
A_{31D → 21P} = 6.369 ×10^{7} s^{−1}.
The more precise (vacuum) wavelengths for these transitions
are 5877.1 Å and 6679.7 Å respectively. The first
ionization energy for helium is E_{1} = 24.587 eV.
The electron energy losses and sources can be used directly
in Eq. (2.48).
The data in he_crII_alad.dat are transformed into netCDF files
by the ratecalc code using the input files
heionII.input and herecII.input (in the data directory).
See the DATA_HISTORY file for additional details.
2.10 Elastic Scattering
2.10.1 NeutralIon Elastic Scattering
Note: this subsection was extracted from Ref. [23].
A few minor changes in notation have been made for consistency.
Only limited experimental data are available for ionneutral elastic
scattering at low energies. The monograph by Janev[16] is one of the
most widely used data sources available, but only contains charge
exchange at moderate to high collision energies and no data for
elastic scattering. As divertor technology has progressed, low
temperature operating environments have become a standard mode of
operation. Under these conditions, both charge exchange and
smallangle elastic scattering must be included in simulations. Below
approximately 2 eV for D^{+} + D collisions, separate treatment of
smallangle elastic scattering and resonant charge exchange becomes
impossible due to quantum mechanical effects. The need for a
comprehensive treatment of ionneutral elastic collisions at low
energies motivated the Controlled Fusion Atomic Data Center at
ORNL to produce a collection of total differential elastic scattering
cross sections for 31 various isotopic combinations of hydrogen atoms
and molecules[24].
DEGAS 2 represents an improvement over the
original DEGAS code in that the pseudocollision algorithm originally
employed in scoring results[3] is replaced with a tracklength type
estimator. This allows efficient determination of momentum and energy
exchange with the background plasma over a wide range of background
densities. The rate expressions used by the tracklength estimator
involve the momentum transfer (or diffusion)
cross section, σ_{d}[25].
A general
expression for the integral cross sections is[25]
σ_{l} = 2 π  ⌠ ⌡

π
0

( 1 − cosθ)^{l}σ(θ, E_{r}) sinθ d θ, 
 (2.60) 
where l=0
corresponds to the total elastic scattering cross section, l=1 to the
diffusion cross section and l=2 to the viscosity cross
section; θ is the centerofmass scattering angle, E_{r}
is the relative
energy of the collision, and σ(θ,E_{r}) is the differential elastic
scattering cross section. The total elastic collision cross section,
σ_{t}, is used in Monte Carlo simulations to determine the time between
elastic collisions. In this manner the size of the total elastic
collision cross section influences the execution time of numerical
simulations.
A simulation technique that reduces the total
collision cross section while keeping the momentum transfer cross
section invariant can be developed. Using the expression for σ_{d} and
dividing Eq. (2.60) into small and large angle components,
σ_{d} (E_{r}) = 2 π  ⌠ ⌡

θ_{min}
0

C(E_{r})( 1 − cosθ) δ(θ− θ_{min}) d θ+ 2 π  ⌠ ⌡

π
θ_{min}

( 1 − cosθ) σ(θ, E_{r}) d θ. 
 (2.61) 
Substituting σ(θ,E_{r}) from the ORNL data set here and directly
computing σ_{d}(E_{r}) [via Eq. (2.60)],
allows Eq. (2.61) to be solved for
C(E_{r}) given a particular value of θ_{min}.
The resulting value of C(E_{r})
is then used in a similar expression for the viscosity cross section.
This provides viscosity cross sections accurate to second order in
θ_{min}. The value of θ_{min}
is chosen so as to reduce the total elastic
cross section as much as possible while retaining physical realism.
For low energy elastic collisions a significant portion of the
collisions involve smallangle deflections. Due to this behavior, the
majority of the reduction in total elastic cross section is obtained
with a fivedegree minimum scattering angle (Fig. 2.5). Minimum
scattering angles greater than five degrees result in little
computational savings. The resulting reduction in σ_{t}
is approximately
a factor of two at relative energies E_{r} > 10^{−1} eV. The values for
σ_{d}
are unchanged during application of a smallangle cutoff.
Figure
Figure 2.5: Comparison of total and diffusion cross sections for D^{+} + D
collisions with and without the smallangle cutoff approximation.
The
use of this smallangle cutoff leads to an efficient algorithm for
simulating the dynamics of elastic collisions in Monte Carlo
simulations. Previous models based on classical collision kinetics
for direct elastic scattering face difficulties in that the cross
section is a rapidly varying function of the impact parameter[25].
This
precludes the use of efficient tablelookup methods and forces a time
consuming evaluation of the collision integral for each individual
interaction.
Utilization of the differential elastic cross
section, σ(θ,E_{r}),
allows development of a tablelookup interpolation
for both direct elastic scattering and resonant charge exchange.
Following the work by AbouGabal and Emmert[26],
we divide the total
elastic cross section into n pieces and interpret the result as a
cumulative probability distribution function,
P_{i} = 2 π  ⌠ ⌡

cosΘ_{i}
−1


σ(θ, E_{r})
σ_{t}(E_{r})

d( cosθ), i = 0, 1, …n. 
 (2.62) 
Equation 2.62
is inverted to obtain the scattering angle as a function of relative
energy and n equally spaced values of the cumulative probability
spacing, Θ_{i}(P_{i},E_{r}), where Θ_{i}
is the scattering angle in the
centerofmass frame. Examples of Θ_{i}(P_{i},E_{r})
for the D^{+} + D and D^{+} + D_{2}
collisions at a centerofmass energy of 1 eV are presented in
Fig. 2.6.
At this energy, inclusion of direct elastic scattering increases the
probability that colliding particles undergo small angle deflections.
Consideration of charge exchange alone would shift the probability
distribution to favor collisions resulting in Θ = π. As evident in
Fig. 2.6, Θ_{i}(P_{i},E_{r}) for likespecies
collisions exhibits a higher
probability of undergoing largeangle collisions compared to the
ionmolecular collision case. This is due to the absence of a charge
exchange contribution to the scattering function for ionmolecular
elastic collisions. Monte Carlo collision processing begins with the
sampling of a suitable ion from a Maxwellian distribution at the local
ion temperature, yielding, E_{r}. The P_{i}
are next sampled from a uniform
distribution. The collision angle follows from a bilinear
interpolation into this data table without evaluation of complex
integral functions.
Figure
Figure 2.6: Scattering function Θ_{i}(P_{i},E_{r}) for
D^{+} + D and D^{+} + D_{2} ionneutral elastic
interactions at an interaction energy of 1
eV.
Previous treatments of charge exchange in
Monte Carlo algorithms employed a 180^{°} collision scattering model in
the centerofmass frame. The approach is simple to implement,
requiring just the exchange of the ion and neutral velocity vectors
and yields a momentum transfer cross section of about 2 σ_{cx},
where σ_{cx}
is the charge exchange cross section. This is, in fact, roughly equal
to the total diffusion cross section in ionneutral
collisions, σ_{d} ∼ 2 σ_{cx} [27,28,29]
The accuracy of this simple
expression is demonstrated in Fig. 2.7, which compares the
diffusion cross section with results for the "spin exchange" cross
section from the ORNL data[24].
Spin exchange is defined as the exchange
of the intrinsic quantum mechanical property of particle spin; for
distinguishable particles this is equivalent to charge exchange. For
likespecies collisions, spin exchange and charge exchange are
equivalent only within classical regimes. The relation
σ_{d} ∼ 2 σ_{cx}
holds within the classical regime above approximately 2 eV for D^{+} + D.
At lower energies, the particles must be treated as indistinguishable
and separation of direct elastic scattering and charge exchange is no
longer possible. At these low energies, spin exchange is no longer
equivalent to charge exchange and the relation
σ_{d} ∼ 2 σ_{cx} no longer
holds. This illustrates the difficulties in treating ionneutral
elastic collisions by charge exchange alone under lowtemperature
detached conditions.
Figure
Figure 2.7: Comparison of diffusion and charge exchange cross sections for
D^{+} + D elastic collisions.
The addition of these data present the user with yet another option for
handling hydrogen charge exchange in DEGAS 2. Figure 2.8
compares the new CFADC data (labeled "Krstic") with previous
versions and with experimental data. The cross sections have been
plotted as a function of the relative velocity so that both the
hydrogen and deuterium data can be included. The Krstic cross sections
are the "spin exchange" cross sections described in Ref. [24].
They are thus smaller than the total elastic scattering
cross sections that are to be used in DEGAS 2. This underscores the point
that these new elastic scattering cross sections are intended to
replace the existing charge exchange data, not used in addition to
them.
The default data used heretofore
in DEGAS 2 have been computed from the "JanevSmith" fit[18].
This fit was derived from older data published by Barnett[30].
Barnett's data were a compilation of various experimental results.
Included was the work by Newman[31]. Why the
Barnett compilation diverges slightly from the
Newman data at the lowest energies is
not clear. The Reviere[32] curve in Fig. 2.8 represents the
charge exchange cross section employed in the original DEGAS code. Not
surprisingly, they are used in DEGAS 2 in the Degas_box_bench
example. Note that Reviere's intended application was neutral beam
penetration. Hence, the lack of agreement at divertorrelevant
velocities is not surprising.
Figure
Figure 2.8: Comparison of various charge exchange cross sections. The curve
labeled "JanevSmith" represents the default data used in DEGAS 2
prior to the addition of the newer CFADC elastic scattering data,
labeled here by "Krstic" (one curve for H, one for D).
2.10.2 NeutralNeutral Elastic Scattering
Low temperature (a few eV or less) plasma conditions may result in
neutral densities that are comparable to the plasma density. In such
cases, collisions amongst the neutral species can be important in
determining the exchange of momentum and energy in the system. For example, hydrogen
molecules generated at a surface typically have energies comparable to
room temperature, say, 0.03 eV. Yet, these molecules may be surrounded
by much more energetic atoms, say 3 to 5 eV, resulting from molecular
dissociation or charge exchange. Elastic collisions between the two
species can significantly increase the molecules' energy, enabling them
to penetrate much further into the plasma before becoming dissociated
or ionized.
When the neutralneutral mean free path is very short,
K_{n} = λ_{mfp} / d << 1 (d is a measure of the
system size or of gradient scale lengths; K_{n} is the
Knudsen number) cases,
the most efficient approach would be to use a fluid model
for the neutral species. This regime is also known as the "viscous
flow" regime[33]. However, K_{n} is usually close to one in most
fusion problems, even in vacuum regions. In that case, a kinetic
treatment is required. Fortunately, such "transition regime" cases
can be handled effectively with a Monte Carlo code. When K_{n} >> 1,
the "molecular flow" regime, neutralneutral collisions can be
neglected. For a practical discussion
of these three regimes as they apply to the tokamak divertor, see
Ref. [34].
The great difficulty with treating neutralneutral collisions is that they
make the problem nonlinear. Namely, the collision kernel is then a
quadratic function of the neutral distribution function rather than a
linear one. The exact problem is very difficult to solve, and one is
forced to resort to approximations. Reiter et al. addressed this issue
first with the EIRENE code and concluded that an iterative "BGK" approach
(referring to the original algorithm developed by Bhatnagar,
Gross, and Krook[35]; proposed independently by
Welander[36]) was sufficient[37,38]. The BGK
approximation replaces the collision integral with a much simpler
expression in which the distribution function relaxes to a Maxwellian
distribution
with a velocity independent collision time. A physically relevant result
is ensured by determining the parameters of the Maxwellian distribution so as
to enforce conservation of mass, momentum, and energy in the collisions.
The collision time remains arbitrary and is chosen so as to reproduce
measured viscosity and diffusion rates; this procedure is discussed in
greater detail in Sec. 2.10.3.
Neutralneutral collisions are included in DEGAS 2 in exactly the same
way as in EIRENE. While more details on nature of the BGK algorithm
are provided in Kanzleiter's thesis[39], the essential elements
and references are covered in the more readily available Ref. [38].
The BGK algorithm approximates the binary collision integral for species i as

∂f_{i}
∂t

≈ 
M_{i} − f_{i}
τ_{i}

+ 
∑
j ≠ i


M_{ij} − f_{i}
τ_{ij}

, 
 (2.63) 
where the M represent Maxwellian distributions,
M_{α} ≡ n_{α}  ⎛ ⎝

m_{α}
2 πk T_{α}
 ⎞ ⎠

3/2

exp  ⎡ ⎢
⎣

− 
m_{α} ( 
→
v

α

− 
→
U

α

)^{2} 
2 k T_{α}
 ⎤ ⎥
⎦

. 
 (2.64) 
Here, the individual particle velocity is →v_{α}; the mass
m_{α} is always that of the species represented on the
lefthand side of Eq. (2.63).
In Eq. (2.63), a single subscript refers to a physical species.
Dual subscripts below refer to fictitious, "mixed" species introduced
for the purpose of handling binary collisions between unlike neutral
species.
Conservation of mass, momentum, and energy in the selfcollision case
yields
 

 ⌠ ⌡

+∞
−∞

d^{3}v f_{i}( 
→
v

), 
  (2.65) 
 

 ⌠ ⌡

+∞
−∞

d^{3}v m_{i} 
→
v

f_{i}( 
→
v

), 
  (2.66) 
 

 ⌠ ⌡

+∞
−∞

d^{3}v 
1
2

m_{i} v^{2} f_{i}( 
→
v

). 
  (2.67) 

Conservation of mass, momentum, and energy are also enforced in the
the mixedspecies case to determine the M_{ij}. Additional constraints
are obtained by requiring that the ratio of the momentum and energy
relaxation rates is the same as for the full collision integral.
The results are[37,38,39] are
 


→
U

ji

= 
m_{i} 
→
U

i

+ m_{j} 
→
U

j

m_{i} + m_{j}

, 
  (2.68) 
 

k T_{i} − 
2 m_{i} m_{j}
(m_{i} + m_{j})^{2}

 ⎡ ⎣

(k T_{i} − k T_{j} ) − 
m_{i}
6

( 
→
U

i

− 
→
u

j

)^{2}  ⎤ ⎦

. 
  (2.69) 

One additional important constraint on the collision times is obtained,
n_{j} τ_{ij} = n_{i} τ_{ji}. 
 (2.70) 
With τ_{ij} ≡ (n_{j} 〈σv 〉_{ij})^{−1}
representing the typical time between physical species i and
the fictitious species represented by M_{ij},
this becomes just
May[37,38] sets the values of the reaction rates for
mixed collisions using empirical diffusion coefficients;
the rates for selfcollisions
come from measured viscosity coefficients. The end result is that
in both cases the reaction rate 〈σv 〉_{i} ∝ T_{i}^{1/4}
or 〈σv 〉_{ij} ∝ T_{ij}^{1/4}. The validity of these expressions is examined in the
next subsection.
The practical implementation of the algorithm is iterative. By default, the
initial iteration of the code proceeds without neutralneutral collisions.
The reason for this is that the neutral "background" with which the
test species collide initially has zero density. A nonzero density
could be specified during the setup of the background, if desired.
The fluid moments of Eqs. (2.66), (2.67), and (2.67)
are computed at the end of the iteration and inserted into the background
density, velocity, and temperature arrays for the neutral background
species. The calculation is then restarted; this time the test neutral
species scatter against the neutral background. At the end of the run,
the fluid moments are again updated. The iterations continue either a
specified number of times or until some convergence criterion is
satisfied.
Accurately simulating neutralneutral scattering requires greater care
in setting up the geometry and deciding on the number of flights than
is exercised in linear runs. First,
the entire problem space must be divided up finely enough to resolve all
anticipated spatial variations of the neutral density
and the other fluid moments. In comparison,
without neutralneutral
scattering, the vacuum regions can be treated as large monolithic zones.
The variation of the neutral parameters between iterations is asymptotically
limited by Monte Carlo noise. In other words, if too few flights are used,
the background densities are not accurately determined, causing the
elastic scattering of the test species to be doubly inaccurate.
The variations in neutral density, etc. in a given simulation will be more
effectively reduced by increasing the number of flights than by
carrying out more iterations.
2.10.3 Validation of BGK Model and Rates
As described by May[37], the reaction rates for the BGK algorithm
can be determined by matching physical quantities, the viscosity and
diffusion cofficients, against measured values. However, since May's thesis
work and the publication of Reiter[38], the ORNL
CFADC has
computed differential scattering cross sections [24], analogous to those
described in Sec. 2.10.1, for these neutralneutral systems.
In this subsection we will show that the two approaches are largely
consistent, at least near room temperature.
We proceed in a manner analogous to May. Namely, we will obtain from the
CFADC differential cross sections values for the viscosity and diffusion
coefficients. These will be transformed
into reaction rates that can be compared with those described in
Sec. 2.10.2.
For selfcollisions, consider the viscosity. The general expression
derived in Chapman and Cowling[40] [their Eq. (10.1,4)] is
μ_{1} = 
5
8


k T
Ω_{1}^{(2)}(2)

, 
 (2.72) 
where Ω_{1}^{(l)}(n) corresponds directly to the Ω^{l,n}
integrals described by Bachmann[25]. The subscript "1" here
is to be associated with the cross section and refers to collisions of species
"1" with itself (as opposed to collisions between species "1" and
"2" which will be labeled with "12" below).
Likewise, we have for the diffusion coefficient [Eq. (9.81,1) in
Ref. [40]],
D_{12} = 
3
16


k T
(n_{1} + n_{2})m_{r} Ω_{12}^{(1)}(1)

, 
 (2.73) 
where n_{1} and n_{2} are the species densities and m_{r} is the
reduced mass for the system, m_{r} = m_{1} m_{2} / (m_{1} + m_{2}).
The Ω^{(l)}(n) integrals represent averaging of the collision
cross sections over Maxwellian distributions for both species. The
ratecalc code is able to integrate these cross sections over
a single Maxwellian distribution,
namely, the I_{l,n} described in Reiter[17]
and the documentation for ratecalc.
Fortunately, the double integral we need here can be found as a limiting case of
the single integral:
〈I_{l,n} 〉(u_{p},u_{n}) = 
lim
u^{′} → 0

I_{l,n}(u_{p}^{′},u^{′}), 
 (2.74) 
where u_{p}^{′2} ≡ u_{p}^{2} + u_{n}^{2}, and
u_{p} = √{2 T_{p} / m_{p}} is the thermal velocity of the
second species. Likewise, u_{n} = √{2 T_{n} / m_{n}} is the
thermal velocity associated with the Maxwellian distribution used in the
averaging over u. The 〈I_{l,n}〉 is then related to
Ω^{(l)}(n),
〈I_{l,n} 〉 = 8  ⎛ ⎝

u_{p}^{′}
u_{p}
 ⎞ ⎠

n

Ω^{(l)}(n/2). 
 (2.75) 
The same result is quoted in Ref. [25] and in [17],
although in the latter case the variable equivalent to u_{p}^{′}
is incorrectly defined.
With these expressions, we can compute the single species viscosity and
two species diffusion coefficients from the CFADC differential scattering
cross sections. The corresponding expressions for the "Maxwell molecule"
interaction [intermolecular force given by κ/ r^{5}; e.g., see
Eq. (12.1,1) in Ref. [40])] are
μ_{1} = 
k T
3 πA_{2}(5)


⎛ √


, 
 (2.76) 
where A_{2}(5) = 0.436 is the result of a dimensionless integral.
And from Eq. (14.2,2) of Ref. [40],
D_{12} = 
k T
2 πA_{1}(5)


1

, 
 (2.77) 
where A_{1}(5) = 0.422 is the result of another dimensionless integration.
The reaction rates needed for the BGK algorithm are related to the constant
in the Maxwell molecule force law by (see Ref. [41] or [37]),
〈σv 〉_{11} = 2 √2 πA_{1}(5) 
⎛ √


, 
 (2.78) 
and
〈σv 〉_{12} = 2 πA_{1}(5) 
⎛ √


. 
 (2.79) 
Now, we can use Eq. (2.76) to eliminate κ in Eq. (2.78)
in favor of μ_{1}. Then, we can insert μ_{1} as given by
Eq. (2.72) to obtain an effective BGK reaction rate that has the
same viscosity as that computed directly from the CFADC differential
cross sections:
〈σv 〉_{11} = 2.065 Ω_{1}^{(2)}(2). 
 (2.80) 
The analogous procedure for the diffusion coefficient yields
〈σv 〉_{12} = 
16
3

Ω_{12}^{(1)}(1). 
 (2.81) 
Finally, we can use Eqs. (2.74) and (2.75) to
replace the Ω integrals with the equivalent integrals obtainable
from ratecalc,
〈σv 〉_{11} = 
2.065
32


lim
u^{′} → 0

I_{2,4} (u_{p}^{′} = √2 u_{p},u^{′}), 
 (2.82) 
and
〈σv 〉_{12} = 
2
3


u_{p}^{2}
u_{p}^{2} + u_{n}^{2}


lim
u^{′} → 0

I_{1,2} (u_{p}^{′} =  √

u_{p}^{2} + u_{n}^{2}

, u^{′}). 
 (2.83) 
In evaluating the latter expression, we will assume that one of the
species (e.g., molecules at the wall temperature) are much cooler
than the other, i.e., u_{n} << u_{p}. Equation (2.83) then
becomes a function of only u_{p}.
The CFADC data exist only for the D + D_{2} and D + D systems
(and isotopic equivalents; molecules are computationally intensive, hence,
the absence of moleculemolecule data). These are compared with the
May data in Fig. 2.9.
One concern with May's treatment of the D_{2} + D_{2} case is that the
empirical viscosity formula quoted in Ref. [37] does not
appear in the literature, nor is its origin clearly stated. The best
guess for its origin is that it was obtained by rescaling the
"Fuller" expression for the diffusion coefficient[42]
into a viscosity using the theoretical expressions to relate the two.
In Fig. 2.10, we compare May's expression with two that do
appear directly in Ref. [42] ("Chung" and "Lewis") and with
individual data points listed in the CRC Handbook[43]. The
agreement of these values near room temperature ( ∼ 0.03 eV) is
expected. However, May's formula has a stronger temperature scaling than
the others and exceeds them by about a factor of 15 at 10 eV. The
relevant temperature range for this process is expected to be between
1 and 10 eV.
Figure 2.9 shows that at least near room temperature,
May's expressions are also consistent with the CFADC data. However,
in the 1 to 10 eV range, they are considerably larger than the
CFADC values,
suggesting that the values used in DEGAS 2 (and EIRENE) may be too
large. However, we will
demonstrate in Sec. 3.7.2 that the temperature dependence
ascribed to these rates by May is inconsistent with Eq. (2.71).
Consequently, DEGAS 2 now uses constant neutralneutral reaction rates,
i.e., independent of the "background"
temperature. Since agreement of the May expressions with the CFADC
and empirical viscosities is good at room temperature,
and since room temperatures are physically relevant (at least for
desorbed D_{2}) we fix the reaction rates at its
room temperature (300 K = 0.02585 eV) value according to May's formulas.
Figure
Figure 2.9: Comparison of reaction rates for the D + D
_{2} and D + D
systems. The original rates suggested by May[
37] and effective
rates computed from the CFADC data are shown.
Figure
Figure 2.10: Comparison of viscosity formulas from May[
37],
Reid[
42] (the "Chung" and "Lewis" expressions), and the CRC
Handbook[
43].
2.11 Recycling
Recycling of both ions and neutral species in DEGAS 2 is controlled by
the set of PMI ("Plasma Material Interactions", e.g., see
the PMI setup file) and one additional
parameter generically referred to as the "recycling coefficient".
The value of the recycling coefficient is a constant associated with
each wall and target sector (see the sector class).
These are prescribed during geometry definition, e.g., via
the recyc_coef keyword in definegeometry2d.
The best way to understand the recycling coefficient, ℜ,
is to note that P_{a} = 1 − ℜ is the probability for
a flight (the "projectile") striking that sector to be absorbed.
Adsorption, is effectively
one of the plasmamaterial interactions that the
code will choose from in deciding the fate of the projectile.
The second PMI that is always present is referred to as the
"default PMI" since the code defaults to that PMI should
none of the other prescribed PMI (including adsorption) be
selected. Since all PMI must have an
associated probability or "yield", a "default PMI" is
assigned a negative (i.e., unphysical) yield.
The most common default PMI is desorption.
The third PMI we consider here is "reflection" or "backscattering"
in which the incident particle returns immediately to the
plasma with most of its incident energy.
The data associated with the reflection PMI
are used to compute a probability or yield, P_{r}, for
the process.
This probability is, in general, a function of the
projectile energy and angle with respect to the surface normal.
Consider now the common case of an incident hydrogen atom for which
both reflection and desorption (default) PMI's have been specified,
as well as a recycling coefficient ℜ. The probability
assigned to the desorption PMI is
P_{d} = 1 − P_{r} − P_{a} = ℜ − P_{r}. In the case
of unity recycling (the surface is in steady state), the physical
interpretation is that atoms not reflected are adsorbed,
eventually recombine with another atom, and are desorbed as
molecules. In a steady state simulation, the time interval
associated with "eventually" is meaningless, so the process might
as well be instantaneous.
In the case of less than unity recycling, i.e. net adsorption,
the value of P_{d} is reduced accordingly. If ℜ < P_{r},
P_{d} is negative and there is no desorption and the
reflection probability is effectively ℜ.
If ℜ = 0, the sector is equivalent to an "exit".
One operational difference is that flights striking an
exit are terminated immediately. Flights striking a
zero recycling sector will still undergo the process of
selecting a PMI, even though the result will always
be adsorption.
Chapter 3
Running DEGAS 2
3.1 Getting DEGAS 2
DEGAS 2 is now govern by a
licensing
procedure developed and managed by the PPPL Theory and Computation
Department. Prospective users of the code should visit this site and
follow the instructions spelled out there.
The DEGAS 2 source code is available in three formats:
 degas2.tar An uncompressed tar file; extract with:
tar xf degas2.tar
 degas2.tar.Z The tar file compressed with the compress
utility. To extract the code,
uncompress degas2.tar.Z
tar xf degas2.tar
 degas2.tar.gz The tar file compressed with the gzip
utility. To extract the code,
gunzip degas2.tar.gz
tar xf degas2.tar
There are a number of older versions of the code in this directory. The most
recent is always linked to the names listed above. Please consult with the
authors if you think you need an earlier version.
3.2 Getting Supporting Tools
DEGAS 2 relies on several pieces of utility software. In some cases (netCDF),
these utilities are absolutely required. These utilities
can be obtained from the indicated sites. You
will probably want to have your system administrator download, compile,
and install these packages systemwide for your convenience.
3.2.1 GNU Software
The DEGAS 2 Makefile automates several chores, including generating and compiling
FORTRAN code, as well as weaving T_{E}X files from from the FWEB source,
The Makefile also controls the application of platformspecific lines of code
and compilation options. Although not completely necessary, using DEGAS 2
without using GNU Make to run this Makefile would be extraordinarily inconvenient.
The gmake (the GNU Make utility) can be found at
ftp://prep.ai.mit.edu/gnu/make/
Likewise, the Emacs editor is not essential, but
its use simplifies a number of tasks.
For instance, Emacs provides a facility which allows the user to search for variable
and routine names across many files. For example, say you're
editing randomtest.web and you'd like to see the source
code for subroutine next_seed called there. All you
have to do is execute the Emacs command M. next_seed (it
will ask for the name of the tags file; the default value of TAGS
is correct and you just need to hit Enter here). All of DEGAS 2
source and header files are included in this list of tags. Additional
facilities such as a tagqueryreplace (e.g., to change all
occurences of a variable name) are available. See the Emacs Info
file under the node "Tags" for more information.
3.2.2 netCDF
netCDF provides a deviceindependent binary file
format. Furthermore, these files are selfdescribing. Most of the effort
required to read and write these files is handled by macros. These are
documented in the file netcdf.hweb. By adding the commands
in src/.emacs to
your .emacs file, Emacs will be able to convert automatically these files
to text and back to binary for perusal and editing.
netCDF is fully documented in the Info facility of Emacs. The netCDF
libraries are required to run DEGAS 2 since virtually all of the input
data is stored in netCDF files. The associated netCDF utilities ncdump
(for translating a netCDF file into humanreadable text) and ncgen
(for generating a netCDF binary file from a text file) are useful for
reading and making minor changes to DEGAS 2 input files.
Two short scripts, ncdumpfilter and ncgenfilter, allow
Emacs to invoke these utilities via the commands in the .emacs file.
These are not part of the netCDF distribution and are included amongst
the DEGAS 2 source files for your convenience. These should be copied
to a directory which is covered by your shell's PATH variable.
An additional note: if you use ncgen, be sure to either specify
the b or o netcdf_filename option to let the utility know
that you want it to generate a netCDF file (see: man ncgen).
netCDF can be obtained from
http://www.unidata.ucar.edu/software/netcdf/.
The version presently being used for DEGAS 2 development
is 4.1.1. Note, however, that the version of ncgen
that comes with this library has a bug that causes problems
when creating large files. This is remedied in netCDF
V. 4.2, but a workaround with V. 4.1.1 is to use instead
ncgen3.
If you are using ncgenfilter, you would need
use ncgen3 there as well.
3.2.3 FWEB
The FWEB package provides a facility for code documentation as well
as a powerful preprocessor. The typesetting abilities of FWEB allow authors
to insert T_{E}X (or L^{A}T_{E}X) comments into their code. In DEGAS 2, the
preprocessor is utilized extensively for writing macros which allow the
code seen by the user to be significantly more powerful than pure FORTRAN.
In fact, the code can be processed into FORTRAN 77 or FORTRAN 90 as
needed. However, FWEB is not optional. The user must have access
to at least its ftangle executable on some system in order to
generate compilable FORTRAN files (once generated, these can be copied
to other systems for compilation).
For examples of the overall structure of FWEB programs, see the *.web
files; most of the macros are located in the header *.hweb files.
FWEB is fully documented in the Info facility of Emacs.
FWEB (version 1.62) is available at
ftp://ftp.pppl.gov/pub/fweb/fweb1.62.tar.gz
To typeset DEGAS 2 source code in "humanreadable" formats such as DVI,
PostScript, or PDF, the user will need FWEB's fweave utility and
a version of T_{E}X or L^{A}T_{E}X.
L^{A}T_{E}X is available in lots of places. To begin with, see the
introductory material and pointers on the
official CTAN T_{E}X network page at
http://www.tex.ac.uk/ctan. The
recommended distribution for Unix (including Linux) systems
is T_{E}X Live (http://www.tug.org/texlive/.
T_{E}X Live is also available for Windows.
3.2.4 Triangle
Jonathan Shewchuck's
Triangle program is used by the definegeometry2d
code to break polygons up into triangles. The home for Triangle
is
http://www.cs.cmu.edu/ quake/triangle.html
The Triangle V. 1.6 was used in the creation of
the definegeometry2d examples described in
Sec. 3.7.3.
Anyone downloading a more recent version should notify the DEGAS 2
authors. The DEGAS 2 Makefile will expect to find
the Triangle files, triangle.c, triangle.h, etc.,
in $HOME/Triangle. If your copy is installed elsewhere, you can
either create an appropriate link or modify Makefile.
3.2.5 Silo and HDF
DEGAS 2 produces no direct graphical output. Instead, the geomtesta
utility (see Secs. 3.4.1 and the
source code documentation.) generates
deviceindependent binary files in either the HDF (more specifically,
HDF Version 4) or Silo formats that can be processed by
external graphical packages, such as VisIt or IDL. The VisIt
package in particular is freely available
(see https://wci.llnl.gov/codes/visit/home.html), can read Silo files directly, and can generate
both 2D and 3D images. A similar alternative to
geomtesta, ucd_plot, is now available. It's
more restrictive in that it requires the Silo library
and works only with 2D geometries created by
definegeometry2d. On the other hand, the data
contained in the Silo file are effectively
represented on the actual DEGAS 2 mesh instead of
a pixelated representation. Visualization of the
ucd_plot Silo has only been performed with VisIt.
The Silo library can be found at https://wci.llnl.gov/codes/silo/index.html. It can use either its own
PDB I/O driver or HDF5. In the latter case, the HDF5 library
must also be installed. For HDF5 and / or HDF4, see
http://www.hdfgroup.org/.
Note that as of HDF version 4, some of the netCDF functionality is
implemented in HDF. To avoid conflicts with the official netCDF libraries,
do not install the netCDF portion of the HDF distribution.
Note that definegeometry2d also produces an HDF or Silo
file, geomtestc.hdf or geomtestc.silo, as part of
its test of the geometry.
3.2.6 Concurrent Version System
The UNIX CVS (Concurrent Version System) utility
is used to coordinate work on a
code by a number of authors.
The central DEGAS 2 source code is maintained
in a CVS repository. Users planning to coordinate with the code authors should
become familiar with CVS. Users who will do extensive code development on their
own will find CVS a valuable tool in tracking their work. The average enduser
should be able to work effectively without ever dealing with CVS.
The CVS repository is a central, safe version
of the code; authors check out copies in their own disk space in order to
work on the code. By doing a cvs update, a given author can update his
local sources to incorporate changes others have saved to the central
repository. The cvs commit command allows him to save his own changes to
the repository.
The Info page for CVS is the best reference for commands. The manual
page is complete, but lengthy. The man page for
rcsintro may also be helpful as an introduction. Many cvs commands
may be executed from within Emacs. See the Emacs Info page on
Pclcvs for more details.
3.3 Structure of the DEGAS 2 File System
Here is an example structure of the directories under the toplevel
degas2 directory. Those labeled with "(D)" represent
source and data directories maintained with CVS. The others are created
by the user as described in the compiling section. Note that the Makefile
assumes that this degas2 directory sits in the directory pointed to
be the $HOME$
environment variable on your workstation. If you have installed
the code elsewhere, see Sec. 3.6.3.
 data
 Input atomic and surface physics data and information (D)
 src
 All source code (D)
 dg_carre
 Source and executable files for DG and Carre (D)
 Doc
 Location of the various formats of this document and ancilliary
files (D)
 examples
 Input and output files from example and benchmark runs (D)
 Aladdin
 Source code and data for use with IAEA Aladdin package (D)
 ALPHA
 Object and executable files for Dec Alpha machines
 LINUX
 Object and executable files for Linux machines
 LINUX64
 Object and executable files for 64bit Linux machines
 SUN
 Object and executable files for Sun Solaris and SunOS machines
 CRAY
 Preprocessed source code for export to Cray computers
 tex
 Intermediate files and final DVI documentation files generated
from source code
If you need to have multiple run directories for a given system (e.g., focussed
on different problems or using different compilers), you can create additional
directories with names like SUNfoo where "foo" is any string you like.
3.4 Components of the Code
DEGAS 2 is not a single code, but a complex package of
setup, test, simulation, and postprocessing tools. Which ones are
needed depend on the user's objectives and familiarity with DEGAS 2.
A few additional "targets" are listed in the Makefile. These are
either sufficiently out of date or infrequently used to warrant
mention here.
3.4.1 Main Code
Virtually every user will be running these executables,
usually in this order.
 problemsetup
 Reads a text file describing the species,
materials, atomic, and surface physics required for the
current problem. The code
accesses the required "reference" data, and compiles them
into a netCDF file (the problem.nc file).
 definegeometry2d
 Generates a DEGAS 2 netCDF geometry file
from from a text
input file. Input files for simple geometries may be
generated manually; more complex geometries can be specified via other computer
generated text files
files referenced by the input file. For example, definegeometry2d is
able to read data produced by the DG package
(formerly part of the B2Eirene or
SOLPS distribution, but now included with
DEGAS 2) or by UEDGE. Of the three tools for setting up
DEGAS 2 geometries definegeometry2d is the most flexible.
 readgeometry
 Generates a DEGAS 2 netCDF geometry file from an
existing DEGAS, UEDGE, or SONNET geometry description. This package has
largely been superseded by definegeometry2d. However, users with legacy
input files for the old DEGAS code will still want to use readgeometry.
A short text input file helps readgeometry
through the transformation process.
 defineback
 Uses a short text input file to control the mapping of plasma
data in external text files onto the DEGAS 2 geometry. Designed to be used in
conjunction with definegeometry2d, defineback can utilize data
produced by other codes, including UEDGE. In particular,
defineback can be run in an iterative mode with UEDGE
analogous to that provided by updatebackground (see below
and in the source file
defineback.web).
 readbackground
 Plasma data in an external file is mapped onto
the zones of the DEGAS 2 geometry. Presently, only DEGAS and UEDGE files
can be used; readgeometry and readbackground read the same data
file.
 updatebackground
 Is similar to readbackground, but is intended
for use in iterations with the UEDGE code only. The source
information generated by the initial run of readbackground is
stored in a separate file (oldsourcefile,
Sec. 3.5.1). On subsequent
iterations, small changes in the source terms are accounted for
by allowing nonunit weighting factors at each source segment.
Once the accumulated changes exceed a certain size (specified by
parameters in the sources class), the DEGAS 2 source sampling
arrays are reinitialized (i.e., to once again use unit weights)
and the oldsourcefile is rewritten. The effectiveness of this
procedure is currently being evaluated and should be considered
as experimental.
 tallysetup
 This code specifies the types of "scores" (the
essential output of the code), along with their
dependencies, which will be computed in the DEGAS 2 run. Although the input
files(s) should not be modified by the novice user, this code needs to be
run to factor the values of physics and geometry parameters for the
problem at hand into the specification of these scores.
 flighttest
 The core of DEGAS 2. This program launches, follows,
and scores the Monte Carlo trajectories which make up a DEGAS 2 simulation.
(The name is a holdover from early versions in which this did tracking
with a minimum of hardwired physics).
 outputbrowser
 Reads the output netCDF file generated by
flighttest (in addition to all other netCDF files for the
problem) and allows the user to interactively browse its contents.
A script facility is provided for batchlike operation.
 geomtesta
 Another inappropriately named routine which serves as
the principal postprocessing code for DEGAS 2. Although originally
designed as a diagnostic for the geometry, it has proven too easy to
continue extending this code to warrant development of an honest
postprocessing tool. The DEGAS 2 input and output files are read in by
the code. Using the geometry information, the zonebased plasma and neutral
data are transcribed onto a high density uniform rectangular mesh and
dumped into Silo or HDF files suitable for further manipulation by
external graphics packages such as VisIt or IDL.
 ucd_plot
 Alternative postprocessing code that uses DEGAS 2's
own mesh for representing output data. Creates a Silo format file
for visualization with VisIt.
 degas2_xgc.a
 This library contains the core routines of DEGAS 2
and an interface for a subroutine based coupling to a plasma code.
The primary quantities exchanged between the codes are moments of the
plasma and neutral distribution functions. An interface into the
DEGAS 2 atomic physics routines is also provided. This library was
designed for and is used with the XGC0 guiding center neoclassical
particle transport code as part of the Center for Plasma Edge
Simulation [44].
3.4.2 Other Setup Routines
These are occasionally used to add new atomic
or surface physics reactions, or to generate new atomic and surface physics
data files.
 datasetup
 This routine possibly belongs in the first list. It
reads text files describing the complete lists of elements, species, reactions,
materials, and plasmamaterial interactions available to DEGAS 2 and generates
the corresponding DEGAS 2 netCDF files. These lists are occasionally referred
to in the code as the "reference" lists; the input file for problemsetup
specifies subsets of these lists. A new reaction is added to DEGAS 2 by
inserting the requisite information into the "reactionfile"
and rerunning datasetup. The new reaction can
then be added to the problem input file.
 ratecalc
 A routine for computing averages of atomic physics
cross sections. Since it is designed to read data from the Aladdin
database, ratecalc is presently limited in its flexibility. The principal
application thus far has been the generation of reaction rates and higher
moments for charge exchange and other interactions between atoms, molecules,
and ions.
 reactionwrite
 A simple routine designed to translate collisional
radiative data for the electron impact ionization of hydrogen from the
format used by the old DEGAS code into netCDF files (ionization and
recombination) for DEGAS 2. This code has been run a handful of times
as needed to keep up with changes in the content and format of the atomic
physics data files.
 pmiwrite
 Is analogous to reactionwrite, but handles almost all
of the plasmamaterial interaction data. Presently, two data files from
the old DEGAS code and one from EIRENE are used as input. A few processes
describable without external data are also included. This code can be used,
in a somewhat clumsy manner, to add new plasmamaterial interaction data
to DEGAS 2 using either one of these three existing data bases or via
a new datafile (with the rest of the code serving as an example).
 boxgen
 Is the third of geometry / background generation tool provided
with DEGAS 2. It is intended to serve as an example of how to generate
geometry and background files from scratch. As the name implies, boxgen
sets up a simple box geometry with a linearly varying plasma. However,
no input file for boxgen has been developed; all modifications have to
be effected through direct source code modification. The typical user
in search of a simple linear geometry to study can be overwhelmed by the
complexity of the code encountered there. For this reason, the
definegeometry2d / defineback pair may be more suitable for
setting up simple problems.
3.4.3 Test Routines
 reactiontest
 Can be run after problemsetup
to check the computation
of the reaction rate, and collision products (including velocities and
scoring data) for a given reaction.
 pmitest
 The analog of reactiontest for plasmamaterial
interactions. Again, this code can be executed for a given problem once
problemsetup has been run. Since the velocity distributions of some
plasmamaterial interactions are described statistically, pmitest
is also set up to average the outcome of a given interaction over a
specified number of trials and print out a list of the product velocities
for external manipulation.
 sourcetest
 Tests the distribution in phase space of the source
terms for a particular problem. Since the source is described in the
"background" netCDF file, defineback,
readbackground, boxgen, or something equivalent
must be run prior to using this routine. Like pmitest, the code is able
to run many trials, compute averages, and printout data from each
individual trial.
 dataexam
 Provides the user with a convenient interface to the
atomic physics data files. Although netCDF files can be converted to
a humanreadable text form, interpreting those data is since DEGAS 2
collapses the multidimensional structure components of the atomic physics data
into a single onedimensional array. This routine reads the "raw" (at
the "reference" level) data files prior to their insertion into the
problem netCDF file and, thus, can be run at any time. An analogous
routine for the plasma material interactions is needed.
 sysdeptest
 Tests several systemdependent features of DEGAS 2.
This would be run only if the code was being ported to a new operating
system or if the operating system on a currently supported architecture
underwent a signficant upgrade.
 randomtest
 Is used to verify that DEGAS 2 will give the same
results on
different architectures (provided the random number seeds are set the
same!).
3.4.4 Miscellaneous Routines
 matchout
 Is used to compare the current DEGAS 2 output netCDF
file with another from an equivalent run (i.e., all of the arrays must
be exactly the same size) on the command line. The most frequent use
of matchout is to verify that two separate runs have yielded
results which are the same to within roundoff error. Because of the size of
the ouptut netCDF file, verifying this by visually comparing the files
is impractical.
 matcheir
 Is specifically designed for comparing with output from
the EIRENE code. The name of the EIRENE file is specified on the command
line. See Sec. 4.3.
 datamatch
 Provides a means for comparing the contents of the two
netCDF data files specified on the command line.
3.4.5 DG and Carre
The graphical interface DG code is used here
in conjunction with definegeometry2d, as
in B2Eirene (its original application), to simplify specification of
2D divertor tokamak geometries . Since DG is also available
separately from DEGAS 2, its manual, dg.pdf is contained in the
dg_carre/DG/DOC directory. Carre[14] generates 2D
quasiorthogonal meshes of the sort used by B2 and UEDGE. This
version of Carre includes pre and postprocessing scripts that
simplify its use with DG. A number of utility programs for
manipulating magnetic equilibrium files also come with this
DG distribution; see the DG manual for details.
3.5 Input Files
The essential inputs to a Monte Carlo neutral transport code are:
 The identity of the "test" species (usually neutral in plasma problems)
that are of interest,
 The identity and phase space distributions of the "background"
(typically electrically charged in plasma problems)
species with which
the test species will interact,
 The interaction processes between the test and background species.
The relative probability of these processes as a function of the
test and local background properties are required as well as a
prescription for specifying the outcome of each interaction.
 The simulation geometry, including a framework for
holding the background species data, as are a means for tracking
the test species through the background.
 The boundaries of the geometry, including a prescription
for specifying the outcome of interactions between the test species and
these boundaries.
DEGAS 2 utilizes several text based input files for controlling the actions
of its component executables (see Sec. 3.4).
In this section we will discuss a few of these input files. The more complex
input files are described in their respective codes; links to the documentation
for those files will be given instead.
More
detailed information on running DEGAS 2 will follow in later sections.
These text input files are, for the most part,
read with the same set of text processing utilities that adhere to the following
guidelines. In some cases, particularly
those involving files generated by other codes, the files are read directly
with formatted FORTRAN read statements.
 Spacing and blank lines are ignored. The user is free to utilize white
space in whatever way facilitates reading and maintaining the input file.
 Comments, either a complete line or at the end of a line, are begun
with a # sign.
 File names can be a relative or full path name.
3.5.1 degas2.in
degas2.in is the input file which controls (most of) the other
input files. Each line consists of DEGAS 2's symbolic name for the
file (this can be changed only by modifying readfilenames.hweb
and readfilenames.web) followed by the path name for the file
to be used in the current problem. The order of the lines in this file
is unimportant.
elementsfile ../data/elements.nc
backgroundfile bk_uers.nc
geometryfile ge_uers.nc
problemfile pr_uers.nc
reactionfile reactions.nc
speciesfile ../data/species.nc
aladinfile ../data/aladinput.nc
aladoutfile ../data/aladoutput.nc
elements_infile ../data/elements.input
problem_infile pr_uers.input
reaction_infile reactions.input
species_infile ../data/species.input
materials_infile ../data/materials.input
materialsfile ../data/materials.nc
pmi_infile ../data/pmi.input
pmifile ../data/pmi.nc
cramdproblemfile cramdprob.nc
tallyfile tally_uers.nc
outputfile degas2_uers_out.nc
oldsourcefile os_uers.nc
snapshotfile sn_uers.nc
Note that some of these are
outdated or rarely used (unused entries may be deleted from the file, if
desired);
see the readfilenames class for
more specific information. Most of the files come in pairs with the name
XXX_infile corresponding to a text input file and
XXXfile being a netCDF file generated by a program
like datasetup, problemsetup, etc.
3.5.2 elements_infile
elements_infile lists all of the elements available to DEGAS 2 for
use in constructing the species (see Sec. 3.5.3).
Additional detail is given at the beginning of the
file elementsetup.web.
3.5.3 species_infile
Inputs for species are contained in the file with the symbolic
name species_infile. Additional detail is given at the beginning
of the file speciesetup.web.
3.5.4 reaction_infile
The file reaction_infile describes all of the reactions available in
DEGAS 2. Additional detail is given at the beginning of the
file reactionsetup.web.
Note that adding a new reaction to this file involves two tasks beyond
inserting the appropriate lines in reaction_infile. First, the
netCDF file for the atomic physics data must be generated (see
Sec. 3.9). The
second task would be to write subroutines for setting up the products and
handling collisions (see Sec. 3.9.3).
This would be necessary only if a new reaction type
were being added.
3.5.5 ratecalc Input
The ratecalc code is used to compute atomic physics reaction
rates from cross section data, as well as a few other tasks.
The code is documented internally. Follow this
link to
the introductory section to learn more.
3.5.6 materials_infile
The materials described in materials_infile are essentially just labels
which will be used in conjunction with the plasma material interactions
(see Sec. 3.5.7) to specify how test
species interact with nontransparent
surfaces. Additional detail is given at the beginning of the
file materialsetup.web.
3.5.7 pmi_infile
The file pmi_infile describes all of the plasmamaterial interactions
(PMI) available in DEGAS 2. Additional detail is given at the beginning
of the file pmisetup.web.
3.5.8 problem_infile
As noted already in the section describing the DEGAS 2 components
(see Sec. 3.4), there are two levels of physics input to DEGAS 2.
The input files noted thus far comprise the reference level: in
practice, the sum total of the data available to the code (although
in principle it could be smaller). The second level is the problem level
of data. This prescribes the species, reactions, materials, and
PMI which will serve as the physical
model to be used in carrying the simulation at hand. In some
parts of the (internal) code, these are also referred to as the
subset data since they represent a subset of the reference data
(see also Sec. 3.9 and Fig. 3.2).
Additional concepts alluded to above in connection with the
reaction input file (see Sec. 3.5.4) are those of
test species
and background species. Most simply, test species are the ones
DEGAS 2 will track as they collide off of background species.
The use of the word species here is important: both of the test
and background lists are subsets of the "species" list
(see Sec. 3.5.3). One more precise distinction between the
two species types is that we assume that we know the distribution
function (in space and velocity) of the background species; in
fact, such information is required input to DEGAS 2. On the other
hand, we are attempting to compute the test species
distribution function, moments of which serve as the primary output
of DEGAS 2.
More information about the input file can be found at the beginning of
the file problemsetup.web.
3.5.9 readgeometry Input
The documentation for readgeometry input is maintained in the
source code, readgeometry.web. For readers of the PDF version
of this manual, here's a
link to the
corresponding PDF file.
3.5.10 definegeometry2d Input
The documentation for definegeometry2d input is maintained in the
source code, definegeometry2d.web. For readers of the PDF version
of this manual, here's a
link to the
corresponding PDF file.
3.5.11 defineback Input
The documentation for defineback input is maintained in the
source code, defineback.web. For readers of the PDF version
of this manual, here's a
link to the
corresponding PDF file.
3.5.12 readbackground Input
Currently, the only possible inputs to readbackground are the
UEDGE (a specifically formatted text file generated during
UEDGE postprocessing) and DEGAS formats (an input file for
the old DEGAS code). This routine makes some specific assumptions about
the contents of these files. Since this situation is unsatisfactory from
a number of viewpoints, a better longterm approach is being contemplated.
For that reason, no additional documentation is provided at this point.
3.5.13 tally_infile
The documentation on the input file for tallysetup is maintained
in the source code, tallysetup.web. This
link will take the reader to the corresponding
PDF file. The tally class contains more
extensive and detailed information.
3.5.14 outputbrowser Input
This postprocessing utility can be run interactively or with an
input file. When run interactively, the session is logged into a
script file called outputscript (in the run directory) which
can be subsequently used as input to outputbrowser. Details on
the use of outputbrowser and on the format of the script
file can be found in the introduction to the
file outputbrowser.web.
3.5.15 geomtesta Input
The documentation for geomtesta is maintained in
the source code, geomtesta.web. For readers of the
PDF version of the manual, here's a
link to the corresponding PDF file.
3.5.16 ucd_plot Input
The documentation for ucd_plot is maintained in
the source code, ucd_plot.web. For readers of the
PDF version of the manual, here's a
link to the corresponding PDF file.
3.6 Compiling DEGAS 2
The Makefile has been written to compile on a variety of architectures
and operating systems. Those that have been used in thus far,
and their corresponding
symbols (in brackets) are:
 Linux 32bit [LINUX]
 Linux 64bit [LINUX64]
 Sun Solaris [SUN],
 IBM AIX [IBM],
 Digital OSF1 [ALPHA],
 Silicon Graphics IRIX [SGI],
 Others (e.g., [CRAY]) are available for crosscompilation.
Note that due to the proliferation of relatively cheap and fast
Linux clusters, usage of DEGAS 2 on other systems has essentially
ceased and support for those systems cannot be guaranteed.
The most frequently used compiler on Linux systems
is Portland Group F90 compiler (PGROUP in
the Makefile), albeit with F77 format code
(FORTRAN90=no). The Pathscale (PATHSCALE in
the Makefile) and gfortran (version 4.3) have been
successfully tested with the F90 formatted code (in
DEBUG mode only). The use of
the F90 option eliminates the need for system dependent code
so that other compilers can be used without modifications to
the sysdep.hweb or sysdep.web files; only the
compiler
related variables in the Makefile or Makefile.local
need be defined. On the other hand, some F9x compilers (e.g.,
LaheyFujitsu F95) may result in code that runs very slowly.
By default, the Makefile assumes that the DEGAS 2 main directory
exists at $HOME/degas2. If you unpacked the tar file
in some other location, you'll need to tell the Makefile where to look
(see Sec. 3.6.3).
Blindly change directory names is not recommended.
Some flexibility is again provided; see below for details.
3.6.1 Basics
Let's start by compiling the randomtest utility (see
Sec. 3.4.3) on a Sun system. To use other systems,
replace SUN with the appropriate symbol from the above list.
First, change to the main DEGAS 2 directory and create the SUN
subdirectory:
cd ~/degas2
mkdir SUN
cd SUN
cp ../src/Makefile .
Once the Makefile is present here, it will "update itself" with respect
to the copy in the src directory when needed. To get an idea of
how ftangle works, just "make" the main FORTRAN file for
randomtest:
gmake randomtest.f
At some point, you may want to visually compare this file with the
original source code randomtest.web (in the src directory, with
the other source files) to gain an appreciation
for the amount of work the FWEB macros do in DEGAS 2. For more information
on these macros, see the documentation in array.hweb and the
other header files.
Now, finish making randomtest:
gmake randomtest
To run randomtest, just type
./randomtest
The output to the right of the "=" sign should match the numbers in
parentheses.
3.6.2 Making Documents
Generating "woven" (i.e., using FWEB's fweave utility)
documentation is similar:
cd ~/degas2
mkdir tex
cd tex
cp ../src/Makefile .
gmake randomtest.dvi
The DVI file can then be printed or viewed (using xdvi), as desired.
In fact, the Makefile also provides additional "targets", e.g.,
randomtest.print (uses lpr to print to your default printer),
randomtest.view (launches xdvi) for you, and
randomtest.ps (a PostScript file generated via dvips). You
can generate woven documents for all of the ".web" source files.
For some of the more useful ones, see Sec. 3.13.
You could also just generate the T_{E}X file (again, from the degas2\tex
directory),
gmake randomtest.tex
You could then use latex, pdflatex, or whatever other T_{E}X application
you wished.
3.6.3 Adjustments to Makefile
A number of the default settings in the Makefile can be overridden by
placing the desired values in a file called Makefile.local in the
working (e.g., SUN) directory. Some of the more frequently
changed variables are:
 DEBUG=no turns on optimization. The default DEBUG=yes
is required for debugging.
 DEGASROOT=somedir will tell Makefile to find the main
DEGAS 2 directory at somedir/degas2. somedir should be an
absolute path name. E.g.,
DEGASROOT = /u/somedir
 FORTRAN90=yes switches from using the default FORTRAN 77
compiler to FORTRAN 90. Be aware that only some of the FORTRAN 90
compilers available work satisfactorily. Some will compile the code,
but run slowly.
 FCF77, FCF90 tell Makefile which "normal" (nonMPI)
compiler to use.
 MPI=yes will enable the MPI commands in the source
code and direct the Makefile to use the MPI compiler flags.
 FCMPI77, FCMPI90 tell Makefile which MPI
compliant compiler to use.
 NETCDFV2=yes tells DEGAS 2 that your system has only the
older version 2 of the netCDF library (by default, DEGAS 2 expects version 3).
 The Makefile will accomodate more than one working directory on
a given machine provided the additional directory has a name like
SUNjunk, where "junk" is some string meaningful to you.
 If you would like to have more than one source directory, we
recommend creating a new degas2 directory heirarchy somewhere else
in your file system and use the DEGASROOT variable to get
the Makefile working there (you can set up symbolic links to common
directories such as data).
3.6.4 CrossCompiling
The mechanism for crosscompilation has been updated to use
ssh and scp and no longer relies on the presence
of AFS. However, this section of the
manual has not been correspondingly updated.
Please contact the code authors if you need this capability.
As an example of how to crosscompile, here are the procedures for
making randomtest on the NERSC Crays. This assumes that you
have access to AFS from your local workstation (and at NERSC, of
course):
mkdir <someafsdirectory>/degas2/CRAY
ln s <sameafsdirectory>/degas2/CRAY ~/degas2/CRAY
cd ~/degas2/CRAY
cp ~/degas2/src/Makefile .
gmake randomtest
This will use ftangle on the local workstation to generate
the FORTRAN source code files from the FWEB files in the src
directory. It is advisable to set up a Makefile.local
(see Sec. 3.6.3) in this
directory that contains the same flags that will be used on the
remote machine, particularly DEBUG, FORTRAN90
and MPI.
Then, log on to the remote machine and do:
ln s <sameafsdirectory>/degas2 ~/degas2
cd <someworkdirectory>/CRAY
cp ~/degas2/CRAY/Makefile .
make randomtest
In addition to the usual variables in Makefile.local, one
needs to add STANDALONE=no. This will tell make that
the FORTRAN (e.g., .f) files have already been generated on
your local workstation, and that it can find them in
~/degas2/CRAY on the remote machine, which should be
linked back to your local workstation via AFS.
The degas2/data directory will likely also be needed on the
remote machine. You can manually copy the directory and its contents
to the remote machine (e.g., via ftp)
or keep a copy in the AFS directory and use
a link to allow the code on the remote machine to find the files.
E.g., for the second option, starting on the local workstation:
cd ~/degas2
cp R data <sameafsdirectory>/degas2
And on the remote machine
cd <sameworkdirectory>
ln s <sameafsdirectory>/degas2/data data
Be aware that in either case, the data directory could eventually
get out of sync with the one on your local workstation. If this
is a likely possibility, one could consider also linking the
data directory on the local workstation to the one in the
AFS directory, making the latter the only real copy.
3.6.5 Miscellaneous Targets
The Makefile provides some other occasionally useful targets. Please
see the Makefile for more details.
 clean removes source, object, and other unessential
files. This is pretty thorough and indescriminate; use with caution.
 TAGS updates the Emacs tags table (see Sec. 3.2.1).
You may need to do this if you subtstantially altered any of the source
files or have other reason to believe that the tags table is outofdate.
You would run gmake TAGS from the src directory only.
 depend updates the Makefile dependencies
(in Makefile.depends). Do this if the header file dependencies
of one or more executables have been changed or if you believe that the
dependencies are out of date. Note that source code (i.e., object files)
dependencies are
explicitly stated in the Makefile and must be updated by hand.
 foof dumps out values of some of the internal Makefile
variables. This is useful for debugging Makefile problems.
3.6.6 CVS Access to DEGAS 2
Properly privileged (in the UNIX sense) authors of DEGAS 2 will be given
access to the
CVS (see Sec. 3.2.6) repository for the code. With this access,
one can obtain the latest version of the source with:
setenv CVSROOT /afs/pppl.gov/common/cvs
cd
cvs co degas2
(this will place the code in ~/degas2; see also Sec. 3.6.3).
Of course, the command setting the CVSROOT variable can be more
conveniently executed from a login (e.g., ~/.login) or shell startup
script (e.g., ~/.cshrc). The user will also have access to
all older versions of the source and data files through the appropriate
CVS commands. More importantly, the user will be able to update
the CVS repository. For this reason, such access is controlled by
the primary authors of DEGAS 2. Contact them if you feel that you can
contribute to the DEGAS 2 source code
(see Sec. 5.1).
3.7 Examples
The examples directory contains several documented DEGAS 2 examples,
with input output, and auxiliary files. They are included not only
for instructive purposes, but also for verifying that successive versions
of the code give the same results. Namely, the user should be able to
download the code and generate exactly (to within roundoff error) the
results contained within each of the example subdirectories. Any
discrepancies should be reported to the primary code authors (see
Sec. 5.1).
3.7.1 Analytic_fluid_bench
This simple run demonstrates that DEGAS 2 matches the results of an analytic
model in the fluid limit. More information on this comparison can be
found in the IAEA proceedings[15] or on the Web as a
PostScript
or
PDF file.
To run this example, follow these steps:
 Set source code switch. The source code file boxgen.web
has within it a few different settings available which can be selected
by changing the value of the flag BOXRUN. With this file open in
an editor, change the line
@m BOXRUN 0
so that it now reads
@m BOXRUN 32
Save the file and exit the editor.
 Switch to the example directory. If the main DEGAS 2
directory sits within your home directory:
cd ~/degas2/examples/Analytic_fluid_bench
 Copy the degas2.in file into your working
directory. E.g., if you are working on a Sun,
cp degas2_boxgen.in ../../SUN/degas2.in
cd ../../SUN
(overwriting any degas2.in file you may have had there already!).
 Prepare reference data. This example uses the standard
reference data files in the data directory. Hence, this step can
be skipped if
those files have not been altered since you downloaded the code.
gmake datasetup
./datasetup
 Prepare problem data. This and subsequent steps are
not optional! Note that the degas2.in file points back to
the examples directory for the problem input file (see
Sec. 3.5.8).
gmake problemsetup
./problemsetup
 Generate geometry and background files. The boxgen
program takes care of both of these:
gmake boxgen
./boxgen
 Set the number of flights. The background file bk_boxgen.nc
contains the specification of the number of flights to be run. Open this
file in an editor (Emacs will be the most convenient) and find the line
source_num_flights = 100 ;
Change the "100" to either "1600" or "6400". The example contains the
output for both cases. Save the file and exit the editor.
 Set up the tallies.
gmake tallysetup
./tallysetup
 Run the code.
gmake flighttest
./flighttest
 Check the text output. The file density.out can be directly
compared with either density_1600.out or density_6400.out,
depending on the number of flights. The numbers in here have only five
digits of precision. If an exact match is not obtained, something is
amiss. Please contact the code authors (see Sec. 5.1) if you feel
that a problem exists with the code as you downloaded it.
 Check the binary output. You can use the matchout
utility (see Sec. 3.4.4) to do this.
 Compare with the analytic solution. The file soln_32
contains the columns:
 x
 Distance along the problem space in meters.
 N1(x)/N(0)
 Relative density (to density at x=0) predicted by
one analytic model (see the above references) solution.
 N2(x)/N(0)
 Relative density predicted by a second analytic model
solution.
 Ti(x)
 Ion temperature profile.
 Ni(x)
 Ion density profile.
The columns of the density.out file contain
 First zone number. These zones correspond directly to the x
values in the soln_32 file.
 Second zone number. This is always 0 since this is a 1D problem.
 Neutral hydrogen density in m^{−3}. To normalize this column,
divide the whole column by the value in the first row. Since this corresponds
to the first value of x in soln_32 and not to x = 0,
multiply the whole column by the relative density from the analytic
solution at the first value of x, normalizing the DEGAS 2 result
to the analytic solution at one point.
 Relative standard deviation of the neutral density (see
Sec. 4.4).
 Neutral hydrogen pressure in pascals.
 Relative standard deviation of the neutral pressure.
 The last two columns are no longer used and should be filled with
zeroes.
3.7.2 NeutralNeutral Scattering Examples
The initial implementation of the BGK algorithm for handling neutralneutral
collisions (see Sec. 2.10.2) in EIRENE
has been tested by comparison
with analytic expressions for the time dependence of the relaxation of
a distribution function (both in selfcollision and
mixedspecies collision systems) and against models describing the
Couette flow problem.[37] The relaxation benchmarks have been
repeated with DEGAS 2, but have not been formally made into an example
run and will not be discussed here.
Further details can be provided upon request. However, the Couette flow
problem can be set up with only a few minor modifications to the code and
does make a suitable example.
Couette Flow
The Couette flow problem of fluid mechanics involves the flow of fluid between
two parallel, sliding plates. The fluid is assumed to have "no slip"
boundary condition at the plates. The viscosity of the fluid drags
along adjacent fluid elements, resulting in a velocity gradient
between the boundary velocities represented by the two plates.
In this example, two semiinfinite plates are a distance d apart
in the x direction.
Periodic boundary conditions are enforced in the y and z directions so that
the problem is effectively infinite in those directions. Particles
are initialized at the x = 0 plate with a thermally distributed
velocity to which a constant velocity in the y direction,
v_{y0},
is added; this represents the velocity of the plate. The second plate
at x = d is treated as stationary. Particles striking it are assumed to be
reemitted with a thermal distribution pointed in the −x direction.
Particles striking the x = 0 plate are absorbed. Hence, in the no
collision limit, particles make only one round trip across the box.
The choice of the distributions used at the plates is crucial
to reproducing the analytic models. In particular, the "Maxwell flux"
distribution [12] must be used. This distribution describes a
recycling thermal flux
of particles moving in a particular direction (in this case, perpendicularly
away from the wall). This is the distribution used in some of the
PMI data files, such as in h2_des_maxw_mo. In this example, however,
the distributions of the particles coming off of the two plates
are enforced through hardwired code,
not through the data files.
The initial velocity is chosen from the distribution
Q( 
→
v

) ∝ v_{x} exp  ⎧ ⎨
⎩

− 
m
2 k T

[v_{x}^{2} + (v_{y} − v_{y0})^{2} + v_{z}^{2}]  ⎫ ⎬
⎭

. 
 (3.1) 
Note the plate velocity v_{y0} in the exponent. For particles
"reflected" at the x = d surface, v_{y0} ≡ 0. This is just a
thermal (Maxwellian) distribution multiplied by v_{x}, the velocity in the
direction normal to the surface and is required for describing sources that
simulate a recycling process.
A brief digression will clarify the relationship between this source
distribution and the thermal distribution that characterizes the
particles in the volume of the fluid. Consider the number of particles
of arbitrary volume distribution f
that recycle at, say, the x = 0 surface in a time interval ∆t
(a source distribution specifies a density in phase space per unit time).
For particles in the velocity interval v_{x} → v_{x} + dv_{x},
this is f dv_{x} A ∆x, where A is the area of the recycling surface
and ∆x = v_{x} ∆t is the distance over which particles with
velocity v_{x} can reach x = 0 in ∆t. So, the phase space
density of the recycling particles per unit time (dividing this number
by ∆t and dv_{x}) is ∝ v_{x} f. Another way of saying this
same thing is that the source needs to emphasize the faster particles
in the distribution since more of them will reach the recycling surface in
a given time interval. One upshot is that the average energy of the
source particles is 〈E 〉_{source} = 2 T while the
average energy in the volume (for v_{y0} = 0) is the familiar 3/2 T.
For clarity, we note that the volume distribution function in the
freemolecular limit is
 

f_{+}( 
→
v

) + f_{−}( 
→
v

), 
  (3.2) 
 


n
2

2  ⎛ ⎝

m
2 πT
 ⎞ ⎠

3/2

exp  ⎧ ⎨
⎩

− 
m
2T

[v_{x}^{2} + (v_{y} − v_{y0})^{2} + v_{z}^{2}]  ⎫ ⎬
⎭

, v_{x} > 0 
  (3.3) 
 


n
2

2  ⎛ ⎝

m
2 πT
 ⎞ ⎠

3/2

exp  ⎡ ⎣

− 
m
2T

(v_{x}^{2} + v_{y}^{2} + v_{z}^{2})  ⎤ ⎦

, v_{x} < 0 
  (3.4) 

With this distribution, one can show that the fluid pressure is
P = 
2
3

〈n E 〉 = nT + 
1
6

mv_{y0}^{2}. 
 (3.5) 
The fluid velocity is
 

〈 
→
v

〉_{+} + 〈 
→
v

〉_{−} 
 
 

  (3.6) 

so that  〈→v 〉 = v_{y0} / 2.
To compute the fluid temperature used in the BGK distributions, we need
to subtract the drift energy,
The primary physical quantity appearing in the analytic
models[45,46,47] of Couette
flow is the xy component of the stress tensor,
Π_{xy} ≡  ⌠ ⌡

d^{3} v m v_{x} v_{y} f( 
→
v

)− n m U_{x} U_{y}, 
 (3.9) 
where →U is the first (velocity) moment of f (the flow
velocity). Note that while we should have U_{y} ≠ 0, U_{x} should
be identically 0 since there is no net flow in the x direction.
In the molecularflow regime, the distribution function is not changed by
collisions as the particles move across the box. Furthermore, particles
coming from the x = d side of the box make no net contribution to
Π_{xy} (since the integrand is an odd function of v_{y}).
Inserting Eq. (3.3), we find
Π_{xy}^{fm} = 
n m
4

 ⎛ ⎝

8 T
πm
 ⎞ ⎠

1/2

v_{y0}. 
 (3.10) 
As collisions become more important,
the plate velocity information is
dissipated into random motions of the fluid to a greater and greater
extent, and Π_{xy} falls below Π_{xy}^{fm}.
The variational theory of Cercignani[46,47] yields a
relatively compact formula for Π_{xy} that matches the
numerical results of Willis[45] to within 0.5 %.
With δ = 1 / K_{n},

Π_{xy}
Π_{xy}^{fm}

= 
a + b δ+ c δ^{2}

, 
 (3.11) 
where for the case of a BGK collision operator,
The xy component of the stress tensor is computed along with other
test particle data during tracking. However, final processing of that
information into a usable form is only done when the COUETTE
macro is enabled in flighttest.web.
The Couette flow problem can be simulated in DEGAS 2 with only a few
code modifications to set up the required boundary conditions (at the moment,
the required particle distributions cannot be specified completely
by input or data files). Because the problem is relatively simple, it
represents a good starting point for learning how to use
neutralneutral collisions in DEGAS 2.
Most of these changes are invoked with a COUETTE
FWEB macro inserted into the files
To run this example, follow these steps:
 Set source code switches. The source code file boxgen.web
has within it a few different settings available which can be selected
by changing the value of the flag BOXRUN. With this file open in
an editor, change the lines
@m DIM 2
@m BOXRUN 0
@m SOLN 1
so that they now read
@m DIM 1
@m BOXRUN 41
@m SOLN 0
Save the file and exit the editor. The DIM flag controls the
dimensionality of the geometry. By setting it to 1, we get a purely
onedimensional geometry, with periodic boundary conditions in the y
and z directions. The SOLN flag controls the analytic
solution used with the "Analytic_fluid_bench" example
(Sec. 3.7.1); since that's not needed here, we can
turn it off.
The files plate.web, sources.web, and flighttest.web each
have the line
@m COUETTE 0
near the top. Edit these files, changing the "0" to a "1" in each case
(be sure to revert to the default setting of "0" before running other
problems).
 Switch to the example directory. If the main DEGAS 2
directory sits within your home directory:
cd ~/degas2/examples/Couette_flow
 Copy the degas2.in file into your working
directory. E.g., if you are working on a Sun,
cp degas2.in ../../SUN/degas2.in
cd ../../SUN
(overwriting any degas2.in file you may have had there already!).
 Prepare reference data. This example uses some nonstandard
data files so, unlike the previous example, datasetup
must be run. Note that the degas2.in file points back to
the examples directory for the input files needed here and in
subsequent steps.
gmake datasetup
./datasetup
 Prepare problem data.
gmake problemsetup
./problemsetup
If you look at the problem input file, you will see that there is only
one reaction, that for the neutralneutral collisions. Note also that
the test species D2 also appears in the list of background species
(the usual background species e and D+ appear here just because
boxgen expects them to be there; they are not used in this simulation).
 Generate geometry and background files. The boxgen
program takes care of both of these (be sure that the BOXRUN macro
has been set to 41):
gmake boxgen
./boxgen
 Set the number of flights. The background file bk_boxgen.nc
contains the specification of the number of flights to be run. Open this
file in an editor (Emacs will be the most convenient) and find the line
source_num_flights = 100 ;
The results in this directory were obtained with 10 BGK iterations of
4000000 flights each. Such a run could easily require
more than a day if run on a single processor (this run was executed
in 15 minutes using 16 Dualcore AMD Opteron 1 GHz processors).
You can probably run a smaller
number of flights (say, 100000) and still get reasonable results.
Save the file and exit the editor.
 Set up the tallies.
gmake tallysetup
./tallysetup
 Run the code.
gmake flighttest
./flighttest
The number of BGK iterations is controlled by the value of the
bgk_max macro in flighttest.web. The value used here is 10. If you just
want to get an idea of how the calculation proceeds, you can choose a
smaller number. Note that the macro parameters controlling the convergence
tests, bgk_cvg_dens and bgk_cvg_pres,
are set to ridiculously small
values (10^{−8}) to ensure that all 10 iterations are completed.
 Check the text output. At each iteration, the code writes out the
test species density and other data in a file with a name of the form
densityxx.out where xx is replaced by the iteration number.
The contents of these files are analogous to those of the other examples
(the density and pressure for each zone are in the third and fifth
columns, respectively;
see the description in Sec. 3.7.1), except for the
last two columns. These are normally not used; here, they contain the
values and standard deviations of the xy stress tensor, Π_{xy}.
The comparison of those values with theory will be discussed below.
At the end of
each iteration, the background netCDF file (bk_bgktest.nc) is
also updated (overwritten) with
the density, velocity, and temperature of the neutral background species.
A third file, cvg_global.txt, will provide information on the
global progress of the BGK iterations towards convergence. The four columns are:
 Iteration number,
 Test species number,
 Global fractional change in density (dimensionless); this is compared
with bgk_cvg_dens to decide whether or not to stop the iterations.
 Global fractional change in pressure. (dimensionless); this is compared
with bgk_cvg_pres to decide whether or not to stop the iterations.
Note that on most machines this file
does not appear until the unit is closed at the end of the run.
These files can be directly
compared with the corresponding files in the examples directory.
If you have repeated the run at the abovedescribed length, you should
be able to match the text files exactly. The level of agreement in the
netCDF file should be several digits (perhaps as many as 10).
If a satisfactory match is not obtained, something is
amiss. Please contact the code authors (see Sec. 5.1) if you feel
that a problem exists with the code as you downloaded it.
 Check the binary output. You can use the matchout
utility (see Sec. 3.4.4) to do this. Again, you should expect
the largest differences to be ∼ 10^{−8} or better.
 Compare with the analytic solution. The rightmost set of data
in the densityxx.out file are the stress tensor values, Π_{xy}.
Ideally, these are constant across the problem space. Some deviation
exists here due to Monte Carlo noise. You
can compute the free molecular value from Eq. (3.10). With both of
these numbers in hand, you can compare with Eq. (3.12). To do that, you
need to know δ. The explicit formula used is
where n = 10^{20} m^{−3} (enforced in flighttest.web by code
enabled with the COUETTE macro), 〈σv 〉 = 1.1×10^{−16} m^{3}/s (from the d2d2_bgktest.nc file), d = 0.1 m
(set in boxgen.web), and m = m_{D2} = 6.689 ×10^{−27} kg
(also from boxgen.web). One would expect the temperature T of the
fluid to be equal to that of the wall by virtue of the initial conditions.
However, in cases where the energy associated with the plate velocity
v_{y0} is significant compared with T_{wall}, the actual temperature
of the fluid is slightly higher than that of the wall as is indicated
by Eq. (3.8).
In fact,
this example with v_{y0} = 1000 m/s
shows better agreement with the theoretical result Eq. (3.12)
if T is set to the simulated fluid temperature (e.g., from the
the bk_bgktest.nc file).
With T = T_{wall} = 0.0258 eV = 4.1336 ×10^{−21} J, we find
δ_{wall} = 0.9894. Plugging this into Eq. (3.12),
 ⎛ ⎝

Π_{xy}
Π_{xy}^{fm}
 ⎞ ⎠

c

(δ_{wall}) = 0.6047. 
 (3.14) 
The freemolecular stress, Eq. (3.10), is
Π_{xy}^{fm}(T_{wall}) = 0.20977 Pa.
If we average the 10 values for Π_{xy} in the density10.out file, we
get Π_{xy, sim} = 0.13045 Pa and
 ⎛ ⎝

Π_{xy}
Π_{xy}^{fm}
 ⎞ ⎠

sim

= 0.6219. 
 (3.15) 
This would qualify as pretty good agreement. However, if we average the fluid
temperature values in the bk_bgktest.nc file, we get T_{sim} = 4.5866 ×10^{−21} J. Then, δ_{sim} = 0.9394, and
 ⎛ ⎝

Π_{xy}
Π_{xy}^{fm}
 ⎞ ⎠

c

(δ_{sim}) = 0.6157, 
 (3.16) 
about a factor of two better agreement. Now, one might argue that we should
also use
T_{sim} in computing Π_{xy}^{fm}. However, the
freemolecular stress, Eq. (3.10) was computed directly
from the general distribution function for arbitrary v_{y0} so that
T = T_{wall} in that formula, without any ambiguity.
Qualitatively speaking, for small
δ the character of the source distribution dominates. As δ
approaches or exceeds 1, the temperature of the fluid in the volume,
T_{sim}
computed from Eq. (3.8), better describes the distribution.
The better agreement is shown over a range of inverse Knudsen
numbers and two values of v_{y0} in Fig. 3.1.
Figure
Figure 3.1: Comparison of 8 DEGAS 2 simulations (symbols)
of the Couette flow problem
with the analytic expression for the normalized shear stress
Eq. (
3.12) obtained by Cercignani[
46,
47]. Two
different plate velocities V = v
_{y0} were used in the simulations.
For the higher velocity, better agreement with the analytic expression
is obtained if the inverse Knudsen parameter is computed from the
actual fluid temperature from the simulation, T
_{sim}, rather
than the wall temperature, T
_{wall}.
CMod 1D: Conservation Checks
This example demonstrates a run with complete hydrogen physics in a simple
geometry. In particular, this run includes both neutralion and
neutralneutral
elastic scattering processes. The complexity of the physics also
makes for a good demonstration of DEGAS 2's conservation checks.
Experiments by Pitcher et al. [54,55] on the Alcator CMod tokamak
indicated that the flow of neutral gas through material structures around
the plasma was being limited by its diffusion rate through the divertor plasma.
Pitcher [55] developed a semianalytic model that reproduced this
behavior. An attempt was made to precisely reproduce the
results of that model with DEGAS 2. The overall effort consisted of several
simulations in which the model assumptions were gradually relaxed towards those
consistent with normal DEGAS 2 operation. This example is from a simulation in
having atomic and surface physics assumptions more like a normal DEGAS 2 run.
The remaining difference from a "full physics" DEGAS 2 run (apart from the
greatly simplified geometry) is that the neutral source is a gas puff rather
than a recycling source. The full physics and geometry simulations are described
in Ref. [56].
To run this example, follow these steps:
 Set source code switch. The source code file boxgen.web
has within it a few different settings available which can be selected
by changing the value of the flag BOXRUN. With this file open in
an editor, change the lines
@m BOXRUN 0
@m SOLN 1
so that they now read
@m BOXRUN 40
@m SOLN 0
Save the file and exit the editor. The second parameter disables the
analytic solutions that are used in the original boxgen
example (Sec. 3.7.1). Note that the
DIM parameter should remain at its default value of
2.
 Switch to the example directory. If the main DEGAS 2
directory sits within your home directory:
cd ~/degas2/examples/CMod1D
 Copy the degas2.in file into your working
directory. E.g., if you are working on a Sun,
cp degas2.in ../../SUN/degas2.in
cd ../../SUN
(overwriting any degas2.in file you may have had there already!).
 Prepare reference data. This example uses some nonstandard
data files so, unlike the previous example, datasetup
must be run. Note that the degas2.in file points back to
the examples directory for the input files needed here and in
subsequent steps.
gmake datasetup
./datasetup
 Prepare problem data.
gmake problemsetup
./problemsetup
 Generate geometry and background files. The boxgen
program takes care of both of these. Be sure the BOXRUN macro was
set to 40:
gmake boxgen
./boxgen
 Set the number of flights. The background file bk_oned.nc
contains the specification of the number of flights to be run. Open this
file in an editor (Emacs will be the most convenient) and find the line
source_num_flights = 40000 ;
Save the file and exit the editor.
 Set up the tallies.
gmake tallysetup
./tallysetup
 Run the code.
gmake flighttest
./flighttest
The code will execute 5 iterations automatically, each with 40,000 flights.
The number of BGK iterations is controlled
by the parameter bgk_max in flighttest.web; the default setting is 5. The
BGK iterations can alternatively be controlled by the convergence parameters
bgk_cvg_dens and bgk_cvg_pres.
The corresponding tests monitor global measures
of the changes in the density and pressure of the BGK species. Note that
these can be satisfied only if there are sufficiently many flights
to achieve a correspondingly high precision in the local density and
pressure values, no matter how many BGK iterations are run. In this example,
the parameters chosen are such that
bgk_max will be governing the number of iterations.
 Check the text output. At each iteration, the code writes out the
test species density and other data in a file with a name of the form
densityxx.out where xx
is replaced by the iteration number.
The contents of these files are analogous to those of the other examples
(the density and pressure for each zone are in the third and fifth
columns, respectively;
see the description in Sec. 3.7.1).
At the end of
each iteration, the background netCDF file (bk_oned.nc) is
also updated with
the density, velocity, and temperature of the neutral background species.
This file can be used to restart the BGK iterations, e.g., with an increased
number of flights.
A third file, cvg_global.txt, providea information on the
global progress of the BGK iterations towards convergence. The four columns are:
 Iteration number,
 Test species number,
 Global fractional change in density (dimensionless); this is compared
with bgk_cvg_dens to decide whether or not to stop the iterations.
 Global fractional change in pressure. (dimensionless); this is compared
with bgk_cvg_pres to decide whether or not to stop the iterations.
Note that on most
machines this file does not appear until the unit is closed at the end of
the run. The output netCDF file from this run is included as well.
These files can be directly
compared with the corresponding files in the examples directory.
The numbers in the text files have only five
digits of precision and should be matched exactly. The level of agreement
in the netCDF file should be several digits (perhaps as many as 10).
If a satisfactory match is not obtained, something is
amiss. Please contact the code authors (see Sec. 5.1) if you feel
that a problem exists with the code as you downloaded it.
 Check the binary output. You can use the matchout
utility (see Sec. 3.4.4) to do this. Again, you should expect
the largest differences to be ∼ 10^{−8} or better.
 Examine test particle balances To simplify examination
of the conservation
checks, an input script for the outputbrowser code is included,
bgkbalances. The results are in bgkbalances.out. The checks can be
performed for the density, momentum, and energy of each of the test
species. Since D2+ does not move, the checks on it are satisfied
trivially and are not considered further.^{2}
For each quantity, the
procedure consists of totaling the negative of the quantity lost to the
walls (the "current into wall" tally), the quantity coming in from the
walls (the "current out of wall" tally), and the ßource rate" tally. The
last item is further
broken down by reaction (or source) to provide an indication of the
relative importance of each process. For this example, we consider only
the Xcomponent of the momentum (dimension #1) since the problem does not vary
in the other two directions.
For this case, the following totals are obtained:
 particles

D_{2} ¯ 3 ×10^{15}
D −2 ×10^{16}
D_{2} 8 ×10^{15}
 momentum 1

D_{2} ¯ < 3 ×10^{15}
D < 1 ×10^{−6}
D_{2} 2 ×10^{−8}
 energy

D_{2} ¯ 3 ×10^{15}
D −5 ×10^{−2}
D_{2} 7 ×10^{−3}
In addition, a global check involving the background particle number can
be performed:
 

Total source of deuterium atoms (the puff source), 
 
 

 
 

 
 

Total number of deuterium atoms lost to the walls, 
 
 

−(3.48326 ×10^{22} + 2. ×4.11582 ×10^{21}) 
 
 

+ (3.23371 ×10^{22} + 2. ×4.76467 ×10^{21}) 
 
 

 
 

Total number of deuterium atoms lost to the plasma (D+ source rate), 
 
 

 
 

 
 

 
  

There will usually be a small remainder in these balances;
its magnitude in this case roughly coincides with the accuracy provided by
outputbrowser. This remainder results from the enforcement of a minimum
weight via Russian roulette in subroutine follow. By default, that
minimum weight is 10^{−3} of the initial weight.
Consequently, particles,
momentum, and energy are conserved
only on the average (i.e., after many flights).
You can demonstrate this effect by reducing the minimum weight (WMIN in
flight.web) to, say, 10^{−7}.
Doing so will make the run substantially longer.
To see the difference in the conservation checks, you may
need to examine the output netCDF file directly (hint: run outputbrowser
first to get "rough" values for these quantities and then use the
search facility in Emacs to locate the corresponding data in the output
netCDF file).
 Transfers between BGK species. In looking at the energy and
momentum
transfers between BGK species, you should notice that no sources
appear for the background partners. The reason for this is that the transfers
to and from these species do not have a clear physical meaning; all
relevant information can be obtained from the corresponding test species.
However, the transfers between BGK test species, say between D and
D2
here, are significant. The process d2d_neut (see the "momentum
source vector by reaction" and "energy source by reaction" sections
of bgkbalances.out under reaction 11) transfers
energy and momentum between the background D2 and the test D
(problem species, "problem_sp", number 6); likewise,
dd2_neut (reaction number 12) connects the
background D and test D2 (problem_sp 7).
One would expect that in a
converged state these two transfer rates would be equal (this actually
goes into the derivation of the algorithm). In this example, they are
(to within the error bars),
mom. 1 to D: ¯ xx −9.39835 ×10^{−3} xx ¯
mom. 1 to D2: ¯ xx +9.19450 ×10^{−3}
mom. 1 to D: −9.88463 ×10^{−3}
mom. 1 to D2: +9.59007 ×10^{−3}
energy to D: −4.46727 ×10^{2}
energy to D2: +4.35459 ×10^{2}
(The molecules come off the wall at low energy and are being heated by the
warmer atoms that arise from elastic scattering with the plasma ions and
from dissociation.)
For this to occur, the reaction rate for these two processes must
be numerically identical. The original implementation of the BGK
algorithm in EIRENE had the same temperature dependent expression for
the two rates. However, because the temperatures of D2 and
D species are different
in general, the resulting rates had different numerical values. The
initial implementation of the BGK algorithm in DEGAS 2
repeated this same mistake. Due to a dearth of applications for
these reactions, this problem has not yet been remedied.
To obtain the nearly equal transfer rates seen above,
the degas2.in file in this example actually points to a
dedicated reactions.input file in the same directory
that uses modified ("hacked") data files for these two elastic
scattering processes. To see the effect of unequal reaction
rates, modify the degas2.in file to instead point
to the reactions.input in the data directory.
Rerun both datasetup and problemsetup, then
run the code again. Note that you should also rerun boxgen
to start the code off in the same initial state. The
result after five iterations can
be compared with
bgkbalances_default.out that comes with the examples:
mom. 1 to D: ¯ xx −9.39835 ×10^{−3} xx ¯
mom. 1 to D2: ¯ xx +9.19450 ×10^{−3}
mom. 1 to D: −5.06271 ×10^{−3}
mom. 1 to D2: +7.51898 ×10^{−3}
energy to D: −2.16540 ×10^{2}
energy to D2: +3.55497 ×10^{2}
 Other results. For reference purposes,
this directory contains output files from a
single iteration (on a Linux machine running
Portland Group Fortran with DEBUG=yes) with 1000 flights. These
output files contain the string 1K.
The primary purpose of these files
is to serve as a known, reproducible point of contact for testing
subsequent code modifications. It also demonstrates that these
conservation checks hold even if the code is run with relatively few
flights. To duplicate these files, you will need to stop the BGK iteration
process after the first (zeroth, actually) iteration. The easiest way
to do this is to change the number of iterations, near the
top of flighttest.web, from:
@m bgk_max 5
to
@m bgk_max 0
You may notice that the tallies used for the conservation checks
predominantly utilize collision estimators while the standard (zoneresolved)
tallies are based on track length estimators. The reason is that the
former will explicitly demonstrate conservation of mass, momentum, and
energy since they are computed using the instantaneous test particle
attributes. The track length estimators will generally provide more
accurate values for the code output, but will exhibit conservation only
in a statistical sense. You can demonstrate this by replacing the
estimator lists in each of the "total" tallies with the lists from the
corresponding zoneresolved tallies.
3.7.3 definegeometry2d and defineback Examples
The definegeometry2d (for setting up geometries, obviously)
and defineback (for specifying the corresponding sources
and plasma background) are sufficiently flexible and general
that specific examples are useful in illustrating the ways that
they can be used.
Gas Puff Imaging Example
This example illustrates the use of the DG and Carre codes
with definegeometry2d, and also demonstrates the 3D
capabilities of definegeometry2d.
It is based on the DEGAS 2 simulation of the NSTX
Gas Puff Imaging (GPI) GPI shot 112811, documented in
[57].
To run the 2D version of this example, follow
these steps:
 Switch to the example directory. If the main
DEGAS 2 directory sits
within your home directory:
cd ~/degas2/examples/He_GPI
 Copy the degas2.in file into your working directory.
E.g., if you are
working on a 64bit Linux machine,
cp degas2.in ../../LINUX64
cd ../../LINUX64
 Prepare reference data.
This example uses the standard reference data
files in the data directory. Hence, this step can be skipped if those files
have not been altered since you downloaded the code.
gmake datasetup
./datasetup
 Prepare problem data.
gmake problemsetup
./problemsetup
 Generate geometry file. Here, we outline
briefly the preliminary steps
leading up to the creation
of the objects referred to in the definegeometry2d input
file. Users interested in learning
more about that process should consult the DG User's Manual (in
the dg_carre/DG/DOC directory;
especially the DEGAS 2 specific sections under
"Customizing DG").
The starting points for the DG code are a specification of the
vacuum vessel shape [(R,Z) coordinates; nstx_walls in this
case] and an equilibrium (values of
the poloidal flux on a rectangular R,Z grid,
g112811.00250.equ, already converted from the EFIT format to DG's
preferred file format). Once these are loaded into DG,
the user proceeds to graphically generate the input for
the Carre code. Carre is run separately outside of DG to
create a nearly orthogonal 2D mesh with one
coordinate that follows the flux surfaces. Although DEGAS 2 doesn't
really need such a mesh, having one facilitates specifying plasma parameters
(typically a function of a flux variable). That mesh is then loaded into DG
and the user proceeds to connect it to the vacuum vessel
points forming the "DG polygons"
(see the entry in the definegeometry2d
manual).
This
process results in three files: the DG file (nstx_25.dg),
the mesh (nstx.carre.021) and
DG's objects (nstx_25.dgo). Only the latter two
are read in by definegeometry2d.
There are two definegeometry2d input
files here, one for the 2D (axisymmetric:
nstx_25_dg2d.in) case and one for 3D (nstx_25_dg3d.in).
The results
contained in the examples/He_GPI directory
are for the much quicker running 2D
case. Note that
many of the 3D lines appear in the 2D file, but are commented
out. Both runs utilize just the portion of the Carre mesh directly
in front of the gas manifold ("puffer" in the definegeometry2d
input file); this is done using the additional arguments
to the sonnet_mesh keyword. One would establish the values of these
indices during the creation of the DG polygons. The 3D run then only
simulates a portion of the torus (arguments to the bounds and
y_values
keywords). These steps keep the problem size and run time down.
The emulation of the GPI camera view is done with the code in
gpicamera.web.
Prior to compiling definegeometry2d, the user should copy this file
to usr2ddetector.web (see the brief explanation at the top of the
default version of the file def2ddetector.web). This is actually
just a placeholder 5 ×5
pixel version of the camera. The full resolution
image is created during postprocessing with the postdetector routine.
cd ../src
cp gpicamera.web usr2ddetector.web
cd ../LINUX64
gmake definegeometry2d
./definegeometry2d ../examples/He_GPI/nstx_25_dg2d.in
 Generate background file. The detailed
specification of the plasma parameters is actually done with
the code in nstxgpi2.web. Like gpicamera.web,
this is a userdefined
subroutine (get_n_t) that needs to be compiled into the driver code.
This routine will read the plasma density and temperature profiles
obtained with the NSTX
Thomson scattering system, ts_112811_base.TXT and map them on to
the entire geometry assuming that the plasma density
and temperature are constant on flux surfaces. On the Carregenerated
mesh, this is straightforward. Outside of it, where there are only
triangles, an approximate procedure is used.
The plasma_file keyword in the defineback input file has two
parameters used by the nstxgpi2.web get_n_t routine.
The second is the name of the Thomson scattering data file mentioned above.
The first argument
provides the subroutine with the location
of the separatrix in the Carre mesh
subset used by definegeometry2d.
Specifically, the separatrix is assumed to be between the cells iz_sep
and iz_sep+1.
If one knows the number of "cells" used in DG's
"Create Surface(s)..." dialog for the core (n_cells)
and the minimum iz
used in the mesh subset (specified via the sonnet_mesh keyword in
definegeometry2d), iz_min,
the value of this is: iz_sep = n_cells  iz_min + 2.
The second piece of input to defineback is the specification of the
neutral sources. In the GPI simulations, this is a gas puff from
the polygon representing the gas manifold, a vertical line
in the 2D case.
The physical locations of the 2D sources are
provided by the "stratum" and "segment" indices
in the file source_nstx_26. The "stratum" number
corresponds to the "stratum" keyword used in definegeometry2d.
In general, the simplest way to determine the segment numbers
is to use the "*" option (see the
defineback documentation) to
select all of the segments associated with a stratum and then
manually identify those desired for the problem at hand.
In this
problem, a constant flux source is used since the sizes of
the segments vary. The overall magnitude of the source is arbitrary
in this problem, but physically realistic. The source
file for the 3D case (source_nstx_26_3d) is similar, but also provides
the toroidal location of the source segments (Segment_iy). These
were chosen so that the source points fall along a nearly straight
line at the same angle with respect to horizontal as the actual gas
manifold.
cd ../src
cp nstxgpi2.web usr2dplasma.web
cd ../LINUX64
gmake defineback
./defineback ../examples/He_GPI/nstx_30_db.in
 Set the number of flights.
The background file bk_nstx_30.nc contains the specification of the
number of flights to be run. Open this file in an editor (Emacs
will be the most convenient) and find the line:
source_num_flights = 100 ;
The results in this directory were obtained with 100000 flights.
However, the 3D production run used 2000000.
 Set up the tallies.
gmake tallysetup
./tallysetup
 Run the code.
gmake flighttest
./flighttest
 Check the text output.
The file density.out can be directly compared with the one in the
examples/He_GPI directory. The numbers in here have only five
digits of precision. If an exact match is not obtained, something
is amiss. These results were obtained with Triangle V. 1.6.
 Check the binary output.
You can use the matchout utility to do this. The largest differences
should be ∼ 10^{−10} or smaller.
 Generate graphical output.
For the 2D simulation,
cp ../examples/He_GPI/geometry.inp .
gmake geomtesta
./geomtesta
The currently recommended approach is to compile the code with the
GRAPH_FILE=SILO or SILO_HDF5
option and use VisIt (built with Silo
support) to view the Silo format file generated by geomtesta.
The primary output variables in this problem are the emission rate
of the helium 5877 Å (he_5877) viewed by the GPI camera
and the neutral helium density, spHeden.
For 3D simulations, two example input files for geomtesta are
provided. Keep in mind that these need to be copied to the
run directory and renamed geometry.inp in order for geomtesta to
find them.
 Generate camera image and associated data.
The camera simulated by the main code uses only a 5 ×5 mesh in
this problem so as to not lengthen the run. The
generation of actual images suitable for comparison with
experiment is deferred to postprocessing with the postdetector
code. In principle, postdetector can be parallelized (although recent
changes to the code have yet to be propagated into the MPIspecific
sections), further mitigating the time expense for producing
the simulated camera image.
While postdetector has been generalized with regard to species
(e.g., He or D_{2}), there are still application specific settings
in the code. These are primarily the camera resolution and
specification of the target plane; note that these same
quantities appear in the preamble of gpicamera.web.
The orientation of the simulated image can also vary with
application; see subroutine wrdatap in postdetector.web.
Open postdetector.web in an editor and find the section of
code:
@#if 0
@m nx 64
@m ny 64
@#else
@m nx 5
@m ny 5
@#endif
Edit the first line to read:
@#if 1
This sets the camera resolution to 64 ×64. The
5 ×5 resolution is used as
the default to facilitate testing.
A few lines below this, check that the application is set to
"NSTX":
@m APP NSTX
As was the case in compiling definegeometry2d, the file
usr2ddetector.web
should be a copy of gpicamera.web. Note that the camera resolution
parameters in usr2ddetector.web
are overridden by those in postdetector.web
and are not needed here.
gmake postdetector
./postdetector
The camera image and other data are written to HDF format files (this
code has not yet been adapted for use with the Silo format). Here
are the contents of the most useful files:
 emtt
 Total emission rate of the 5877 Å line
[in W/(m^{2} ster)]. This is the simulated camera image.
 tane
 Electron density (in m^{−3}) at the
target plane (intersection of the camera
view and the idealized gas puff sheet).
 tate
 Electron temperature (in eV) at the target plane.
 tarr
 Major radius (in m) at the target plane.
 tarz
 Vertical coordinate (in m) at the target plane.
 tns1
 Neutral helium density (in m^{−3}) at the target plane.
Note, however, that the identification of the "target plane zones" is
an automated optimization process that assumes a relatively fine toroidal
discretization. As such, its results are rather difficult to
interpret in the 2D case.
The quantities in the other files require much more space to explain.
The
interested reader can learn more in [58]
 Threedimensional simulation.
The process for producing the analogous threedimensional simulation is
the same as above. The only modifications are:
 The input file for definegeometry2d is
nstx_25_dg3d.in. Note
that definegeometry2d will require much more time (on the order of an
hour) to run.
 The pointer to the source file in the defineback input
file nstx_30_db.in must be changed to point to
source_nstx_26_3d.
 Many more flights are required for the main code,
e.g., 2000000, and
the run time is correspondingly longer.
 One of the two 3D input files should be used for
geomtesta.
The resulting camera image corresponds to Fig. 4(a) in
[57]. However,
the latter also incorporated the effects of camera vignetting, which
are not accounted for by the postdetector routine.
NCSX Example
This example illustrates the use of definegeometry2d with
input files generated with a relatively simple
standalone code. The same concepts could be used
to manually create definegeometry2d input files.
This particular simulation provided
neutral penetraton estimates used
in the design of plasma facing components for the
National Compact Stellarator Experiment [59].
Since the bulk of the recycled neutrals leave
the surface with a cosine distribution (i.e., peaked about the normal
to the material surface), the
greatest penetration will be achieved in the same poloidal plane as
the recycling surface.
Consequently, toroidally axisymmetric
(i.e., a single plasma cross section at a
particular toroidal angle) neutral transport calculations
can place an upper bound
on the neutral penetration distance.
In scoping problems such
as this, satisfactory results can be obtained by assuming that
the plasma parameters are constant on flux surfaces (of course,
if one does know the plasma variation along a flux surface, this
assumption can be relaxed). The DEGAS 2
geometry can then be constructed out of flux surface shapes, e.g.,
from an equilibrium calculation of some sort. The basic approach is
to define a "wall" corresponding to each flux surface and then
connect adjacent walls to form the "polygons" required
by definegeometry2d (see
its documentation for
detailed descriptions of these objects). More specfically, each pair of
adjacent flux surfaces is used to construct two polygons (e.g.,
a lower and upper half) so that
the interior of both is topologically well defined.
This objective can be in principle be achieved by manual manipulation
of the coordinate data. However, the process is sufficiently
straightforward, not to mention tedious, that a code can be constructed that
generates the wallfile needed by
definegeometry2d, as well as the
polygons created out of those walls that form the bulk of
the definegeometry2d input file.
The starting point in this case was a 3D equilibrium
file generated by the VMEC code, containing Fourier representations
of the NCSX flux surfaces. An example code
distributed with the VMEC file was first modified to evaluate these
surfaces at a set of input toroidal angles (actually just two, with
one represented here). A uniform poloidal angle mesh (an
arbitrary choice) was then defined so that a given flux
surface could be represented by a set of discrete points.
There were many more flux surfaces in the VMEC equilibrium than were
needed for this application. The set of surfaces was
reduced (again arbitrarily) so that the maximum distance
between two adjacent surfaces was 5 cm. A "wall" was
created out of each of these surfaces, with the
constituent points being the R, Z coordinates of the
surface at each of the points along the poloidal angle mesh.
Because the relative location of these walls was known,
the generation of the corresponding polygons was a
simple matter.
The code also defined walls and polygons filling
the gap between the VMEC surfaces and the vacuum vessel
and connecting the vacuum vessel to the universal cell.
Details such as these depend on the particular
application and are, thus,
not worth describing here.
Another crucial task that must be begun at this stage is
establishing the mechanism that will allow defineback
to map the plasma parameters onto these
polygons. The first, and generally essential, step to doing
this is to label each polygon with a "stratum" number.
Note that DEGAS 2's internal geometry specification,
contained in the geometry
netCDF file, has no record of the polygons
used by definegeometry2d; it only knows about zones.
The stratum labels are used in conjunction with the
geometry's "sectors", but they are of little use within
the plasma volume. For this reason, definegeometry2d also
generates its polygon_nc_file based on the
2D geometry class (see the
geometry2d class documentation).
Via the g2_polygon_stratum and g2_polygon_zone
arrays in this class,
the user can connect these
zones back to the polygons
specified on input to definegeometry2d.
In this example, this objective is accomplished
by the subroutine in the
ncsxplasma.web file compiled into
defineback.
Be aware that the polygons in the polygon_nc_file
file are derived from, but not identical to, those
input to definegeometry2d. Specifically,
polygons broken up with definegeometry2d's
breakup_polygon will remain as individual
polygons in the polygon_nc_file. But,
those processed with trianglulate_polygon or
triangulate_to_zones are represented by the
resulting triangles within the polygon_nc_file and, thus,
correspond to multiple polygons there. The stratum
labels are propagated to the triangle polygons in the
process, however, permitting them and the zone numbers
to be related to the original polygons.
In this example, the plasma profiles are specified as
a function of the square root of the normalized toroidal flux.
Consequently, the geometry setup code also writes out a file,
radiiz1li383,
containing the toroidal flux corresponding to each stratum.
Because each pair of adjacent flux surfaces is used to construct
two polygons, these two polygons are assigned the same stratum
label in the definegeometry2d input file.
The names of the files containing the plasma profiles, the
effective stratum radii, and the polygon_nc_file
are input to defineback and passed to the routine
in ncsxplasma.web.
One will also need to be able to tell
defineback where to place the neutral sources. This usually
amounts to manually identifying the "strata" corresponding to
those sources within the definegeometry2d input file, and
then following a simple procedure for deducing the specific "segments"
of those strata that are to be used.
To run this example, follow these steps:
 Switch to the example directory. If the main
DEGAS 2 directory sits
within your home directory:
cd ~/degas2/examples/NCSX
 Copy the degas2.in file into your working directory.
E.g., if you are
working on a 64bit Linux machine,
cp degas2.in ../../LINUX64
cd ../../LINUX64
 Prepare reference data.
This example uses the standard reference data files.
gmake datasetup
./datasetup
 Prepare problem data.
gmake problemsetup
./problemsetup
 Generate geometry file.
gmake definegeometry2d
./definegeometry2d ../examples/NCSX/dg2dinz1li383
 Generate background file. As was described above,
the detailed specification of the plasma parameters is done
by the get_n_t routine in ncsxplasma.web. It requires
as input the three files described previously (the effective radii for
each stratum, the plasma profiles, and the polygon_nc_file),
and the plasma density and temperature scrapeoff lengths (the plasma
parameters are only given on closed flux surfaces).
The remainder of the defineback input file
contains the specification of the
neutral sources. In this example, this is a recycling
source at a limiter. The corresponding "stratum" was
established as part of the process of creating the input
files for definegeometry2d.
A general, and relatively quick, technique for
determining the source segment numbers
is to use the "*" option (see the
defineback documentation) to
select all of the segments associated with a stratum and then
manually identify those desired for the problem at hand.
The source strength is arbitrary in this problem, although
a physically realistic value was selected.
cd ../src
cp ncsxplasma.web usr2dplasma.web
cd ../LINUX64
gmake defineback
./defineback ../examples/NCSX/z1li383_db.in
 Set the number of flights.
The background file bk_ncsx.nc contains the specification of the
number of flights to be run. Open this file in an editor (Emacs
will be the most convenient) and find the line:
source_num_flights = 100 ;
The results in this directory were obtained with 10000 flights.
 Set up the tallies.
gmake tallysetup
./tallysetup
 Run the code.
gmake flighttest
./flighttest
 Check the text output.
The file density.out can be directly compared with the one in the
examples/NCSX directory. The numbers in here have only five
digits of precision. If an exact match is not obtained, something
is amiss. These results were obtained with Triangle V. 1.6.
 Check the binary output.
You can use the matchout utility to do this. The largest differences
should be ∼ 10^{−10} or smaller.
 Generate graphical output.
For the 2D simulation,
cp ../examples/NCSX/geometry.inp .
gmake geomtesta
./geomtesta
3.8 Run Control Parameters
Several run control parameters are contained in the
sources class. Ideally, these would be
set by the user in a dialog box just prior to the run. For now,
they must be changed manually by editing the background netCDF file.
The following subsections describe the features, roughly in order of usefulness
to the average user:
3.8.1 Time Dependence
The use of time dependence in DEGAS 2 is described above
in Sec. 2.8. For completeness, we only note here
the principal controlling parameters in the
sources class. The switch
so_time_dependent will be FALSE (0) for steadystate
runs and TRUE (1) for the time dependent case.
For the latter, the time interval for the run will
be from so_time_initial to so_time_final;
both have units of seconds.
There is also a sourcegroup dependent parameter
so_time_varn that controls the time
variation of tahat group.
3.8.2 Checkpoints
By default, the DEGAS 2 output netCDF file is written only at the end of
the run. Each source group can be broken up into an integral number
of pieces with an intermediate output file being written after
each such checkpoint. The array so_chkpt(grp) contains the number
of checkpoints for source group grp. If so_chkpt(grp) is 0,
no file is written. If 1, an intermediate output file will be written
at the end of the source group, and so on.
The intermediate output file does not contain postprocessed results and
cannot be used with postprocessing utilities. The raw data there are
suitable only for a restart of the code.
3.8.3 Restarting
The two main uses of the restarting capability are
 To complete an interrupted (and checkpointed) run without having
to restart from the beginning,
 To increase the number of flights in a run to improve statistics.
In both cases, the user needs to change the value of the flag
so_restart from its default value of FALSE (represented
by 0 in the background netCDF file) to TRUE (1). The code will
expect to find an existing output file, with the name specified in
the degas2.in file. This can be either an intermediate
output file from a checkpoint or a completed one generated at the
end of a run.
The restarted run will extend the number of flights to the value specified
in the background netCDF file (by so_nflights). A few different
situations may arise:
 If the number of flights is the same as those in the output file,
no additional flights are run. This feature permits the data in an
intermediate output file to be postprocessed into a "completed"
output file.
 If the previous checkpointed run was interrupted, the restarted
run resumes from the point of the last checkpoint and completes the
specified number of flights. The results should match those that would
have been obtained if the run were to have run to completion on the first
try.
 If the user wishes to increase the number of flights for
the last source group of a completed run, the code will compute the
required number of flights. The results will be the same as if the
flights were all done at once.
 If the user wishes to increase the number of flights in more
than one source group, the code will extend the run accordingly. However,
the results will not exactly match those obtained in a single run of
the same length. The reason is that by default a run uses a continuous
random number chain for all source groups. In a restart, that chain is
resumed at its end. Hence, the first source groups will be using
random numbers for these new flights different from those they would have
had in a single run. Note that the resulting run should still be
statistically equivalent to the single run. The "spaced seeds"
option has been added to provide a workaround to this minor shortcoming
of the restart procedure.
3.8.4 Seed String
Prior to version 2.6 of DEGAS 2, the initial random number seed was
hardwired in flighttest.web. The user can now change the initial
seed (and, hence, the entire random number chain) via
so_seed_string in the background netCDF file.
The preferred mechanism for specifying the initial seed is as a
character string. This is transformed into the appropriate integer
representation using routines in random.web
(see Sec. 2.5). The default value is "12". Reasons
to use different values include
 An unexpected, but perhaps not statistically significant, result
has been obtained. A repeat of the run with a different random number
seed will aid in establishing the validity of the result,
 The DEGAS 2 code package contains several tools for
statistically analyzing results (e.g., matchout and
matcheir). By changing the random number seed, the effectiveness of
these statistical comparisons can be evaluated.
3.8.5 Spaced Seeds
If the user anticipates needing to extend a multiple source group run to
a larger number of flights (e.g., to achieve a desired variance) and
insists that the results match those of a single run, the flag
so_spaced_seeds should be set to TRUE (1) in the initial
run. With this option enabled, the code will set the initial seed
for each source group to be 100000000 (set by so_seed_spacing)
flights apart. Each source group
in a restarted run will then be able to "pick up where it left off",
up to a total of 100000000 flights.
To minimize the chances of compromising the integrity of the random
number chain, so_spaced_seeds should be left at its default
FALSE (0) value, unless the user truly needs to be able to extend
the run and match the single run results.
3.8.6 Direct Sampling
The default (random) method of sampling the initial particle source in
DEGAS 2 leads to a source distribution that differs from the ideal,
input distribution by a fraction of order 1/√N (N being the
number of flights). Yet, if we were to take those same N flights
and divide them up by hand amongst the source segments, we would
obtain an error of order 1/N. Such an approach is impractical in
general as it requires specifying the number of flights once and for
all at the beginning of the run (a restart would not be possible) and
the possibility of overlap in the initial source positions would arise
if the number of flights were large enough.
An equally effective alternative can be developed based on hashing
algorithms described by Knuth[53]. The idea is to replace the
pseudorandom number ξ usually used in the source sampling process with
where x is a single, fixed random number, i is the flight number,
and
is the golden ratio. This value yields the desired distribution properties
(as can be verified using sourcetest), but minimizes the possibility
of overlap by ensuring that ξ is far from low order rational numbers.
The sampling method used is governed by the flag so_sampling. The
default value of so_random yields the usual random sampling procedure.
This is the recommended value for most applications. The new, direct
sampling method is selected by setting so_sampling to
so_direct. This method was installed to see if it would lead
to more efficient iterations in coupling DEGAS 2 to plasma transport
codes.
3.9 Adding Reactions and PMI
3.9.1 Overview of Data Format
The same, basic approach is used to store all reaction and PMIrelated data
in DEGAS 2. First, the data for each reaction or PMI is stored in its own
netCDF file. That file contains:
 Name of reaction,
 Number of dependent variables,
 Organization of data for each dependent variable ("table" or
"fit"),
 Rank of each dependent variable (number of independent variables),
 Character name for each dependent and independent variable,
 Number of values of each independent variable (for tabular
data),
 Data table, a single onedimensional array. Indices into this
array are computed on the fly based on the number of values used for
each independent variable.
For more details, see the description of the
cross section class. The corresponding information for PMI is in the description of the
PMI format class.
The only general tool available for creating new data files is ratecalc
(see Sec. 3.4.2 and the ratecalc
file itself). Example input files for ratecalc can be found in the
data directory. Two other routines reactionwrite and
pmiwrite (again see Sec. 3.4.2) read existing text files
and write the data out into netCDF files. These two routines can be
used as examples for generating specific data files.
At run time, all of the reaction data files specified in the problem
input file (see Sec. 3.5.8) (and only those) are read
in. The reaction rates (for PMI, the corresponding quantity is termed
the "yield") and the descriptions of their independent
variables and other characteristics are compiled into one set of arrays.
The information is organized in a similar manner to that of the
original data files, but with an additional index corresponding to
the problem reaction number. Additional processing is done to replace
character strings with integer indices, to make for more efficient searches.
Again, all of the actual reaction rate
data values are stored as a single 1D array with macros used to provide
simple access to individual entries for a particular reaction.
The reaction rate is used in computing the current mean free path of the
flight and in choosing amongst the flight's possible reactions when a
collision is indicated.
All other data present in the reaction data files is compiled into a set
of "handling" arrays, indicating that their primary purpose is to
control the specification of collision product velocities for that
reaction. In addition to an index for problem reaction number, there is
one for dependent variable number. The collision routines will
search these data arrays for required information and will stop code
execution (in DEBUG mode, anyway) if they are not found. The flow
of data through the various DEGAS 2 classes is shown in
Fig. 3.2.
Figure
Figure 3.2: Flow of atomic and surface physics data through various DEGAS 2
classes during preprocessing.
The data format, evaluation routines, and scoring mechanisms are designed so
that specific quantities can be scored just by adding the appropriate
data to the data file and creating a corresponding tally with that
quantity as the dependent variable.
One aspect unique to the PMI is that specification of the outgoing velocity
distribution necessarily involves a "fit" since standard FORTRAN will
not permit a function (in this case, the function which interpolates
the data tables) to return a vector.
3.9.2 An Example: Bateman Format Data
A particularly involved example of a PMI is the specification of
"Bateman format" data[3,50]. A few uses of this for
reflection processes are currently in the code. A large set of data
for PMI in an improved version of this format will soon be added.
For each incident
energy and polar angle pair, (E_{in}, θ_{in}), we have a
conditional distribution for the product atom
P_{Ein,θin}(v,α,ϕ) v^{2} dv sinαdα d ϕ, 
 (3.19) 
where α and ϕ are the outgoing polar and azimuthal angles,
respectively.
This distribution is sampled from three separate 1D distributions. First,
the outgoing velocity v_{0} is specified by
f^{1}_{Ein,θin}(v) =  ⌠ ⌡

 ⌠ ⌡

P_{Ein,θin}(v,α,ϕ)sinαdα d ϕ. 
 (3.20) 
The outgoing polar angle α_{0} obeys the distribution
f^{2}_{Ein,θin}(α) =  ⌠ ⌡

P_{Ein,θin}(v_{0},α,ϕ) d ϕ. 
 (3.21) 
And, finally, the outgoing azimuthal angle ϕ_{0} is taken from
f^{3}_{Ein,θin}(ϕ) = P_{Ein,θin}(v_{0},α_{0},ϕ). 
 (3.22) 
The inverse cumulative distribution, F(ω) = G^{−1}(ω), with
is specified, say, 5 values of 0 < ω < 1. For example,
for a normally incident θ_{in}) = 0,
(E_{in} = 1 eV H atom being reflected off
of Fe, the data are
R_{N}(E_{in},θ_{in}) =
8.27750E01
F^{1}(ξ) =
Footnotes:
^{1}This file is written in T_{E}X.
A hyperlinked PDF
file is generated from that source using pdflatex. This document
is distributed as part of the DEGAS 2 code. The user can be best assured
of concurrence between code and manual by referring to the
version of the manual
that came with the code. Copyright © 2013 PPPL.
Permission is granted to make and distribute verbatim copies of this manual
provided the copyright and permission notice is preserved on all copies.
Permission is granted to copy and distribute modified versions of this
manual under the conditions for verbatim copying, provided also that
the entire resulting derived work is distributed under the terms of a
permission notice identical to this one.
Permission is granted to copy and distribute translations of this manual
into another language, under the above conditions for modified versions.
^{2}In fact, the
values of the "total" tallies
for D2+ are all determined by roundoff errors. For this reason,
you will not likely match these values exactly. You should convince
yourself that they are numerically much smaller than the other totals. If you
use matchout to compare with the reference output netCDF file, these
D2+ totals will also stand out as apparently significant discrepancies.
File translated from
T_{E}X
by
T_{T}H,
version 4.03.
On 20 Mar 2013, 17:27.