Figure User's Guide for DEGAS 2
CVS Revision: 1.4 ¹

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 cover all aspects of DEGAS 2 from the user's point-of-view: 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. Information for programmers interested in modifying DEGAS 2 can be found within the source code itself. Links to some of the more useful ones 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 multi-species, 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 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 DEGAS-B2 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₂.

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 1992-1993:

High speed,
Flexibility,
Ease-of-use,
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 quasi-object 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 pre-processor macros) in DEGAS 2 is limited. Furthermore, there is no provision for inheritance within this programming environment.
Provide an easy-to-use and well documented code by writing it with the FWEB[7] package. FWEB provides the powerful preprocessing and macro facility which underlies the "quasi-object oriented" implementation and internal T_EX 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 self-describing, platform-independent binary files. Presently, the main code uses exclusively netCDF[8] (Network Common Data Format) files for input and output. However, a simple post-processing utility is included which generates graphics and data in HDF[9] (Hierarchial Data Format) 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 pseudo-random 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 three-dimensional 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 S_t. We write the probability density function for having the next collision between x and x + dx, p(x) dx, as:

p(x) dx = S_t exp(-S_t x) dx.

(2.1)

Then, the cumulative probability distribution P(y),

P(y) =

ó
õ

y

0

p(x) dx,

(2.2)

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 x 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 x, with 0 £ x < 1. For this example, P(y) is easily computed and inverted so that

y = -

lnx

S_t

(2.3)

where we have surreptitiously replaced 1 - x with x since both are uniform over the unit interval. Again the logic of this process, as depicted in the figure, is that if x is uniformly distributed, then numbers between x to x+ Dx are chosen with the same frequency as particles having collisions between y and y + Dy. 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) =

v_th

exp

é
ë

(v -v₀)²

v_th²

ù
û

(2.4)

For simplicity, limit the v values to |v - v₀| < 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,x) pairs and tossing out ones with x > p(x).

Now, we choose a pair of random numbers (x, x) with x uniformly distributed over the interval v₀ - M v_th £ x < v₀ + M v_th and x uniform over 0 £ x < p_max, where p_max is the maximum value of p(v). If x < p(x), accept this x as the sampled velocity. Otherwise, if x > p(x), reject these (x, 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 well-known Box-Mueller 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 x < 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 non-analog techniques (the rejection technique described in the previous subsection is a non-analog 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,10] can be readily digested.

From an analog point-of-view, 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, c(x) as the sum of two pieces:

Arrivals directly from source: Q(x).
Particles transported from elsewhere:

ó
õ c(y) L(y ® x) dy
(2.5)

And we have

c(x) = Q(x) +

ó
õ

c(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\vec],[v\vec]) and a species index i. The kernel L can then be broken into two pieces:

, i^¢ ®

, i) = T(

;

, i) C(

, i^¢ ®

, i;

(2.7)

where T([r\vec]^¢ ® [r\vec]; [v\vec], i) is the (probability) density of particles leaving a collision at [r\vec]^¢ and having a subsequent collision at [r\vec]. 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\vec]^¢, i^¢ ® [v\vec], i; [r\vec]) which provides the probability that a particle of velocity [v\vec]^¢, species i^¢ will have a collision resulting in a particle of velocity [v\vec], 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 c(x). Furthermore, there is no well-established 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 c(x) is called the post-collision density. The origin of this designation should be clear from the previous discussion. Equivalently, the integral equation can be written in terms of the pre-collision density y(x). The two are related by

,i) =

ó
õ

d³r^¢ c(

,i) T(

;

,i).

(2.8)

It is this pre-collision density which is most closely related to the more familiar particle flux F and particle distribution function f by the definitions:

,i) = S_t(

,i) F(

,i) = S_t(

,i) |

| f(

,i).

(2.9)

The integral equation for y bears a strong resemblance to that for c. 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:

,i)

ó
õ

¥

0

ó
õ

d³v C(

v^¢

,i^¢ ®

,i;

×T(

R ®

;

,i)y(

v^¢

,i^¢),

(2.10)

where

,i) =

ó
õ

¥

0

dR Q(

, i)T(

R ®

;

,i),

(2.11)

and

R ®

;

,i) = S_t(

,i) exp

é
ë

ó
õ

R

0

dR^¢ S_t(

R^¢,

,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 phase-space volume integrals with the path enforced via delta-functions[10] was copied from Case and Zweifel[11] 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³r d³v f(

)g(

(2.13)

For example, if we wish to estimate a reaction rate r º S|[v\vec]|,

I(r) =

ó
õ

d³r d³v SF(

(2.14)

Values of r are accumulated along the random walk of the neutral flight and used to estimate I(r). The mathematical description of this accumulation process is called the estimator of I.

2.3.1 Simple Estimators

Consider a random walk g which consists of k steps, (x₁, x₂, ¼, x_k). I.e., the particle is absorbed or otherwise terminated at the kth step. The simplest estimator involves writing down the value of r at the point of absorption. Thus, the absorption estimator is

x_a(g) =

r_a

w c_V,

(2.15)

where r_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 c_V

c_V (x_m) =

ì
ï
í
ï
î

if x_m Î V

if x_m Ï V

(2.16)

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 g of the probability of each walk P(g) times the corresponding estimator,

E(x_a) =

å
g

P(g) x_a(g).

(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.

x_c1(g) =

k
å
m=1

é
ë

r_a (x_m)

r_t (x_m)

ù
û

wc_V (x_m),

(2.18)

where r_t = r_a + r_s is the total reaction rate for the particle, with r_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(x_c1) =

å
g

P(g) x_c1(g).

(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

r_s (x_j)

r_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

x_c2(g) =

k
å
m=1

w_m

é
ë

r_a (x_m)

r_t (x_m)

ù
û

c_V (x_m).

(2.21)

For formal proofs that these expectation values do indeed yield I(r), 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[47] (who had improved upon an earlier approach described by Spanier[48]). To make the connection to DEGAS 2 clearer, we replace the S and d variables of Macmillan with r and t.

Macmillan's paper, like the earlier one by Spanier, aims to derive the tracklength estimator. Collisions are described as either scatterings, rate r_s, or absorptions r_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 (r_s + r_a) dt. Again, the weight reduction factor applied at each scattering is r_s / (r_s + r_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 a with the delta scatterings. To obtain the estimators used in DEGAS 2, we work with a total fictitious scattering reaction rate

r_t,d º r_s,d + ar_a,

(2.22)

where r_s,d is equivalent to Macmillan's S_s,d. Then, the probability of a delta scattering in a time interval dt is r_t,d dt. At each delta scattering, the weight is reduced by a factor (r_t,d - ar_a) / r_t,d. The probability of a real scattering in the same interval is (r_s + (1 - a) r_a) dt; at each such event, the weight is reduced by r_s / [r_s + (1 - a) r_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) =

( r_t,d t)ⁿ

exp(- r_t,d 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):

x =

r_t^*

(2.24)

where w is the current weight of the particle and

r_t^* = r_t,d + r_s + (1 - a) r_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(x | t) =

¥
å
n=0

P(n | t)

r_t^*

n
å
i=0

æ
è

r_t,d - ar_a

r_t,d

ö
ø

(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(x | t), we find that we can write

E(x | t) = c^¢

r_s + (1 - a) r_a

+ (1 - c^¢) r

1 - exp(- ar_a t)

ar_a

(2.27)

where

c^¢ º

r_s + (1 - a) r_a

r_t,d + r_s + (1 - a) r_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( - ar_a t)

r_s

r_s + (1 - a) r_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( - ar_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 r_t,d ® 0, we should get no score at all. In contrast, MacMillan's pure "collision estimator" still compiles a nonzero score as S_s,d ® 0 through the aS_a contribution to the pseudocollision cross section. Thus, MacMillan's "collision estimator" behaves somewhat like a track-length 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(x | t) =

r_t^*

r_t,d

1 - exp(- ar_a t)

ar_a

(2.30)

Clearly, E(x | t) ® 0 as r_t,d ® 0. The limiting expressions above for r_t,d ® 0 can be obtained from MacMillan's with S_s,d ®- aS_a.

We now consider the expressions used in DEGAS 2. We start by using a = 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 r_t,d = 0,

x_c º

r_s

exp(- r_a t) w c_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 c_V for completeness. At the end of the timestep (after the collision is scored) the weight is reduced by

w^¢ = w exp( - r_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 track-length estimator is obtained with r_t,d ®¥,

x_t º r

1 - exp(- r_a t)

r_a

w c_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 a = 0 in the above expressions:

x_c º

r_s + r_a

w c_V,

(2.34)

and

x_t º rt w c_V.

(2.35)

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),

x_g^TLE = w

1 - exp(-s_i t_V)

s_i

(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 reaction-related, e.g., ion momentum source due to charge exchange. The latter g quantity enter the track-length estimator as averages over the background distribution[14].

The second type of estimator is the collision estimator,

x_g^CE =

å
m_s

w_{m_s}

s_s

(2.37)

where the sum is over m_s, all scattering collisions within V. The quantity w_{m_s} is the weight at m_sth 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,

x_g^CE =

å
m_j

w_{m_j}

s_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₂ 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³r g(

)

ó
õ

d³v f(

) =

ó
õ

d³ r g(

) n(

(2.39)

This leads to the post-processing estimator. The requirement is that g not depend upon the test or background velocities. An example of this is the H_a 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., chord-integrated)],
The dependent variable (e.g., energy, background momentum, H_a 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 (track-length, collision, post-processed, 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\vec] 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: 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 post-collision 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 c₀([r\vec]^¢,[v\vec]^¢,i) using the source term Q([r\vec]^¢,[v\vec]^¢). The subscript 0 denotes the initial track for the flight.

The next stage of the algorithm involves computing the reaction rates r_s, r_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 post-collision density c_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 non-absorbing reaction rates using a random number z.

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\vec]^¢ to [r\vec]. At the same time, this step represents the conceptual transition from post-collision density to pre-collision density y, 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 y_n+1 to c_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 geometry-related 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 3-D 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 2-D 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 3-D, but the output results are averaged over the third dimension. As will be discussed in more detail below, 1-D and 3-D "symmetries" are also permitted. Most geometry related variables refer to Cartesian coordinates. However, in cylindrically symmetric or nearly symmetric (3-D) 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 2-D plasma code like UEDGE. In these cases, the plasma flow velocities are effectively specified at a toroidal angle f = 0, corresponding to y = 0. To use these velocities in a collision at some particular location for which f ¹ 0 (equivalently, y ¹ 0), the velocities need to be rotated through that angle f. Likewise, a code like UEDGE expects to get from DEGAS 2 toroidally symmetric ion momentum sources specified at f = 0. So, the scores for those sources need to be rotated from the local toroidal angle f ¹ 0 to f = 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 2-D geometry in which the y coordinate is ignored.
Cylindrical: A 2-D geometry in which the toroidal angle f is ignored.
1-D: A 1-D (plane) geometry in which the y and z coordinates are ignored.
3-D Plane: This geometry is close to be planar symmetric; all of the geometrical objects are defined using 2-D 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 3-D geometry.
3-D Cylinder: This is the cylindrical equivalent of the 3-D Plane symmetry. The full range of toroidal angles (360 degrees or 2 p radians) is included in this case. Obviously, the problem is divided in the f direction by f = constant surfaces. The analogous simplified characterization of this approach is the "pie slice" technique.
3-D Cylindrical Section: This is also a cylindrical equivalent of the 3-D plane symmetry. The difference is that in this case only a finite range of toroidal angles is considered. The problem is bounded by f = 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 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[5,6]. The physical processes included in such a collisional-radiative model are:

Electron collisional excitation K_e,mn and de-excitation K_d,mn,
Spontaneous radiative transitions A_mn,
Electron collisional ionization K_i,mn,
Three-body recombination K_r,mn,
Radiative recombination b_rad,m,
Dielectronic recombination (for multi-electron atoms) b_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 m-th excited state of the atom, N_m,

d N_m

d t

- N_m

é
ë

n_e

æ
è

K_i,m +

å
n > m

K_e,nm +

å
n < m

K_d,nm

ö
ø

å
n < m

A_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 ( b_rad,m + b_dia,m + n_eK_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 collisional-radiative model[5,6] 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₁ 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⁽ⁱ⁾_m

N₁

ö
ø

N₁ +

æ
è

N⁽ⁱⁱ⁾_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, b_rad,m, b_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 and n_i = 0, providing the "coupling to the ground state",
Solving the equations again with N₁ = 0 and n_i = 1, yielding the "coupling to the continuum".

Tables of the population coefficients N⁽ⁱ⁾_m/N₁ and N⁽ⁱⁱ⁾_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₁

d t

= - N₁ n_e S_eff + n_e n_i R_eff + transport, etc.,

(2.42)

with

S_eff =

å
m ³ 1

æ
è

N⁽ⁱ⁾_m

N₁

ö
ø

K_i,m,

(2.43)

and

R_eff = -

å
m ³ 2

æ
è

N⁽ⁱⁱ⁾_m

n_i

ö
ø

K_i,m+

å
m ³ 1

( b_rad,m + b_dia,m + n_e K_r,m ).

(2.44)

S_eff is thus the effective ionization rate of the ground state neutral, incorporating the multi-step 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⁽ⁱ⁾_m/N₁ and N⁽ⁱⁱ⁾_m/n_i may also be tabulated or used to compute photon emission rates for specific transitions of interest (e.g., Balmer-a).

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¹S and 2³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):

d E_e

d t

E₁ ( - N₁ 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

b_rad,m ( áE_rad,m ñ+ E_m)

- n_e n_i

å
m ³ 1

b_dia,m ( áE_dia,m ñ+ E_m),

(2.45)

where

E₁ 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[14],

áE_x,m ñ = T_e

é
ë

T_e

b_x,m

d b_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₁ ( N₁ n_e S_eff - n_e n_i R_eff ) - N₁ E_loss⁽ⁱ⁾ - n_i E_loss⁽ⁱⁱ⁾,

(2.47)

where the terms E_loss⁽ⁱ⁾, E_loss⁽ⁱⁱ⁾ 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₁ 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_loss^(i)¢

E₁ S_eff + E_loss⁽ⁱ⁾/ n_e,

(2.49)

E_source^(ii)¢

E₁ R_eff,

(2.50)

E_loss^(ii)¢

E_loss⁽ⁱⁱ⁾ / n_e.

(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.8.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[17]. 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 Janev-Smith database[15]. The resulting collisional-radiative 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. 6.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³ s^-1,
Recombination rate R_eff in cm³ s^-1,
Neutral electron losses E_loss⁽ⁱ⁾ in erg s^-1,
Continuum electron losses E_loss⁽ⁱⁱ⁾ in erg s^-1,
Neutral "n=3 / n=1", N₃⁽ⁱ⁾ / N₁,
Continuum "n=3 / n=1", N₃⁽ⁱⁱ⁾ / n_i,
Neutral "n=2 / n=1", N₂⁽ⁱ⁾ / N₁,
Continuum "n=2 / n=1", N₂⁽ⁱⁱ⁾ / 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 row-by-row. There are 10 rows of 6 entries each.

The density of the n = 3 excited state can be found with

N₃ =

æ
è

N⁽ⁱ⁾₃

N₁

ö
ø

N₁ +

æ
è

N⁽ⁱⁱ⁾₃

n_i

ö
ø

n_i.

(2.54)

The corresponding emission rate for the Balmer-a line (a.k.a. "H_a", n = 3 ® 2) is given by N₃ A₂₃ = 4.41 ×10⁷ N₃ cm^-3 s^-1. Note that the energy of this transition is E₂₃ = 13.595 (2^-2 - 3^-2) eV. Similar expressions give the n=2 density and the Lyman-a (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-a emission rate as a power, per atom (for ionization) or ion (for recombination) and per electron:

P_a⁽ⁱ⁾

E₂₃ A₂₃

æ
è

N₃⁽ⁱ⁾

N₁

ö
ø

/ n_e,

(2.55)

P_a⁽ⁱⁱ⁾

E₂₃ A₂₃

æ
è

N₃⁽ⁱⁱ⁾

n_i

Figure User's Guide for DEGAS 2 CVS Revision: 1.4 1

Daren Stotler and Charles Karney Princeton Plasma Physics Laboratory Randall Kanzleiter Rensselaer Polytechnic Institute Jaishankar S. Institute for Plasma Research, India

Abstract

Chapter 1 Introduction

1.1 Purpose of Monte Carlo Neutral Transport

1.2 Historical Background

1.3 Need for DEGAS 2

1.4 DEGAS 2 Features

Chapter 2 Background

2.1 Generic Monte Carlo Algorithms

2.1.1 Sampling a Distribution

2.1.2 Sampling Distance to Next Collision

2.1.3 Sampling More Complicated Distributions

2.2 Particle Transport

2.2.1 Motivating the Integral Equation

2.2.2 Transport and Collision Kernels

2.2.3 Integral Equation Details

2.3 Estimators

2.3.1 Simple Estimators

2.3.2 Generalized Collision and Tracklength Estimators

2.3.3 Specification of Tallies

2.4 Tracking Procedure

2.5 Random Numbers

2.6 Geometry

2.7 Symmetry and Coordinate Systems

2.8 Atomic Physics

2.8.1 Hydrogen Collisional Radiative Model

Figure User's Guide for DEGAS 2
CVS Revision: 1.4 ¹

Daren Stotler and Charles Karney
Princeton Plasma Physics Laboratory
Randall Kanzleiter
Rensselaer Polytechnic Institute
Jaishankar S.
Institute for Plasma Research, India

Chapter 1
Introduction

Chapter 2
Background