The EIRENE Code User Manual

including: B2-EIRENE interface

Version: 11/2009

D. Reiter

EIRENE

THE GREEK GODDESS OF PEACE

The Greek author Hesiodos wrote a genealogy of the gods, The Theogony, in the 8th century B.C.. According to him, EIRENE and her sisters, Eunomia and Dike, were the daughters of Zeus and Themis. These sisters were called the horae, the Greek word for the right time, hora.

Abstract

The EIRENE neutral gas transport code is described. This code resorts to a combinatorial discretization of general 3 dimensional computational domains. It is a multi-species code solving simultaneously a system of time dependent (optional) or stationary (default) linear kinetic transport equations of almost arbitrary complexity. A crude model for transport of ionized particles along magnetic field lines is also included. EIRENE is coupled to external databases for atomic and molecular data and for surface reflection data, and it calls various user supplied routines, e.g. for exchange of data with other (fluid-) transport codes. The main goal of code development was to provide a tool to investigate neutral gas transport in magnetically confined plasmas. But, due to its flexibility, it also can be used to solve more general linear kinetic transport equations, by applying a stochastic rather than a numerical or analytical method of solution. In particular, options are retained to reduce the model equations to the theoretically important case of the one speed transport problem.
Major applications of EIRENE are in connection with plasma fluid codes, in particular with the various versions of the B2 code [1]. The semi-implicit iterative coupling method of B2-EIRENE [2], [3] and it's implementation (code segment: EIRCOP) are also described.

Foreword to 2nd edition

The first edition of this EIRENE code user manual was published as KFA report JÜL-2599 in March 1992, [4]. It basically was a collection of my notes on the meaning of the various input data options for EIRENE. Over the years too many such input flags had accumulated to memorize the meaning of each individual one.

Since the EIRENE code became a quite popular tool also for many other users, it was decided to provide these notes as a kind of user's manual in a somewhat completed and edited form.

In the meantime the distribution of the EIRENE code has become even wider, and it seems timely to update the previous manual, although most of it's content is still a relevant source of information for a user of the code. Apart from several minor corrections, e.g., of spelling errors and unclear language, the major new features as compared to the previous edition are the following:

2.13

An increasing number of EIRENE users is running cases in which EIRENE is coupled to a fluid plasma model to provide sources and sinks. The widely used B2-EIRENE code system [2] is one such example. A new section describing this part of the EIRENE code has been added. It describes the "sandwich" file EIRCOP which was written to permit linkage between EIRENE and a plasma fluid model (one, two, or three dimensional. Initial and boundary value problems can be treated in a rather self-consistent manner.
This code segment was developed largely under support of a KFA-EURATOM contract ([5]), and first applications have been published for ITER configurations, see [2] and [6]. Subsequently a large number of application runs, in particular those carried out at IPP-Garching (Ralf Schneider et al.) and AEA-Culham Laboratories (Geoff Maddison) have lead to identification of many critical issues and limitations of the code and to permanent improvements until these days. Typical convergence behaviour of the combined code system ("saturated residuals") is described in the new section 1.8, written in collaboration with G.P.Maddison, AEA Culham Laboratories.
Detlev Reiter, Spring 1998

Foreword to 3rd edition

During the year 2001 a major revision of the EIRENE code has been carried out, largely in order to implement a somewhat more modern FORTRAN structure into the code. EIRENE code versions since then have sometimes been referred to as "EIRENE-FACELIFT".

In particular a dynamical allocation of storage (rather than the hitherto necessary pre-assignment of storage in the PARAMETER statements collected in the previous PARMUSR file) is now in place. This dynamic memory managment has led to a replacement of the previous "Common-Blocks" now by "Modules", and the entire elimination of the "Equivalence" statements used for storage economy in earlier versions.

Further main upgrades are related to:

incident-species dependent surface interaction models, pumping speeds, sputter models etc., can now be controlled by input flags, rather than the previous quite difficult implementation via the USR-routines.
integration volumes of volume averaged tallies can now be larger than the grid cell volumes, i.e., the plasma background and geometrical discretisation, which is not very storage sensitive, can be made much finer than the neutral particle volume averaged tally profile output.
discretisation of general multiply connected 3D volumes by tetrahedron-grids.
photons as new type of test particle species. The photon_dummy module allows simple photon transport (one-speed Boltzmann-equation) with spatial profiles of spontaneous emission, Doppler line broadening, for optically thin lines only, with inclusion of (multiple) surface reflection. The full photon module additionally allows for various other line broadening mechanisms and photon re-absorption, stimulated emission, radiation trapping (excitation-hopping), iteration with neutral gas excitation population kinetics, etc...

EIRENE versions 2003 and younger have been compiled and successfully tested at FZ-Jülich on the following systems and compilers:

Linux: Suse 11.1 and all predecessors down to Suse 6.***:

pgf90 Portland Compiler Version 5.2-4 - 10.1
ifort Intel Fortran Compiler Version 8.1 - 11.1
lf95 Lahey Fortran Compiler Version 6.2
nagf95 NAG Fortran Compiler Version 5.0 - 5.2

AIX: AIX 4.3 und 5.2

xlf95 XL Fortran for AIX Version 8.1

Windows: Windows 2000 und Windows XP

Compaq Visual Fortran Version 6.6

Detlev Reiter, Spring 2010

1 The neutral gas transport equation;
Monte-Carlo terminology
    1.1 The linear Boltzmann equation for the distribution function f
    1.2 The linear integral equation for the collision density Ψ
        1.2.1 The Green function concept
    1.3 Monte Carlo solution of equation 1
        1.3.1 Unbiased estimators
        1.3.2 Scaling of tallies
        1.3.3 Efficiency (FOM)
        1.3.4 Source sampling
        1.3.5 Sampling from the transport kernel T
        1.3.6 Sampling from the collision kernel C
        1.3.7 elastic collisions
        1.3.8 charge exchange
        1.3.9 electron-impact collisions
        1.3.10 general heavy-particle-impact collisions
        1.3.11 photon processes (emission, absorption, scattering
    1.4 Surface Reflection Models
        1.4.1 The Behrisch matrix reflection model
        1.4.2 TRIM code database reflection models
    1.5 Surface sources
        1.5.1 The electrostatic sheath
    1.6 Combinatorial description of geometry
    1.7 time dependent mode of operation
    1.8 non-linear effects: coupling to plasma fluid models
    1.9 non-linear effects: neutral-neutral collisions
        1.9.1 BGK approximation
        1.9.2 Direct Simulation (DMCS) of self-collisions
    1.10 Radiation transport, photon gas simulations
        1.10.1 Line shape options
    1.11 Charged Particle Transport
        1.11.1 Orbit integration
        1.11.2 Coulomb Collision Models
        1.11.3 The full Fokker-Planck Collision Model, A
        1.11.4 The full Fokker-Planck Collision Model, B: binary Monte Carlo collision procedure
    1.12 EIRENE flow charts
2 Description of formatted input file
    2.1 Input data for operating mode
        2.1.1 The NLERG option for cell volumes
        2.1.2 The NLMOVIE option for making movies of trajectories
    2.2 Input for Standard Mesh
        2.2.1 Mesh Parameters
    2.3 The Input Block for Surfaces
        2.3.1 The Input Block for "Non-default Standard Surfaces"
        2.3.2 Input Data for "Additional Surfaces"
    2.4 Input Data for Species Specification and Atomic Physics Module
        2.4.1 Default atomic and molecular data
        2.4.2 Neutral-Neutral collisions in BGK approximation
        2.4.3 Fitting expressions (IFTFLG)
    2.5 Input for Plasma Background
        2.5.1 Derived Background Data
    2.6 Input Data for Surface Interaction Models
    2.7 Input data for Initial Distribution of Test Particles
        2.7.1 Piecewise constant "Step-functions" for sampling
    2.8 Additional Data for some Specific Zones
    2.9 Data for Statistics and non-analog Methods
    2.10 Data for additional volume and surface averaged tallies
        2.10.1 Additional volume averaged tallies, tracklength estimator
        2.10.2 Additional volume averaged tallies, collision estimator
        2.10.3 Algebraic expression in volume averaged tallies
        2.10.4 Algebraic expression in surface averaged tallies
        2.10.5 Algebraic expression in surface averaged tallies
        2.10.6 Energy spectra in selected cells or surfaces
    2.11 Data for numerical and graphical Output
    2.12 Data for Diagnostic Module
        2.12.1 Line of sight: charge exchange spectrum
        2.12.2 Line of sight: line emissivity
        2.12.3 Line of sight: line shape
        2.12.4 Line of sight: user defined integral
    2.13 Data for nonlinear and time dependent Options
    2.14 Data for interfacing Subroutine "INFCOP" (example)
        2.14.1 Version B2-EIRENE-1999 and older
        2.14.2 Version B2-EIRENE-2000 and younger
3 Problem specific Routines
    3.1 Parameter Statements (for EIRENE-2001 or older)
    3.2 The "Additional Tally" routines UPTUSR, UPCUSR, UPSUSR and UPNUSR
        3.2.1 Tracklength estimated volume tallies, UPTUSR
        3.2.2 Collision estimated volume tallies, UPCUSR
        3.2.3 Snapshot estimated volume averaged tallies, UPNUSR
        3.2.4 Surface averaged tallies, UPSUSR
    3.3 The user surface reflection model REFUSR
    3.4 The user source sampling routine SAMUSR
    3.5 The user routines to overrule input data
        3.5.1 The user geometry data routine GEOUSR
        3.5.2 User supplied background data routine PLAUSR
    3.6 The user routines for profiles PROUSR
    3.7 User supplied post-processed tally routine TALUSR
    3.8 User supplied "general geometry block"
        3.8.1 Subroutine INIUSR
        3.8.2 Subroutine LEAUSR
        3.8.3 Subroutine TIMUSR
        3.8.4 Subroutine VOLUSR
        3.8.5 Subroutine NORUSR
4 Routines for interfacing with other codes:
EIRCOP
    4.1 Routine for interfacing INFCOP
        4.1.1 entry IF0COP
        4.1.2 entry IF1COP
        4.1.3 entry IF2COP(ISTRA)
        4.1.4 entry IF3COP(ISTRAA,ISTRAE)
        4.1.5 entry IF4COP
    4.2 Routines for cycling of EIRENE with external codes: EIRSRT
    4.3 Routines for special tallies needed for code coupling: UPTCOP
    4.4 Statistical noise in Monte Carlo terms for external code, noise-residuals: STATIS_COP
5 Default EIRENE tallies, and selected Modules
    5.1 Tables of EIRENE tallies
        5.1.1 Current status, incl. photon gas tallies (Eirene-02 and younger)
        5.1.2 old version, w/o photon gas tallies (Eirene-01 and older)
References

Introduction and General Information

This manual describes the input required by the EIRENE code to run a Monte Carlo study for a fully 3 dimensional simulation of linear transport (i.e., of test particles) in a prescribed background medium. Although the geometry of the problem and the interaction between test particle species and the background are in principle not subject to any restrictions, the aim of code development was to provide a tool for investigating neutral gas transport in tokamak plasmas. Consequently the choice of preprogrammed options has been made mainly with this application in mind. A large variety of problems in this field can be run without having to resort to any problem specific routines but instead by an appropriate setting of logical and numerical input flags.
Any user of EIRENE should be aware that this code is a moving target as is this manual.
Therefore it is possible that there are some inconsistencies between this description and a particular version of the code. The user should always first check subroutine INPUT, where most of the data are read in using hard wired formats, as most input errors will lead to a rapid exit in the initialization phase of a run.
This manual was written by the author of the code who often may not have been able to anticipate difficulties in understanding the use of EIRENE. He therefore gratefully acknowledges any suggestions to make this manual more informative and clear than it might be at the present status.
This code description consists of the following parts:

In the first part an introduction is given to the general linear transport problem and its solution by Monte Carlo methods. Most of the terminology used in later sections is introduced there.
In part two a description of the formatted input file required by EIRENE to run on a specific problem is given. It mainly consists of explanations of the meanings of the various input flags.
In part three the most important problem specific routines are explained. At present we have restricted this part to routines for evaluation of "user requested tallies", namely the subroutines UPTUSR, UPCUSR and UPSUSR. Other often applied routines such as SAMUSR (user supplied source sampling distribution) or REFUSR (user supplied surface interaction model) are briefly described.
In part four the package EIRCOP for interfacing EIRENE with other codes (e.g., the B2-EIRENE package) is described. Here mainly the location of the storage on the EIRENE work array RWK for plasma data and geometrical information is given. Also the implementation of the method of (semi-implicit) corrections (see section 1.8) at each plasma code time-step to the terms transferred from EIRENE to plasma fluid transport codes is described here.

Chapter 1
The neutral gas transport equation;
Monte-Carlo terminology

General Remarks

To introduce the terminology used throughout this report, we briefly recall the basic definitions and principles of a Monte Carlo linear transport model, following the lead of many textbooks on Monte Carlo methods for computing neutron transport (see e.g., [7]). We begin with the linear transport equation for the pre-collision density, written as integral equation (linear non-homogeneous Fredholm integral equation of 2nd kind). Distinct from standard terminology in the (analytic) transport theory we do not discuss analytic properties of the various terms in the equation, but, instead, point out their probabilistic interpretation, as needed for a Monte-Carlo solution of that equation. Next (section 1.3) we sketch the Monte Carlo procedure to solve such equations, by referring to the two most often applied techniques: "track-length - and collision based estimators". In the third subsection we briefly describe the treatment of boundary conditions (models for interaction of the particles with surrounding surfaces) and discuss some special models which are in use for the neutral gas transport in fusion plasma devices. Then (section 1.5) we discuss the most important source function (non-homogeneous part of the integral equation) and its implementation in a Monte Carlo algorithm, namely the surface source of neutral particles due to recombination of ions incident on solid surfaces at the boundary or inside of the computational area. In section 1.6 we comment on the implementation of geometry within the framework of a Monte Carlo code in general terms and for the EIRENE code in particular. In section 1.7 the time dependent mode of operation of the EIRENE code is described. It merely amounts just to an increase in the dimensionality by one, by adding a time co-ordinate and treating it, formally, in a rather symmetric fashion with the other spatial coordinates. See also reference [8].
Two kinds of non-linearities may be accounted for:
In section 1.8 the nonlinear behavior resulting from background data, which depend on neutral particle transport (sources and sinks), is discussed. The algorithm of the B2-EIRENE code system [3] is described. More details on this can be found in the report [9].
In the final section of this introductory chapter 1.9 the non-linear BGK formalism and the direct simulation method (DMCS) for self collisions between neutral particles, as implemented in the EIRENE code, is described.

1.1 The linear Boltzmann equation for the distribution function f

The term μ-space in this report refers to the phase space of a single particle. The quantity of interest is then the one particle distribution function f

f(r, v, i, t) or f(r, E, Ω, i, t) or f(x) ,

where the state x of the μ-space is characterized by a position vector r, a velocity vector v, a species index i (i stands for, e.g., H,D,T,D₂,DT,He, CH_n,...) and the time t.

The number density n_i(r) for species i then reads:

n_i(r,t) =

⌠
⌡

d³v f(r,v,i,t)

Instead of v we sometimes utilize the kinetic energy E, E=m/2 ·v², and the unit (speed) vector Ω = v/|v| in the direction of particle motion. Hence:

f(r, v, i ,t) =

⎛
⎝

⎞
⎠

(r, E, Ω, i, t)

We start by assuming that collisions are discontinuous events.

Further: Lets consider only one specific species i₀, now omitting this species index. We assume that there are only collisions of this species i₀ with only one further species (labelled b), and that exactly one particle of each of these species will also be present after the collision event (i.e., chemical reactions excluded for the time being, but see below). The familiar Boltzmann equation for the distribution function f for this species i₀ reads

[

∂

∂t

+ v·∇_r +

F(r,v,t)

·∇_v] f(r,v,t)

⌠
⌡

σ(v′,V′;v,V)|v′−V′| f(v′)f_b(V′)

−

⌠
⌡

σ(v,V;v′,V′)|v−V| f(v) f_b(V)

Q(r,v,t)

(1.1a)

Q is any external source (particles per unit time injected per unit volume in phase space). Integrations are over the velocities v′,V and V′. σ(v′,V′;v,V) is the cross section for a binary particle collision process. It is equal to zero unless the conservation laws for total energy and momentum are fulfilled.

The first two arguments in σ, namely the velocities v′,V′ in the first integral, correspond to the species i₀ and b, respectively, prior to a collision. These are turned into the post collision velocities v, V, again for species i₀ and b, respectively. The first integral, therefore, describes transitions (v′,V′→ v,V) into the velocity space interval [v,v+dv] for species i₀, and the second integral describes loss from that interval for this species.

Furthermore, m is the particle mass and F(r,v,t) is the volume force field.

The right hand side is the collision integral [(δf)/(δt)]|_{_b}. If there are more than just one possible type of collision partners, then the collision integral has to be replaced by a sum of collision integrals over all collision partners b, including, possibly, b=i₀ (self collisions)

δf

δt

|_{_coll} =

δf

δt

|_{_gain}−

δf

δt

|_{_loss} =

∑
b

δf

δt

|_{_{coll_b}}.

(1.2)

This is readily generalized to the semi-classical Boltzmann equation for chemical reactions (including, for example, vibrational relaxation or exchange of internal energy as special cases) symbolized as i₀ + j₀ \rightleftarrows i₁ + j₁. These species indices label both the chemical species and/or the internal quantum state. In this case the sum in the collision integral is over j₀, i₁ and j₁ and the cross sections in the corresponding collision integrals σ_i₀j₀^i₁j₁(v,V,v′,V′) are differential for scattering at a certain solid angle and post collision energies with simultaneous transition from (i₀, j₀) to (i₁, j₁). Further generalizations to include particle splitting, absorption or fragmentation into more than 2 post collision products are straight forward, but can more conveniently be formulated in the C-collision kernel formulation used below.

All these collision operators are bi-linear in the distribution functions. The first term on the right hand side is due to scattering into the element dv of velocity space and we shall abbreviate it by defining the collision kernel ("redistribution function") C as a proper integral over pre- and post collision velocities of species b-particles:

δf

δt

|_{_gain} =

⌠
⌡

d³v′ C(v′→ v)|v′| f(v′)

(1.3)

The kernel C can be a quite complicated integral, as is involves not only multiple differential cross sections, but also, possibly, particle multiplication factors, e.g. in case of fission by neutron impact, dissociation of molecules by electron impact, or stimulated photon emission from excited atoms. It can also include absorption, in which case the post collision state must be an extra "limbo" state outside the phase-space considered. Due to both particle multiplication and/or absorption the collision kernel C is not normalized to one, generally.

If the distributions f_b are assumed to be given, then the kernel C is linear and the expression above becomes a linear integral operator.

The second term on the right hand side is much simpler, because the function f(v) can be taken out of the integral. We even take the product |v| ·f before the integral. The remaining integral is then just the total macroscopic cross section Σ_t, i.e., the inverse local mean free path (dimension: 1/length). It is solely defined by total cross sections and independent of particle multiplication factors, since we only consider binary collisions (exactly two pre-collision partners always).

This term is then, often, taken on the left hand side of the Boltzmann equation with a positive sign, in the form:

δf

δt

|_{_loss} = Σ_t(r,v) |v|f(v)

(1.4)

For the linear case (f_b given) this "extinction coefficient" Σ_t is independent of the dependent variable f = f_i₀, and this term (out-scattering) just describes the loss of particle flux of i₀ particles due to any kind of interaction of them with the host medium.

With these formal substitutions the Boltzmann equation takes a form which is often more convenient, in particular in linear transport theory:

[

∂

∂t

+ v·∇_r +

F(r,v,t)

·∇_v] f(r,v,t) +Σ_t(r,v) |v|f(v)

⌠
⌡

d³v′ C(v′→ v)|v′|f(v′)

Q(r,v,t)

(1.1b)

Often the characteristic time constants for neutral transport phenomena are very short (μs), compared to those for plasma transport (ms). We can, therefore, often neglect explicit time dependence in the equations describing the neutral particles. This is done in most applications. The extension to time-dependent problems is rather straight forward and the procedure in the EIRENE code for such cases is described in 1.7 below.
A computationally crucial simplification is the neglect of all neutral - neutral interactions, rendering the equation for f(x) linear. For the status of options in EIRENE to deal with "self-collisions" see below, section 1.9.
For stationary (time-independent) problems the scalar transport flux ("angular flux") Φ, where

Φ(x) : = | v| ·f( r, v, i) ,

(1.5)

is sometimes used in preference to the distribution f(x) as dependent variable.

In particular for stationary (∂/∂t = 0) and force free (F = 0) problems, as e.g. often encountered in linear transport theory such as neutronics, neutral particle transport in plasmas, etc., the transport equation then reduces to the more compact form:

∇_r Φ(r,v,t) +Σ_t(r,v)Φ(r,v) = Q(r,v,t) +

⌠
⌡

d³v′ C(r;v′→ v)Φ(r,v′)

(1.1c)

Alternatively, in computational domains with non-vanishing collisionality (i.e., if Σ_t(x) ≠ 0 everywhere) the (pre-) collision density Ψ is used, i.e.,

Ψ(x) = Σ_t (x) ·Φ(x) = ν_t(x) ·f(x) , Σ_t = Σ_t(r,E,Ω,i) = ν_t(r,E,Ω,i)/|v|

(1.6)

where, again, the "macroscopic cross section" Σ_t is the total inverse local mean free path (dimension: 1/length), and ν_t is the collision frequency (dimension: 1/time). This cross section can be written as a sum Σ_t = ∑Σ_k over macroscopic cross sections for the different types (identified by the index k) of collision processes. Further details about this "macroscopic cross section" and its relation to the conventional ("microscopic") cross sections are given below, section 1.3.5.

1.2 The linear integral equation for the collision density Ψ

By formally integrating the characteristics for (1.1c) the same transport equation can also be written in integral form. This formal procedure is outlined below in paragraph 1.2.1.

The resulting integral equation is often most conveniently written for the collision density Ψ (1.2) rather than for transport flux Φ:

Ψ(x) = S(x) +

⌠
⌡

dx ′Ψ(x ′) ·K(x ′→ x) .

(1.1d)

This equation has the general form of the backward integral equation of a Markovian jump-process and it is therefore particularly well suited for a Monte Carlo method of solution. The formal relation between the integro-differential form (1.1c) and this integral form is very useful to generalize the Monte Carlo procedure, e.g., to time-dependent equations, and to Boltzmann-Fokker-Planck equations (which contain diffusive contributions or diffusive approximations for some processes, in addition to the jump processes described by the Boltzmann collision integral). It also allows to make connection to the so called "Greens-functions Monte Carlo" concept (originally developed for quantum mechanical problems involving solutions to the Schrödinger equation). A corresponding discussion is postponed to section 1.2.1. A direct intuitive interpretation of the integral equation is already sufficient to understand the Monte Carlo method of solution and shall be given first.

In (1.1d) x′ and x are the states at two successive collisions (jumps). The integral ∫dx′ in (1.1d) is to be understood as an integral over all initial coordinates, i.e. over the entire physical space, the full velocity space and a summation over all species indices. The transition kernel K is usually decomposed, in our context, into a collision- and a transport kernel, i.e., C and T , where

K(r′, v′, i ′→r, v, i) = C( r′;v′, i ′→ v, i) ·T(v, i; r′→ r) .

(2.7)

The kernel C is (excluding normalization) the conditional distribution for new co-ordinates (v,i) given that a particle of species i′ and with velocity v′ has undergone a collision at position r′. This kernel can further be decomposed into:

C( r′, v′, i ′→v, i) =

∑
k

p_k C_k ( r′; v′, i ′→ v, i) , p_k =

Σ_k

Σ_t

(2.8)

with summation over the index k for the different types of collision processes under consideration and p_k defined as the (conditional) probability for a collision to be of type k. The normalizing factor

c_k (x ′) =

∑
i

⌠
⌡

d³v C_k (r′, v′, i ′→v, i ) , C_k =

c_k

C_k

(2.9)

gives the mean number of secondaries for this collision process. The function C_k then is a conditional probability density. The particle absorption process can conveniently be described by adding an "absorbing state" x_a to the μ-space (generally referred to as "one-point compactification of this space" in the language of mathematical topology). This "limbo state", once it is reached, is never left again if the kernels T or C are employed as transition probabilities. Formally, an additional collision kernel C_a (x → x_a ) and an absorption probability p_a = Σ_a / Σ_t must be included in the collision kernel. The quantity Σ_a comprises all collision processes with no next generation particles.
The kernel T describes the motion of the test particles between the collision events. Let, again, Ω′ denote the unit vector in the direction of particle flight v′/|v′|, and let Ω₂′ and Ω₃′ be two further unit vectors such that these three vectors form an orthonormal basis at the point r′. Neither velocity nor species change along the transition described by T, i.e. v′=v and i′−i. Omitting the corresponding delta functions in velocity space and the Kronecker delta δ_i′i the transport kernel T then reads as follows:

(T(l)=) T ( v′, i′; r′→ r )

Σ_t ( r,v, i) ·exp

⎡
⎣

−

⌠
⌡

l=Ω′( r−r′)

0

ds Σ_t ( r′+ sΩ′, v′, i′)

⎤
⎦

·δ(Ω₂′( r−r′)) ·δ( Ω₃′(r− r′)) ·H (Ω′( r− r′))

(2.10)

Σ_t ( r,v, i) ·F( v′, i′; r′→ r ) 0 ≤ l ≤ ∞

(2.11)

with H(x) = 0 if x ≤ 0 , and H(x) = 1 if x > 0, the Heaviside step function. The two remaining delta functions restrict the motion to a path in the direction of the initial velocity v′.

Thus, although strictly being a conditional (on x′) distribution in phase space, for an infinite medium T can be interpreted as the distribution density for the distance l for a free flight starting from r′ to the next point of collision r = r′+ l ·Ω′.

For a finite medium this distribution can be generalized to (writing shorter l for (r′+lΩ′,v′,i′):

T(l)

⎧
⎪
⎨
⎪
⎩

Σ_t ( l)

·exp

⎡
⎣

−

⌠
⌡

l

0

ds Σ_t ( s )

⎤
⎦

l < l_max

δ(l−l_max)

·exp

⎡
⎣

−

⌠
⌡

l_max

0

ds Σ_t ( s )

⎤
⎦

l ≥ l_max

(2.12)

Here l_max denotes the distance along the flight direction from r′ to the boundary for the computational domain, to any internal surface at which the test flight shall be stopped, e.g. for scoring surface fluxes there, or even to the cell boundary. The "Transport Kernel" T has dimension: [1/Dimension of phase space] (T(x′→ x)dx is a probability). The function F (Dimension: [length times Dimension of T]) defined in expression (2.9b), and analogously from (2.10) for a finite medium, will turn out to be the relevant Green function for the transport problem, see Section 1.2.1. The integral

α(r′,r) =

⌠
⌡

Ω( r− r′)

0

ds Σ_t ( r′+ s Ω)

in equation (2.9a) is well known as "optical thickness of the medium" in linear transport theory.
The inhomogeneity S in equation (1.1d) is, excluding normalization, the distribution density of first collisions, whereas the integral term in equation (1.1d) describes the contribution to Ψ from all higher generations. The quantity S can be written as:

S(x) =

⌠
⌡

dx ′Q(x ′) ·T(x ′→ x) ,

(1.6a)

with a source density Q. As the problem is linear, Q can be normalized to 1 and, thus, Q can be considered a distribution density in μ -space for the "primary" birth points of particles, as, e.g., opposed to the "secondary" birth point distribution (or "post collision density") χ of particles after a collision event

χ(x) =

⌠
⌡

dx ′Ψ(x ′) ·C(x ′→ x).

(1.6b)

It can be shown that a unique solution Ψ(x) exists subject to appropriate boundary conditions and under only mild restrictions (basically on the constants c_k and p_a ) to ensure that the particle generation process stays sub-critical.
Usually, a detailed knowledge of Φ or Ψ is not required, but only a set of "responses", R, defined by

R = < Ψ| g_c > =

⌠
⌡

dx Ψ(x) ·g_c (x)

⎛
⎝

= < Φ| g_t > =

⌠
⌡

dx Φ(x) ·g_t (x)

⎞
⎠

(2.13)

where g_c (x), g_t (x) are given "detector functions".

For example all terms in the plasma fluid equations resulting from neutral plasma interaction can be written in this way. This can be seen by considering a numerical grid, composed of M mesh cells (spatial and/or temporal), for the numerical solution of the fluid equations. The detector functions for many responses needed for fusion plasma applications are hard wired in EIRENE, generalization to any arbitrary response by resorting to "user defined detector functions" is described in Section 3.2.1 for tracklength estimates, and in Section 3.2.2 for collision estimates.

Lets therefore define an entire set of detector functions g_m, one for each mesh cell of an external code, each including a characteristic function

g_m=g ×ch_m(r,t) , m=1,2,...,M,

(2.14)

i.e., ch_m (r,t)=1 inside the numerical mesh cell (or time interval) labeled with the cell index m, and ch_m (r,t) = 0 outside this cell. Thus profiles of cell volume averaged responses are readily obtained in a single Monte Carlo run.

Estimates of surface fluxes, or point estimates e.g. in time, are also included in this concept if proper use is made of delta functions to reduce dimensionality of the response:

g_m,α=g ×ch_m(r,t) ×δ^α(r,t) , m=1,2,...,M,

(2.15)

Here, e.g., a surface cell m would discretise a surface S characterized by a surface delta function δ^S.

Monte Carlo estimates of the volumetric source terms in the fluid equations due to the trace particles (neutral particles, but also trace ions described by the Monte Carlo procedure) are therefore cell averages, surface averages, time averages or point estimates (e.g. in time), averaged over a sub-manifold with reduced dimensionality, however, only if the probability for particle histories crossing this sub-manifold is not zero. Surface averages in EIRENE are described in Section 3.2.4, point averages in time are described in Section 3.2.3. This option to reduce dimensionality of responses usually excludes point estimates in real space, for which special estimators not mentioned here would be required.

Nevertheless, in concluding, depending upon the numerical algorithm in the fluid code, one may then have to interpret these Monte Carlo estimates properly, e.g. they may have to be interpolated to the grid points, or be properly rescaled in case of time-dependent applications.

1.2.1 The Green function concept

We return to Equation (2.9b), where the function F(v,i;r′→ r)) was defined. For a finite domain, see (2.10), and introducing again the shortcut l for (r′+lΩ,v,i) this becomes:

F(l) =

⎧
⎪
⎨
⎪
⎩

exp

⎡
⎣

−

⌠
⌡

l

0

ds Σ_t ( s )

⎤
⎦

l < l_max

l ≥ l_max

(2.16)

to be written

1.3 Monte Carlo solution of equation 1

A statistical solution to equation (1.1b) is straight forward, because it is formulated in probabilistic terms as follows.

A discrete Markoff chain is defined using Q as an initial distribution and L = T ·C (order of C and T reversed compared to K) as a transition probability. Histories ωⁿ from this stochastic process are generated according to ωⁿ = (x₀, x₁, x₂, ... , x_n), (where x_j = x_a for all j ≥ n and x_i ≠ x_a for all i < n), with x_n being the first state after transition into the absorbing state x_a. x₀ denotes the initial state distributed as described by Q. Thus, the length n of the chain ωⁿ is itself a random variable. A random sampling procedure to generate such chains is carried out in Monte Carlo codes by converting machine generated (pseudo-) random numbers ξ_i₁ , ξ_i₂, ... into random numbers with the distributions Q, T and C. Having computed N (several thousand) chains ω_i, i=1,2,...,N, the responses R are estimated as the arithmetic mean of functions ("statistics", or "estimators") X(ω), i.e.

R ≈

N
∑
i=1

X(ω_i) .

(3.17)

1.3.1 Unbiased estimators

One possible choice for X(ω) is the so called "collision estimator" X^c,

X^c_g ( ωⁿ_i) =

n
∑
l=1

g_c (x_l ) ·

l−1
∏
j=1

c(x_j )

(1 − p_a(x_j))

(3.18)

This estimator is, for example, used in the DEGAS code [10]. It can be shown that the statistical expectation E(X^c) produces:

R = E(X^c ) =

⌠
⌡

d( ω) X^c_g ( ω) h( ω)

(3.19)

with h(ω) being the probability density for finding a chain ω from the Markoff process defined above, i.e. X^c is an unbiased estimator for R.
Other estimators ("track-length type estimators") are employed frequently. These estimators are unbiased as well but have higher moments different from those of X^c. Instead of evaluating the detector function g_c(x) at the points of collisions, x_l , (as X^c does,) they involve line integrals of g_t(x) along the trajectories, e.g.,

X^t_g ( ωⁿ_i ) =

n−1
∑
l=0

⎧
⎨
⎩

⌠
⌡

x_l+1

x_l

ds g_t (s)

⎫
⎬
⎭

l−1
∏
j=1

c(x_j )

(1 − p_a(x_j))

(3.20)

with R = E(X^t ) = E(X^c). It has been shown, e.g., in [7], that the collision estimator, derived not for the pre- but for the post-collision density integral equation results in a track-length type "conditional expectation estimator" X^e, which, together with X^c and the "traditional" track-length estimator X^t, may be used as one further option in the EIRENE code.
This estimator X^e is obtained from X^t by extending the line integration, which is restricted to the path from x_l to x_l+1 in formula (3.17), to the line segment from x_l to x_end. Here x_end is the nearest point on a boundary along the test flight originating in x_l. I.e., the line integration may be extended into a region beyond the next point of collision, into which the generated history would not necessarily reach. This "conditional expectation estimator" reads:

X^e_g ( ωⁿ_i ) =

n−1
∑
l=0

⎧
⎨
⎩

⌠
⌡

x_end

x_l

ds g_t (s) ·exp

⎛
⎝

−

⌠
⌡

s

0

ds ′Σ_t ( s ′)

⎞
⎠

⎫
⎬
⎭

l−1
∏
j=1

c(x_j )

(1 − p_a(x_j))

(3.21)

If x_end is taken to be the nearest point on each mesh cell boundary, then the estimator (3.18) reduces, as a special case, to the method employed by the NIMBUS code [11]. However, the estimator (3.18) is more general, as the integration may be extended over arbitrarily many cells. The length of the integration path is controlled in EIRENE by an input flag WMINC, see section 2.10.

1.3.2 Scaling of tallies

With increasing N, the number of Monte Carlo histories, the unbiased estimators X given in previous section 1.3.1 provide arbitrarily precise approximations for the responses R_g=〈Ψ|g〉, for detector function g and dependent variable Ψ (equation 1.1d). The precise meaning and physical dimensions of g and Ψ depend upon the problem at hand, e.g. also on symmetry (ignorable spatial coordinates, or time).

Because the inhomogeneous part of the transport equation (1.1d) S was assumed to be normalized to one (for source sampling),

S (r,v,t)=

S(r,v,t), s =

⌠
⌡

d³rd³vdt S(r,v,t)

(3.22)

this normalization factor s has to be multiplied to estimators (3.14) to turn responses R_g (2.11) into estimates with correct dimensionality and in absolute units.

Mostly the responses of interest are intensive quantities, such as density, pressure, collision density (number of collision per unit time and volume) rather than the extensive quantities obtained by the phase space integration in (2.11), such as volume, total energy or number of particles. Thus by dividing the estimate by the appropriate cell-volume V_m (in space-time), see (2.12), finally the profiles of cell averaged intensive quantities (e.g. density profiles) are obtained. The final unbiased estimate, in absolute units, for any of the unbiased estimators X_g for detector function g discussed above, then is:

R_g = 〈Ψ|g〉 ≈

(N) =

V_m×N

N
∑
i=1

X_g(ω_i) , N >> 1

(3.23)

with N being the number of Monte Carlo histories. Correct interpretation of the results of an EIRENE run hence requires knowledge of the ratio s/V_m. The total source strength scaling factor s, i.e. the integral of the inhomogeneous part S of the governing Fredholm integral equation, is input to EIRENE, variable "FLUX" in input block 7 (see 2.7). The space-time cell volumes V_m are either automatically calculated by EIRENE in subroutine "VOLUME", for cells belonging to standard grid (input block 2, Section 2.2), or are provided from external considerations by a number of options, e.g. input blocks 3b, 5, or 8. In this latter case great care is needed, however, that the cell volumes are proportional to the volumes as seen by the test flights for otherwise not only unscaled profiles, but even wrong (biased) profile shapes will be obtained.

Scaling in problems with ignorable coordinates

As seen in the previous paragraph, the absolute values of estimates are scaled with the ratio s/V of source strength to cell volume. Depending on ignorable coordinates in any particular problem, or on whether stationary (i.e.: time is an ignorable coordinate) or time dependent transport problems are considered , the dependent variable Ψ in the governing integral equation (1.1d), the source strength s and "cell volume" V may have different interpretations:

stationary (time independent) problems

one ignorable spatial coordinate In stationary problems with one ignorable coordinate, say, the z-coordinate, the source strength s (input flag "FLUX") is the flux, particles per unit length dz in direction of z. Likewise, the volume V is per unit (mostly "toroidal") length, i.e. if dz=1 then V is the cell area in the two remaining coordinates x,y. In

1.3.3 Efficiency (FOM)

The efficiency for Monte Carlo Codes is the inverse of the figure of merit (FOM) of the calculation, defined as

FOM = statistical variance ·computing cost.

(3.24)

Note that FOM should be approximately independent of the running time, because the number of histories generated is (excluding overhead) proportional to the CPU costs and inversely proportional to the statistical variance σ².

It is one of the major advantages of Monte Carlo methods over other numerical schemes that the error estimates, empirical variances [(σ)\tilde]², are directly provided by the method itself, not requiring any further considerations. The options to activate evaluation of statistical variance in an EIRENE run for any computed quantity (tally) are described in input block 9 (section 2.9). EIRENE provides numerical and graphical output for the "empirical relative standard deviation", in %.

g,rel

(N) =

(N)/

(N) ×100

(3.25)

and [(R)\tilde]_g^N is the Monte Carlo estimate for tally R_g based on N Monte Carlo histories (3.20). Here R_g may be any volume or surface averaged tally for detector function g, either tracklength- collision- or snapshot estimated. The variance σ²₁ per history is obtained as the (unbissed) estimate ("empirical variance") as:

σ²₁ ≈

2
1

(N) =

N−1

N
∑
i=1

(X_i−

)²=

N−1

⎡
⎣

N
∑
i=1

X_i² −

⎛
⎝

N
∑
i=1

X_i

⎞
⎠

⎤
⎦

(3.26)

where X_i = X_g(ω_i) is the contribution of Monte Carlo history ω_i to estimator X_g for tally g, and [ˉ(X)] is the arithmetic mean [ˉ(X)]=1/N∑X_i over all histories. N is the number of (statistically independent) Monte Carlo histories. The subscript g is omitted here and from now on. The variance per history is turned into the final estimate for the variance by the "law of large numbers" (and the "central limit theorem" of probablity theory):

(N)=

2
1

(N)

(3.27)

For large N the variance per history [(σ)\tilde]₁² converges to a constant (namely to σ²₁), and the final Monte Carlo estimate of variance [(σ)\tilde]² (3.24) has the expected probabilistic 1/N scaling, i.e. the empirical standard deviation scales as: [(σ)\tilde] ∼ 1/√N.

Multiple cell crossings One should note that except in simple 1D cases in general one history ω_i can contribute more than once to the estimate for a particular cell of the grid. E.g. test particles can cross one cell m (see (2.12b) more than once, with each cell crossing "j" leaving a score X_i,j. Then: X_i=∑_j X_i,j. The summation in (3.14) and hence also the second sum on right hand side of (3.23) can still be accumulated "on the fly" while generating the Monte Carlo histories.

Sampling, Non-analog sampling

One often encounters in literature the description of special, "new" Monte Carlo techniques, that are greatly superior to so called "standard methods".

It is true that one can devise very intelligent methods to optimize performance, but each optimization always only works well for one very particular problem, it can equally well entirely wreck the performance of an only slightly different case.

For a general purpose Monte Carlo solver for transfer problems, as EIRENE, therefore, no general advise can be given, although the performance often could greatly be improved by adapting the method to a particular problem. In particular, non of these "intelligent methods" have been (and will ever be) hard wired into EIRENE.

EIRENE contains a set of such methods, referred to as "non-analog" methods and controlled by the flags in input block 9, see section 2.9. Activating those must be accompanied by a very careful statistical analysis of results, not just the run time per particle, nor even the standard deviation alone, suffices.

We first note here that any random distribution function f(x) arising during the course of generating the chains ω ("particle histories") can be replaced by another "non-analog" function, g(x), if a weighting factor,

w(y) = f(y) / g(y)

(3.28)

is included in the formula for the estimator X(ω), y being an actual random number generated from g . This choice can in some cases increase the efficiency of the algorithm.
The only restrictions on the choice of non-analog functions (besides practical ones) are:

if g(x) = 0, then f(x) = 0;

(3.29)

and conditions to ensure that the non-analog process remains sub-critical as well.

The condition 3.26 is checked in an EIRENE run whenever non-analog distributions are applied, and, if violated, an error exit with the message: "violation of Radon-Nikodym condition" occurs.

Note: a violation of condition 3.26 can not be detected otherwise, e.g., by monitoring overflow in the weighting factor w. Of course, values y resulting in a zero in the denominator of w (equation 3.25) have sampling probability zero according to distribution g(x), and, hence, are never sampled. Still the results would be biased.

Note: when inspecting an EIRENE geometry plot with particle trajectories (or even a movie: NLMOVIE=TRUE, input block 1) in order to get an intuitive feeling for the particular transfer process considered, then all non-analog sampling options must first be turned off, (NLANA flag in input block 1), for otherwise the pictures may be grossly misleading.

We consider the procedures for random sampling from univariate distributions as known, and refer to the many textbooks on that, in particular to the "random sampling library" [12]. If none of the direct methods apply, then still either the non-analog method mentioned above, or the "rejection method" can be used.

Sampling from a multivariate distribution f(x₁,x₂,...,x_n) can always be reduced to a sequence of samplings from univariate distributions, by noting that:

f(x₁,x₂,...,x_n)=f₁(x₁) ·f₂(x₂|x₁) ·f₃(x₃|x₁,x₂) ·...

(3.30)

Here, f₁ is the marginal distribution obtained from f by integrating over all but the first independent variables. It is a univariate distribution of x₁.

f₂ is a conditional marginal distribution obtained from f firstly by integrating f over all but the first 2 variables x₁, x₂, and secondly then taking the conditional distribution, conditional on x₁ (i.e., the univariate distribution of x₂ for given values of x₁).

Likewise, f₃ is a conditional marginal distribution, n-3 fold integration of f, and then distribution conditional on x₁, x₂, and so on.

Random sampling from f(x₁,x₂,...,x₃) then proceeds by sampling first x₁ from the (univariate) first factor in equation 3.27, then x₂ by sampling from the (univariate) second factor, and so on.

One particular example of this scheme is the "TRIM database surface reflection model" (see below, section 1.4.2). There, essentially, tables of the conditional marginal distributions mentioned here are pre-computed with the TRIM code. For convenient random sampling, these tables are stored for the inverted cumulative distribution, i.e., as conditional quantile functions.

Stratified Source Sampling

The EIRENE code resorts to a "stratified sampling" technique. This means that the primary source distribution (see Equation 1.1a) can be decomposed into a sum of independent sources:

Q =

∑
i

Q_i

and the solution is obtained by linear superposition of the solutions for each sub-source ("stratum") Q_i.

The stratified sampling technique is a well known statistical procedure, see any textbook on statistics or Monte-Carlo methods.

It sometimes can significantly affect the efficiency of a Monte-Carlo run. This can go both ways, depending upon the particular problem and the particular stratification used. Therefore, as with the non-analog options, no general recommendation can be made.

The effects of stratification, however, are usually more easy to assess (predict) than with general non-analog methods.

In order to make connection with the textbooks we need to note here only the following:

EIRENE provides estimates of linear functionals ("moments", "responses")

< g,φ_Q > = < Q,φ^*_g >

(3.31)

As described above, g is the detector function ("weighting function") φ_Q is the solution to the kinetic transfer equation (Boltzmann eq.) in which Q was used as (inhomogeneous) source term. φ^*_g is the solution of the corresponding adjoint equation, with g as source term.

The above equation means that EIRENE estimates can simply be regarded as Monte-Carlo estimates of multi-dimensional integrals, namely of the adjoint function φ^*_g integrated over phase space with Q defining a distribution (measure) in phase space "μ_Q".

It is then obvious that one can split that integral into a sum of integrals by arbitrarily decomposing the domain of integration, i.e., by decomposing the source Q into sub-sources ("strata") Q_i with

R_g =

⌠
⌡

dxφ_Q ·g = < Q,φ^*_g > =

∑
i

< Q_i,φ^*_g > =

∑
i

⌠
⌡

dμ_{Q_i} φ^*_g

(3.32)

This is the concept of stratified sampling, which, for Monte Carlo particle transfer procedures in particular actually turns out to be a "stratified source sampling".

Further details about this in our particular context are given in the description of input block 7, section 2.7 below.

1.3.4 Source sampling

to be written:

and see special section on "plasma surface recycling source sampling".

1.3.5 Sampling from the transport kernel T

The conditional sampling distribution T(r′→ r) for the next point of collision (event), given a test flight is at position r, traveling with velocity v, is given by (2.9a). Let us introduce the distance of the flight l in direction Ω =v/|v|, hence r =r′+lΩ, then the sampling distribution for l is simply

T(l)=Σ_t(r′+lΩ)exp

⎡
⎣

−

⌠
⌡

l

0

ds Σ_t(r′+sΩ)

⎤
⎦

, Σ_t=1/λ_t

(3.33)

Distances l are hence readily sampled from this exponential distribution using the inversion method and the cumulative distribution G_T of T

ξ = G_T(l)=1−F(l)=1−exp

⎡
⎣

−

⌠
⌡

l

0

ds Σ_t(r′+sΩ)

⎤
⎦

(3.34)

with ξ (and hence also 1−ξ) being a uniformly distributed random variable. Therefore, cumulating exp[−∫₀^l ds Σ_t(r′+sΩ)] along a flight and comparing it with ln(ξ) with a uniform random number ξ provides a random sample of the free flight distance l.

Sampling the distance l of a free flight, for test particle species "i" with velocity v, from the transport kernel T requires knowledge of the ("monochromatic") mean free path λ = λ(r,|v|) along the trajectory, and, hence, of the "macroscopic cross section" Σ = Σ(r,|v|), which also may vary along the trajectory.

The "macroscopic cross sections" Σ (by abuse of language, dimension: 1/length) with or without subscript k (the label for the type of the collision process) or subscript t (for "total", the sum over k), are defined as:

Σ =

|v|

, λ =

|v|

(3.35)

with the mean free path λ and the collision frequency ν:

ν(r,v,i) = n_b 〈σ·v_rel〉_b, v_rel = |v−v_b|

(3.36)

Here n_b and v_b are the target particle density and velocity, respectively, and the brackets 〈...〉_b denote averaging with the velocity distribution f_b(r,v_b) of the background medium species b.

σ is the total scattering cross section, σ = σ(|v_rel|).

Hence, e.g., for all isotropic target distributions f_b one has:

ν = ν(r,E,\not Ω,i) and Σ = Σ(r,E,\notΩ,i).

Furthermore, let f_b be a drifting Maxwellian, with temperature T_b and drift V_b: f_b = f_b(T_b,V_b).

Then:

ν = ν(T_b,|v−V_b|,i)

Finally, if the thermal speed v_b^therm in this Maxwellian: v_b^therm = √{T_b/m_b} is very large compared to typical values of |v−V_b| (i.e.: nearly incompressible flow conditions), then:

ν ≈ ν(T_b,0,i) = n_b 〈σ·v〉_b

and here 〈σ·v〉_b is the usual "Maxwellian rate coefficient" for the collision process, taken at temperature T_b of the background particles (species b) and at zero velocity of the test particles (species i). For example, in case of neutral atom or molecule collisions with electrons this approximation is usually valid.

1.3.6 Sampling from the collision kernel C

The general rules for sampling from multivariate distributions outlined in section 1.3.3, equation 3.27 also apply here.

It is not necessary, but usually most convenient and intuitively most clear, to decompose the kernel C according to the physical nature of the various collision processes taken into account.

The collision kernel C was written as (see equation 2.7)

C( r′, v′, i ′→v, i) =

∑
k

p_k C_k ( r′; v′, i ′→ v, i) , p_k =

Σ_k

Σ_t

with summation over the index k for the different types of collision processes under consideration and p_k defined as the probability for a collision to be of type k. The normalizing factor

c_k (x ′) =

∑
i

⌠
⌡

d v C_k (r′, v′, i ′→v, i ) , C_k =

c_k

C_k

was the mean number of secondaries for this collision process. The function C_k then is a conditional probability density for (v,i), given the pre-collision coordinates (r′,v′,i′).

Sampling from C (given that a collision of some kind has occurred) proceeds most conveniently by firstly sampling the type k of the process from the discrete distribution

(p_k, k=1,..,K),

and then by sampling the post-collision state of the test-particle from C_k. Finally the weight of the test-flight after collision is increased (or reduced) by multiplying it with the normalization constant c_k.

The decomposition of C into sub-kernels C_k is somewhat arbitrary, and the criteria may also be computational aspects rather than the different physical nature of the various processes taken into consideration.

Sampling (v,i) from C_k can be done along the same idea: first factoring out the discrete distribution for the species index i, and then by sampling v from C_k,i, the conditional distribution for v given the type of the process is k and the species of the post collision trajectory is i. Let (n_k,i′,i′=1,...,I_k) be the number of post collision particles of species i′ from process k, and n_k = ∑_i′ n_k,i′. Then

(

k,i′

i′=1,...I_k) = (n_k,i′/n_k, i′=1,...I_k)

is the discrete distribution for the post-collision particle species i′, given a collision process of type k.

1.3.7 elastic collisions

1.3.8 charge exchange

1.3.9 electron-impact collisions

1.3.10 general heavy-particle-impact collisions

1.3.11 photon processes (emission, absorption, scattering

1.4 Surface Reflection Models

Within the terminology fixed in section 1.2 the effects of the interaction of neutral particles (and escaping ions) with surfaces surrounding or inside the computational volume can be described by additional boundary conditions. It is, however, much more convenient for Monte Carlo applications to describe this interaction as just one special type of collision process in the collision kernel C of equation (2.7), i.e. as

C = p_w ·C_w +

∑
k

p_k C_k ,

(4.37)

with p_w (r) = 1, if r is a point on a surface, and p_w (r) = 0 elsewhere. The transport kernel T can be modified such that the maximum length l_max for free flight to the next point of collision is the distance to the nearest surface along the track of the particle. The reflection kernel C_w is further decomposed it into a kernel C_wf (fast particle reflection), C_wt (thermal particle re-emission) and C_wa (particle absorption) with respective probabilities, p_f , p_t and p_a, such that

C_w = C_wf + C_wt + C_wa = p_f C_wf + p_t C_wt + p_aC_wa , p_f + p_t + p_a = 1 .

(4.38)

Here C are, as above, the normalized versions of the collision kernels C.

We discuss in detail only the fast particle reflection model C_wf since C_wt is generally a mono-energetic and cosine angular distribution or a Maxwellian flux distribution at wall temperature, and C_wa only describes the transition into the "limbo state" x_a "absorbed particle".
It is usual in neutral gas models to replace C_wf with numerical or analytical fits to moments of C_wf, e.g., such as the particle and energy reflection coefficients C_P and C_E. The particle - wall interaction is often supplemented by simple assumptions on the angular distribution (cosine or specular) of the particles leaving the surface.

One has for the particle reflection coefficient:

C_P =

∑
i

⌠
⌡

d v C_wf (r; v′, i ′→ v ,i)

(4.39)

Note that C_P = p_f = c_wf in the terminology of equation 2.8.

The energy reflection coefficient is defined as:

C_E =

E ′

∑
i

⌠
⌡

d v E ·C_wf (r; v′, i ′→ v,i) =

E ′

(4.40)

Here E = m v² / 2 is the energy of the reflected particle. Hence C_E = [ˉ(E)]/E′·C_P, with [ˉ(E)] denoting the mean energy of the reflected particles.

1.4.1 The Behrisch matrix reflection model

One model, the matrix model due to Behrisch [13], is hardwired in EIRENE (and in many other codes such as AURORA, BALDUR, TRANSP). This model has been used as a standard option in many benchmarks [14]. The model is based upon a table of values for C_{P_⊥} (E ′) and a stochastic matrix B_⊥ for the transition E ′→ E, for H atoms striking a stainless steel target at normal incidence. For other incident angles ϑ_in the moments are modified in EIRENE by the relations:

C_P ( ϑ_in ) = 1 − ( 1 − C_{P_⊥}) ·cos^e1 ( ϑ_in )

(4.41)

and

C_E ( ϑ_in ) = 1 − ( 1 − C_{E_⊥}) ·cos^e2 ( ϑ_in ) .

(4.42)

This latter moment equation can be satisfied by modifying the tabulated stochastic matrix B_⊥( E ′→ E) in many different ways. In EIRENE, we first sample E from B_⊥ (given E ′) and then set

E = E ′− ( E ′− E ) ·cos^e2 ( ϑ_in ) .

(4.43)

The linearity of E in E gives an expectation value [ˉ(E)] with [ˉ(E)] / E ′·C_P( ϑ_in ) = C_E (ϑ_in ), with C_P and C_E as given by equations (4.38) and (4.39), respectively.

To describe the angular distribution of re-emitted particles, let e_x , e_y , e_z define an orthonormal basis at the strike-point of a test-particle on a surface. Here let e_x be the unit vector parallel to the outer surface normal and let the incident particle travel in the local xz-plane. Furthermore, let ϑ_r and ϕ_r be the polar and azimuthal angles sampled from a cosine distribution around − e_x. The speed unit vector Ω of the reflected particle in EIRENE is then given by

Ω_x

− cos(

) f₁ ( ϑ_in) + sin(

)sin( ϕ_r) f₂ ( ϑ_in ) ,

Ω_y

sin(

)cos( ϕ_r ) ,

Ω_z

sin(

) sin( ϕ_r )f₁( ϑ_in ) + cos(

) f₂(ϑ_in ) ,

(4.44)

with

sin(

)

f₁( ϑ_in) sin( ϑ_r ) ,

f₁( ϑ_in )

cos^e3 ( ϑ_in ) ,

f₂( ϑ_in )

( 1 − f₁ ·f₁ )^[1/2] .

In total we have defined three parameters e1, e2, e3. The first two parameters e1 and e2 account for particle and energy reflection coefficients that continuously increase with incident angle (recommended: e1=1, e2=0.5 from comparison with TRIM code calculations, see below).
The third parameter e3 can be used to take into account an increasing contribution of specular reflection at near grazing incidence. e3=0 gives a pure cosine distribution, e3=1 reproduces the distribution given in [10]. In general this model describes a fraction of specular versus cosine reflection which increases with ϑ_in, and in the (purely mathematical) limit of ϑ_in = 90^°, pure specular reflection. Furthermore, the parameter e3 controls the speed at which the specular part a_spec (Ω) of the angular distribution a(Ω) = a_spec ( Ω) + a_cosin ( Ω) increases with ϑ_in .
The interaction of neutrals and ions with a solid surface for target-projectile combinations other than H onto stainless steel are modelled using a reduced energy scaling.

1.4.2 TRIM code database reflection models

Complementary to the surface model described above, EIRENE is coupled to computer generated reflection databases, produced with the TRIM code (see references [15] and [16]). At present a discretized form of the reflection kernel C_wf, including all correlations between the vector components of v′ and v, is available for the target-projectile combinations listed in table 1.1. ZNML is a target material flag set for each surface in the input file, blocks 3a and 3b, and is described below in block 6b, "Data for local reflection models". These databases are, essentially, tables of conditional quantile functions obtained from large random samples of test particle trajectories in the solid.

Table 1.1: Surface reflection database (TRIM Code)

No	Projectile	Target	ZNML

1	H	Fe	5626.
2	D	Fe	5626.
3	H	C	1206.
4	D	C	1206.
5	He	Fe	5626.
6	He	C	1206.
7	T	Fe	5626.
8	T	C	1206.
9	D	W	18474.
10	He	W	18474.
11	H	W	18474.
12	T	W	18474.
13	D	Be	904.
14	D	Mo	9642.
15	Ne	C	1206.
16	Ne	Be	904.
17	H	Cu	6429.
18	H	Mo	9642.
19	T	Mo	9642.
20	He	Mo	9642.

We consider these "database -" reflection models to be the most complete option available at present and recommend to further extend the databases to additional projectile - target combinations, and in particular to sputtering and self-sputtering data (such as, e.g., C onto C, or Be, etc.).
For an INTOR benchmark case for neutral particle Monte Carlo Codes [14] the choice of a particular surface reflection model had no strong impact on the results. This was shown by running the EIRENE code with all three models ("Behrisch Matrix", TRIM database and MARLOWE database (at that time available as a third option. Meanwhile the MARLOWE database is not available anymore).

It was found that this had little influence on global parameters as, e.g., the pumping efficiency for the particular geometry (5.4 %, 5.3 % and 5.2 % respectively.). On the other hand, it can easily be imagined that in other cases (different plasma conditions, more detailed geometries) the choice of the reflection kernel C_wf can be more important for the result. Certainly in all cases, in which the results are sensitive to the neutral particle velocity distribution near surfaces, details of the reflection kernel matter. This is, for example, the case for computed line shapes of certain electronic transitions (H_α lines etc.), see e.g., reference [17].

1.5 Surface sources

It has been found as a result of benchmark calculations with several different neutral gas transport Monte Carlo codes on the same input data that the main divergence in the results from the different codes tested was due to different primary source terms Q in the transport equation [equations (1.1a,1.1b) and (1.6a)] , although target fluxes and temperatures had been prescribed to be identical. Here we describe one particular EIRENE surface source option in detail. The way the source is modelled by EIRENE allows consistency with kinetic boundary conditions for plasma fluid equations at plasma facing surfaces. To describe this, we first write the volumetric source Q (particles per unit time and volume) in terms of a surface (wall) flux Γ_w and surface delta-function:

Q( r, v, i ) = Γ_w (r, v, i ) ·δ( π) ,

(5.45)

where π is a coordinate normal to the surface S_w: π(r)=0 , i.e.: e_π = ∇_rπ(r), pointing away from the computational volume. Γ_w is a surface source (flux-) density distribution, backward directed (into the computational volume, i.e., only v_π < 0).

A typical example is a sonic parallel plasma flow (Mach number M_||=1 along the magnetic field lines) entering the magnetic presheath. In this magnetic presheath the ion trajectories are bend over such that at the entrance of the electrostatic sheath the sound speed is already reached with respect to the surface normal: M_π=1.

To convert from a forward directed ion flux distribution Γ_Sh,ion at the plasma-sheath interface S_Sh to a neutral flux distribution back from the wall surface S_w, we fold Γ_Sh,ion with a (magnetic plus electrostatic) sheath transmission kernel T_Sh and with the surface reflection kernel C_w, equation 4.34:

Γ_w ( r, v, i )

∑
i ′

⌠
⌡

d³ r′d³ v′d³v" Γ_Sh,ion ( r′,v" , i ′)

·T_Sh (i ′; r′,v" → r, v′) ·C_w ( r;v′, i ′→ v, i )

(5.46)

The summation over i ′ is performed over all ion species entering the sheath region.

A typical example for Γ_Sh,ion is a sonic parallel plasma flow (Mach number M_||=1 along the magnetic field lines) entering the magnetic pre-sheath. In this magnetic pre-sheath the ion trajectories are bend over such that at the entrance of the electrostatic sheath the sound speed is already reached with respect to the surface normal: M_π=1.

A commonly adopted simplification is made by ignoring the typically small spatial volume between the plasma sheath interface and the wall surface. Hence: T_Sh factorizes into a 3D delta function δ³(r−r′) and a deterministic transition kernel for the velocities v" → v′.

Random numbers from Γ_w are generated conveniently by sampling, firstly, (r′, v" , i ′) from Γ_Sh,ion, next generating (r,v′) from T_Sh and then, finally, (v, i) from the conditional distribution C_w .

A frequently employed assumption for Γ_Sh,ion is a forward drifting Maxwellian ion flux distribution at the interface between plasma and electrostatic sheath S_Sh. In this case the bending of trajectories from parallel to normal (to the surface) sonic flow in the magnetic pre-sheath is already included in Γ_Sh,ion , e.g. via boundary conditions in a CFD model for plasma ion flow, and the action of operator T_Sh is solely the acceleration in the electrostatic sheath. Such a forward drifting Maxwelllian ion flux distribution is given as

Γ⁰_Sh,ion =

α²

·v_π · exp

⎛
⎝

−

2 β²

(v− V_d )²

⎞
⎠

, v_π > 0

(5.47)

with β = √{ [(kT)/(m)] }, normalizing constant α² and with the drift velocity V_d. The superscript 0 indicates normalization of the flux to one.

Consider the three orthogonal velocity components v₁,v₂,v_π, then Γ⁰_Sh,ion(v) factorizes in the three uni-variate distributions

Γ⁰_Sh,ion(v)=Γ⁰_Sh,ion,1(v₁)·Γ⁰_Sh,ion,2(v₂) ·Γ⁰_Sh,ion,π(v_π)

(5.48)

The first two are simply shifted Gaussians. However, sampling the π-component v_π of v from Γ_Sh,ion,π is not trivial because the cumulative distribution function F of a 1D forward drifting Maxwellian (with cut-off at v_π=0):

F ( v_π ) =

⌠
⌡

v_π

0

dv_π ′ Γ⁰_Sh,ion ( v_π ′)

(5.49)

cannot be inverted explicitly. Instead we sample this component [(v)\tilde] _π from a "non-analog" truncated Maxwellian flux distribution F (see section 1.3 above) with zero drift in π-direction but with different temperature [(T)\tilde]_π ,([(T)\tilde]_π ≠ T , [(V)\tilde] _d,π = 0 ). This can be done by setting, e.g.,

√

−ln(ξ) ·2

, ( note:

since

d,π

= 0 ),

(5.50)

ξ being a uniformly distributed random number between 0 and 1.

In order to find the weight correction factor associated with this procedure, we first note that the normalizing factor α_π² for the correct v_π distribution reads:

α_π² = β² ·exp(−M²_π ) g(M_π ) ,

(5.51)

with M_π = V_d,π/ (√{ 2 } ·β) and g(x) = 1 + √{π} x (1 + erf(x)) ·exp(x²).

The conditions

= 0 ,

then give the following result for the weight correction factor w ( [(v)\tilde]_π ) to be applied for each sampled non-analog velocity component [(v)\tilde]_π:

w (

) =

Γ⁰_Sh,ion (

)

0
Sh,ion

(

)

β²

exp

⎛
⎜
⎝

V_{d, π}

β²

−

2
π

2 β²

⎛
⎝

1 −

β²

[(β²)\tilde]

⎞
⎠

⎞
⎟
⎠

g(M_π)

(5.52)

This non-analog sampling of the source term introduces additional statistical fluctuations in the results, since the initial test particle weights fluctuate with variance

σ² =

⌠
⌡

∞

0

w² (

) ·

0
Sh,ion

(

) − 1

⎛
⎝

2M_π

√b

⎞
⎠

/g(M_π)² ·

(2−b) b

− 1

(5.53)

where b = 2 − β² / [(β²)\tilde]. We use the ansatz b = C ·M_π + 1 and compute numerically that value of C which minimizes σ² in the range 0 ≤ M_π ≤ 1. We find C = 0.6026. This means that the ratio of the "real" temperature T and the non-analog temperature [(T)\tilde] should be chosen to be T / [(T)\tilde] = 1 − 0.6026 ·M_π. This choice optimizes the statistical performance. In principle, of course, any other [(T)\tilde] could have been used without introducing bias into the algorithm.
The mean energy of ions generated by this sampling procedure is (note: Γ⁰ is the distribution already normalized to flux one):

(M_π)

m_i

⌠
⌡

∞

−∞

⌠
⌡

∞

−∞

⌠
⌡

∞

0

d³v′ v′² ·Γ⁰_Sh,ion (v′) =

m_i

⌠
⌡

∞

−∞

dv₁

⌠
⌡

∞

−∞

dv₂

⌠
⌡

∞

0

dv_π (v₁²+v₂²+v_π²) ·Γ⁰_Sh,ion (v′)

(5.54)

T ·γ_E ,

with the half sided integration for the π-component of v′, and:

γ_E = ( M_σ² + 2 ) +

⎡
⎣

M_π² + 0.5·

g(M_π ) − 1

g(M_π )

⎤
⎦

(5.55)

We have decomposed the drift velocity V_d into a normal (to the surface) π-component V_d,π and a σ-component V_d,σ parallel to the wall surface.

M_σ =

| V_d − V_d,π |

√2 ·β

| V_d,σ |

√2 ·β

(5.56)

Note that for M_π = 0 this reduces to the expected formula γ_E = 2 + M_σ². Furthermore, if M_σ = 0 (normal incidence) and M_π → 1 (Bohm sheath condition for T_e = T_i), then γ_E for our forward sided fluxes approaches the value for total Maxwellian heat fluxes: γ_E = 2.449 + M_π² ≅ [5/2] +M_π², as it must.

We conclude that the above scheme is both an accurate and efficient random number generator for truncated shifted Maxwellian flux distributions of particles incident onto or emitted from surfaces. Of course, after having generated [(v)\tilde] with the weight w ([(v)\tilde] _π) and before reflecting the particle via the kernel C_w , the acceleration of an ion in a sheath potential in direction π (normal to the target surface), as formalized in the sheath transmission kernel T_Sh mentioned above in Eq. (5.43), can easily be included.

1.5.1 The electrostatic sheath

The electrostatic sheath potential in front of a target surface is derived from assuming a given net electrical current j_pl and a known ion particle current Γ⁺ (i.e., j⁺ = q_iΓ⁺ electrical ion current, with q_i=eZ_i being the ion charge, e the elementary charge and Z_i the ion charge number) flowing over the target sheath edge. Both these currents, and the electron current j⁻=q_eΓ⁻ = −eΓ⁻, are taken to flow parallel to the B-field at the plasma-sheath edge.

Often, but not necessarily, it is assumed: j_pl = j⁺ + j⁻ = 0, i.e., locally ambipolar flow conditions.

We start with a Maxwell-Boltzmann density distribution n_e(x) for the electrons, with electron temperature T_e, in the electron retarding electric field with potential Φ(x). Here x is a coordinate parallel to the magnetic field B. The electron current density along the field is then (after reduction, possibly, due to secondary electron emission):

j⁻_T = −e ·1/4 ·n_e,T ·v_e^av ·(1−γ_s.e.e.) = −e ·1/4 ·n_e,S ·exp[e∆Φ_S,T/kT_e] ·

⎛
√

8 kT_e

πm_e

·(1−γ_s.e.e.)

Here γ_s.e.e. is the secondary electron emission coefficient and v_e^av is the average electron velocity. The subscripts S and T stand for the plasma-sheath interface and for the target surface, respectively and ∆Φ_S,T is the (typically negative) voltage drop from the plasma sheath edge to the target.

In order to determine the sheath voltage drop ∆Φ_S,T we equate the component normal to the target of this current density with the corresponding component of the net ion current density:

−j⁻_T ·cos(Ψ) = (j⁺_T − j_pl) ·cos(Ψ)

Hence the cosine of the angle Ψ of field line inclination against the target surface cancels out, as does the elementary charge e. The ion electrical current j⁺_T ( = j⁺_S, no ion sources within the sheath region) is given as

j⁺_T = e

∑
i

Z_i ·n_i,S ·V_||,i,S = e

∑
i

Z_i ·n_i,S ·M_i,S ·c_i,S.

The sum is over all ion species, Z_i is the charge state of species i, n_i,S is the corresponding ion density and V_||,i,S the parallel (to the B-field) ion flow velocity of ion species i. M_i is the Mach number, and c_i is the ion acoustic speed for species i.

Hence, the sheath potential difference ∆Φ in terms of given j_pl, T_e, γ_s.e.e. and ions flows Γ_i, reads:

− e ∆Φ/kT_e

⎡
⎢
⎣

√

2π

(1−γ_s.e.e.)

⎛
⎝

∑
i

⎛
⎝

Z_i

n_i,S

n_e,S

·V_||,i,S

√

[(m_e)/(kT_e)]

⎞
⎠

−

j_pl

e ·n_e,S ·

√

[ˉ([(kT_e)/(m_e)])]

⎞
⎠

⎤
⎥
⎦

(5.57)

⎡
⎢
⎣

√

2π

(1−γ_s.e.e.)

⎛
⎝

∑
i

⎛
⎝

Z_i

n_i,S

n_e,S

·M_i,S ·c_i,S

√

[(m_e)/(kT_e)]

⎞
⎠

−

j_pl

e ·n_e,S ·

√

[ˉ([(kT_e)/(m_e)])]

⎞
⎠

⎤
⎥
⎦

(5.58)

This expression (5.54) is evaluated by the function SHEATH of the EIRENE code, taking local plasma data for T_e and n_i and V_||,i (n_e from quasi-neutrality) at the point of impact of an ion at a wall boundary (e.g., the divertor target plate).

Special forms of (5.55) are often quoted in the literature. Writing for the ion acoustic speed c_i explicitly

c_i =

⎛
√

⎛
⎝

γkT_i + kT_e

m_i

⎞
⎠

(γ denoting the adiabatic coefficient) and replacing this in (5.55):

c_i ·

⎛
√

m_e

kT_e

⎡
⎣

⎛
⎝

γkT_i

kT_e

+ 1

⎞
⎠

m_e

m_i

⎤
⎦

1/2

we recover well known expressions for simpler cases. Consider, e.g., a single ion fluid case, n_e = n_i, with Z_i = 1, and M_i = 1 (Bohm sheath condition for isothermal flow, γ = 1), furthermore zero net electrical current: j_pl = 0 and zero secondary electron emission γ_s.e.e = 0. Then we find:

− e ∆Φ/kT_e =

⎛
⎝

2πm_e

m_i

⎛
⎝

T_i

T_e

+ 1

⎞
⎠

Hence, e.g., for a pure hydrogen plasma, and for T_e ≈ T_i, we have −e ∆Φ/kT_e ≈ 2.5 and ≈ 2.84 for a pure deuteron plasma.

1.6 Combinatorial description of geometry

EIRENE is a Monte Carlo Code. Defining geometry in a Monte Carlo code is complex and laborious.

Of the three typical options for definition of the target geometry in Monte Carlo Codes:

"Constructive solid geometry (CSG)",
"boundary representation (BREP)"
"Voxel Geometry"

EIRENE uses a combination of the latter two concepts.

In CSG a solid object is "constructed" from the intersection, union or other boolean operations of several "half-spaces". Each half space in specified by a relation F(x,y,z) ≤ 0, typically F is a low order polynomial. For example the neutronics code MCNP used e.g. for ITER nuclear safety studies, is based on CSG.

In BREP geometry, surfaces and volumes are build up from constituent parts. These parts may be related with algebraic functions (elementary BREPs) or complex functions such as Bezier functions or B-splines (advanced BREPS). Most CAD software uses BREPs in their native definition of geometry.

In the EIRENE BREP options surfaces are build from (bounded) algebraic functions. These surfaces are referred to as "additional surfaces", see section 2.3.2. EIRENE uses only algebraic functions, up to second order (some extensions to raise the order to fourth order have been started recently but are not ready to use).

In Voxel geometry, a 2D or a 3D object is constructed from identically shaped volume elements ("voxels"). In EIRENE these may be triangles or quadrangles in 2D, or tetrahedrons in 3D. New options to accommodate trilinear hexaedra as 3D voxels are currently being developed in connection with 3D edge codes such as EMC3. The bounding surfaces of these voxels are referred to as "standard grid surfaces" in EIRENE, see section 2.2. These volume discretisation meshes (VOXELS) can be read into EIRENE from external files. For example the mesh in B2 CFD plasma fluid code is read into EIRENE directly from the CFD grid generator in coupled B2-EIRENE applications. Wall meshes are build from these by assigning a surface reflection model to some of the (voxel-) surfaces, whereas internal (grid-) surfaces remain transparent for flow of test-particles.

Note added in 2008:

An attempt has been made in 2008 to use CAD files directly for 3D wall surfaces. This was done via the already existing CAD-MCNP interface. We had then tried to turn the CSG input for MCNP (see above) into a BREP form for EIRENE. In principle this works. But in practice the de-tour via the CSG concept leads to prohibitively complex BREP input. A direct CAD-BREP interface seems now much more efficient to us and more natural anyway, but this has not been started yet in EIRENE due to lack of manpower.

A complete set of surfaces in an EIRENE geometry may consist of a combination of both: standard surfaces (hence: defining standard grid cells, (voxels), and additional surfaces, defining a general BREP discretisation (additional cells). A more complicated situation arises when additional surfaces intersect standard grid voxels. If these surfaces are transparent, they can be used to conveniently "measure" fluxes at any position. But if they are reflecting, or absorbing, e.g. to define a complex shaped bounding surface in an otherwise simple grid, then some care is needed because the automatic cell volume evaluation, which is available in EIRENE for standard grid cells (voxels), may not correspond anymore to the cell volume accessible for test particles. Cell averaged quantities (divided by cell volume) may then not necessarily be scaled properly. For a possible work-around see NLERG option in input block 1, section 2.1.

From the above it follows than in an EIRENE model surface elements building up the wall can be defined as bounded algebraic surfaces, with triangles as special case, i.e. as "additional surfaces, BREP", or can be bounding surfaces of the voxels, if the "grid of voxels" extends up to the real wall, or a combination of both.

As can be seen from the flowcharts, all the geometrical calculations needed for tracing test flights and for estimating volume- and surface averaged tallies are compiled in a "geometry block", which can easily be exchanged. In most cases (e.g., all installations of EIRENE outside FZ-Jülich) we load the most complete, fully 3D block "GEO3D", which allows fully 3 dimensional spatial resolution. Lower dimensional (0D, 1D or 2D) are contained as special cases, e.g. by choosing only grid surfaces with proper symmetry.

There is, in general, not a very large CPU benefit from running EIRENE with simpler specialized and geometry optimized blocks such as GEO2D or GEO1D for 2D or 1D models, respectively, because the saving in CPU costs in most cases is only around 10 - 20 %, so these options have not been supported since the late eighties of the last century anymore.
The following operations have to be executed by the block:
Suppose the three dimensional computational volume is discretized by a mesh of n surfaces S_i, i=1,2,...,n . Each of these surfaces is defined by a set of Λ_i equations (unbounded surfaces):

f_λⁱ (x,y,z) = 0 , λ = 1,2,..., Λ_i ,

(6.59)

and, optional, for each unbounded surface λ a subset of Π_λ inequalities (e.g. for turning the unbounded surfaces into finite segments):

g_{λ, π}ⁱ (x,y,z) ≤ 0 , π = 1,2,..., Π_λ .

(6.60)

This means that for all unbounded surfaces S_i,λ only the points, which fulfill all corresponding inequalities (6.57) are regarded by EIRENE as valid parts of the finite surface element.

Assume a test particle at position r₀ = (X0,Y0,Z0) moving with speed unit vector Ω₀ = (VELX,VELY,VELZ). The block then has to provide the nearest intersection r₁ of the straight line ("sample trajectory")

G(t): r₀ + t Ω₀

(6.61) amongst all possible intersections with the surfaces S_i. Moving the "particle" from r₀ to r₁ and then repeating the block until the total track length l as sampled from the transport kernel T(l) (2.9a) is reached, allows to generate histories with all information available along the track to evaluate the estimators mentioned in section 1.2.
In EIRENE f and g are general 2nd order (with linear surfaces as special case) algebraic equations of the form

h(x,y,z) = a₀ + a₁ ·x + a₂ ·y +a₃ ·z

a₄ ·x² + a₅ ·y² +a₆ ·z² +

a₇ ·xy + a₈ ·xz +a₉ ·yz

(6.62)

Hence in the geometrical block nothing else but second order equations
t² + p ·t + q = 0
have to be solved accurately and as efficiently as possible. We note in passing that each of this surfaces may be transparent, purely absorbing, semi transparent or (partially) reflecting. Surface sources may be defined on each surface. These surfaces may be "invisible" from one side and transparent, absorbing or reflecting from the other and each surface can "operate switches" if intersected by a test flight, such as abandoning the evaluation of atomic or molecular reaction rates (entrance into vacuum regions) along the flight, abandoning evaluation of volume averaged responses (e.g., in some uninteresting regions) or abandoning geometrical calculations (for those surfaces, for which it can be known a priori that they will not intersect a flight path originating at a particular surface or in a particular cell).

1.7 time dependent mode of operation

The section is currently rewritten. Please refer to the J.Nucl.Mater. paper: D.Reiter et al., Vol 220-222 (1995), 987-992.

The following issues have to be considered:

The coupling is explicit, hence CFL-type conditions for the time-step have to be obeyed. Here this means: the time-step must be shorter than the time between re-ionization and subsequent neutralization (in the volume or at the targets), because the fate of a particle during its intermediate ionized phase is not included in the explicit Monte Carlo procedure.

The noise level for collision estimators increases with smaller time-steps. Hence: avoid all collision estimators in t-dependent mode.

Please contact us for an update on this situation, if you wish to use that option.

1.8 non-linear effects: coupling to plasma fluid models

The section is currently rewritten. Please refer to the report Jül-2872, Feb 1994, by G.Maddison and D.Reiter for details (ref. [9]).

Figure 1.1: Structure: EIRENE code

Figure 1.2: Structure: EIRENE code

Figure 1.3: Structure: EIRENE code

Figure 1.4: Convergence of B2/EIRENE, measured by B2-residuals. Method: implicit coupling, method of "saturated residuals", ref.[9]

1.9 non-linear effects: neutral-neutral collisions

Neutral-neutral collisions make the Boltzmann collision integral non-linear (bi-linear). More generally, inclusion of any binary collisions between two test-particle species i₀, j₀, elastic or not, have this effect. Important examples are, of course, elastic self collisions (i₀,i₀)_el, e.g., H₂−H₂ collisions leading to viscous effects in the gas, elastic cross collisions (i₀,j₀)_el, such as H−H₂ collisions, often leading to a heating of the molecular gas due to hot charge exchange atoms. Furthermore there may be relevant inelastic binary particle collisions (i₀,j₀)_inel, such as (H₂−H⁺₂→ H+H₃⁺) which is an important channel in hydrogen plasma chemistry, and (hν+ H(1s) → H(n=2)), the absorption of resonance Lyman-alpha photons under non-LTE conditions. In large dense divertor plasmas the ratio H(1s) to H(n=2) may have to be calculated consistently with the radiation field (photon Boltzmann equation), competing electron impact collisions and neutral atom transport.

Two methods of solution for this non-linear problem have been tested in EIRENE, 1) a parametrization of the single particle distribution functions involved in non-linear collision terms and an iterative scheme (see next subsection: BGK model collision terms), and 2) a Direct Monte Carlo Simulation (DMCS) technique based upon a swarm simulation, an operator splitting ("Streaming" and "Collisions") and a random simulation of binary collisions in the swarm. Currently only the first of these two options is still being used because, despite its approximate character, it proved to be much more effective than the full DMCS method, for the cases relevant to fusion research applications studied so far. A simple estimate for the relevance of neutral-neutral gas collisions is given by the so called Knudsen number Kn:

Kn =

λ_el

√2 σ_el n_N l

(9.63)

Here λ_el is the mean free path for elastic neutral neutral collisions, σ_el is the typical cross section, n_N the gas density (cm⁻3) and l a characteristic dimension (cm). Fixing the hard sphere molecular diameter at d=4 ·10⁻⁸ cm results in a typical cross section of 1−2 ·10⁻¹⁵ cm² to be used in this Knudsen number. For Kn << 1 a neutral fluid model (Euler, or Navier Stokes) is appropriate, while Kn >> 1 leads to the free molecular flow regime, in which only the linear collisions terms (collisions of neutrals with a given plasma host medium) are retained.

1.9.1 BGK approximation

currently under development. Please contact us for details and status.

See flowchart 9.1 (figure 1.12).

Further details can be found if ref. [18].

Consider first the elastic self collisions amonst particle species A. The BGK approximation to the Boltzmann collision integral reads:

∂f_A

∂t

(9.64)

1.9.2 Direct Simulation (DMCS) of self-collisions

currently under development. Please contact us for details and status. Further details can be found in ref. [19]

1.10 Radiation transport, photon gas simulations

The theory of light and its interaction with matter is one of the best that physics can offer, and can predict virtually every observed phenomenon with great accuracy [Feynman, 1985, QED: The strange Theory of Light and Matter, Princeton University Press, Princeton]. In our simulations of radiation transport with the EIRENE code we assume a geometric optics model (i.e., essentially, "particle theory of light"), which is sufficient for all our purposes. Light is emitted, scattered, reflected and absorbed, but moves on straight lines between these events. I.e., a continuously varying index of refraction is not allowed, and we ignore most properties of light that depend on a wave or quantum model for their explanation (e.g. diffraction, or interference). However, the linear wave-optics phenomenon of polarization and the linear quantum optics phenomenon of fluorescence [absorption of a certain wavelength (line), and re-emission at another wavelength (line)] can be added to our linear geometric optics model.
We are concerned with the kinetic theory of interacting particles and photons, where the behaviour of each species is governed by a kinetic equation: the kinetic equation for photons being nothing else but the well known equation of radiation transfer, i.e., the "strangely normalized photon-Boltzmann equation". In contrast to plasma radiation (e.g., Bremsstrahlung) the radiation processes studied here are due to single, uncorrelated particles such that emission and absorption of the plasma is obtained by simple summation of all contributions of all individual particles. We therefore consider all particles as uncorrelated and in well defined one-particle states.

For continuum radiation (bound-free or free-free) the frequency of the photon in the atoms rest frame may be replaced by the frequency in the lab. frame, because all cross sections are slowly varying functions of frequency. By contrast, for line-radiation due to bound-bound transitions, the Doppler effect must be taken into account in all quantities that vary rapidly as function of frequency, such as cross sections, line-profile, or redistribution functions

The key atomic parameters for radiation transport are the Einstein (or Einstein-Milne) coefficients A_mn, B_mn and B_nm of spontaneous emission, induced (stimulated) emission and absorption, respectively. We will discuss them here first for a gas of stationary atoms and include Doppler- and other line broadening effects later.

For a thermal radiation field (temperature T) in a gas of atoms with density n₁ (lower state) and n₂ (upper state) the rate equation expressing balance between emission and absorption reads:

n₁ B₁₂ R_ν(T) = n₂[A₂₁+B₂₁ R_ν(T)]

R_ν(T) is a specific function characterizing the wavelength and temperature dependence of the thermal radiation field (i.e., the Planck-function). For example the specific intensity I_ν (energy/area/time/solid-angle/frequency-interval) is often used, or the (angle)-integrated intensity J_ν = ∫I_ν dΩ, or the specific energy density u_ν = 1/c J_ν (energy/volume/frequency-interval), etc. Also the wavelength dependence may be given either using frequency ν (Hz), ω (rad/s), wavelength λ or energy E = h ν = (^h/_2p) ω.

Depending upon all these choices the numerical factors in the definitions of the B_nm coefficients (sometimes even in the A_nm) differ in the literature. For our purpose (Monte Carlo solution of the full radiation transfer equation) we need to define monocromatic mean free pathes for photons, rates and rate coefficients and cross sections for the individual processes.

In order to take advantage of the analogy between radiation and particle transport equations as much as possible, we will choose the numerical constants in the Einstein coefficients such that these parameters have the same meanings and even the same units in both cases of radiation and particle transport. I.e., we will use the "absorption coefficient" (better: "extinction coefficient", if scattering is allowed) α_ν = n σ ≡ Σ_t, (invers monocromatic mean free path), which directly corresponds to the "total macroscopic cross section" (dimension: 1/length) defined earlier for particle transport problems.

1.10.1 Line shape options

Emission of photons from a particular line is simulated by sampling firstly the location of emission, and secondly the frequency and direction of emission. Spatial distribution of photon emission is treated in EIRENE in exactly the same way as for any other particle type (neutrals, ions), and actually using the same modules of code for that. Emission from the transition i→ k from an upper state i=nl of an atom is treated, by abuse of language, as "volume recombination source": The position is sampled from the spatial distribution A_ikP^nl([(r)\vec]), which, because the Einstein emission coefficient A_ik is constant, is the same as the spatial distribution P^nl([(r)\vec]) of upper state particle density. ..... Given the birth point [(r)\vec] of the photon, and hence all local parameters of the host medium at this point, the frequency ν (in fact: the energy E=hν) is sampled from the normalized line shape function Φ(ν), .... The following line shape options are currently available:

ISHAPE=0 δ-profile: mono-energetic emission (not yet available for absorption)
ISHAPE=1 pure Doppler profile
ISHAPE=2 Lorentz profile, truncated at zero
ISHAPE=3 Doppler convoluted Lorentz profile, i.e. Voight profile
ISHAPE=4 LORVDW: convoluted Lorentz and van der Waals pressure broadened profile
ISHAPE=5 Doppler broadened LORVDW profile
ISHAPE=6 normal Zeeman triplet, δ-profile in each component
ISHAPE=7 Doppler broadened Zeeman-δ-profile in each component
ISHAPE=8 Doppler broadened Zeeman-Stark-profile, 10 Lorentzian model, Rosato 2006

1.11 Charged Particle Transport

EIRENE also follows charged particles since EIRENE₁₉₈₇. The particle tracing is carried out in Subroutine FOLION. The differences of this module to the module FOLNEUT for neutral particles (and photons) are:

Motion is not along straight lines between events, but instead determined by the electro-magnetic fields.
Time stepping is not of Poisson type, but with constant increments, because the equation of motion has to be solved numerically.
after each time-step a Coulomb collision kernel is modelled, i.e. a diffusive step in velocity space.
after each time-step a diffusion-step is modelled, i.e. a diffusive step in physical space (not yet written).
ionized test particles which hit a surface see an electrostatic sheath potential, by which they are either reflected, slowed down or accelerated, depending on the sign of their charge relative to the sheath potential.

1.11.1 Orbit integration

Particles following the B-field

A: no electric field

Guiding centre approximation

1.11.2 Coulomb Collision Models

Simple Coulomb Collision Models

Strongly simplified Coulomb collision models are usually sufficient for molecular ions, because of their very short lifetime until further fragmentation (dissociative excitation, dissociative recombination, etc...) to neutral products. In order to roughly account for thermalization effects and background plasma cooling or momentum transfer the energy, parallel and perpendicular velocity of such ions can be modified "continuously" during particle orbit integration, hence approximating the velocity space diffusion (Fokker Planck equation) by mean values of energy and momentum of the test particle. There are two classes of simple collision models: deterministic and stochastic.

A: Deterministic relaxation in velocity space The energy and/or velocity of charged test particles are modified in a fully deterministic way along the particle trajectory. I.e.: along with orbit integration also energy, or certain velocity components are evolved according to ordinary (non-stochastic) differential equations, using appropriate relaxation times.
A1: Only energy relaxation
The simplest Coulomb Collision model is an averaged BGK-type relaxation of energy of test particles to the mean energy of the Maxwellian background ions. In the earliest versions of the EIRENE code (mid eighties) it was originally based upon the hydrocarbon ion equilibration model given by Langer, [20] and it is a slight generalization thereof. The equation for the energy ε_a of a test-ion a in a background of ions b with mean energy ε_b = 3/2 T_b is given as:

dε_a

=−ν_ε^a,b ·ε_a , ε_a(t=0)=ε_a⁰

(11.65)

This is a non-linear equation because the energy exchange collision rate (energy transfer rate) ν_ε^a,b depends on ε_a itself. The full expression for ν_ε^a,b is given below, it involves the Maxwell integral Ψ and its derivative Ψ′, see below.

In the limit of small test particle velocities (compared to thermal velocities of field particles b), i.e.: v_a << √{(2T_b/m_b)}, i.e. x^a/b = μ_bv_a²/2kT_b << 1, as often appropriate for molecular ions, one has [21]:

ν_ε^a,b = 2 ν_s^a,b − ν_⊥^a,b − ν_||^a,b.

(11.66)

Here the energy transfer rate ν_ε^a,b is given via the "slowing down rate" ν_s^a,b, the transverse diffusion rate ν_⊥^b and the parallel diffusion rate ν_||^b. In this limit on has for the basic "Spitzer rates" (loc. cit.):

ν_s^a,b

n_b Z_a² Z_b²λ_a,b·6.8 ×10⁻⁸

μ_b^1/2

μ_a

⎛
⎝

μ_b

μ_a

⎞
⎠

T_b^−3/2

(11.67)

ν_⊥^a,b

n_b Z_a² Z_b²λ_a,b·1.4 ×10⁻⁷

μ_b^1/2

μ_a

T_b^−1/2ε⁻¹

(11.68)

ν_||^a,b

n_b Z_a² Z_b²λ_a,b·6.8 ×10⁻⁸

μ_b^1/2

μ_a

T_b^−1/2ε⁻¹

(11.69)

and hence, using Equation (11.62) :

ν_ε^a,b

n_b Z_a² Z_b²λ_a,b·6.8 ×10⁻⁸

μ_b^1/2

μ_a

T_b^1/2

⎡
⎣

T_b

⎛
⎝

μ_b

μ_a

⎞
⎠

−2

−

⎤
⎦

(11.70)

ν^c_ε +ν^v_ε

(11.71)

with the "constant" and "variable" contribution ν^c_ε and ν^v_ε, respectively, defined as

ν^c_ε

2 n_b Z_a² Z_b²λ_a,b·6.8 ×10⁻⁸

μ_b^1/2

μ_a

T_b^3/2

⎡
⎣

μ_b

μ_a

⎤
⎦

⎡
⎣

μ_b

μ_a

⎤
⎦

(11.72)

ν^v_ε

2 n_b Z_a² Z_b²λ_a,b·6.8 ×10⁻⁸

μ_b^1/2

μ_a

T_b^3/2

⎡
⎣

−

3T_b/2

ε_a

⎤
⎦

⎡
⎣

−

3T_b/2

ε_a

⎤
⎦

(11.73)

Using this decomposition and definition of [(ν)\tilde]_ε in equation (11.61) we find:

dε_a

=−

⎡
⎣

μ_b

μ_a

⎤
⎦

·ε_a +

3T_b

(11.74)

The solution of this linear differential equation is

ε_a(t)=ε_a⁰ ·exp(−

[1+μ_b/μ_a] t)+

3T_b

⎡
⎣

1−exp(−

⎤
⎦

(11.75)

With μ_b=1 for protons, μ_a=16 for CH₄, λ_a,b ≈ 10, and [1+μ_b/μ_a] ≈ 1, hence [(ν)\tilde]_ε ≈ 8.8 ×10⁻⁸ n_pT_p^−3/2 this is also the solution Eq. (11) given in [20] for thermalization of hydrocarbons in a Maxwellian hydrogen plasma. And this is exactly the approximation which was also used in EIRENE since about 1987 for all test ions traveling in a bath of background ions.

To be done: in case of drifting Maxwellian background: first transform into frame such that stationary Maxw. background, then relax test particle energy in this frame, then transform back...

The full expression for the energy relaxation frequency ν_ε^a,b, without the approximation of slow particles a, i.e. for all values of x^a/b is:

ν_ε^a,b(x^a/b) = 2 [μ_a/μ_b Ψ(x^a/b) − Ψ′(x^a/b)] ν₀^a/b

(11.76)

where x^a/b has been defined above, ν₀^a/b = ... and

Ψ(x) =

√

⌠
⌡

x

0

dy e^−y √y

= erf(√x)−

√

e^−x√x

(11.77)

Ψ′(x) =

dΨ/dx

√

e^−x√x

(11.78)

With this general expression for the collision frequency Equation 11.61 can then only be integrated numerically along the particle trajectory.

1.11.3 The full Fokker-Planck Collision Model, A

The full Fokker-Planck Coulomb Collision Model, which was originally developed for a stand alone trace impurity ion transport code DORIS [22], is accommodated now as collision kernel in EIRENE: Thesis J. Seebacher, Univ. Innsbruck, soon to come

1.11.4 The full Fokker-Planck Collision Model, B: binary Monte Carlo collision procedure

1.12 EIRENE flow charts

Figure 1.5: Structure: EIRENE code

Figure 1.6: Flowchart: Program EIRENE

Figure 1.7: Flowchart 2: Subroutine INPUT

Figure 1.8: Flowchart 6: Subroutine MCARLO

Figure 1.9: Flowchart 6.3: Particle loop in subr. MCARLO

Figure 1.10: Flowchart 9.1: Iterative mode for neutral-neutral interactions

Chapter 2
Description of formatted input file

General Remarks

All input data from unit IUNIN are read in subroutine INPUT using formatted I/O. There are 14 "blocks", each one starting with a card

*** n. TEXT

where n is the number of the input block, n=1,2,...,14. Some blocks are subdivided into so called ßub-blocks", each one starting with a card

** TEXT .

Each input block is described in one of the following 14 sections 2.1 to 2.14:

2.1 Input data for operating mode
2.2 Input data for standard mesh
2.3A Input data for "Non-default Standard Surfaces"
2.3B Input data for Ädditional Surfaces"
2.4 Input data for species specification and atomic physics module
2.5 Input data for plasma background
2.6 Input data for surface interaction models
2.7 Input data for initial distribution of test particles
2.8 Additional data for some specific mesh zones
2.9 Data for statistics and non-analog methods
2.10 Data for additional volumetric and surface averaged tallies
2.11 Data for numerical and graphical output
2.12 Data for plasma diagnostic module DIAGNO
2.13 Data for nonlinear and time-dependent mode
2.14 Data for interfacing with external codes

This last block 14 may be read in from the user supplied part (with any format specified there) from this same unit IUNIN, (but equally well from any other unit specified there), e.g., from the subroutine INFCOP for interfacing with other codes.

Units

cgs units are used, with two exceptions:

All temperatures and energies are in eV.

All particle fluxes are in Ampere (i.e. #/s ·1.602210⁻¹⁹), even neutral particle fluxes.

E.g., to convert a gas flow rate of x sl/min (standard liters per minute) first divide x by 22.4 (then: mol/min) next multiply by Avogadro's no. N_A = 6.0221 ·10²³, (i.e., then particles per minute), then divide by 60 (then: particles per second), and finally multiply by the elementary charge 1.6022 ·10⁻¹⁹ (then Ampere, as used for fluxes in EIRENE).

All densities are in cm⁻³.

All velocities are in either cm/s or in (isothermal) Mach-number units.

All geometrical data are in cm.

Note: consequently, e.g., energy densities are in eV/cm³, and energy fluxes are given in eV ·Amp, i.e., in Watt.

Hence: In order to convert from energy density units into pressure units in N/m² (i.e., Pascal: Pa) one must multiply by 1.6022 ·10⁻¹⁹ (eV to Nm) and multiply by 10⁶. Check: 1 Pa gas at 273 K corresponds to about 2.7 ·10¹⁴ ×0.026 eV/cm³ ≈ 7 ·10¹² eV/cm³

For example, energy densities as computed by EIRENE defaults tallies 4, 5, 6 (see table 5.2) must be converted into "temperature densities" (a factor 2/3, if one neglects the energy of the directed motion) and then be re-scaled to pressure units as described above.

For further scaling to other pressure units note, e.g.:

101.3 Pa = 760 mTorr

Further: momentum is in g cm/s, and a momentum flux or momentum rate of change is in g cm/s ·Amp (per cell) and the momentum source density in the EIRENE tallies for interfacing with plasma codes is then in g cm/s Amp /cm³.

Hence, conversion from g cm/s Amp into Pa ·m² (momentum source per cell in plasma momentum balance equations in SI units) is done by firstly multiplying with 10⁻⁵ (then: kg ·m/s ·Amp) and then dividing by 1.6022 ·10⁻¹⁹ (hence: kg ·m/s²). After dividing by the cell volume (m³) the momentum source is in Pa/m, or N/m³ (force density).

Input format

The format statement numbers in subroutine INPUT start with 666... for reading and with 777... for writing. The following conventions are used:

: Standard Fortran naming conventions, i.e., variable names starting with the letters I,J,...,N are Integer, all other variables are Real.

In addition to these standard convention the following rules are employed:

: Character strings start with the letter C... or with T...
: Logical variables start with NL..., LG... or with TRC...

The format for a card containing only real variables: (6E12.4)

The format for a card containing only integer variables: (12I6)

The format for a card containing only logical variables: (12(5L1,1X))

The format for a card containing text: (A72)

The format for a card containing K (K < 6) integer variables first and up to 6-K real variables next: (K(I6,6x),(6-K)E12.4)

An input card starting with '*' is recognized as a comment line, i.e., it does not contain input data. These comment lines should not start with '**' or '***', because in this latter case EIRENE would assume that the input for the current block is completed. It would expect the first card of the following block (or '* comment...' lines belonging to this next block) as next card. Any data not read from the input file, before the following block is read, are set to default values.

An arbitrary number of such comment lines may be included in the input file, but each block of comment lines must always be put at the beginning of a block or sub-block. I.e., comment lines are only identified by EIRENE if preceded by a card starting with '**' or '***' indicating a new sub-block or block respectively.

Example (see section 2.7 below):

*** 7. Data for initial distribution of test particles

* this is a comment

* there are NSTRAI=2 strata in this run

* this is a comment

** Data for first stratum

* this is a surface recycling source of helium ions, IPLS=2

* this is a comment

FFFF

** Data for second stratum

* this is a hydrogen ion volume re-combination source, IPLS=1

FFFF

*** 8. Data for some specific zones

and so on.

In very few cases formats other than these are used. They are then explicitly mentioned in this manual and in case of doubt the user should check in subroutine INPUT.

2.1 Input data for operating mode

General Remarks

The variables in this block control some of the more general options in EIRENE such as overall running time, usage of dump files but also a few parameters depending on the simulation model (drifts included or not, etc.).

The Input Block

TXTRUN

*** 1. Data for operating mode

first card: integer flags controlling global options for the entire run.

NMACH NMODE NTCPU NFILE NITER0 NITER NTIME0 NTIME

next input line is only for EIRENE₂₀₀₂ and younger, and optional in these versions (defaults given below).

NOPTIM, NOPTM1, NGEOM_USR, NCOUP_INPUT, NSMSTRA, NSTORAM, NGSTAL, NRTAL, NREAC_ADD

next line (second mandatory line): logical flags for global options.

NLSCL NLTEST NLANA NLDRFT NLCRR NLERG NLIDENT NLONE LTSTV

next input lines are only for EIRENE₂₀₀₅ and younger, and optional in these versions. There can be any number of these cards starting with character string CFILE in the input file.

"CFILE" DBHANDLE DBFNAME

Meaning of the input flags for operating mode

TXTRUN

At least one line of text to identify the run on printout and plot files.

Each text card must start with *. The first of these text lines will appear on all plots, the full text will be echoed at the beginning of the printout file

NMACH

Code-number for the computer used

= 1: CRAY
= 2: IBM
= 3: FACOM
= 4: VAX

This flag is currently not in use. It only affects the round-off error margin used for geometrical computations. In case NMACH=1, the error tolerated in the geometry routines is EPSGEO=1.D-10, and EPSGEO=1.D-6 for all other values of NMACH.

NMODE

Operating mode

= 0: EIRENE run as stand-alone-code. Code coupling segment couple_Dummy.f may be used. Input block 14 has a fixed format, see section 2.14.
≠ 0: EIRENE calls the code coupling subroutine INFCOP in the code coupling segment couple_Name.f where 'Name' is a character string identifying the particular external code, to which EIRENE is to be coupled (e.g.: Name=B2, Name=DIVIMP, Name=U-file, Name=FIDAP, Name=EMC3, Name=OSM, etc.).
Hence: the routine INFCOP is called by EIRENE for communication with external data sources, e.g., other codes. The calling program for the first 3 entries of INFCOP is subroutine INPUT, and the call to subroutine INFCOP is after reading 13 blocks from the formatted input file (READ (IUNIN,...)). A 14^th block of the input file is read from subroutine INFCOP with the format (if any) specified there. At entry IF0COP geometrical data are provided (overruling corresponding data in input block 2).
At entry IF1COP plasma profiles are defined (overruling the background data specified in block 5).
Any other data (e.g., surface- or volume source distributions overruling the data in input block 7) are expected from entry IF2COP.
> 0: If NMODE > 0, subroutine INFCOP is called once again after each completed stratum (at the entry IF3COP(ISTRAA,ISTRAE)) to return data to another code (post-processing). The call to IF3COP is from subroutine MCARLO, at the end of the DO 1000 -loop over the strata. One final call to the interfacing routine (entry IF4COP) is implemented after the calls to the EIRENE printout- and plot routines, e.g. to perform global balances etc. after summation of the contributions from the individual strata.

NTCPU

Maximum number of CPU seconds allowed for this run. NTCPU must be less than or equal to the time parameter in the job-card (if any).

If more than one iteration is carried out (NITER, see below, this section) or more than one time cycle (NTIME, see below, this section), then each iteration or cycle can take up to NTCPU seconds.

In 2008 the definition of NTCPU was slightly changed: the initialisation overhead is not any longer included in NTCPU. I.e., this variable NTCPU now is close to the true Monte Carlo sampling time. The total cpu-time, including initialisation overhead as well as post-processing is printed at the end of an EIRENE run (see: "total cpu-time of this run" in printout file).

NFILE

Flag for the use of dump files FT10, FT11, FT12, FT13, FT14 and FT15.

= 0: Neither reading nor saving of data is done.

= JKLMN:

N = 1: EIRENE writes output data from this run into files FT10 and FT11 to save them for other plot or printout options in later 'read only' runs with the same input file
N = 6: same as N=1, but only the data for the sum over all strata are written (sufficient, e.g., for BGK-iterations, or for post processing, (subroutine DIAGNO, see input block 12) as long as this post processing only is to be done on total results (not on contributions from individual strata).
N = 2: EIRENE reads output data from an earlier run from files FT10 and FT11. No new particle histories are computed, only the requested output is printed and plots are produced. Iteration cycles or time cycles (if any) are abandoned.
N = 7: same as N=2, but only the data for the sum over strata are read. See N=6 option described above.
M = 1: EIRENE writes geometry data into file FT12. These data are the output from subroutines GRID and VOLUME in the initialization phase.
M = 2: EIRENE reads geometry data from file FT12. The geometry subroutines GRID and VOLUME are not called. This option should be used if the geometry has not changed as compared to an earlier run. It reduces the CPU costs for the overhead.
L = 1: EIRENE writes plasma data and atomic and molecular data ("A&M data") into file FT13. These data are output from subroutines PLASMA, XSECTA, XSECTM, XSECTI, XSECTP in the initialization phase, or at the end of a run, if some data have been modified in subroutine MODUSR, see NITER flag below. Also the primary source parameters (input block 7) are saved, for later iterations or time-steps.
L = 2: EIRENE reads plasma data, A&M data from file FT13. It also reads the entire input block 7 from FT13, and overwrites the input read from this block 7 by those data. Routines PLASMA, XSECTA, XSECTM, XSECTI, XSECTP are not called. This option should be used if neither the plasma background nor the selection of atomic processes to be used, nor the primary source model (block 7) has changed as compared to an earlier run. It reduces the CPU costs for the overhead.
L = 3: Acts as if both NFILE-L=1 and NFILE-L=2. For continuation of iterative calculations, for example.
L = 4: Same as NFILE-L=3, except: primary source data (input block 7) are not read from fort.13, but are taken as specified in input block 7. This allows to modify these primary sources parameters during iterations or time-steps not only via module MODUSR, but directly via input file.
L = 6,7,8,9: Same as L=1,2,3,4, respectively, but XDR file format is used.
K = 1: EIRENE saves some data for optimizing non-analog sampling and stratified source sampling on file FT14.
K = 2: EIRENE reads from file FT14 and optimizes operation. Currently only allocation of CPU time in stratified source sampling is optimized.
K = 3: Acts as if both NFILE-K=1 and NFILE-K=2.
J = 1: Only for time-dependent option (see block 13). EIRENE writes "census data" at the end of last time-step onto file FT15.
J = 2: EIRENE reads "census data" from file FT15, and uses it as initial distribution for the coming next time-step.
J = 3: Acts as if both NFILE-J=1 and NFILE-J=2.

NITER0

Initial iteration number: (Default: NITER0=1) Irrelevant parameter. Only needed for book keeping and printout. EIRENE labels the iterations from NITER0 to NITER.

NITER

Number of iterations, if EIRENE runs in "iterative mode".

> 0: EIRENE calls user supplied subroutine MODUSR after completing the run; some model parameters may be modified here for the next iteration step, and some results from the previous step may be saved on a file.
> 1: EIRENE recalls itself but does not read from the formatted input file again. This recalling is repeated NITER times (including the first iteration). The CPU time NTCPU is used for each iteration. Hence the true CPU time then is NTCPU ·NITER .

NTIME0

Initial time-cycle number: (default: NTIME0=1). Irrelevant. Only needed for book keeping and printout. EIRENE labels the time steps from NTIME0 to NTIME.

NTIME

Total number of iterations ("cycles") in time carried out in one single run. The total time per cycle is defined as NTMSTP · DTIMV (see below, input block 13).
After each time-cycle the subroutine TMSTEP is called. In this routine the "census arrays", which store the test particle population at time t_i:

t_i = t_i−1 + ∆t = t₀ + i ∆t = TIME0 + ITIME·[NTMSTP ·DTIMV]

are filled and prepared for the next time-cycle. The census arrays from the previous time cycle (if any), i.e., at t = t_i−1, are overwritten here.

After the last time-cycle, the census arrays are written on file fort.15, in order to permit continuation in time in a next run.

The background conditions (and source distributions or any other input parameters) for the next time cycle can be modified in subroutine TMSUSR(t_i), which is called from subroutine TMSTEP.

If the population on the census array is not empty or known from a previous run (NFILE-J flag, see above), then this census population defines one additional stratum for the current cycle. I.e., the census population then determines the initial condition for the distribution function f(r,v,i,t=t₀). The source strength FLUX is computed from that initial condition.

NOPTIM

Default: 1

NOPTM1

Default: 1

NGEOM_USR

for value =1: user defined (external) geometry package (LEVGEO=10), then no grid storage is provided in EIRENE. Default: 0

NCOUP_INPUT

=1: Storage for data transfer via coupling routines, =0: no such storage. Default: 1

NSMSTRA

Sum over strata disabled for value 0, enabled for value 1. Default: 1

NSTORAM

Storage vs. speed in atomic data evaluation. Maximum storage, fastest computation: =9. Minimum storage, maximum work (slowest option) =0. Default: 9

NGSTAL

Storage for spatial distribution of surface tallies on non-default standard surfaces for value =1. For value =0: only total (spatially integrated) surfaces fluxes. Default: 0

NRTAL

Condensing mesh cells into fewer larger cells, for output volume tallies. The underlying fine mesh has NRAD cells (see input block 2). The coarser mesh, obtained from condensing cells into one larger cell, has NRTAL cells. Default: 0: Then internally: NRTAL=NRAD, and no condensation is carried out.

NCLTAL(IRAD) = IRTAL: grid cell IRAD is condensed into the larger cell IRTAL. I.e.: output is average over a larger cell IRTAL, which is comprised of all cells IRAD such that NCLTAL(IRAD) = IRTAL, IRAD=1, NRAD

The input data in input block 5 (background medium) are always given on the fine mesh (size: NRAD)

For output tallies (cell averaging) several cells IRAD1, IRAD2,... can be condensed into one larger cell IRTAL. Cell volumes, scoring, statistics are automatically done on the coarser grid (size: NRTAL).

The index array NCLTAL(NRAD) can be defined in the problem specific "user" routines (see section 3), e.g.: INIUSR, GEOUSR, etc..

Default: NCLTAL(IRAD)=IRAD for IRAD=1,NRAD

NREAC_ADD

Storage for additional reaction decks read onto EIRENE arrays in USR-routines, post-processing, etc. Default: 0

NLSCL

Some volume averaged tallies are re-scaled in order to exactly preserve the total number of particles, which otherwise would be the case only up to statistical precision (due to the use of track-length estimators). EIRENE computes three factors FATM, FMOL and FION such that particle balances for atoms, molecules and test ions, respectively, are accurately observed, if NLSCL = TRUE.

NLTEST

Tests for consistency between cell numbers and geometrical data along the particle tracks are carried out at each point of collision. If inconsistencies are detected, the history is stopped and an error message is printed. The contribution of these particles to the particle- and energy balances is stored in the bins "PTRASH" and "ETRASH" respectively.

NLANA

De-activates (NLANA=.TRUE.) all non-analog sampling distributions, such as biased source sampling, splitting, etc.. Select NLANA=.TRUE., if particle trajectory plots are used to get an intuitive picture of what is going on physically.

NLDRFT

Drift component in the bulk ions velocity distribution is included. (see input block 5). Otherwise an isotropic Maxwellian distribution is assumed for the background particles.

NLCRR

Correlation sampling is used.

NLERG

The case is automatically reduced to a case for estimating cell volumes by utilizing an "ergodic property".

NLIDENT

In multi-processor calculation mode NLIDENT forces all processors working on the same stratum to use the same sequence of random numbers and hence to carry out exactly identical work. This flag can be used to test the parallelized code version.

NLONE

The case is automatically reduced to a single species and "one speed transport" problem. (not ready, don't use)

LTSTV

free flag used for testing purposes

Optional input cards, for pathways to external databases. Their presence is identified by the starting character string CFILE of these cards.

DBHANDLE: Name of the database which is specified in this "CFILE" card. The format of a data file is identified in EIRENE by its name. Possible database names currently recognized are:
AMJUEL, HYDHEL, METHAN, H2VIBR, HYDRTC, ADAS, SPECTR, PHOTON, POLARI,

SPUTER, TRIM, TRMSGL,
with volumetric collision databases in the first block, and surface interaction databases in the second.
DBFNAME: Absolute or relative path and name of file to be used for this database.

Example:

CFILE AMJUEL /home/path/AMdata/amjuel.tex
CFILE TRIM /home/path/Surfacedata/TRIM/trim.dat

Default: If the specification of an external database name and path, which is later referred to in the run, is not given here, then these data files are expected, usually with the same name, in the directory in which the calculation takes place. Due to historical reasons and for backward compatibility for some data files also different default names are recognized.

Default file names to be used in this directory currently are

AMJUEL,
HYDHEL,
METHAN or "METHANE", for database name METHAN
H2VIBR,
HYDRTC,
SPECTR,
PHOTON,
POLARI,
SPUTER,
TRIM or "fort.21", for database name TRIM

For ADAS (AMdata) and TRMSGL (Surfacedata) there is an exception from this scheme. In these cases no default file names are available. And input variable DBFNAME specifies only the path to the directory which is the root of the ADAS and TRIM data tree, respectively.
Example:

CFILE ADAS /home/path/AMdata/Adas_Eirene_2010/adf11/
CFILE TRMSGL /home/path/Surfacedata/TRIM/

The particular data file from these trees (e.g.: SCD96 in case of ADAS, or A_on_B in case of TRMSGL), the chemical element and the charge stage of a particular ion Z=1,... Z+1 (e.g.: C, 4), are detailed later in input block 4 (for volumetric collision processes based on ADAS adf11 format), and input block 6 (TRMSGL, for the surface reflection database based on TRIM) see below, sections 2.4, 2.6, respectively.

2.1.1 The NLERG option for cell volumes

All non-transparent surfaces are made perfectly reflecting (mono-energetic, cosine distribution). All volumetric collision processes are de-activated. Except, possibly, for the very first flight, there is only either one atomic species (IATM=1) or one molecule species (IMOL=1). Unless otherwise specified in input block 13, a single time-step, with time limit DTIMV=1.0 (s) is set and up to NPRNLI=10 test flights are launched. Under such conditions, the particle density should be constant throughout the grid. Statistically significant deviations from that constancy indicate wrong cell volumes as compared to the volumes seen by the test particles (e.g. due to additional surfaces intersecting the grids, or complicated shapes of some additional cells (see input blocks 2 and 3). In a second run the volumes of these cells can then be set explicitly, e.g. in input block 8.

2.1.2 The NLMOVIE option for making movies of trajectories

documentation to be written

2.2 Input for Standard Mesh

General Remarks

EIRENE follows atoms, molecules or test ions in a 3-dimensional computational box. This box is discretized by zones (mesh cells), the boundaries of which are defined either by regular "standard mesh surfaces", i.e., co-ordinate surfaces, which are described here, and/or by zones whose boundaries are defined more generally by "additional surfaces" (see below block 3B).

As long as test particles travel inside the "standard mesh", all cell indexing computations are done automatically. At transition into the more general "additional mesh" and inside these additional cells, this indexing has to be specified by the flags ILSWCH, ILCELL etc. by the user, in input block 3B.

There are up to 3 sets of standard surfaces defined by the arrays RSURF, PSURF and TSURF.

Here RSURF(1), RSURF(2), ..., RSURF(NR1ST) are grid-points in the radial or x direction (depending on the geometry level, see below) defined by the input data in sub-block 2A. I.e., there are NR1STM= NR1ST-1 cells in radial- or x direction. With regard to the terminology in 1.6, eq. 6.56, a corresponding co-ordinate surface is labelled with a double index 1,IR and is defined by the equation

f₁^1,IR (x,y,z) =constant=RSURF(IR).

PSURF(1), PSURF(2), ..., PSURF(NP2ND) are grid-points in the poloidal or y direction (sub-block 2B), and there are NP2NDM=NP2ND-1 cells in this direction. A corresponding co-ordinate surface is labelled 2,IP and is defined by the equation

f₁^2,IP (x,y,z) =constant=PSURF(IP).

In case of the NLPLG option ( (mesh in x-y plane generated by a set of polygons) for simplicity and by an abuse of language we sometimes refer to the polygons themselves as radial surfaces or "radial polygons" RSURF and refer to the polygons obtained by connecting all points with a particular index along the various radial polygons as "poloidal polygons" PSURF. No use is made in EIRENE of the arrays RSURF and PSURF in connection with the polygon mesh option. The NLPLG option permits to run EIRENE on 2D computer generated meshes with grid cuts in one direction.

Furthermore, TSURF(1), TSURF(2), ..., TSURF(NT3RD) are grid-points in toroidal or z direction (sub-block 2C). A co-ordinate surface is labeled 3,IT and is defined by the equation f₁^3,IT (x,y,z) =constant=TSURF(IT).

In addition to these grids and, depending on the geometry options chosen, EIRENE then computes further "derived" grid arrays, for example a "(radial) surface labeling grid" RHOSRF, which may or may not be identical to the RSURF grid, and a "zone centered (radial) surface labeling grid" RHOZNE(IR) defined by:

RHOZNE(IR) = 0.5 · (RHOSRF(IR)+RHOSRF(IR+1))

These grids may be used for the input of closed form radial- or x direction plasma profiles or as abscissae in plots of e.g. poloidally (or y) and toroidally (or z) averaged tallies.

We have found it convenient to allow for several (NBMLT) identical copies of such "standard meshes", separated and bounded by arbitrary "additional surfaces". Input flags for this "mesh multiplication" option are comprised in sub-block 2D.

As mentioned above, in addition to these regular grids one can define cells of almost arbitrary complexity by using the "additional surfaces" defined in sub-block 3B. Such surfaces are labeled 0,IA in EIRENE. Appropriate cell number switching flags must be set on transparent additional surfaces separating these general "additional cells". Data relevant for the additional cells (e.g. the volume) are read in sub-block 2E. The number of additional cells is NRADD.

Therefore, the total number of cells in an EIRENE run is NSBOX, with

NSBOX = (NR1ST · NP2ND · NT3RD) · NBMLT + NRADD .

Cell volumes in the standard mesh region, as well as the center of mass vector in each of these cells are computed automatically. The corresponding data in the additional cell region must be specified by the user, e.g. in the input block 2E or 8, or in the user supplied routine GEOUSR (section 3.6), which is called by EIRENE in the initialization phase.

The Input Block

*** 2. Data for Standard Mesh

 (INDGRD(J), J=1,3)
 ** 2A. x- or radial co-ordinate surfaces
 NLRAD
 IF (NLRAD) then
    NLSLB  NLCRC  NLELL  NLTRI  NLPLG  NLFEM  NLTET NLGEN
    NR1ST  NRSEP  NRPLG  NPPLG  NRKNOT
    IF (INDGRD(1).LE.5) THEN
       IF (NLSLB.OR.NLCRC.OR.NLELL.OR.NLTRI) THEN
          RIA  RGA  RAA  RRA
          IF (NLELL) THEN
             EPIN  EPOT  EPCH  EXEP
             ELIN  ELOT  ELCH  EXEL
             IF (NLTRI) THEN
               TRIN  TROT  TRCH  EXTR
            ENDIF
          ENDIF
       ELSEIF (NLPLG) THEN
          XPCOR  YPCOR  ZPCOR  PLREFL
         (NPOINT(1,K) NPOINT(2,K), K=1,NPPLG)
          DO 21 I=1,NR1ST
             (XPOL(I,J) YPOL(I,J),  J=1,NRPLG)
 21       CONTINUE
       ELSEIF (NLFEM) THEN
         XPCOR  YPCOR  ZPCOR
         NRKNOT
         (XTRIAN(I),I=1,NRKNOT)
         (YTRIAN(I),I=1,NRKNOT)
         DO ITRI=1,NR1ST
           I,NVERT(1,ITRI),NVERT(2,ITRI),NVERT(3,ITRI),
             NEIGH(1,ITRI),NSIDE(1,ITRI),IPROP(1,ITRI),
             NEIGH(2,ITRI),NSIDE(2,ITRI),IPROP(2,ITRI),
             NEIGH(3,ITRI),NSIDE(3,ITRI),IPROP(3,ITRI)
         ENDDO
       ENDIF
    ENDIF
 ENDIF
 ** 2B. y- or poloidal grid surfaces
 NLPOL
 NLPLY  NLPLA  NLPLP
 NP2ND  NPSEP  NPPLA  NPPER
 IF (INDGRD(2).LE.5) THEN
    YIA  YGA  YAA  YYA
 ENDIF
 ** 2C. z- or toroidal grid surfaces
 NLTOR
 NLTRZ  NLTRA  NLTRT
 NT3RD  NTSEP  NTTRA  NTPER
 IF (INDGRD(3).LE.5) THEN
    ZIA  ZGA  ZAA  ZZA  ROA
 ENDIF
 ** 2D. mesh multiplication
 NLMLT
 NBMLT
 IF (NBMLT GT 1) (VOLCOR(NM), NM = 1,NBMLT)
 ** 2E. additional cells outside standard mesh
 NLADD
 NRADD
 IF (NRADD.GT.0) (VOLADD(NM), NM = 1,NRADD)

Meaning of the Input Variables

INDGRD: This index controls the meaning of input variables of different standard grid options.
INDGRD(1) for the radial (or x-direction) grid
INDGRD(2) for the poloidal (or y-direction) grid
INDGRD(3) for the toroidal (or z-direction) grid

Data for first standard mesh: radial or x grid RSURF

NLRAD= .TRUE.

A radial or x grid is defined. Otherwise the complete sub-block 2A may be omitted. Depending upon the logical parameters in the next input card the "geometry - level" variable LEVGEO is set internally.

LEVGEO=1: cartesian coordinates x (and y)
LEVGEO=2: polar coordinates r (and θ)
LEVGEO=3: general curvilinear coordinates: a full 2D mesh (polygonal coordinate lines) is used in the x - y plane. Grid cuts are permitted in the y-direction.
LEVGEO=4: a 2D "finite element" mesh of triangles is used in the x - y plane
LEVGEO=5: a 3D "finite volume" mesh of tetrahedrons used
LEVGEO=6: a general, user defined geometry block is used. All geometrical calculations are performed in problem specific routine VOLUSR, TIMUSR, ...etc.

If NLRAD=.FALSE., then no spatial grid is defined and the default geometry level LEVGEO = 1 is used. Volume discretisation may still be achieved "by hand" by defining "additional surfaces" (input block 3b) and appropriate cell number switching.

INDGRD(1)

= 1: standard grid option; the input parameters are used as described below.
= 2: As for INDGRD(1)=1, but in case of LEVGEO = 2 by this option a radial grid is defined such that the spacing of radial surfaces is not equidistant at some poloidal position, but instead such that the area enclosed between two neighboring surfaces is kept constant. The input variable RGA is irrelevant in this case.
= 3,4,5: not in use
= 6: Data for NR1ST radial surfaces are read from the work array RWK. These data must have been written onto RWK in the user supplied subroutine INFCOP. By this option grid parameters may directly be transferred into EIRENE from other files, e.g. from plasma transport codes (see below: section 2.14. Data for interfacing routine INFCOP).

One and only one of the next 6 logical variables must be .TRUE.

NLSLB = .TRUE.: Geometry level: LEVGEO = 1
Cartesian geometry, the x co-ordinate is discretized by setting
x_I=RSURF(I) I=1,NR1ST .
Furthermore, the flux-surface labeling grid RHOSRF(I) is identical with the grid RSURF(I).
NLCRC = .TRUE.: Geometry level: LEVGEO = 2
Cylindrical or toroidal geometry, mesh of concentric, circular surfaces. Polar coordinates are used in the x-y plane. The radial surfaces are given by r² = x² + y² = const. and r is discretized by setting r_I=RSURF(I) I=1,NR1ST. Furthermore RHOSRF(I) = √{ [(AREA)/(π)] } where AREA is the area inside surface number I. Thus for this option one has again: RHOSRF(I) = RSURF(I), I=1,NR1ST .
NLELL = .TRUE.: Geometry level: LEVGEO = 2
Mesh of nested, but not necessarily concentric or confocal elliptical flux surfaces. The equation for the "radial" surface is (x − EP)² + ( y/EL )² = r². r is discretized by setting r_I=RSURF(I) I=1,NR1ST.
EP and EL may vary with r. These parameters are stored in the arrays EP(I),EL(I), I=1,NR1ST which now are used in addition to RSURF to define one surface.
RHOSRF: as in NLCRC option.
Note: RHOSRF and RSURF may differ in this case.
NLTRI = .TRUE.: Geometry level: LEVGEO = 2
to be written: triangularity in mesh of nested closed algebraic surfaces
NLPLG = .TRUE.: Geometry level: LEVGEO = 3
The mesh in the x-y plane is described by NR1ST polygonal arcs of length NRPLG each. A polygon may consist of several "valid" and "invalid" parts (to account for "grid cuts" in CFD meshes). The "invalid" parts of a polygon are not seen by test particles and are allowed for in EIRENE only in order to facilitate index mapping in case of runs coupled to plasma transport models, which resort to computer generated meshes including grid cuts.
The polygons must not intersect each other.
In this case RHOSRF(1)=0., and RHOSRF(I) is the area enclosed by polygon number 1 and polygon number I.
NLFEM = .TRUE.: Geometry level: LEVGEO = 4
The mesh in the x-y plane consists of NR1ST triangles, composed from NRKNOT knots.
In this case a flux surface labeling gird RHOSRF is not defined.
NLTET = .TRUE.: Geometry level: LEVGEO = 5
3D discretisation of volume by tetrahedrons. For this grid option please make contact to the authors.
NLGEN = .TRUE.: Geometry level: LEVGEO = 6
Arbitrary geometrical configuration. Mesh consists of NR1ST arbitrarily shapes cells (in any dimension). Particle tracing routines must be provided by user (VOLUSR, SAMUSR, TIMUSR, LEAUSR)
NR1ST: Number of grid-points in the radial (or x-direction) standard mesh
if NR1ST ≤ 1, no radial (or x-direction) standard mesh is defined.
NRSEP: This flag is active for LEVGEO = 1 or LEVGEO = 2. Otherwise it is irrelevant.
The standard mesh is composed by two equidistant x- or radial grids of co-ordinate surfaces with different grid density. There are NR1ST-NRSEP+1 grid-points in the first, and NRSEP grid-points in the second part. The grid-point RSURF(NR1ST-NRSEP+1) belongs to both parts.
RIA: left endpoint of standard grid (internally set ≥ 0 if LEVGEO =2); RSURF(1)=RIA
RGA: boundary separating first and second part of standard grid with different grid-point densities; RSURF(NR1ST-NRSEP+1)=RGA
RAA: right endpoint of standard grid: RSURF(NR1ST)=RAA
RRA: if RRA > RAA, one additional zone is defined
RSURF(NR1ST+1)=RRA, and NR1ST is increased by one.
(irrelevant, if RRA ≤ RAA)

if NLELL = .TRUE. :

EPIN: Value of EP(r) for cylindrical co-ordinate surface number 1 with r₁=RIA (see: NLCRC = .TRUE. option above)
EPOT: Value of EP(r) for cylindrical surface number NR1ST with
r_NR1ST=RAA
EPCH: Value of EP(r) for cylindrical surface number NR1ST+1 with
r_NR1ST+1=RRA
(irrelevant if RRA ≤ RAA)
EXEP: The variation of the "shift function" EP(r) with r is given by
EP(r) = EPIN + ( [(r − RIA)/(RAA −RIA)] ) ^EXEP·(EPOT − EPIN)
ELIN: Value of EL(r) for cylindrical surface number 1 with
r₁=RIA (see: NLELL = .TRUE. option)
ELOT: Value of EL(r) for cylindrical surface number NR1ST
with r_NR1ST=RAA
ELCH: Value of EL(r) for cylindrical surface number NR1ST+1 with
r_NR1ST+1=RRA
(irrelevant for RRA ≤ RAA)
EXEL: The variation of the "ellipticity function" EL(r) with r is given by
EL(r) = ELIN + ( [(r − RIA)/(RAA −RIA)] ) ^EXEL·(ELOT − ELIN)

if NLTRI = .TRUE. :

to be written

if NLPLG = .TRUE. :

NR1ST: number of polygons for discretisation in "radial" or x direction
NRPLG: Number of points per polygon
NPPLG: Number of valid parts on each polygon. Each polygon is described by the x and y co-ordinates of NRPLG points. It is not necessary that all this points are used for the polygon. One can cut the polygon into several valid parts interrupted by parts which are not seen by the test particles. ( Default : NPPLG = 1 ). This option facilitates the use of 2-d computer generated meshes which contain topological grid cuts.
XPCOR,YPCOR: shift whole mesh by that vector in x,y-plane
RFPOL: if RFPOL > 0., one additional polygon zone is defined, at a distance RFPOL outside the polygon NR1ST, and then NR1ST is increased by one.
(irrelevant, if RFPOL ≤ 0.)
NPOINT(1,J): Index of the first point of the valid part number J (same for each radial polygon)
( Default : NPOINT(1,1) = 1 )
NPOINT(2,J): Index of the last point of the valid part number J (same for each radial polygon)
( Default : NPOINT(2,1) = NRPLG )
XPOL(K,I): x-co-ordinate of the polygon point number K on polygon number I
YPOL(K,I): y-co-ordinate of the polygon point number K on polygon number I

if NLFEM = .TRUE. :

NR1ST: number of triangles for discretisation in x-y plane
XPCOR,YPCOR: shift whole mesh by that vector in x,y-plane
NRKNOT: There are NRKNOT knots, by which the triangles are defined
XTRIAN,YTRIAN: x and y co-ordinates of the knots, respectively.
NVERT(I,ITRI): Each triangle ITRI is defined by 3 points P₁, P₂ and P₃ from the set of NRKNOT knots.
NVERT(I,ITRI) is the number of point P_I (I=1,2,3) in the set of knots for triangle ITRI.
NEIGH(I,ITRI): The three sides of each triangle S₁, S₂, S₃ are defined such that S₁ connects P₁ and P₂, S₂ connects P₂ and P₃ and S₃ connects P₁ and P₃.
NEIGH(I,ITRI) is the number of the neighboring triangle at side S_I for I=1,2,3.
NSIDE(I,ITRI): NSIDE(I,ITRI) is the number (1,2 or 3) of the side of the neighboring triangle, which corresponds to side S_I of the triangle ITRI to be written
IPROP(I,ITRI): ISTS = ABS(IPROP) is the integer, by which a particular surface property is assigned to side I of the triangle ITRI. ISTS=0 stands for default grid option (transparent surface, cell indexing is done automatically). Otherwise ISTS is the number of an additional (ISTS < NLIMI) or non-default standard (NLIM < ISTS < NLIM+NSTSI) surface option as read in sub-blocks 3a or 3b, respectively. By default the normal vector of each side of a triangle points out of the triangle. In case IPROP < 0, this vector points into the triangle. This is relevant at surface options with ILSIDE ≠ 0.

if NLTET = .TRUE. :

Documentation not available here. Make contact with the authors.

if NLGEN = .TRUE. :

NR1ST: number of cells in otherwise arbitrary mesh option NLGEN.

Data for second standard mesh: poloidal or y grid PSURF

NLPOL: A poloidal or y grid is defined. Otherwise the complete block 2B may be omitted and the volume averaged tallies are then automatically integrated over this co-ordinate.
INDGRD(2) = 1: standard grid option
INDGRD(2) = 2,3,4,5,6: not in use (default: INDGRD(2)=1)
NP2ND: Number of grid-points in y- or poloidal direction
YIA,YGA,YAA,YYA: for the LEVGEO=1 option: the y grid PSURF is defined in the same way as the x grid, using the parameters YIA,YGA,.... (cm) instead of RIA,RGA,....

for the LEVGEO=2 option: the θ grid PSURF is defined in the same way as the r grid, using the parameters YIA,YGA,...(poloidal angle in degree) instead of RIA,RGA,...

for all options LEVGEO > 2: this input card is irrelevant.

Defaults: (needed for scaling, i.e., cell volumes):

YIA=0., YGA=0., YAA=1., YYA=1. in case of LEVGEO=1,

and YIA=0., YGA=0., YAA=360., YYA=360. in case of LEVGEO=2.

Data for third standard mesh: toroidal or z grid TSURF

NLTOR: A toroidal or z grid is defined. Otherwise the complete block 2C may be omitted and the volume averaged tallies are then automatically integrated over this co-ordinate.
In case NLTOR = TRUE, sub-block 2C must be read.
INDGRD(3) = 1: standard grid option
INDGRD(3) = 2,3,4,5,6: not in use (default: INDGRD(3)=1)

One and only one of the next four boolean variables must be TRUE

NLTRZ=TRUE: cylindrical approximation is used, i.e., TSURF is a grid in z direction. The co-ordinate surfaces are given by z=TSURF(L)
(Default: NLTRZ = TRUE)
NLTRA = TRUE: toroidal approximation is used. The coordinate line is a polygonal approximation of a circle, or a segment thereof.
In case NLTOR = TRUE, the 3rd grid TSURF is a grid of toroidal angles. The toroidal segment (or the full torus) is approximated by NT3RD-1 straight cylinders.
In case NLTOR = FALSE, there are NTTRA toroidal periodicity boundaries, such that the toroidal segment is approximated, again, by NTTRAM=NTTRA-1 straight cylinders, but without toroidal resolution in the results. (Note: If NLTOR, NTTRA=NT3RD, internally.)
The torus axis of the entire mesh can be shifted in radial direction by adding a radial offset ROA (see below) to the x coordinates of the poloidal mesh (RSURF,PSURF).
The radial shift of the poloidal mesh RMTOR of the approximated torus is computed from ROA such that the volume inside radial surface NR1ST is exactly equal to the volume of an exact torus with poloidal cross section defined by the shifted radial grid (first standard mesh: RSURF).
Due to the approximations made by defining a torus by NTTRAM = NTTRA - 1 straight cylinders, this condition is fulfilled only approximately for the other radial surfaces. RMTOR converges to ROA with increasing NTTRA. NTTRA ≅ 30 is already a very good approximation.
(Default: NLTRA = FALSE).
NLTRT = TRUE: torus co-ordinates R,PHI,THETA. Presently being developed for NLSLB,NLPLG and NLTRI options. Not ready for use.
NT3RD: Number of grid-points in z- or toroidal direction
(default: NT3RD = 1, i.e. no grid is defined)
NTTRA: only needed in case NLTRA and .NOT.NLTOR. See above.
ZIA,ZGA,ZAA,ZZA,ROA: The 3rd grid ZSURF is defined in the same way as the x grid, using the parameters ZIA,ZGA,.... (cm) instead of RIA,RGA,....
In case of NLTRZ = TRUE , a z-grid is defined. ROA is irrelevant.
In case of NLTRA = TRUE , ZIA and ZAA are toroidal angles (in degrees). A grid of toroidal angles is defined. For example use ZIA=0.0 and ZAA=360.0 for a full torus. In case NLTOR this grid also defines the toroidal resolution. If .NOT.NLTOR, then periodicity at the endpoints ZIA and ZAA is automatically enforced.
(Defaults: NLTRA = FALSE, ZIA = 0 , ZAA = 1)

Data for mesh multiplication

NLMLT: the complete Standard mesh data are copied NBMLT times
NBMLT: Number of identical copies of the standard mesh. Transition from one such mesh (called "block" in EIRENE) into another one has to be defined by transparent additional surfaces (see block 3B)
VOLCOR: The volumes of all cells of the standard mesh as computed by EIRENE (in subroutine. VOLUME) are multiplied by one common factor VOLCOR for each "block".

Data for additional zones outside the Standard mesh

NLADD: There are cells in the computational volume, which are defined through "additional surfaces" as cell boundaries. E.g. a standard mesh (is there is one) is augmented by "additional cells" in this case. If NLADD = FALSE , the complete block 2D may be omitted and the defaults are used.
NRADD: Number of additional zones. The cell indexing along test flights in these zones has to be specified explicitly by making use of the ILSWCH, ILCELL parameters (block 3B)
(Default: NRADD = 0 ).
VOLADD: Volume (cm⁻³) of each additional zone as seen by the test-particles.

2.2.1 Mesh Parameters

As stated above, the computational volume in an EIRENE run is subdivided into NSBOX cells,

NSBOX = (NR1ST · NP2ND · NT3RD) · NBMLT + NRADD .

The can mesh consist two different parts.

The first is the so called "standard mesh" part, defined by the grids RSURF, PSURF and TSURF. There may be NBMLT copies of that mesh, separated by additional surfaces (see below, section 2.3.2 (default: NBMLT=1).

The second part is the so called "additional cell region". Up to NRADD additional cells, of arbitrary shape, are defined by their boundaries: the "additional surfaces" defined in section 2.3.2.

NR1ST is the number of surfaces in the 1st (radial -or x-) grid. There are NR1STM = NR1ST - 1 cells.

The cell index of a particular cell is called NRCELL.

1 ≤ NRCELL ≤ NR1STM

NP2ND is the number of surfaces in the 2ND (poloidal -or y-) grid. There are NP2NDM = NP2ND - 1 cells.

The cell index of a particular cell is called NPCELL.

1 ≤ NPCELL ≤ NP2NDM

NT3RD is the number of surfaces in the 2RD (toroidal -or z-) grid. There are NT3RDM = NT3RD - 1 cells.

The cell index of a particular cell is called NTCELL.

1 ≤ NTCELL ≤ NT3RDM

For the cells of the standard mesh, NACELL = 0

For the cells in the "additional cell region" :

1 ≤ NACELL ≤ NRADD

and

NRCELL = 0, NPCELL = 1, NTCELL = 1

The last radial (or x-) cell (number NR1ST), for any NPCELL, NTCELL inside the standard mesh, is not a real geometrical cell. It is used to store the radial (or x-) average of the volume averaged tallies.

Likewise, the last poloidal (or y-) cell (number NP2ND) is not a real geometrical cell, but used to store the poloidal (or y-) average of the volume averaged tallies (for any NRCELL, NTCELL inside the standard mesh)

NTCELL = NT3RD: analogically.

A particular cell can be specified in one of two possible ways:

Either by giving the 5 cell numbers:

NRCELL, NPCELL, NTCELL, NBLOCK, NACELL

or, alternatively, by giving the position in the 1-dimensional cell array

NCELL

The relation between these two options is:

NCELL = NRCELL + ((NPECLL-1)*NR1P2) + (NTCELL-1))*NP2T3 + NBLCKA

with

NBLCKA = (NBLOCK -1) * NSTRDT +NRADD

2.3 The Input Block for Surfaces

Grid surfaces may be assigned special properties (e.g.: reflecting, absorbing, periodicity, modified cell index switching, etc...). Their definition is described in subsection 2.3.1.

In addition to these grid surfaces also additional surfaces can be seen by the histories (vacuum boundaries, special surfaces for scoring fluxes, diagnostic surfaces, ....). Such surfaces can be general linear or second order surfaces in 3D space. Their definition and properties are described in subsection 2.3.2.

2.3.1 The Input Block for "Non-default Standard Surfaces"

General remarks

Some of the "standard grid surfaces" RSURF, PSURF and TSURF may be defined as reflecting or absorbing (rather than the default transparent character of co-ordinate surfaces), or as transparent but with other switching features (cell indexing, transition into additional cell region etc.) than the default settings. E.g., surface averaged tallies (see block 2.11, end of section 3.2 and Table 5.3) are only evaluated on either "Additional Surfaces" or on these "Non default Standard Surfaces", whereas they are not computed for those standard surfaces which are not explicitly identified in this block.

The geometrical properties of these "Non default Standard Surfaces" are specified here. This is described below.

Particle-surface interaction models

Although the (local) data for particle-surface interaction models (PSI-models) for each specific surface can be read in this input block, their meaning is described in block 2.6 together with the globally valid input data for particle-surface interactions.

*** 3A. Data for "Non-default Standard Surfaces"

 NSTSI
 DO 31 ISTSI=1,NSTSI

unless otherwise stated, N = NLIM+ISTSI is the index of the following data:

    TXTSFL
    ISTS IDIMP INUMP(ISTSI,IDIMP) IRPTA1 IRPTE1 IRPTA2 IRPTE2 IRPTA3 IRPTE3
    ILIIN ILSIDE ILSWCH ILEQUI ILTOR ILCOL ILFIT ILCELL ILBOX ILPLG

optional for non transparent surfaces (ILIIN > 0):

The next 3 (or 4) lines comprise the block for local particle-surface interaction data (in analytic terminology: the boundary condition at this surface element). If they are omitted, the default particle-surface model is activated for this particular surface element, see section 2.6.

    ILREF  ILSPT ISRS  ISRC
    ZNML EWALL EWBIN TRANSP(1,N) TRANSP(2,N) FSHEAT
    RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
    RECYCS  RECYCC  SPTPRM  (this line may be omitted,
                             then: default sputter model, see \ref{sec2.6})
 31 CONTINUE

also, optional for non-transparent surfaces (ILIIN > 0),

read a surface interaction model identifier to link one of the surface "local reflection models" defined in block 2.6 below to this current surface. The presence of such a link is identified by the string ′SURFMOD′, followed by a name ′modname′. Later, in input block 2.6, a "local reflection model" block with that name ′modname′ must be included. This option allows quick changes of particle-surface interaction parameters, affecting many surfaces at once, by changes in just one location of the input file.

   SURFMOD_modname

Meaning of the Input Variables for "non-default standard surfaces"

NSTSI

Total number of non-default standard surfaces that do not act as prescribed by the default transparent standard co-ordinate surface model.

TXTSFL

Text to characterize a surface ("name of the surface") on the printout file.

ISTS

irrelevant; labelling index

IDIMP

flag to identify mesh from which this particular surface is chosen.

= 1: surface from the x- (radial) standard mesh (RSURF) Note that for the unstructured grid options NLTRI and NLTET (i.e. for the 2D triangular grid option or for the general 3D grids of tetrahedra) all surfaces are referred to as 1st grid (x or radial) surfaces, by abuse of language.
= 2: surface from the y- (poloidal) standard mesh (PSURF)
= 3: surface from the z- (toroidal) standard mesh (TSURF)

INUMP

Number of the surface in mesh RSURF, PSURF or TSURF respectively

IRPTA,IRPTE

Only a subregion of the surface acts by the "non-default options" specified for this particular surface. This subregion is defined by these flags.

If JMP is a surface from the first mesh, then IRPTA2→IRPTE2 and IRPTA3→IRPTE3 are the surface index ranges of the 2nd and 3rd mesh, respectively, for which this surface acts as non-default surface. IRPTA1 and IRPTE1 are irrelevant.

If JMP is a surface from the 2nd mesh, then IRPTA1→IRPTE1 and IRPTA3→IRPTE3 are the surface index ranges of the 1st and 3rd mesh, respectively, for which this surface acts as non-default surface. IRPTA2 and IRPTE2 are irrelevant.

If JMP is a surface from the 3rd mesh, then IRPTA1→ IRPTE1 and IRPTA2→ IRPTE2 are the surface index ranges of the 1st and 2rd mesh, respectively, for which this surface acts as non-default surface. IRPTA3 and IRPTE3 are irrelevant.

The third card ILIIN, ..., ILPLG is identical to the corresponding card for "additional surfaces", see block 2.3.2 below, "Input Data for Additional Surfaces". One exception is the flag ILSIDE, which controls the sign of the surface normal vector (hence the orientation of a surface). In case of unstructured "standard grids" NLTRI or NLTET (triangles in 2D and tetrahedrons in 3D) there is no well defined default surface orientation and the flag ILSIDE is irrelevant in such cases.

The input data in the block for the local reflection model are described below, see 2.6 : "Input Data for Particle-Surface Interaction Models".

2.3.2 Input Data for "Additional Surfaces"

General remarks

Internally each additional surface is defined by an algebraic equation and some algebraic inequalities specifying the boundary of that surface, i.e., that part of the surface which is seen by the test particles. By changing the sign of all coefficients in the algebraic equation the orientation of the surface normal vectors are reversed.

The intersection with the nearest surface along a trajectory is found by checking surfaces in the sequence of their input (in subroutine TIMEA of the geometry block GEO3D). This must be taken into consideration, if there are two parts of one surface specified by different boundaries and with non-empty intersection. In this case always the surface later in the input sequence is seen by the test particles.

One can define two distinct plane surfaces as one surface of second order, provided that this is compatible with the options available for surface boundaries.

We will refer to positive and negative directions of flights for particles intersecting a surface. By "positive" it is meant that the angle between the trajectory of the particle and the surface normal at the point of intersection is less than 90 degrees, "negative" has the opposite meaning.

speeding up geometry calculations:

Checking surfaces along a test flight can be abandoned depending upon the current position of a test particle.

For particles in cell no. NCELL all surfaces MSURF with IGJUM3(NCELL,MSURF)=1 are not checked. For particles on surface MS all other surfaces MSURF with IGJUM1(MS,MSURF)=1 are not checked. These matrices IGJUM1 and IGJUM3 can be set either in this block, see CH1 cards and CH3 cards described below, or in some user-specified problem specific routines (e.g. GEOUSR).

ignorable coordinates:

Surfaces with ignorable coordinates are specified either by setting the corresponding coefficients in the surface equation to zero, or, in case of 2-point, 3-point, 4-point of 5-point options, by setting the corresponding coordinates of the points to -1.D20 and +1.D20, respectively.

coordinate systems, transformations

All surfaces are specified in cartesian coordiantes.

In case of NLTRA (toroidal geometry approximation), surfaces must be defined in the local coordinate system of toroidal cells centered at the "magnetic axis", (xa,ya) = (RMTOR,0.), i.e., excluding the major radius of the torus.

To facilitate input of the geometrical data, each single (or a set of) surface(s) can be transformed by a Cartesian mapping after specification of the surface coefficients and boundaries in some "convenient" coordinate frame, by including "TRANSFORM cards" at the end of an input block for a particular surface. One such "TRANSFORM-deck" can act on an entire range of surfaces, but only on those which have been read previous to the "TRANSFORM-deck". Hence, in order to transform the entire set of additional surfaces by one deck, this deck must come after the last surface (no. NLIMI, see below).

Particle-surface interaction models

Although the (local) data for particle surface interaction models for each specific additional surface can be read in this input block, their meaning is described in block 2.6 together with the globally valid input data for particle surface interactions.

The Input Block

*** 3B. Data for Additional Surfaces

 NLIMI
 ``CH0-cards" (may be omitted)
 (format: 3A,69A)

arbitrary number of strings ±n/m, separated by blanks. n and m must be integer variables 1 ≤ n,m ≤ NLIMI.

 DO  ILIMI=1,NLIMI
   TXTSFL
   ``CH1-cards" (may be omitted)
   (format: 3A,69A)

arbitrary number of strings ±n/m, separated by blanks. n and m must be integer variables 1 ≤ n,m ≤ NLIMI.

   ``CH2-cards" (may be omitted)
   (format: 3A,69A)

arbitrary number of strings ± n/m, separated by blanks. n and m must be integer variables 1 ≤ n,m ≤ NLIMI.

   RLBND  RLARE  RLWMN  RLWMX
   ILIIN ILSIDE ILSWCH ILEQUI ILTOR ILCOL ILFIT ILCELL ILBOX ILPLG
   IF (RLBND.LT.0.) THEN
     A0LM  A1LM  A2LM  A3LM  A4LM  A5LM
     A6LM  A7LM  A8LM  A9LM

let RLBND = -KL, then first read L cards:

     ALIMS  XLIMS  YLIMS  ZLIMS

and then read K blocks next:

     ALIMS0  XLIMS1  YLIMS1  ZLIMS1  XLIMS2  YLIMS2
     ZLIMS2  XLIMS3  YLIMS3  ZLIMS3
   ENDIF
   IF (RLBND.EQ.0.) THEN
     A0LM  A1LM  A2LM  A3LM  A4LM  A5LM
     A6LM  A7LM  A8LM  A9LM
   ENDIF
   IF (RLBND.GT.0..AND.RLBND.LT.2.) THEN
     A0LM  A1LM  A2LM  A3LM  A4LM  A5LM
     A6LM  A7LM  A8LM  A9LM
     XLIMS1  YLIMS1  ZLIMS1  XLIMS2  YLIMS2  ZLIMS2
   ENDIF
   IF (RLBND.GE.2.) THEN
     P1(1,..)  P1(2,..)  P1(3,..)  P2(1,..)  P2(2,..)  P2(3,..)
     IF (K .GT. 2) THEN
         P3(1,..)  P3(2,..)  P3(3,..)  P4(1,..)  P4(2,..)  P4(3,..)
         IF (K .GT. 4) THEN
             P5(1,..)  P5(2,..)  P5(3,..)
         ENDIF
     ENDIF
   ENDIF

optional for non transparent surfaces (ILIIN > 0):

    ILREF  ILSPT ISRS  ISRC
    ZNML EWALL EWBIN TRANSP(1,N) TRANSP(2,N) FSHEAT
    RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
    RECYCS  RECYCC  SPTPRM  (this line may be omitted,
                             then: default sputter model,
                             see \ref{sec2.6})

also optional for non transparent surfaces (ILIIN > 0):

   SURFMOD_modname

optional: read one or several blocks of five cards each for orthogonal mapping. The presence of each such block is identified by the first card of that block containing the string ′TRANSFORM′ followed by the other four cards:

   TRANSFORM
   ITINI,ITEND
   XLCOR,YLCOR,ZLCOR
   XLREF,YLREF,ZLREF
   XLROT,YLROT,ZLROT,ALROT

 ENDDO   (end of ILIMI loop)

Meaning of the Input Variables for additional surfaces

CH0 n1/m1 n2/m2 ...: surfaces from the range n1 to m1, n2 to m2, ..., are ignored by EIRENE. Specifying a surface in such a CH0-card is identical to taking it out from the input file. It may be more convenient in some cases, however, to use the CH0 option, because the labelling index of the remaining valid surfaces is not altered then. No input is read for these surfaces, and the input segment for the next valid surface (identified by the string '*text') is read next.
NLIMI: number of surfaces in the input block
TXTSUR: Text to identify a surface ("name of the surface") on the printout file
CH1(ILIMI) n1/m1 n2/m2 ...: surfaces from the range n1 to m1, n2 to m2, ..., are considered invisible for a particle located on this current surface ILIMI. Intersection of trajectories starting from surface ILIMI with those "invisible" surfaces is not checked.
CH2(ILIMI) n1/m1 n2/m2 ...: Only for second order surfaces ILIMI, n1-m1, ... The first of the two possible intersections is ignored for particles located on surface no. ILIMI.

General Data for Surfaces

RLBND

flag for different options to define the boundary of the surface

RLARE

Area (in cm²) of the surface element which is seen by the test particles. (Default: 666.0) (needed only for scaling of non default surface averaged tallies) If RLARE is not specified here, (i.e., if a value less than or equal to zero is read) then EIRENE tries to evaluate this area itself. For some surfaces this is still not possible automatically.

RLWMN

lower weight limit for space weight window for particles crossing the surface in positive direction. (not in use)

RLWMX

upper weight limit for space weight window for particles crossing the surface in positive direction. (not in use)

ILIIN

defines the type of surface

> 0

non-transparent surface

= 1: reflecting, partly or purely absorbing surface.
local reflection model has to be specified unless default model is to be used; all surface tallies (see Table 5.3) are updated and a switch can be operated.
= 2: purely absorbing surface ;. surface tallies for incident fluxes (i.e., POT... and EOT... tallies in Table 5.3) are updated and the particle history is stopped then.
= 3: mirror for incident test particles. I.e., specular reflection for neutral test particles, and for charged test particles the sign of the velocity component parallel to the B-field is reversed.
=m4: periodicity surface, with regard to x, y, or z coordinate, depending upon whether this surface is a standard x, y, or z grid surface, respectively. Move particle to x / radial grid surface no. m, m integer (or to y / poloidal or to z / toroidal surface no. m, respectively) and continue track from there with otherwise identical particle parameters.
This option is currently implemented only for Cartesian grids (NLSLB and NLTRZ) and for non-default standard grid surfaces only.

The periodicity options are not fully implemented yet. Please contact the author for the current status of your particular version.

≤ 0

transparent surface (for example: hole in one of the other "additional surfaces"). Particle and energy fluxes onto and from these surfaces do not contribute to global balances.

= 0: Particle history is not interrupted
No surface tallies are updated, no switches can be operated. Fastest option.
= -1: Particle history will be stopped and restarted. A switch can be operated. I.e., this surface is used only for switching (see below: ILSWCH) or re-initializing the particle's track at the point of intersection. No surface tallies are updated.
= -2: as -1, and, additionally:
if a particle is crossing the surface in the positive direction, (one sided-) surface tallies are updated, e.g., by default: partial particle and energy currents J⁺ (Amp) and K⁺ (Watt). These are stored in the POT... and EOT... tallies of Table 5.3. If the particle crosses the surface in the direction opposite to the surface normal, then negative partial particle and energy currents J⁻ (Amp) and K⁻ (Watt) are updated. These are stored in the PRF... and ERF... tallies of Table 5.3.
= -3: Net currents (e.g. J⁺ − J⁻), are evaluated, and stored on the POT... and EOT... tallies (see Table 5.3). The PRF... and ERF... tallies are empty for these surfaces.
≤ -4: Not in use. Currently: same as ILIIN=-2 option.

ILSIDE

= 0: both sides of the surface act as described by ILIIN option (default).
= 1: particles incident on the surface in the negative direction will be absorbed (i.e., ILIIN = 2 option from that side).
= 2: particles incident on the surface in the negative direction will be killed and the message
"ERROR IN ADDCOL"
or
"ERROR IN STDCOL"
will be printed. The contribution of these particles to the particle- and energy flux balances will be called PTRASH and ETRASH respectively. This option should be used for geometry testing whenever the user expects particles incident only from one side.
= 3: particles incident on the surface in the negative direction will not see the surface, i.e., this surface acts like a (semi) transparent surface (ILIIN = 0 option) from that side.
= -1: as 1, but with the opposite direction of the surface normal
= -2: as 2, but with the opposite direction of the surface normal
= -3: as 3, but with the opposite direction of the surface normal

ILSWCH

= IJKLMN, i.e. six digits I, J, K, L, M and N

= 0

no switch is operated

N

EIRENE flag ITIME

N = 1: The calculation of the step sizes in the standard mesh is abandoned for a particle which crosses the surface in the positive direction, and is reactivated, if the particle strikes in the negative direction
N = 2: as 1, but with the direction of the surface normal reversed for this option.

M

EIRENE flag IFPATH

M = 1: Abandon the calculation of the collision rates (entry into the vacuum) for a particle which is striking the surface in the direction of the surface normal. For particles incident from the other direction, evaluation of collision rates is reactivated.
M = 2: as 1, but with the direction of the surface normal reversed.

L

EIRENE flag IUPDTE

L = 1: Abandon the updating of volume-averaged tallies for a particle which is striking the surface in the direction of the surface normal. For particles incident from the other direction, updating of volume averaged tallies is reactivated.
L = 2: as 1, but with the direction of the surface normal reversed for this option.

I,J,K flags for switching cell numbers at transition into a different mesh cell.

K

for particles in an additional cell, i.e., not in one of the "standard mesh" blocks:

K = 1: Increase the actual additional cell number NACELL for a particle striking the surface in the direction of the surface normal by ILACLL. Decrease NACELL by ILACLL if the particle is striking in the negative direction. Specification of ILACLL is via the input variable ILCELL, see below.
K = 2: as K = 1, but with the direction of the surface normal reversed for this option.

for particles inside the "standard mesh ", i.e., not in the "additional cell region"

K = 1: Increase the standard mesh block number NBLOCK for a particle striking the surface in the direction of the surface normal by ILBLCK. Decrease NBLOCK by ILBLCK if the particle is striking in the negative direction. Specification of ILBLCK is via the input variable ILCELL, see below.
K = 2: as K = 1, but with the direction of the surface normal reversed.

J

for particles at the boundary between "additional" and "standard" mesh regions.

J = 1: entrance into standard mesh, block no.
NBLOCK = ILBLCK
or exit from standard mesh into additional cell
NACELL = ILACLL.
Specification of ILACLL and ILBLCK is via the input variable ILCELL, see below. If ILACLL = 0, then no switch to additional cell is operated. (E.g.: for surfaces which are reflecting from this side).
J = 2: as J = 1. The direction of the surface normal does not matter here.

I

similar to J-flag, i.e., for transitions between standard and additional meshes, but different cell number switching.

I = 1: Entrance into standard mesh, block no.
NBLOCK = NACELL+ILBLCK,
if the particle is striking in the positive direction, or
NBLOCK = NACELL-ILBLCK,
if the particle is striking in the negative direction. Exit from standard mesh into additional cell
NACELL = NBLOCK+ILACLL,
for a particle striking the surface in the positive direction, or
NACELL = NBLOCK-ILACLL,
if the particle is striking in the negative direction.
Specification of ILACLL and ILBLCK is via the input variable ILCELL, see below.
I = 2: as I = 1, but with the direction of the surface normal reversed.

If a test particle history starts from a surface (NLSRF option), then ILSWCH acts as if this particle had struck the surface prior to the birth process in the positive direction. This default setting is only available for ILSIDE ≠ 0 and can (must) be overruled by the SORIFL flag (see section 2.7), e.g. if a surface source needs to be defined on a surface with ILSIDE = 0.

ILEQUI

The algebraic equations for the surfaces J and IABS(ILEQUI(J)) will be described by exactly the same coefficients (up to a common sign, if ILEQUI(J) .lt. 0). For example a triangle can be specified by the three corners and another part of the same plane surface can be specified directly by its algebraic coefficients. To avoid round-off errors one should use the ILEQUI option in such cases, in particular if surface J is a transparent "hole" in surface ILEQUI(J), or vice versa.

Default: ILEQUI = 0

ILTOR

For NLTRA option only (see block 2c):

if ILTOR > 0 :

the surface is defined with respect to the local coordinate system of the toroidal segment with "cell-number" ILTOR, hence: 1 ≤ ILTOR ≤ NTTRAM.

if ILTOR=0 :

the surface is defined with respect to any local coordinate system. I.e., the surface equations are taken to be the same in each local system.

If the surface equations are z-independent, then this surface is toroidally symmetric (within the NLTRA-approximation).

Otherwise the surface has NTTRAM-fold periodicity.

This flag is irrelevant for NLTRZ (i.e., if cylindrical coordinates are used) or for NLTRT.

Default: ILTOR = 0

ILCOL

Flag for the color that is used for plotting this surface on 2d or 3d geometry plots. If ILCOL ≤ 0 than -ILCOL is used and the surface area is filled in by that color on the 3d geometry plots.

Default: ILCOL=1

ILFIT

= MN

This option is relevant only for surfaces with one ignorable coordinate, i.e. it only works for the 2. ≤ RLBND ≤ 3. surface boundary options.

It is a tool to facilitate a neat fitting of surfaces, in particular for connecting curved and plane surfaces (i.e. to avoid particle leakage due to numerical round-off errors in the algebraic surface coefficients.

M and N (3 digits each, M may be omitted if not needed) are the numbers of the surfaces (which must have the same ignorable coordinate) the boundaries of which should match to those of the actual surface J. The boundaries of these neighboring surfaces M and N must be specified by the RLBND = 1 or RLBND = 1.5 option.

The fitting is achieved by a small automatic internal modification of the surface data P1 and/or P2 of surface number J in subroutine SETFIT. Printout of the modifications made there is activated with the TRCSUR flag (input block 11).

ILCELL

Parameter ILBLCK and ILACLL for the ILSWCH flags described above. Let ILCELL = NM, with N and M being integers with 3 digits each. Then N = ILBLCK and M = ILACLL.

ILBOX

ILPLG EIRENE can write out information for a finite element mesh generator to produce a grid of triangles for a muliply connected 2D domain with cracks and holes. The various (inner and outer) boundaries are given as polygonal lines, which are composed of selected standard grid surface segments (NLPLG option) and/or additional surfaces (2 ≤ RLBND < 3 option).

This flag identifies closed polygonal lines composed of additional surfaces given by the 2-point option and/or of standard surfaces in the x-y-plane. For example if ILPLG(I)=NN, for surfaces I = I1, I2, ...IN, (NN a positive integer) then these N surfaces form a closed polygonal line in the x-y-plane. The region inside this closed line is part of the computational domain. By a negative integer value of NN a closed polygonal region can be excluded from the computational domain, i.e., a hole in the domain is specified by these surfaces. EIRENE writes an output file appropriate for a finite element mesh generator (available from FZ-Juelich) to produce a triangular discretization of the resulting (possibly multiply connected) domain. This option can be used to discretize arbitrarily complex 2D domains with internal and external boundaries given by the additional or non-default standard surfaces.

These finite element grids can be combined with the regular grids by using the problem specific geometry routines (see section 3) or the code interfacing routines INFCOP (see section 4).

Surface coefficients and boundaries of surfaces

The surface equation:

A0LM

A1LM ·x + A2LM ·y + A3LM ·z +

A4LM ·x² + A5LM ·y² + A6LM ·z² +

A7LM ·xy + A8LM ·xz + A9LM ·yz = 0

(3.1)

The positive surface normal (n_x , n_y , n_z) depends upon the point of impact (x,y,z) and is defined by the vector

n_x

A1LM + 2 ·A4LM ·x +A7LM ·y + A8LM ·z

n_y

A2LM +A7LM ·x + 2 ·A5LM ·y +A9LM ·z

n_z

A3LM +A8LM ·x + A9LM ·y +2 ·A6LM ·z

For a plane surface the following reduction is valid:

(n_x , n_y , n_z) = (A1LM,A2LM,A3LM)

The boundary of the valid part of the surface may be described by 4 different options, depending upon the value of the flag RLBND.

RLBND = 0

No boundary inequalities specified, i.e. the whole surface is seen by the test particles.

0 < RLBND < 2

= 1: Only that part of the surface, which lies inside the right parallelepiped defined by the two vectors (XLIMS1, YLIMS1, ZLIMS1) and (XLIMS2, YLIMS2, ZLIMS2), is seen by the particles. I.e. the three inequalities
XLIMS1 ≤ x ≤ XLIMS2
YLIMS1 ≤ y ≤ YLIMS2
ZLIMS1 ≤ z ≤ ZLIMS2
are checked at the point of intersection (x,y,z).
= 1.5: Complement to RLBND = 1.
Only the surface element outside the parallelepiped is seen by the particles.

RLBND ≥ 2

In this case the surface will be defined by the input of the coordinates of at least 2 and at highest 5 points on a plane surface. If there are only 2 points, the surface is parallel to one axis. If there are 3 or more points, then the boundary of this plane surface is a closed polygon (P₁ ,..., P_n , P₁). Therefore, the correct order of points at input is relevant. The orientation of the positive surface normal vector is defined by the first 3 points, and it is given by the vector product (P3 − P1) ×(P3 − P2). Thus, the orientation can be reversed e.g. by interchanging P₂ and P₃.

2.1: plane surface parallel to z axis. The surface equation of this plane reads ax+by+c=0 with the coefficients a,b and c such that the points P₁, P₂ lie on this surface and the valid part of that surface ranges from P₁ to P₂ in the xy-plane. The z-coordinates of these two points define the boundaries in z direction
= 2.2: Complement to RLBND = 2.1
= 2.4: as RLBND=2.1 option, but with z and y exchanged. I.e., now the y coordinates of the points P₁, P₂ are the boundaries of the surface ax+bz+c=0 in y direction.
= 2.5: Complement to RLBND = 2.4
= 2.7: as RLBND=2.1 option, but with z and x exchanged. I.e., now the x coordinates of the points P₁, P₂ are the boundaries of the surface ay+bz+c+0 in x direction.
= 2.8: Complement to RLBND = 2.7
= 3: plane triangle defined by the corners P₁ ,P₂ , P₃
P₁=(P1(1),P1(2),P1(3))
P₂=(P2(1),P2(2),P2(3))
P₃=(P3(1),P3(2),P3(3))
= 3.5: complement to RLBND = 3; only The plane surface outside the triangle is seen by the test particles.
= 4: plane quadrangle; surface inside the polygon
(P₁ , P₂ , P₄ , P₃ , P₁).
Here P₁, P₂, P₃ are as in the RLBND=3 option, and P₄ = (P4(1), P4(2), P4(3))
Thus this surface is the union of the triangles with vertices P₁, P₂, P₃ and P₂, P₄, P₃ respectively.
= 4.5: complement to RLBND = 4; only the part of the plane surface outside the quadrangle is seen by the test particles
= 5: plane quint-angle; surface inside the polygon (P₁, P₂, P₄, P₅, P₃, P₁) P₁, P₂, P₃, P₄ as RLBND=4, and P₅ = (P5(1), P5(2), P5(3))
= 5.5: complement to RLBND = 5; only the part of the plane surface outside the quint-angle is seen by the test particles.

RLBND < 0

-KL

The surface is bounded by L linear inequalities and by K second order inequalities.

ALIMS + XLIMS * x + YLIMS * y + ZLIMS * z ≤ 0

(L inequalities)

ALIMS0

XLIMS1 ·x + YLIMS1 ·y + ZLIMS1 ·z

XLIMS2 ·x² + YLIMS2 ·y² + ZLIMS2 ·z²

XLIMS3 ·xy + YLIMS3 ·xz + ZLIMS3 ·yz ≤ 0

(K inequalities)

The meanings of the input variables for the local reflection model are described below (see: 2.6: "Input Data for Surface Interaction Models").

The meanings of the input variables for a transformation block are:

ITINI, ITEND: The transformation defined by the next 3 cards is carried out for additional surfaces from number ITINI through ITEND The transformation is carried out as soon as this transformation deck is found. Hence, if such a transformation deck is found after "additional surface" no. IS, then one must guarantee: ITINI ≤ IS and ITEND ≤ IS.
XLCOR, YLCOR, ZLCOR: Translation:
The origin of the coordinate system is shifted by the vector
XLCOR,YLCOR,ZLCOR
XLREF, YLREF, ZLREF: Reflection:
to be written
No transformation if XLREF,YLREF,ZLREF = (0,0,0)
XLREF, YLREF, ZLREF: Rotation:
The vector XLROT,YLROT,ZLROT defines the axis of rotation.
ALROT is the angle of rotation (in degrees)
No transformation if ALROT=0 or axis of rotation= (0,0,0)

2.4 Input Data for Species Specification and Atomic Physics Module

General remarks

EIRENE can handle up to NATM "atom" species, NMOL "molecule" species, NION "test ion" species, NPHOT "photon" species (lines) and NREAC different atomic, molecular or photonic reactions between these "test particles" and the "bulk ions" or electrons. There may be up to NPLS "bulk ion" species, and one electron gas derived internally from the assumption of local charge neutrality. Amongst the heavy background particles (the "bulk ions") may be species with charge state zero, i.e., neutral particles. By abuse of language, we will refer to them as "bulk ions" as well, but meaning in this case: heavy background particles.

NATM, NMOL, NION, NPHOT, NPLS and NREAC are specified in the parameter block "PARMUSR", see: "Parameter Statements" (section 3.1 below, for versions EIRENE₂₀₀₀ and older), or are determined automatically from the input file.

Masses in EIRENE (RMASSA, RMASSM, RMASSI, RMASSP for atoms, molecules, test ions and bulk ions, respectively) are in units of the proton mass, which is taken to be 1.0073 in atomic mass units.

First EIRENE reads "non-default" atomic data (if any, i.e., if NREAC > 0) specified in so called "reaction cards", IR ≤ NREACI ≤ NREAC


IR  FILNAM  H123  REAC  CRC  MASSP  MASST  DP RMN  RMX

(format: I3,1X,A6,1X,A4,Axxx,A3,2I3,3E12.4)

from external A&M data files such as FILNAM = HYDHEL, = METHAN, or = AMJUEL, or other data files of appropriate format.

For collision reaction rates, electron cooling rates, etc., from the "ADAS database format AFD11", or any other properly tabulated data source, i.e. for FILNAM = ADAS, an additional card is read for each IR entry, see below, to identify the species and ion charge state within the ADF11 file.


IR  ADAS  H123  REAC  CRC  MASSP  MASST  DP

(format: I3,1X,A6,1X,A4,Axxx,A3,2I3,E12.4)
ELNAME, IZ

(format: 4X,A2,1X,I3)

For photonic emission and absorption data (line shapes), i.e. for FILNAM = PHOTON, the line is shortened to


IR  FILNAM  H123  REAC

(format: I3,1X,A6,1X,A4,Axxx)

for reading from the spectroscopic line-shape database PHOTON, but additional lines are needed to specify the line broadening options. See below.

There also exist an option for specifying constant (or very simple functional dependencies) cross sections, reaction rate coefficients, reaction rates, etc. directly via the input file, without resorting to an external database. This has proven useful for testing purposes, e.g. comparison with simple analytic or numerical solutions. In this case: FILNAM = CONST, further details: see below.

Later EIRENE will assign the relevant atomic molecular of photonic processes to each test particle (and also bulk particle) from either this set of NREACI processes, or from a small set of "default processes", which are hard wired in EIRENE.

In versions younger than 2002 each reaction in these databases can contain a string 'fit-flag = value', with IFTFLG=value. IFTFLG is used in EIRENE to identify the type of fitting expression to be evaluated with the fitting coefficients from the database.

By default: IFTFLG=0 for all data, i.e., single or double polynomial fits.

All atomic rate coefficients and cross sections can be scaled with a constant factor (FREAC, see below), for sensitivity studies.

Excluded from this scaling option are all photonic processes and those elastic collisions for which EIRENE uses an interaction potential (H.0-type data), because here the scaling would not have the expected effect on transport, but rather it would only modify the effective small angle scattering cut-off in the binary collisions.

Atomic and molecular processes are always of the following type:

a ·A + b ·B → m ·M + n ·N + ∆E

Here A,B,M and N are labels for the type of pre- and post collision particles. a,b,m and n are the stoichiometric coefficients, and ∆E is the amount of internal energy transferred into (or at the expense of) the kinetic energy of the collision products M and N.

The following conventions are always in use:

At least one of A, B, M or N must stand for a "test particle", for otherwise this process is not relevant for the transport problem solved by EIRENE.

particles A are always "bulk particles", i.e., from the list of "bulk ions" (input block 5) or electrons.
particle B may be one "test particle". Then this process A + B→ ... must be included in the list of reaction decks for this test particle B. Hence: b = 1 then.
particle B may be one "bulk particle". Then this process must be included in input block 7 as a volume-source for particles M and/or N,
particles M and/or N can be either secondary "test particles" or secondary "bulk particles".
non-linear processes, i.e., processes in which both A and B should, in principle, stand for test particles, can be included by specifying particle A as bulk particle in input block 5 and, simultaneously, as test particle in input block 4. The parameters for the "artificial bulk species" A are iterated in a sequence of EIRENE runs. In cases in which species M and/or N are also available as ("real" or "artificial") bulk species, there is the choice to specify secondaries from a reaction either as test-particles and to continue the history with those, or, alternatively, to specify them as bulk particles, and stop the trajectory. In this latter case one must add a volume source to launch these "secondaries" via a proper spatially distributed source in the next iteration. By using this option carefully, and combining it with the stratified source sampling technique, coupling to external codes/models can be made either more or less implicit.
Some care is needed here to avoid double-counting.

Next EIRENE expects one so called "species block" for each test particle species. Later, in input block 5, it will also expect one species block for each of the heavy background particles, i.e., for the "bulk ions".

A "species block" has the format:

ISPZ$ TEXT$ NMASS$ NCHAR$ NPRT$ NCHRG$ ISRF$ ISRT$ ID1$
NRC$  NFOL$ NGEN$  NHSTS$ ID3$
 (format: I2,1X,A8,1X,12(I2,1X))

In case ID1 ≠ 3, i.e., two or less post collision particle (not counting electrons)

 DO  K= 1, NRC$
   IREAC$  IBULK$  ISCD1$  ISCD2$  ISCDE$  IESTM$  IBGK$
   EELEC$  EBULK$  ESCD1$  EDUMMY  FREAC$  FLDLM$
 ENDDO

New option in 2006: In case ID1 = 3, i.e., three post collision particles (not counting electrons)

 DO  K= 1, NRC$
   IREAC$  IBULK$  ISCD1$  ISCD2$  ISCD3$ ISCDE$  IESTM$  IBGK$
   EELEC$  EBULK$  ESCD1$  EDUMMY FREAC$  FLDLM$
 ENDDO

New option in 2006:

In case ID1 = 4, i.e., four post collision particles (not counting electrons)

 DO  K= 1, NRC$
   IREAC$  IBULK$  ISCD1$  ISCD2$  ISCD3$ ISCD4$ ISCDE$  IESTM$  IBGK$
   EELEC$  EBULK$  ESCD1$  EDUMMY FREAC$  FLDLM$
 ENDDO

where $ stands for either A (atoms), M (molecules), I (test ions), PH (photons) or P (bulk ions). I.e., a species block consists of one species specification card ISPZ$ ..., and NRC$ "reaction decks" of two cards each.

For some particle species (in particular for hydrogenic atoms, molecules and molecular ions, and for helium atoms) EIRENE has default A&M data, to which it resorts, if no reaction cards are in the input deck (i.e., if NRC$=0, see below) for a particular species. These default models are described at the end of this section. They are overruled by the reaction cards.

In order to de-activate all possible reactions (also the default reactions), e.g., to simulate the collision-less "Knudsen flow" for all or a particular test particle species, one must set NRC$ =−1.

Usually, a particular particle type and species in the species blocks is identified by an integer IPART$=LMN (3 digits). Here N fixes the type of the particle, M is the number of particles characterized by IPART$ and L is the species index within the specified group (type) of particles.

The following values of N are currently in use:

N=0: photons (versions 2004 and younger)
N=1: atoms
N=2: molecules
N=3: test ions
N=4: bulk ions
N=5: electrons

The Input block

*** 4. Data for Species Specification and Atomic Physics Module

Note: in EIRENE₂₀₀₃ or later the format of input flag REAC described below was generalized from A9 to Axxx, with xxx ≤ 50, to accommodate the longer reaction information for photonic processes into the same format. The old format A9 is of course automatically still recognized as special case.

NREACI
DO

for "real particles", databases: FILNAM = HYDHEL, AMJUEL, METHAN, H2VIBR

   IR  FILNAM  H123  REAC  CRC  MASSP  MASST  DP RMN  RMX
   (format: I3,1X,A6,1X,A4,Axxx,A3,2I3,3E12.4)
     IF (RMN.GT.0.) READ IFEXMN, FPARM(J),J=1,3
     IF (RMX.GT.0.) READ IFEXMX, FPARM(J),J=4,6

for "real particles", database: FILNAM = ADAS (or any equivalent ADF11 format)

   IR  ADAS  H123  REAC  CRC  MASSP  MASST  DP
   (format: I3,1X,A6,1X,A4,Axxx,A3,2I3,3E12.4)
   ELNAME, IZ
   (format: 4X,A2,1X,I3)

i.e.: one additional input card is read to identify the chemical element and charge state. Extrapolation flags RMN and RMX are not in use in this case and may be omitted. Default (hard wired) extrapolation schemes are used if the plasma density and temperature is outside the tabulated range.

for "real particles", database: FILNAM = CONST. Fit parameters are directly read from input file, up to nine coefficients. In EIRENE₂₀₀₅ or later: Optional : Further parameter FTFLAG, length: not more than 9 characters.

   IR  CONST  H123  FTFLAG(optional) CRC  MASSP  MASST  DP
   (format: I3,1X,A6,1X,A4,(A9),xxxX,A3,2I3,3E12.4)

i.e.: Reading of FTFLAG is optional. Default for FTFLAG: =0. Parameter REAC does not exist. Extrapolation flags RMN and RMX are not in use in this case and may be omitted.

for "photons", database FILNAM = PHOTON

   IR  PHOTON  H123  REAC
   (format: I3,1X,A6,1X,A4,Axxx)
   IPRFTYPE IPLSC3 IMESS IFREMD NRJPRT
   DO IFREMD
     II KENN  IK6
    (format: I6,1X,A2,3X,I6)
   ENDDO

I.e.: after the first line there is a special format for this set of input cards in case of photonic processes, read from the spectral database PHOTON. In case IFREMD > 0 the additional cards specify options for "foreign pressure line broadening".

Then the subroutine SLREAC is called, which picks the atomic data (fit coefficients) for reaction "IR" from the file FILNAM. It then stores this information on the arrays in Module COMAMF with label IR.


 * 4A. atoms species cards
 NATMI
 DO 44  IATM= 1, NATMI

* read NATMI species blocks with $ = A

 44 CONTINUE
 * 4B. molecules species cards
 NMOLI
 DO 46  IMOL= 1, NMOLI

* read NMOLI species blocks with $ = M

 46 CONTINUE
 * 4C. test ions species cards
 NIONI
 DO 48  IION= 1, NIONI

* read NIONI species blocks with $ = I

 48 CONTINUE

and, for versions younger than 2003:

 * 4D. test ions species cards
 NPHOTI
 DO 49  IPHOT= 1, NPHOTI

* read NPHOTI species blocks with $ = Ph

 49 CONTINUE

Meaning of the input variables

NREACI

Total number of different reactions to be read.

The next block has different meanings for "real particles" and "photons". Cross section and rate coefficients are specified for particles, but emission and absorbtion line shapes are specified for photons.

"real particles"

An interaction potential V(r), cross section σ plus a reaction rate coefficient 〈σv〉 for one process counts as one reaction, but higher order rate coefficients such as energy or momentum weighted rate coefficients for the same process count as new reaction and must be labelled by a different index IR (see below).

Storage is provided for up to NREAC different additional reactions (see "Parameter Statements", section 3.1), i.e., one must guarantee NREACI.LE.NREAC

IR

Labeling index of the reaction that is read in. The condition

IR.LE.NREACI

must be fulfilled (otherwise: error exit).

FILNAM

Name of the data-set that contains the required reaction.

HYDHEL: Hydrogen and Helium data [23]
METHAN: Hydrocarbons data [24]
H2VIBR: Data for H₂ molecules and their isotopomeres, including effects of vibrational and electronic excitation, as well as rates needed for radiation trapping calculations.
AMJUEL: supplement to HYDHEL and METHAN for neutral gas transport Monte Carlo codes, e.g. multi-step reaction rates, etc.
CONST: reaction IR is a collision process with explicitly specified fit coefficients for cross sections or reaction rate coefficients, (depending upon input flag "H123"), e.g., constant, power law, etc. These fitting coefficients are directly read from the formatted input file rather than from an external file.
ADAS: Collisional radiative rate coefficients from ADAS ADF11 files, which must be located in a directory, the path to which is specified in input block 1, section 2.1, by one of the "CFILE" cards described there: CFILE ADAS path/adf11/
PHOTON

H123

Identification flag for the type of data: interaction potential parameters, cross section, rate coefficient, momentum loss rate coefficient, energy loss rate coefficient, etc.

case 0: FILNAM = CONST

currently available only for: H123 = H.1, = H.2, = H.5, or = H.8, i.e., only for single parameter fits.

: in case IFTFLG = 10, 110, 210,... , "single parameter", constant collision parameter. One more line
F(0), REAL
is read.
For all other values of IFTFLG, including the default IFTFLG=0, two more input lines with 9 fitting coefficients
F(I), I=0,5 (REAL)
F(I), I=6,8 (REAL)
are read, in subroutine SLREAC from the formatted input file.
The natural logarithm of cross-section: ln(σ) (H123 = H.1) or of a rate coefficient: ln(R) (H123 = H.2, = H.5, or = H.8) is computed as:
ln(σ) = ∑_I=0⁸ F(I) ·ln(E_lab)^I
and, likewise, a rate coefficient R is evaluated as:
ln(R) = ∑_I=0⁸ F(I) ·ln(T)^I
Here E_lab is the relative energy of impact, but with the mass of the charged particle (i.e., for a target of neutral particles at rest, see below), and T is the electron or ion temperature, depending on the type of impacting plasma particle.
This option permits specification of constant cross sections or rate coefficients (only the first of the nine parameters nonzero), or, e.g., specific power law energy dependencies or temperature dependencies (second of the nine parameters).

case 1: FILNAM = HYDHEL, METHAN, AMJUEL, H2VIBR

H.1: single parameter fit for cross section σ(cm²) as function of energy ELAB (eV), with ELAB = m_Lab / 2 ·v_rel² (For the definition of m_Lab = MASSP : see below.
H.2: single parameter fit for rate coefficient 〈σv 〉(cm³/s) as function of target temperature (eV), assuming "zero velocity" projectile
H.3: double parameter fit of rate coefficient 〈σv〉(cm³/s) as function of projectile energy (eV) and target temperature (eV)
H.4: double parameter fit of rate coefficient 〈σv 〉(cm³/s) as function of target density (cm⁻³) and target temperature (eV)
H.5: single parameter fit of target particle momentum weighted rate coefficient 〈 σ v · m_p · |v_p| 〉 (cm³/s ·AMU ·cm/s) as function of target temperature (eV)
H.6: double parameter fit of target particle velocity weighted rate coefficient 〈σv ·m_p ·| v_p | 〉(cm³/s ·AMU ·cm/s) as function of projectile energy (eV) and target temperature (eV)
H.7: double parameter fit of target particle velocity weighted rate coefficient 〈σv ·m_p ·| v_p | 〉(cm³/s ·AMU·cm/s) as function of target density (cm³) and target temperature (eV)
H.8: single parameter fit of target particle energy weighted rate coefficient 〈σv ·m_p /2 ·v_p² 〉(cm³/s ·eV) as function of target temperature (eV)
H.9: double parameter fit of target particle energy weighted rate coefficient 〈σv ·m_p /2 ·v_p²〉(cm³/s ·eV) as function of projectile energy (eV) and target temperature (eV)
H.10: double parameter fit of target particle energy weighted rate coefficient 〈σ v · m_p/2 · v_p²〉 (cm³/s ·eV) as function of target density (cm³) and target temperature (eV)
H.11: single parameter fit for any other data, e.g. to be used in special user supplied programs.
H.12: double parameter fit for any other data, e.g. to be used in special user supplied programs, (i.e. not understood by EIRENE)

case 2: FILNAM = ADAS

H.4: double parameter table of rate coefficient 〈σv 〉(cm³/s) as function of target density (cm⁻³) and target temperature (eV) read from ADAS adf11 files. Filenames starting with scd... (ionization rate coefficients contain ionization rate coefficients to charge state Z, Z=1,... Z_max) and filenames starting with acd...... contain recombination rate coefficients, from charge state Z, Z=1, Z_max. The particular value of Z is specified in the next input card.

case 3: FILNAM = PHOTON to be written

REAC

Name or number of the reaction as used in the file FILNAM (e.g.: REAC = 2.15.2 in file METHAN for the reaction e + CH₄⁺→ CH₃⁺ + H + e ). If FILNAM = CONST, then REAC can be used to input the flag IFTFLG. Whether or not REAC is to be interpreted as IFTFLG is identified by presence of the string FT (e.g.: 'REAC' = 'FT 110' would lead to IFTFLG = 110 for this reaction).

CRC

Identification for the type of the reaction process. Depending upon the value of this flag various assumptions are automatically made concerning the atomic data in the initialization phase.

CRC=EI electron impact collision, i.e., ionization, excitation or de-excitation or dissociation. (In older EIRENE version CRC=DS was also in use. This is now automatically identified with CRC=EI)

CRC=CX (quasi-) resonant charge-exchange

CRC=EL elastic collision

CRC=PI in-elastic ion impact collision (not fully implemented)

CRC=RC re-combination

CRC=OT other

MASSP

Mass number (AMU) of the first particle involved in the collision, for the original data as being read from the data-file. EIRENE automatically re-scales data according to the mass number of the actual first particle involved in this particular collision process (e.g. if data are given for H atoms in the data-file, but are used for D or T atoms in an EIRENE run).

MASST

Mass number (AMU) of the second particle involved in the collision, for the original data as being read from the data-file. EIRENE automatically re-scales data according to the mass number of the actual second particle involved in this particular collision process.

The next input deck is read only in case FILNAM = ADAS:

ELNAME: Name of element in ADF11 file, e.g. Fe for iron, C for carbon, W for tungsten, etc....
IZ: Charge state in ADF11 file. This value is always between 1 and Z_max. For ionization rate coefficients IZ specifies the ionization rate coefficient for the ionization process from IZ−1 to IZ. For recombination process IZ specifies the recombination rate coefficient for the recombination process from IZ to IZ−1.

The next input decks are read only in case FILNAM = PHOTON:

IPRFTYPE
IPLSC3
IMESS
IFREMD
NRJPRT: to be written

Then IFREMD cards are read

II, KENN, IK6: II: labelling index, irrelevant
KENN: to be written
IR: to be written

Note:

The nomenclature used in reference [23] has lead to a certain confusion with regard to the masses of the particles involved in a particular collision process. A series of test calculations has revealed the following definitions implicit in the tables of references [23] and [24]:

For the cross section energy scale, the relevant (laboratory frame) mass is that of the charged particle, assuming an ion- or electron beam incident onto a cold neutral gas at rest (i.e., the mass in the energy scale is neither the reduced mass nor is the energy given in units energy/AMU).

For rate coefficients, depending upon both the beam energy and the target temperature, however, the neutral particles are considered to be the beam particles, i.e., their mass is the relevant mass for the energy scale now, whereas the charged particles are considered as target, with their mass being the relevant mass for the temperature scale.

Following these conventions, EIRENE assumes that cross sections have been measured with the charged particle as projectile and the neutral target at rest. Thus MASSP is assumed to be the relevant mass number of the original data for the energy scale in cross sections (H.1) and for the temperature scale in the rate coefficients (H.2, H.3, H.4). If MASSP = 0 , MASSP is defaulted to the electron mass (AMU). I.e., MASSP is the mass number of those collision partners in the database, which would play the role of bulk particles in an EIRENE run, unless the automatic re-scaling to the true EIRENE bulk ions mass RMASSP(IPLS) of the bulk particle is carried out.

MASST is assumed to be the relevant mass number for the projectile energy scale in H.3, H.6 and H.9, i.e., for the test particles in an EIRENE run.

EIRENE then automatically converts the energy/temperature scales in the data fits to the correct scales for the particular particles masses (isotopes) involved in a collision event during the simulation process, if these do not happen to coincide with MASSP and MASST.

DP: (only for H.8, H.9, H.10) additional (e.g., potential) energy lost or gained, which is not already included in the rates. This may be needed due to the logarithmic fits used and if changes in sign arise. For example, in case of recombination of hydrogen ions, AMJUEL H.10, 2.1.8, DP=+13.6 must be specified.
In reaction decks H.8, H.9 or H.10 in AMJUEL, this value of DP is now automatically read (and overwrites the value specified here). Hence: in the EIRENE versions later than June96 this input flag is redundant.

The next flags control the options to extrapolate atomic data fits read from files. The default data in EIRENE and some reactions in the files AMJUEL already contain that information. EIRENE searches for the string ELABMIN and ELABMAX for cross section data. If it finds them, it uses them as RMN and RMX, respectively. In such cases the next card IFEXMN,.... and IFEXMX,.... may also be omitted, because the corresponding information is found from the atomic data file as well. The extrapolation expression is that of option IFEXM. = 5 (see below).

Likewise, in case of reaction rates, EIRENE searches for TMIN, TMAX (for single parameter fits H.2, H.5, H.8 and H.11) and additionally for PMIN and PMAX (if H.3, H.4, H.6, H.7, H.9, H.10 or H.12) and the extrapolation parameters again are expected in the atomic data file.

RMN

minimum energy or temperature for data as read from file. Below RMN the data are extrapolated using the input card

IFEXMN,FPARM(J),J=1,3.

See below, and e.g., subroutine CROSS for the various options.

Some cross section data in AMJUEL already contain the extrapolation information, and hence the default RNM=0 can be used there.

RMX

same as above, but for high parameter range extrapolation.

IFEXMN, FPARM(J), J=1,3

various extrapolation expressions ex_low(x) at the low end (x < RMN, x = E_lab or x = T) are available:

IFEXMN = 1: ex_low(x) = 0
(e.g., cross section with nonzero threshold at RMN)
IFEXMN = 2: y=x/FPARM(1)
ex_low(x) = FPARM(2) ·y^FPARM(3) ·log(y)
recommended for cross section extrapolation at high energy end, for reactions with nonzero threshold at FPARM(1) [eV].
IFEXMN = 3: y=log(x)
ln(ex_low(x)) = FPARM(1) + y ·FPARM(2)
IFEXMN = 4: not in use
IFEXMN = 5: y=log(x)
ln(ex_low(x)) = ∑_I=1³ FPARM(I) ·y^I−1
IFEXMN < 0: linear extrapolation in log-log scale. Coefficients FPARM(J), J=1,2 are automatically determined, e.g., in function CROSS, FPARM(3) = 0.D0, and then extrapolation option IFEXMN = 5 is used with these coefficients (same as IFEXMN = 3 in this case).

IFEXMX, FPARM(J), J=4,6

same as above, but for high parameter range extrapolation ex_high(x).

*4A. Atoms Species Cards

NATMI: Total number of atomic species blocks
IATM: irrelevant; labelling index

*4B. Molecules Species Cards

NMOLI: Total number of molecule species blocks
IMOL: irrelevant; labelling index

*4C. Test Ions Species Cards

NIONI: Total number of test ion species blocks
IION: irrelevant; labelling index

*4D. Photon Species Cards

NPHOTI: Total number of photon species blocks
IPHOT: irrelevant; labelling index

The meaning of the variables in a species block is as follows:

TEXT$

Name of the species on printout and plots

NMASS$

Mass number

NCHAR$

Nuclear charge number

NPRT$

Flux units carried by one particle of this type and species. The "WEIGHT" of a test flight has dimensions: "atomic flux", i.e., "equivalent atoms per second". For example, a methane molecule CH₄ should have NPRT_CH₄ = 5, if H atoms and C atoms have NPRT_H = NPRT_C = 1.

NPRT is used to convert fluxes into equivalent atomic fluxes, for scaling, (see input flag FLUX in block 7), and for flux conservation in non-diagonal (species changing) events at surfaces and in the volume. (NPRT is irrelevant for atoms, and always set to one for them.)

NCHRG$

Charge state (irrelevant for atoms and molecules, always set to zero for them)

[ISRF$ ] Species index flag for "fast particle reflection model" (see block 6A), irrelevant for molecules

[ISRT$ ] Species index flag for "thermal particle re-emission model" (see block 6A)

ID1$

not in use in versions 98 to 2002. Was "sputtered species index" in versions older than "98", but this species index has then become a surface property, see input block 3. In versions younger than 2002 this flag can be used to increase the number of secondary test particle groups from 2 or less (default) to 3, if ID1=3. This option then allows more complex fragmentation processes of large molecules as compared to the default.

NRC$

> 0: total number of reaction decks to be read for this species
= 0: default model is to be used.
= -1: no collision processes for this particle species at all

NFOL$

controls the motion of test particles.

≥ 0: default test particle tracing model is to be used.
= -1: motion of test particle is not followed (static approximation). The test particle is destroyed immediately at it's point of birth by a collision. A collision estimator is used rather than a track-length estimator for all default volumetric tallies, as described in section 3.2.

NGEN$

Maximum number of generations for this test particle species.

> 0: For each particle history the generation counter XGENER is reset to zero at birth, after surface events and after a non-diagonal event (modification of type and/or species of a particle in a collision event). It is increased at entropy producing events (elastic and diagonal charge exchange collisions). If XGENER > NGEN.. after a collision, no test particle secondaries are generated. The particle flux, parallel momentum flux and energy flux is put into the tallies PGENA, VGENA and EGENA for atoms, resp., PGENM, VGENM and EGENM for molecules, resp., PGENI, VGENI and EGENI for test ions, respectively and PGENPH, VGENPH and EGENPH for photons, respectively.
= 0: In case NGEN, (default) no generation limit is activated for atoms (NGENA=0), molecules (NGENM=0) or test ions (NGENI=0).

NHSTS$

(only in versions 2004 and younger): NHSTS=-1 turns off the trajectory plot for this particular species. (default: NHSTS=0)

ID3$

not in use

IREAC$

Identification flag for collision data. IREAC$ = IR, and IR is the labelling index from the reaction card (see above).

The next three flags control number, species and type of particles involved in this particular collision process.

IBULK$

pre collision bulk particle identifier

= NLM

M: type flag ITYP for impacting bulk particle (electron (ITYP=5) or bulk ion (ITYP=4),
irrelevant, defaulted to M = 4 for bulk ion impact collision and defaulted to M = 5 for electron impact collisions
L: irrelevant, defaulted to L = 1, because number of pre collision bulk particles is known from the type CRC of each collision process.
N: Species index of pre collision bulk ion (corresponds to IPLS in bulk ion species cards, see block 5). Irrelevant for electron impact collisions, CRC=EI.

ISCD1$

first heavy secondary particle group identifier

= IJKLM

M

type flag ITYP for these first secondaries group (atoms (ITYP=1), molecules (ITYP=2), test ions (ITYP=3) or bulk ions (ITYP=4))

L

number of secondaries of this type and species

JK

Species index of secondary (IATM, IMOL, IION, or IPLS), two digits

I

relative importance of this first secondary group as compared to the second secondary group. If I ≤ 0, I is defaulted to I = 1

(currently not in use)

ISCD2$

second heavy secondary particle group identifier

= IJKLM.

Same as ISCD1$, but for second secondary group.

Note:

The number of secondary electrons need not be specified, it is computed from charge conservation.

In case of collisions with only one heavy particle secondary (electron impact collision CRC=EI, or re-combination CRC=RC) ISCD1 is used and ISCD2 is irrelevant.

In case of charge exchange (CRC=CX) ISCD1 is used for the new state of the impacting bulk particle after the collision, whereas ISCD2 controls the options for the post collision state of that particle which was the test particle prior to the collision. Consistency of the particle masses of pre and post collision particles with this convention is checked in the initialization phase of each run.

In case of elastic collision (CRC=EL), the ISCD1 and ISCD2 flags are are irrelevant, because the colliding particles retain their identity by default.

In case of electron impact collision (CRC=EI), the ISCD1 and ISCD2 flags are symmetric, i.e. they can be interchanged without any effect on the calculation.

ISCDE$

flags for post collision velocity space distributions for all particles involved in the collision event. Each of the five digits controls a different type of particle.

Let us write JKLMN = ISCDE

Let then furthermore the character # stand for J, K, L, M, or N respectively. Then J controls electrons, K controls pre collision (plasma) bulk ions (if any involved), L the secondary particles (if any). Each of these single integer flags controls the meaning of the energy parameters

EP=EELEC (# [^=] J),

EP=EBULK (# [^=] K),

EP=ESCD (# [^=] L),

for the electrons, for the pre collision bulk ions, and for the heavy secondary particles, respectively.

# = 0: constant loss or gain of EP (eV) per collision event.
# = 1: bulk temperature dependent loss -(1.5 ·T + E_Drift) per collision, where
T = T_e for # [^=] J and T = T_i for # [^=] K;
T_e and T_i are the local electron and ion temperatures (eV) resp. at the point of collision, and E_Drift is the kinetic energy of the drift motion (VXIN,VYIN,VZIN) of the background ions (only in case NLDRFT=TRUE), see input block 5.
If # stands for a heavy particle secondary, i.e. # [^=] M and/or # [^=] N, then a fraction EP of the energy loss of the impacting bulk particle is distributed equally over all (if more than one) heavy particle secondaries of group 1 and 2 resp.
# = 2: Let # stand for a pre collision bulk particle: # [^=] K or # [^=] L. The energy loss of the impacting bulk particle is calculated from the velocity and energy weighted rate coefficients. In this case EP must be the labelling index IR of the reaction card by which these moments have been read from the data-file.
Let now # stand for a post collision heavy particle group, i.e # [^=] M and/or # [^=] N. Then the velocity vector of the secondary particle is sampled from a cross section weighted (and optional, if NLDRFT, shifted by the local drift velocity of the pre collision bulk particle) Maxwellian distribution at the local ion temperature T_i. The cross section is taken at the relative velocity of the impacting test particle and a particle sampled from the above mentioned (shifted) Maxwellian. This option permits to identify one particular collision partner from the background population, correctly accounting for the energy dependence of the cross section.
# = 3: to be written
# = 4: to be written
# = 5: to be written

IESTM

Three digits IESTM=LMN. N, M and L are flags for the choice of estimators for particle, momentum and energy sink and source term contributions from the collision process. (Default: =0 track-length estimators.)

IBGK

Three digits IBGK=NML. Particle identifier for BGK self and cross collisions. Format as for IBULK, ISCD1 and ISCD2 species identifier flags. I.e., it has three digits, pointing to the type and the species index within that type-group. The second digit is irrelevant and defaulted to one (only binary collisions are considered). This flag controls options for non-linear test-particle - test-particle collision.

Let ISPZ be the species index of a test-particle. If the collision is a self-collisions, this flag must point to the species ISPZ itself. Otherwise it must point to the 2nd test particle species ISPZ2, which is involved in the BGK-cross collision term.

IBGK must be consistent with IBULK, which in this case is the "artificial background species" (= result from previous iteration) in the linearization of the collision operator.

I.e.: if test particle species A collides with test particle species B (of same or of different species and type), then a bulk species AB_bulk must be present in input block 5, such that species A and AB_bulk have the same mass, nuclear charge and charge numbers. The density, flow field and temperature of collision partner AB_bulk are iterated, using either the corresponding parameters estimated for species A (in case of self collisions) or using parameters obtained from the formulas for T_AB, V_AB (see section 1.9.1) for cross collisions between species A and B.

The flag IBULK must point to the corresponding species IPLS=AB_BULK in block 5.

An example for input specifications of non-linear BGK collisions is given below in section 2.4.2. (Default: = 0, no iteration for non-linearity in collision kernel, no "artificial background species" involved in this collision, and no BGK-tallies are updated for iterating the artificial background species.)

EELEC

parameter EP for electron energy loss per collision

EBULK

parameter EP for pre collision bulk ion energy loss per collision

ESCD1

parameter EP for velocity space distribution of secondaries

EDUMMY

(was earlier: ESCD2) in versions 1999 and older: parameter EP for velocity space distribution of second secondaries group. Now: this parameter is not in use anymore. Distribution of energy EP over secondaries is now by default inversely proportional to the secondary particle masses (from momentum conservation in collision). For backward compatibility of input files now: EP = ESCD1 + EDUMMY.

FREAC

multiplier to scale the cross section (and rate coefficients) of this process. Default: FREAC = 1.0. If FREAC ≤ 0, it is also defaulted to 1.0. To turn off this process, set FREAC to a very small number: e.g.: FREAC = 1.0000 E-30

FLDLM

Flag for "generation limit" or "fluid-limit", Default: FLDLM = 0.0. This flag allows to cut trajectories after a certain maximum number of CX or EL collisions (generations) (FLDLM > 0.0) or if the mean free path compared to a geometric dimension falls below a certain value (FLDLM < 0.0) (fluid limit). The precise status of this option should be inquired from the code authors.

2.4.1 Default atomic and molecular data

If default atomic data are requested for some test particle species, these are activated by setting the flag NRCM(IATM)=0 (or NRCM(IMOL)=0, NRCI(IION)=0, resp.) for a particular atom IATM (molecule IMOL or test ion IION), and by omitting the corresponding cards
IREAC$,... and
EELEC$,...
Then EIRENE tries to find, automatically the corresponding species types and indices of the collision products in the species blocks. It then tries to define all variables in the cards
IREAC$,... and
EELEC$,...
for collision rates and collision kinetics, using default assumptions.

Assume that the bulk ion H⁺ is specified as bulk ion species no. "a", D⁺ as bulk species "b", T⁺ as bulk species "c" and He⁺ as bulk ion species no. "d", then the default atomic collision models for H, H₂, H₂⁺ (and isotopes) and He are identical to a model that would result from the following specifications in input block 4.

For atoms D (likewise for H or T), e.g., IATM=1, and atoms He (e.g., IATM=2):

* ATOMIC REACTION CARDS  NREACI=
5
  1 HYDHEL H.2 2.1.5     EI  0  1
  2 HYDHEL H.1 3.1.8     CX  1  1
  3 HYDHEL H.2 2.3.9     EI  0  4
  4 HYDHEL H.1 5.3.1     CX  4  4
  5 HYDHEL H.1 6.3.1     CX  4  4
*NEUTRAL ATOMS SPECIES CARDS:  NATMI SPECIES ARE CONSIDERED,
NATMI= 2
 1 D         2  1  1  0  1 -1  0  2
     1   115   *14     0 00000
 -1.3600E 01  0.0000E 00
     2   *14   111   *14 00000
  0.0000E 00  0.0000E 00  0.0000E 00
 2 He        4  2  1  0  2  2  0  3
     3   115  **14     0 00000
 -2.4588E 01  0.0000E 00
     4  **14   211  **14 00000
  0.0000E 00  0.0000E 00  0.0000E 00
     5 ***14   211 ***14 00000
  0.0000E 00  0.0000E 00  0.0000E 00

Here ∗ is the species index of the D⁺ bulk ion in input block 5, ∗∗ is the species index of the He⁺ bulk ion in input block 5, and ∗∗∗ is the species index of the He⁺⁺ bulk ion in input block 5. If a corresponding background species is not found in block 5, then the corresponding default process is de-activated.

Hence: D and He atoms are ionized by electron impact, using the Corona data 2.1.5 and 2.3.9, respectively, from reference [23], and with constant electron energy loss per ionization of 13.6 and 24.588 eV, respectively. D atoms also undergo resonant charge exchange with D⁺ ions, utilizing the cross section data for reaction 3.1.8 (loc.cit.).

The ISCDE parameters for CX processes are zero: = 00000 (or: = 00001, which has the same meaning because the fifth digit is presently irrelevant). The second digit (i.e., here = 0) controls the options for the velocity of the background ion going into the collision. For this choice the following model is activated:

With respect to ions going into a CX collision the ion velocity distribution f_i(v_i) is assumed to be an isotropic, mono-energetic distribution , with v_i⁰ = √{1.5 kT_i}, then shifted by v_Drift,i.

The mean energy of ions undergoing charge exchange is taken to be 1.5·T_i eV plus the kinetic energy in their drift motion m_i/2 (v_Drift,i)². This is used for the corresponding part in the CX momentum- and energy exchange rate coefficient. The reaction rate coefficient itself is evaluated as

〈σ·v_rel〉 = σ(v_eff) ·v_eff

with an effective relative velocity

v_eff =

√

3/2 kT_i/m_i + (v₀ − v_Drift,i)²

This approximation to the rate coefficient, for this particular choice of f_i, is almost exact. It is assuming that the effect of the scalar product −v_i·v₀ in the relative velocity between ion- and neutral velocity becomes negligible after integrating σv_rel over solid angle, which it would indeed exactly if σ ∝ 1/v_rel.

The velocity of the impacting ion is sampled, consequently, from a drifting, isotropic and mono-energetic distribution weighted by σ·|v_rel|, with v_rel the relative velocity between the two collision partners.

Alternatively to the weighting mentioned above, a rejection method can be used to avoid the weighting and still simulate the same collision integral. Check in subr. VELOCX to find out which of the two physically equivalent methods are in use.

Default reactions 5.3.1 and 6.3.1 (added in April 2006, to simplify input for pure He plasma simulations) are the two resonant charge exchange processes

He + He⁺ → He⁺ + He

He + He⁺⁺ → He⁺⁺ + He

The weighting (or rejection-technique) was absent in versions younger than 99 for this particular type of (default) charge exchange collision integral approximation, which, strictly speaking, had therefore led to a slight violation of the second law of thermodynamics (H-Theorem), due to an inconsistency between the in-scattering term (determined by σ_CX and the out-scattering term (determined by 〈σv_rel〉) in the Boltzmann CX collision integral, although mass and energy had been strictly conserved in each collision.

Alternative to the choice ISCDE = 00000 made in the default model one may also use ISCDE = 01000 (with all the rest kept identical). In this case the rate coefficients are evaluated in the same way as described above, but neutrals emerging from CX collisions are sampled from the local drifting Maxwellian distribution of ion collision partners rather than from a drifting mono-energetic distribution.

For hydrogenic molecules (H₂, D₂, HT,...) or molecular ions (H₂⁺, D₂⁺, HT⁺,....) default dissociation, ionization and dissociative recombination data are available.

EIRENE uses the 6 reaction rates 2.2.5, 2.2.9, 2.2.10, 2.2.11, 2.2.12 and 2.2.14 from the data-set HYDHEL (loc.cit.) and again tries to identify the reaction products from the mass and nuclear charge number of the respective molecule IMOL (or test ion IION). In order to distinguish D₂ from HT, is also uses the name TEXTM (or TEXTI), by looking for the appearance of the character "D" in the name (then: D₂, or D₂⁺) or for "H" or "T" (then: HT, or HT⁺). The default reaction kinetics are set as if they would have been specified by the following species blocks:

* ATOMIC REACTION CARDS  NREACI=
6
  1 HYDHEL H.2 2.2.9     EI  0  2
  2 HYDHEL H.2 2.2.5     EI  0  2
  3 HYDHEL H.2 2.2.10    EI  0  2
  4 HYDHEL H.2 2.2.12    EI  0  2
  5 HYDHEL H.2 2.2.11    EI  0  2
  6 HYDHEL H.2 2.2.14    EI  0  2
  7 AMJUEL H.8 2.2.14    EI  0  2
       .
       .
       .
* NEUTRAL MOLECULES SPECIES CARDS:  NMOLI SPECIES, NMOLI=
1
 1 D2        4  2  2  0  0  1  0  3
     1   115   113     0 00000
 -1.5400E 01  0.0000E 00  0.0000E 00  0.0000E 00
     2   115   121   000 00000
 -1.0500E 01  0.0000E 00  3.0000E 00  3.0000E 00
     3   115   111   *14 00000
 -2.5000E 01  0.0000E 00  5.0000E 00  5.0000E 00
* TEST ION SPECIES CARDS:  NIONI ION SPECIES, NIONI=
1
 1 D2+       4  2  2  1  0  1  0  3 -1
     4   115   111   *14 00000
 -1.0500E 01  0.0000E 00  4.3000E 00  4.3000E 00
     5   115   *24       00000
 -1.5500E 01  0.0000E 00  0.2500E 00  0.2500E 00
     6   115   121   000 30300
  7.0000E 00  0.0000E 00  7.0000E 00  0.0000E 00

As above, ∗ stands for the species index of the D⁺ bulk ion. Note that for reaction 6, the energy dependence of the cross section leads to a mean electron energy loss per event [ˉ(E)]_el of about 0.88 ·T_e, see reference [25]. This approximation is used in the default model for the electron energy loss E_el per collision and for the total kinetic energy release E_K shared by the two product atoms. The reaction no. 7 above: AMJUEL H.8 2.1.14, provides values very close to this, by an independent fit to the correct electron energy loss for this process

〈σv_e E_el〉 =

E_el

·〈σv_e〉 = kT_e 〈σv_e〉·

⎛
⎝

3/2 +

d ln〈σv_e〉

d ln(kT_e)

⎞
⎠

(4.2)

The velocity distribution of the dissociation products is isotropic in the center of mass system (taken to be the system moving with the incident molecule), and the dissociation energy is shared by the dissociation products so as to preserve momentum. I.e., in case of unsymmetrical hydrogenic molecules the two dissociation products to not receive the same share of the dissociation kinetic energy release, but instead this is distributed inversely proportional to their masses.

In case of mixed molecules, such as DT, HT or HD in case of default reaction models some rate coefficients are automatically splitt into two halfs and assigned to the proper product atomic particles.

For hydrogenic bulk ions (H⁺, D⁺, T⁺) default volume recombination rates are available [26].

Further examples

to be written

2.4.2 Neutral-Neutral collisions in BGK approximation

In the example considered here we have the 4 neutral-neutral collision processes: D on D2, D on D and D2 on D2, D2 on D. In input block 4 there are three BGK relaxation rate coefficients (note: the rate coefficient for D2 on D is equal to the one for D on D2, hence 3 rate coefficients are sufficient).

     .
     .   some 20 other reaction cards in this example
     .

 21 CONST  H.2           EL  1  1
 -2.1091E+01  0.2500E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00
  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00
 22 CONST  H.2           EL  2  2
 -2.0589E+01  0.2500E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00
  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00
 23 CONST  H.2           EL  1  2
 -2.0357E+01  0.2500E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00
  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00  0.0000E 00

In input block 4a (for atoms), we specify 4 collision processes for D atoms. Process no. 3 and 4 are neutral-neutral collisions.

1 D         2  1  1  0  1 -1  1  4  (i.e.: 4 processes to be specified for this D atom)
     .
     .  2 other process decks, e.g. CX and electron impact ionization
     .
    21  2014   (=IBULK)   1001     0   111  (=IBGK , self collision)
 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    23  2214   (=IBULK)   1001     0   112  (=IBGK, cross collision with IMOL=1)
 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00

In input block 4b (for molecules), we specify 7 collision processes for D2. Process no. 6 and 7 are neutral-neutral collisions.

1 D2       4  2  2  0  0  1  0  7   (i.e.: 7 processes to be specified for this D2 molecule)
     .
     .
     .   5 other process decks.
     .
     .
    22  2114  (=IBULK)    1001     0   112   (=IBGK, self collision)
 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
    23  2314  (=IBULK)    1001     0   111   (=IBGK, cross collision with IATM=1)
 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00

In input block 5 there are at least 23 bulk ion species in this example. Of those species no. 20 to 23 are "artificial" BGK species for neutral-neutral self collisions specified above by IBULK = 2014, 2114, 2214 and 2314 in blocks 4a and 4b above.

    .
    .
    .  19 other bulk ion species decks.
    .
    .
20 D(B)      2  1  1  0  1  1  0  0
21 D2(B)     4  2  2  0  1  1  0  0
22 DD2(B)    2  1  1  0  1  1  0  0
23 D2D(B)    4  2  2  0  1  1  0  0

2.4.3 Fitting expressions (IFTFLG)

The default fitting expressions (IFTFLG=0) are single or double polynomial expressions, as described in [23] or [24] for the databases HYDHEL and METHANE, respectively.

The following options are available in recent versions of the code:

a) for cross sections, rate coefficients, rates:

IFTFLG < 100 : rate coeff. or cross section, [cm³/s] or [cm²].

IFTFLG > 100 : rate = density * rate coeff, [s] (also used for spontaneous decay, e.g. radiative decay.

mod(IFTFLG,100) = 10: only one fitting coefficient is read, i.e. the cross section, rate coefficient or rate is constant.

b) for interaction potentials (elastic processes):

to be written

2.5 Input for Plasma Background

General remarks

The background medium (mostly plasma) consists of NPLS (sometimes in the code: NPLSI) different so called "bulk particle" species, also referred to as "bulk ions", by abuse of language,. It is described by several blocks of NSBOX data each, (i.e. one datum per grid cell), namely

one array for the electron temperature T_e (eV) (TEIN),
one common or NPLS (one for each bulk ion species) arrays for the ion temperature(s) T_i (eV) (TIIN),
NPLSI arrays (one for each bulk ion species), ion densities n_i ( cm⁻³ ) (DIIN),
one common or NPLS (one for each bulk ion species) arrays for the (cartesian) drift velocities V = ( V_x , V_y , V_z ) (VXIN, VYIN, VZIN), either in cm/s or Mach number units.

From the assumption of quasi-neutrality, the specified charge numbers of each ion species, and the data for n_i, an array of NSBOX data for the electron density n_e (cm⁻³ ) (DEIN) is computed (see next subsection: "derived background data").

The unit vector b = ( b₁ , b₂ , b₃ ) (BXIN, BYIN, BZIN), parallel to the magnetic field is set to (0.,0.,1.) by default, i.e. pointing in z- (or toroidal-) direction. The default magnetic field strength is |B|=1 T in all grid cells. This default setting may be overwritten by using the INDPRO(5) flag and the corresponding input parameters (see below). Here either a (radial or x-) B-field "pitch" profile can be selected (INDPRO(5)=1,...4, and geometry level LEVGEO=1,2,3), with pitch = B₂/B_tot = b₂. From this a cell centered B-field unit vector is internally constructed, by assuming that the first set of coordinate surfaces (x- or radial surfaces) are flux-surfaces (B-field aligned and constant pitch). The magnetic field direction then has components (0., pitch,√{1−pitch²}) in the three coordinate directions.

Currently only INDPRO(5)=3 profile option allows to define also the B-field strength (BFIN) [T], solely by the input flags described below.

Alternatively and in all other options an external data source for the B-field vector in cartesian coordinates (INDPRO(5)=5,...,7, all geometry levels LEVGEO) is used directly. In this latter case the transfer of magnetic field data (BXIN,BYIN,BZIN and BFIN) for each cell into EIRENE proceeds from user supplied profile subroutine PROUSR (INDPRO(5)=5) or, from external file (INDPRO(5)=6,7), via work array RWK, and then from subroutine PROFR).

The same applies to the "additional plasma profile array" ADIN, INDPRO(6), to permit usage of EIRENE output routines for plasma profiles such as pressure, energy fluxes, etc., which are not directly needed in an EIRENE run but which may be available from an external plasma code, here: INDPRO(6)=5,...,7 only.

It also applies to the space and species dependent weight function WGHT (to be used for non-analog sampling techniques, INDPRO(7)), which is defaulted to 1.0D0 for all cells and all NSPZ species, and it applies to the zone volume array VOL cm⁻³, INDPRO(12). Defaults are the cell volumes computed in EIRENE subroutine VOLUME from the standard mesh data.

Plasma data in the additional zones IADD, IADD=NSURF+1,NSURF+NADD outside the standard mesh are defaulted to the "EIRENE vacuum data" given below. Plasma data other than these in these zones have to be specified either explicitly in input block 8 (see below), or in the user supplied subroutine PLAUSR or with the INDPRO=7 option (see below, this section).

Cell volumes in the additional cell region are defaulted to 1.0 cm⁻³, unless INDPRO(12)=4, 5, 6, 7.

A standard mesh zone ICELL is automatically identified as a vacuum zone with regard to plasma species IPL (no collisions with the bulk particles IPL there, LGVAC(ICELL,IPL)=.TRUE., for electrons: IPL=NPLS+1)

if at least one of the following conditions is met:

for electrons: IPL=NPLS+1:

TEIN(ICELL) ≤ TVAC ,

DEIN(ICELL) ≤ DVAC

for bulk ion species IPL=IPLS:

TIIN(ICELL,IPL) ≤ TVAC,

DIIN(ICELL,IPL) ≤ DVAC

and for all background species, IPL=0:

LGVAC(ICELL,I)=TRUE for I=1,NPLSI and for I=NPLS+1.

The EIRENE vacuum data are defaulted (in subroutine PLASMA) to:

DVAC = 1.0000E 02 cm⁻³

TVAC = 2.0000E-02 eV

Furthermore, zero plasma drift velocities are set in cells which are treated as vacuum zones.

VVAC = 0.0000E 00 cm/s

Note: A cell can be considered a vacuum cell with respect to electrons, (LGVAC(...,NPLS+1) = TRUE), without being a vacuum cell for all background ions. This is because, by abuse of language, also neutral particle species may be used as background "ion" species (NCHARP = 0), to include neutral-neutral collisions, e.g., by the iterative option NITER (input block 1).

The Input Block

*** 5. Data for Plasma background

 NPLSI
 DO 51  J= 1, NPLSI

* read NPLSI species blocks with $ = P

 51 CONTINUE
 (INDPRO(J), J=1,12)
 IF (INDPRO(1).LE.5) THEN
 TE0  TE1  TE2  TE3  TE4  TE5
 ENDIF
 IF (INDPRO(2).LE.5) THEN
 (TI0(I)  TI1(I)  TI2(I)  TI3(I)  TI4(I)  TI5(I)  I=1,NPLSI)
 ENDIF
 IF (INDPRO(3).LE.5) THEN
 (DI0(I)  DI1(I)  DI2(I)  DI3(I)  DI4(I)  DI5(I)  I=1,NPLSI)
 ENDIF
 IF (INDPRO(4).LE.5) THEN
 (VX0(I)  VX1(I)  VX2(I)  VX3(I)  VX4(I)  VX5(I)  I=1,NPLSI)
 (VY0(I)  VY1(I)  VY2(I)  VY3(I)  VY4(I)  VY5(I)  I=1,NPLSI)
 (VZ0(I)  VZ1(I)  VZ2(I)  VZ3(I)  VZ4(I)  VZ5(I)  I=1,NPLSI)
 ENDIF
 IF (INDPRO(5).LE.5) THEN
 B0  B1  B2  B3  B4  B5
 ENDIF
 IF (INDPRO(12).LE.5) THEN
 VL0  VL1  VL2  VL3  VL4  VL5
 currently option INDPRO=6 and INDPRO=7 are not available for volumes.
 ENDIF

Meaning of the Input Variables for Plasma Parameters

The meaning of the variables in the bulk ion species cards is as in block 2.4 for the species cards specified there.

New options since 2003:

However, after the ID3 flag, there may be 2 additional input flags: (only for Eirene₂₀₀₃ and younger):

..... CDENMODEL NRE

Format: ....,1X,A10,1X,I2.

These flags, if included, control additional options to set the density, temperature and flow field of the selected species IPLS. If the string CDENMODEL is found here for a species IPLS, then, after reading the species and reaction cards for IPLS, one (default) or NRE further input cards are expected.

CDENMODEL = 'Multiply'

docu to be written

CDENMODEL = 'Constant'

docu to be written

CDENMODEL = 'fort.13'

Read data for this species from file fort.13.

The additional card (only one) reads (format 2I6):

ISP, ITP

ISP is the number of a bulk species on the file fort.13 (e.g. written in an earlier EIRENE run, see NFILEL option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are set from those of species ISP on fort.13

ITP is the type of the particle, i.e., always: ITP=4 for this option. Internally set: ITP=4, bulk particle type, independent of input value for ITP.

CDENMODEL = 'fort.10'

Read data for this species from file fort.10.

The additional card (only one) reads (format 3I6):

ISP, ITP, ISTR

ISP, ITP is the number and type of a test particle species on the file fort.10, respectively, from stratum no. ISTR (i.e. written onto fort.10 either in the present run or an earlier EIRENE run, see NFILEN option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are derived from from those of species ISP, ITP, ISTR on fort.10.

Note: at the end of a full EIRENE run the corresponding bulk particle profiles are reset to those evaluated in this run, for further use in post processing routines such as the diagnostic module (see Section 2.11), printout, plotting or iterative mode (input flags NITERI, NITERE, and user subroutine MODUSR, see Section 2.1)

CDENMODEL = 'Saha'

Set parameters from Saha-equilibrium of ionization states. The additional card (only one) reads:

ISP, ITP, ISTR, ....

to be written

CDENMODEL = 'Boltzmann'

Set parameters from Boltzmann-equilibrium of excited states. The additional card (only one) reads (format: 3I6,6x,2E12.4)

ISP, ITP, ISTR, G_BOLTZMANN, DELTA_E

docu to be written

CDENMODEL = 'Planck'

(only for background photons) to be written

CDENMODEL = 'Corona'

Set parameters from Corona-equilibrium. The additional card (only one) reads (format: 3I6,1X,A6,1X,A4,A9,A3,E12.4):

ISP, ITP, ISTR, FILNAM , H123 , REACTION, CR, A_CORONA

ISP, ITP is the number and type of a particle species, either one of the already defined bulk species IPLS: ISP=IPLS, (for ITP=4), or from the file fort.10 (for ITP=0,1,2 or 3) from stratum no. ISTR (i.e. written onto fort.10 either in the present run or an earlier EIRENE run, see NFILEN option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are set from those of species ISP, ITP, ISTR.

This model then constructs a background species IPLS which is in corona-equilibrium with these ISP,ITP, ISTR particles.

CDENMODEL = 'Colrad'

set parameters from collisional radiative equilibrium, using reduced population coefficients from atomic database. The NRE additional cards read:

ISP, ITP, ISTR, FILNAM  H123 ....

This model then constructs a background species IPLS which is the sum of NRE species in collision-radiative-equilibrium with these ISP,ITP, ISTR particles.

INDPRO

Flag-array for the type of profile. The last digit (between 1 and 9) controls the type of profile. A second digit and/or the sign control further options, as described below. INDPRO is an array of length 12. Each element in this array controls one particular input tally, namely:

INDPRO(1): for TEIN
INDPRO(2): for TIIN By the default (0 < INDPRO(2) < 10): one common T_i profile for all "bulk ion" species is set. I.e., read only one profile card.
New options since 2001:
Values of INDPRO(2) larger than 10: use only the last digit, and one T_i profile card must be read for each bulk ion species IPLS. For example: INDPRO(2)=15 or =25, means: one separate ion temperature must be specified for each bulk ion species, and the profile type is 5. (INDPRO(2)=-5 would do the same.)
INDPRO(3): read NPLSI cards, for DIIN , IPLS = 1,NPLSI
INDPRO(4): read NPLSI cards, for VXIN, VYIN, VZIN , IPLS = 1,NPLSI
New options since 2001:
By the default (0 < INDPRO(4) < 10): one separate flow field for each bulk species is set, velocity is given in cm/sec.
Negative value of INDPRO(4) means: Mach number units instead. The sound speed cs is taken to be the isothermal ion acoustic speed of species IPLS.
ABS(INDPRO(4) is then used as flag for the choice of profile type.
Values of INDPRO(4) larger than 10: use only the last digit, and only one common flow field is set for all bulk ion species. Note the difference to the T_i options: there the meaning of INDPRO larger than 10 was exactly opposite to the meaning here for the flow fields (due to historical reasons and for backward compatibility of input files. EIRENE had originally by default one single common ion temperature, but one flow field for each bulk ion species).
INDPRO(5): for PITCH (later, in initialization phase, converted into cartesian unit B-field vector BXIN,BYIN,BZIN)
Pitch is defined as B_y/B_tot (LEVGEO=1), B_θ/B_tot (LEVGEO=2), or as B_pol/B_tot (LEVGEO=3), where B_pol is the direction along the polygons.
New options since 2001: In case INDPRO(5)=3 (flat profile) the two redundant input parameters B2, B3 are used to define a constant B-field strength [T], see below under profile type INDPRO=3.
INDPRO(6): for ADIN
INDPRO(7): for WGHT (to be written)
INDPRO(12): for VOL.

For the profiles no. 6 and 7 only the options INDPRO=5 or INDPRO=6 exist.

For the profile no. 12 (cell volumes) the options INDPRO(12) = 4, 5, 6, or INDPRO(12) = 7 are active options. For all other values of INDPRO for these latter four profiles the default profiles described above ("General remarks") are set.

For each profile up to 6 parameters P0 , ..., P5 are read, e.g. TE0, ... , TE5 for the T_e -profile, TI0, ... , TI5 for the T_i -profile, and so on.

Depending upon the value of INDPRO one of the profile routines PROFN, PROFE, PROFS, ... etc. is called from subroutine PLASMA.

INDPRO = 1-4

NR1STM plasma data are defined on the one- dimensional zone-centered grid

RHOZNE(J), J=1, ...,NR1STM (NR1STM = NR1ST - 1 ) .

Vacuum data are specified in zone NR1ST.

These NR1ST data are copied NBLCKS times to define

NBLCKS identical profiles, totally : NBLCKS · NR1ST = NSURF data (subroutine MULTI).

Here the number of copies is NBLCKS = NMULT · NP2ND · NT3RD. see "Input Data for Standard Mesh" (block 2).

INDPRO = 1

(see subroutine PROFN)

RHOSRF(1) ≤ x ≤ P5:

P(x) = P1+( P0 − P1 ) · ( 1 −( [(x − RHOSRF(1))/(P5 − RHOSRF(1))] ) ^P2) ^P3

in particular :

x = RHOSRF(1) → P(x) = P0

x = P5 → P(x) = P1

P5 ≤ x ≤ RHOSRF(NR1ST) :

P(x) = P1 · exp(- (x - P5) / P4)

INDPRO = 2

(see subroutine PROFE)

RHOSRF(1) ≤ x ≤ P5 :

P(x) = P0 ·exp(- (x - P1) / P2)

in particular : x = P1 → P(x) = P0

P5 ≤ x ≤ RHOSRF(NR1ST) :

P(x) = P(P5) ·exp(- (x - P5) / P4)

(P3 : no meaning)

INDPRO = 3

(see subroutine PROFS, PROFS(P0,P1,P5,PVAC)

RHOSRF(1) ≤ x ≤ P5 : P(x) = P0

P5 ≤ x ≤ RHOSRF(NR1ST) : P(x) = P1

Parameters P2, P3, P4 have no meaning, except for the magnetic-field (INDPRO(5)=3): Here PROFS(P0,P1,P5,PVAC) is used to define the pitch angle profile, and it is also used to set the B-field strength (Tesla) by calling PROFS(P2,P3,P5,PVAC)

INDPRO = 4

read from input steam fort.P0

INDPRO = 5

EIRENE calls a user supplied plasma profile routine:

Subroutine PROUSR (RHO, INDEX, P0, P1, P2, P3, P4, P5, PVAC, NDAT).

NDAT data must be defined on

RHO(J), J=1,NDAT e.g. by using the profile parameters P0,...,P5 and any user supplied profile function. See section 3.6.
(NDAT=NSURF)

The flag INDEX is as for the "INDPRO = 6" option, see below, and sections 3.6 and 4 . It determines, which of the background medium profiles has to be provided at a particular call to subroutine PROUSR.

PVAC is the EIRENE vacuum value for the profile in question.

INDPRO = 6

EIRENE calls a the routine:

Subroutine PROFR (RHO,INDEX,N1I,N1,NDAT).

N1I (e.g., NPLSI) profiles, of NDAT= NSURF data each, are read onto RHO(I1I,J), J=1,NDAT, I1I=1,N1I from a work array RWK. N1 (e.g., NPLS) is the leading dimension of the array PRO as specified in the calling program. The data must be written onto array RWK in the initialization phase, e.g., in subroutine INFCOP or MODUSR. The location of a particular profile on this work array is determined by the value of INDEX, see section 4 below. By this option plasma parameters may directly be transferred into EIRENE from other files, e.g. from plasma transport codes.

NDAT (= NSURF) is the number of cells in the "Standard Mesh". Parameters in the additional cell region

ICELL=NSURF+1,...NSURF+NRADD

are set to the "EIRENE vacuum values"

For further details of subroutine INFCOP see 2.1 "Input Data for operating mode", input variable NMODE and sections 2.14 and 4 for an example.

INDPRO = 7

Same as INDPRO=6 option, except that NDAT=NSBOX rather than NDAT=NSURF. Hence in this case background data are read from RWK also for the additional cell region

ICELL=NSURF+1,NSURF+NRADD.

2.5.1 Derived Background Data

Given the data describing the background medium (plasma) in each cell (on common block COMUSR) a set of further, "derived background data" profiles is computed and stored also on COMUSR. This is done by a call to subroutine PLASMA_DERIV. The following quantities are defined there:

DEIN

Electron density (cm⁻3), derived from all NPLSI ion densities and the charge states NCHRGP(IPLS) of each ion species.

EDRIFT

Kinetic energy (eV) in drift motion, for each background species IPLS, derived from the mass MASSP(IPLS) and the flow velocity VXIN,VYIN,VZIN.

LGVAC

Vacuum flag (logical) to identify regions, in which no collision data for all or certain background species are computed, i.e., the collisionalities are set to zero in these cells for those background species. See also explanation above, with respect to the EIRENE vacuum data TVAC,DVAC,VVAC.

NSTGRD

flag for special treatment of some selected cells

=0: default
=1: This cell is a dead cell, not to be seen by the particles (isolated from the computational domain. NSTGRD(ICELL)=1 can be specified in problem specific routines (...USR) or in routines interfacing EIRENE with external codes: INFCOP (code segment: COUPLE....)
=2: indicates that this cell belongs to a grid cut, i.e., is not a valid cell for particle tracing.
=3: indicates that this cell is not a real cell but only the storage position for spatially averaged tallies. E.g.: cell no. NR1ST, 2*NR1ST,....etc. .

2.6 Input Data for Surface Interaction Models

General remarks

An outline of surface interaction models in general terms was given in section 1.4. As for the implementation of such models in EIRENE, there are two parts to the surface interaction data block. The first part contains data which are general to the EIRENE reflection model ("Block for General Reflection Data"). The second part may be different for each surface element and thus must be specified for each reflecting "non-default" surface of the "standard mesh" and for each reflecting "additional surface". It is referred to as "Block for Local Reflection Data". If this block is missing and if the surface MSURF is neither transparent (ILIIN < 0) nor purely absorbing (ILIIN=2) nor a "mirror surface" (ILIIN=3) nor a "periodicity surface", (ILIIN ≥ 4) then the default reflection model (Fe−Target, ... ) is activated for this surface.

The surface interaction model comprises the following information:

p_f (RPROBF): probability for reflection of a "fast" particle
(= 0. for incident molecules),
p_t (RPROBT): probability for re-emission as "thermal particle"
(with surface temperature, see the flags EWALL)
p_a (RPROBA): probability for absorption
(surface sticking, pumping, etc...)
p_s : probability for (physical) sputtering
p_c : probability for (chemical) sputtering
For the reflected "fast particles" and for the sputtered particles the type, species and the distribution of energy, polar and azimuthal angle must be given.
For the re-emitted thermal particles the type, species and only the distribution of energy must be given; the angular distribution follows the cosine law by default.

All these data may be functions of incident energy E_in, of incident polar angle θ_in against the surface normal and of projectile and target species. The three probabilities p_f, p_t and p_a must add up to 1; p_s and p_c are sputtering yields per incident particle and, thus, may occasionally be larger than 1. They are called a probability here only by abuse of language.

All this information is contained in each of three so called "Reflection Models", namely the "Database Reflection model" (see references [15], [16]), the "Modified Behrisch Matrix model" (references [13], [27]) or the "User Supplied Reflection model" (see section 3.3), and in several so called "Sputter Models" ([28], [29], [30].

These models may, to some extend, be modified by the input flags described here. For example, a surface may be purely absorbing (for all test particles incident onto it, ILIIN = 2 option), acting like a mirror for the neutral test particles, by the ILIIN = 3 option (see section 2.3A and 2.3B), or enforce periodicity, e.g., because of symmetry conditions (ILIIN ≥ 4 option). In these two latter cases all test particles of all species are reflected with the species unchanged, with probability one and elastically: E_out = E_in. In case of the "mirror reflection option" the normal component v_n of the incident velocity vector v is reversed: v_n → −v_n. In case of a "periodicity surface" both the position and the speed unit vector are altered according to the specific periodicity of the configuration.

If, instead, a surface reflection model is chosen (ILIIN = 1), then some of the input flags described here (RECYCF, RECYCT, RECYCS, RECYCC) as well as the flags (RINTEG, EINTEG, AINTEG) may be used to modify conveniently the otherwise rather unhandy reflection data and formulas.

The Input Block for General Reflection Data

*** 6A. Data for General Reflection Models

 NLTRIM
``Path-Card" for TRIM surface reflection files"
 (machine specific, may be omitted)

If a "Path Card" is included, then read an arbitrary number of "Target-Projectile Specification Cards" of the format

 A_on_B

Up to NHD6 (see PARAMETER statements, section 3.1) data files for different Target-Projectile combinations can be read.

If no "Path Card" has been specified, or after the "Target-Projectile Specification Cards", continue reading:

 (DATD(I) I=1,NATM)
 (DMLD(I) I=1,NMOL)
 (DIOD(I) I=1,NION)
 (DPLD(I) I=1,NPLS)
 ERMIN  ERCUT  RPROB0  RINTEG EINTEG AINTEG

optionally, an arbitrary number of surface models, labelled by the character string ′modname′, may be read next.

SURFMOD_modname
    ILREF  ILSPT ISRS  ISRC
    ZNML EWALL EWBIN TRANSP(1,N) TRANSP(2,N) FSHEAT
    RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
    RECYCS  RECYCC  SPTPRM  (this line may be omitted,
                             then: default sputter model, see below)
next: within each SURFMOD-sub-block there may be an
arbitrary number of lines
VARNAME  SPZNAME   VALUE
(format: CHARACTER*8,X,CHARACTER*8,*)

Meaning of the Input Variables

NLTRIM: TRIM database is used, if "Database Reflection Model" is specified in at least one block for local reflection data. Data are read from data-set FT21 (no "Path Card" specified, old option), or from the domain specified by the "Path Card" (new option, see next card). If the old option is used, then the complete TRIM file is read, containing the first 12 TRIM target-projectile combination data-sets listed in section 1.4, i.e., the files H_on_Fe to T_on_W.
If the parameter NHD6 < 12, (section 3.1) then only the first NHD6 files are read from FT21.
"Path Card": This card is machine specific. If EIRENE finds a card containing the string 'PATH' or 'path', it assumes that this card specifies the path to the domain containing the TRIM surface reflection data files.
A_on_B: Name of a particular TRIM data file in the domain specified by the path card. E.g., H_on_Fe would include the data file for hydrogen onto iron into the EIRENE run. Up to NHD6 such "Target-Projectile Specification Cards" may be included. For a complete list of such files currently available see again section 1.4.

Note: if during a Monte Carlo simulation a projectile A hits a target B, for which no TRIM data file has been specified, but still NLTRIM=TRUE, then EIRENE searches the "closest" of all its TRIM data files (with respect to a reduced mass argument) and uses this target-projectile combination together with a reduced mass scaling of incident particle energy. Hence, a TRIM data set is chosen such that the reduced mass scaling factor is as close to one as possible amongst the files available.

DATD: distribution for sampling the species index of reflected or otherwise emitted atoms. The NATMI relative frequencies
DATD(IATM) IATM = 1,NATMI
are used to produce the corresponding cumulative distribution DATM in order to facilitate sampling (inversion method).
DMLD: as for DATD, but for molecules. Cumulative distribution is DMOL.
DIOD: as for DATD, but for test ions. Cumulative distribution is DION.
DPLD: as for DATD, but for bulk ions. Cumulative distribution is DPLS.

The activation of these distributions at surface events during the particle history generation process is controlled by the surface species index flags ISRF$, ISRT$ read in the species sub-blocks of section 2.4 and 2.5. By default the following conventions are used (and can be overruled by calls to subroutines REFUSR, SPTUSR only, see section 3.3):

ISRF$

fast particle reflection species flag:

> 0: If ISRF$ ≤ NATMI, then ISRF$=IATM, the species labelling index for the reflected fast atom.
If ISRF$ > NATMI not in use, warning and error exit.
= 0: no fast particle reflection for this species: p_f = 0
< 0: not in use, warning and error exit.

ISRT$

thermal particle re-emission species flag:

incident atoms, test ions and bulk ions:

> 0: If ISRT$ ≤ NATMI, then ISRT$ = IATM, the species labelling index for the re-emitted thermal atom.
If ISRT$ > NATMI, then IATM is sampled from the distribution DATM(IATM).
= 0: no thermal particle re-emission for this species: p_t = 0
< 0: molecule is re-emitted, and -ISRT$ is used to identify the species labelling index IMOL for the re-emitted molecule, exactly as described above for re-emitted atoms. The relevant sampling distribution for species index IMOL in case -ISRT$ > NMOLI is DMOL(IMOL).

incident molecules:

> 0: If ISRT$ ≤ NMOLI, then ISRT$ = IMOL, the species labelling index for the re-emitted thermal molecule.
If ISRT$ > NMOLI, then IMOL is sampled from the distribution DMOL(IMOL).
= 0: no thermal particle re-emission for this species: p_t = 0
< 0: not in use, warning and error exit

These latter two surface-species index flags are considered to be particle properties, hence they are read in the particle specification blocks 2.4 and 2.5. The flags in the next card can be used to modify the pre-programmed fast particle reflection models, to some extend at least.

ERMIN: For incident particle energies below ERMIN, the "fast" particle reflection model is switched off. Only the "thermal" particle model is used.
ERCUT,RPROBF: These variables may be used to modify the default "Behrisch Matrix" reflection coefficients for particles incident on a surface at low energies E_in. The original data [13] are used only for E_in > ERCUT and for normal incidence ϑ_in = 0 .
In the range ERMIN < E_in < ERCUT the particle reflection coefficient p_f (E_in ,ϑ_in = 0) is replaced by a smooth cubic interpolation curve p_f (E_in) such that p_f (0) = RPROBF .
The original "Behrisch Matrix" is recovered by setting ERCUT ≤ 0.

The next three flags modify the particle, energy and momentum reflection coefficients, and/or the resulting distributions in post-reflection energy and angle, respectively. These flags only apply to incident particles of type 1, 3, 4 and 5, i.e., not for incident molecules because for those RPROBF = 0 by default.

RINTEG

> 0: Fixed (independent of energy and angle of incidence) particle reflection coefficient.
The fast particle reflection probabilities RPROBF are set to p_f=MIN(1−p_a,RINTEG), regardless of the reflection model selected by the flag ILREF in block 6B. p_a is kept as specified, and p_t is then recomputed as p_t = 1 − p_f − p_a.
RINTEG ≥ 1 − p_a enforces the fast particle reflection model for all un-pumped incident particles, i.e., RINTEG is internally reset to 1 − p_a.
= 0: Default: fast particle reflection probability as defined by the reflection model chosen.
< 0: The fast particle reflection probabilities RPROBF are set to 1.0, i.e., even pumping is turned off (as distinct from the choice RINTEG=1.0, which would preserve the pumping speed at a surface).

EINTEG

> 0: Fixed (independent of energy and angle of incidence) energy reflection coefficient.
This choice replaces the energy random sampling procedure in the fast particle reflection model by an reflection assumption: E_out = E_in ·EINTEG.
= 0: Default: no modification of energy distribution in the reflection model for fast particle reflection.
< 0: elastic (E_out = E_in) reflection for all particles reflected according to the fast particle reflection model. Hence: same as EINTEG=1.0

AINTEG

(not ready to use)

> 0: Fixed (independent of energy and angle of incidence) momentum reflection coefficient.
This choice replaces the angle random sampling procedure in the fast particle reflection model by a momentum reflection assumption such that, on average,: v_out = v_in ·AINTEG.
= 0: Default: no modification of angular distributions in reflection model for fast particle reflection.
< 0: not in use

Note:

Some "General Reflection Model" data of input block 6A are copied identically NLIMI + NSTSI times onto appropriately dimensioned arrays (with the same names) in order to make them dependent upon the surface labelling index as well.

The first index is then the species index, the second index is the surface number.

Presently this "localization option" is available for the flags

to be written

RF0USR or SP0USR of the user supplied reflection routine REFUSR; RF0USR and SP0USR are called in the initialization phase of an EIRENE run (see section 3.3).

By a similar strategy, some of the species independent flags in the next sub-block 6B can be made species-dependent. See, for example, RECYCF, RECYCT, RECYCS and RECYCC below.

*** 6B. Data for Local Reflection and Sputtering Models

The next 3 (or 4) lines comprise the sub-block for local reflection data. Such sub-blocks can be included in each "surface deck" (block 3A, 3B) to overwrite the default reflection model for each individual surface.

ILREF ILSPT ISRS ISRC

ZNML EWALL EWBIN TRANSP(1) TRANSP(2) FSHEAT

RECYCF RECYCT RECPRM EXPPL EXPEL EXPIL

RECYCS RECYCC SPTPRM (this line may be omitted
then: default sputter model, see below)

An arbitrary number of such sub-blocks of 3 (or 4) lines each may also be defined under a certain label "SURFMOD_modname" and be included in the input file in block 6 directly after the ERMIN,... deck

In blocks 3a and 3b surfaces may be assigned a particular local reflection model by a card reading SURFMOD_modname. Many surfaces may then be linked to the same surface reflection sub-block.

Example:

SURFMOD_BERYL_SPT_300K
     1     2     0     0
 9.04000E+02-2.60000E-02 0.00000E+00 0.00000E+00 0.00000E+00 2.80000E+00
 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00
 1.00000E+00 1.00000E+00 1.00000E+00
SURFMOD_CARB_SPT_1153K
     1    12     0     0
 1.20600E+03-1.00000E-01 0.00000E+00 0.00000E+00 0.00000E+00 2.80000E+00
 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00
 1.00000E+00 1.00000E+00 1.00000E+00
TRANSP1  D2      0.5000E 00
TRANSP2  T2      0.2500E 00
RECYCF   D       0.0000E 00
RECYCT   D2      0.9500E 00
SURFMOD_CARB_SPT_812K
     1     2     0     0
 1.20600E+03-7.00000E-02 0.00000E+00 0.00000E+00 0.00000E+00 2.80000E+00
 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00
 1.00000E+00 1.00000E+00 1.00000E+00

Meaning of the Input Variables

ILREF

Flag for choice of local reflection model

= 1: TRIM database reflection model is used. NLTRIM must be .TRUE..
= 2: "modified Behrisch Matrix model" is used
= 3: user supplied reflection model (see section 3.3: Subroutine REFUSR)

Default: ILREF = 2

ILSPT

Flag for choice of local sputtering model

Let ILSPT = MN, with M and N single digit integers each. Then N controls the options for physical sputtering, and M controls chemical sputtering. See subroutine SPUTER.

N = 0: no physical sputtering at this surface
N = 1: constant physical sputtering rate (see parameter RECYCS below)
N = 2: modified Roth-Bogdansky formula for sputter yield, Thompson energy distribution and cosine angular distribution for emitted particles (see references [29] and [30]).
N = 9: (was option N=3 in Eirene-2004 and older)
user supplied sputtering model (see section 3.3: Entry SPTUSR to subroutine REFUSR)
M = 0: no chemical sputtering at this surface
M = 1: constant chemical sputtering rate (see parameter RECYCC below)
M = 2: "Roth formula" for chemical sputter yield, thermal distribution for emitted particles (see reference [31]), "weak flux dependence option A6".
M = 3: "Roth formula" for chemical sputter yield, thermal distribution for emitted particles (see reference [31]), "strong flux dependence option A7".
M = 4: "Roth formula" for chemical sputter yield, thermal distribution for emitted particles (see reference [31]), "new flux dependence option A8 (2004)".
M = 5: not in use
M = 6: "Haasz-Davis 1998 formula" for chemical sputter yield
M = 7: "Haasz-Davis 1998 formula" for chemical sputter yield, and multiplicative factor for flux dependence (Roth, 2004).; M = 9 (was option N=3 in Eirene-2004 and older)
user supplied sputtering model (see section 3.3: Entry SPTUSR to subroutine REFUSR)

Default: ILSPT=0

The next two surface-species index flags ISRS$ and ISRC$ control the species of sputtered particles. Hence they are considered surface properties and are read in the surface specification decks in blocks 2.3.1, 2.3.2 and/or 2.6.

ISRS$

sputtered particle species flag (physical sputtering).

> 0: both the sputtered particle and the reflected particle (if any) will be followed. Their contribution to surface particle and energy fluxes is stored in surface averaged tallies 1 to 24, i.e., sputtered particles are not explicitly distinguished from reflected particles in the balances. Furthermore the "sputtered flux surface tallies" 25 to 28 are updated. The species index of the sputtered particle (atom) is IATM=ISRS. Hence, on input, ISRS ≤ NATMI.
= 0: There is no physical sputtering for the particular surface element and incident species (note: ISRS=ISRS(ISPZ,MSURF), i.e., p_s = 0 here.
Only the reflected particles are followed. Their contribution to surface particle and energy fluxes is stored in surface averaged tallies 1 to 24.
< 0: Same as ISRS > 0, however, the species index of the sputtered particle is determined automatically from comparing the charge and mass numbers of the atomic species (input block 4a) with the corresponding data of the surface element.

ISRC$

sputtered particle species flag (chemical sputtering).

> 0: both the sputtered particle and the reflected particle (if any) will be followed. Their contribution to surface particle and energy fluxes is stored in surface averaged tallies 1 to 24, i.e., sputtered particles are not explicitly distinguished from reflected particles in the balances. Furthermore the "sputtered flux surface tallies" 25 to 28 are updated. The species index of the sputtered particle (atom) is IATM=ISRC, if ISRC ≤ NATMI, or (molecules) IMOL, if ISRC = NATMI+IMOL and NATMI < ISRC ≤ NATMI+NMOLI. Hence, on input, ISRC ≤ NATMI+NMOLI.
= 0: There is no chemical sputtering for the particular surface element and incident species (note: ISRC=ISRC(ISPZ,MSURF), i.e., p_c = 0 here.
Only the reflected particles are followed. Their contribution to surface particle and energy fluxes is stored in surface averaged tallies 1 to 24.
< 0: Same as ISRC > 0, however, the species index of the sputtered particle is determined automatically from comparing the charge and mass numbers of the atomic species (input block 4a) with the corresponding data of the surface element. I.e., in case of Carbon surfaces the sputtered particle is a C-atom, if such an atom has been specified in input block 4a

ZNML

= KLMN (4 digits)

KL: atomic weight of wall material. Note: the nearest integer of the mass number in the TRIM runs is used. For example, a copper target is specified in the TRIM files with an atomic weight of 63.54, and the corresponding TRIM file is used for surfaces with KL=64.
MN: nuclear charge number of wall material

Example: Carbon: ZNML=1206.

Example: Molybdenum: ZNML=9642.

Example: Copper: ZNML=6429.

Default: ZNML = 5.626E3 (stands for Fe).

EWALL

< 0: -EWALL = TW is a (surface-) temperature (eV) in a Maxwellian flux distribution for the thermal particle energy. The resulting mean energy is E_mean = 2 ·TW.
> 0: +EWALL = Energy of mono-energetic (thermal) particles.
The relation between surface temperature TW and the mean energy of particles then reads EWALL = E_mean = 1.5 ·TW.
= 0: Energy is sampled from a Thompson distribution, using the flag EWBIN (see below) as parameter for the surface binding energy

Default: EWALL = +0.0388 ( ≅ TW = 0.026 eV ≅ 300 K)

Note that the EWALL > 0 option enables EIRENE to include boundary conditions in "one speed transport equation" approximations, which often are of great interest in general linear transport theory.

EWBIN

see above, EWALL = 0 option.

Default: EWBIN = 0.0 (irrelevant for "Default Model")

TRANSP(1)

Semi-transparency for particles incident from the positive side. Renders a non-transparent surface (ILIIN > 0) semi-transparent. The probability for passing through the surface is TRANSP. Hence: the probability for reflection/re-emission etc. is 1-TRANSP.

Default: 0.0 (i.e., fully reflecting surface).

Irrelevant for transparent surfaces.

TRANSP(2)

Semi-transparency for particles incident from the negative side on a non-transparent surface.

Default: 0.0 (i.e., fully reflecting surface).

Irrelevant for transparent surfaces.

FSHEAT

surface sheath potential factor. The sheath potential is FSHEAT·T_e, with T_e the electron temperature at the point of incidence. This sheath potential is applied if ions (test ions or bulk ions) hit a non-transparent surface.

If FSHEAT ≤ 0.0, then a sheath potential computed from the local background plasma flow conditions is used (function SHEATH), assuming ambi-polar flow, a Boltzmann distribution for electrons and zero secondary electron emission. See section 1.5. In case of zero background plasma flow velocity at the place of incidence, a default of FSHEAT = 2.8 is used (corresponding to T_e = T_i, and a single ion species D⁺ plasma flowing at ion acoustic speed parallel to the [(B)\vec]-field into the sheath.

Default: FSHEAT = 0.0

RECYCF

Multiplier for reflection probability RPROBF: The probability

p_f = RPROBF for the fast particle reflection model is modified to

RPROBF = AMIN( RECYCF ·RPROBF, RECYCT) ,

where RPROBF was evaluated from the reflection model specified by ILREF.

The total recycling coefficient RECYCT = p_f + p_t is unchanged. Hence, by the use of RECYCF the probability for thermal particle emission p_t is altered.

Default: RECYCF = 1 for incident atoms, test ions and bulk ions.

Default: RECYCF = 0 for incident molecules.

RECYCT

Recycling coefficient (must not be negative):

A fraction p_f = RPROBF of the incident particles is reflected as described by the "fast particle reflection model".

The fraction p_t = RPROBT (RECYCT - RPROBF) will be re-emitted by the "thermal particle reflection model".

The fraction p_a = (1 - RECYCT) will be absorbed at the surface.

The flag RECYCT is used to specify surface pumping in the following way:

Let A be the surface area [cm²], then the pumping speed L [liters/s] for particles with temperature T [K] and mass m [AMU] is

L = A ·(1 − RECYCT) ·3.638 ·

√

T/m

If the surface area is given in [m²], L in [m³/s], then the numerical factor becomes 36.38 .

Note: RPROBF = 0 for incident molecules (ITYP=2) by default.

Default: RECYCT = 1. , i.e. p_a = 0

RECPRM

free model parameter for user supplied recycling models ILREF=9.

Default: RECPRM = 0.

EXPPL

(only for ILREF = 2 option)

angular dependence of fast particle reflection coefficient. The formula

R( Φ) = 1 − (1 − RPROBF) ·cos^EXPPL (Φ)

is used, e1 = EXPPL, see equation 4.38, where:

RPROBF: Reflection probability from "Behrisch Matrix" model, which is valid only for normal incidence.; Φ Angle of incidence against surface normal; R(Φ) Reflection probability for particles incident with angle Φ
note: R(Φ) = 1 for Φ = 90^° and EXPPL > 0.

Default: EXPPL = 1. (recommended from a comparison with the TRIM database)

EXPEL

as EXPPL, but for the energy reflection coefficient.

Default: EXPEL = 0.5 (recommended from a comparison with the TRIM database). e2 = EXPEL, see equations 4.39 and 4.40

EXPIL

Index for angular distribution of re-emitted atom or molecule; affects both models ILREF = 1 and ILREF = 2

EXPIL = 0: cosine distribution, independent of choice of fast particle reflection model
EXPIL > 0: mixed cosine-specular model.
In case of ILREF=1 the angular distribution given by the database is used.
In case of ILREF=2 the specular contribution increases according to equations 4.41 with angle of incidence Φ and with e3 = EXPIL (recommended: EXPIL ≤ 1).

Default: EXPIL = 0.

RECYCS

The meaning of this flag for the "physical sputtering" options depends upon the value of the first digit N of ILSPT:

N=0: no physical sputtering, YIELD1 = 0. RECYCS is irrelevant.
N=1: constant physical sputtering yield, YIELD1 = RECYCS
N=2: RECYCS is a multiplier for the sputtered particle flux YIELD1. YIELD1 is computed from the incident species, energy, angle and surface parameters by the sputter model N=2. Hence: the sputtered particle yield YIELD1 as computed from subroutine SPUTER is modified to
YIELD1 = RECYCS ·YIELD1 .
N=9: RECYCS is a free model parameter, which can be used in the user supplied sputter model for any particular surface element.
Default: RECYCS = 1.

RECYCC

The meaning of this flag for the "chemical sputtering" options depends upon the value of the second digit M of ILSPT:

M=0: no chemical sputtering, YIELD2 = 0. RECYCC is irrelevant.
M=1: constant chemical sputtering yield, YIELD2 = RECYCC.
M=2: RECYCC is a multiplier for the sputtered particle flux YIELD2. YIELD2 is computed from the incident species, energy, angle and surface parameters by the sputter model M=2. Hence: the sputtered particle flux YIELD2 as computed from subroutine SPUTER is modified to
YIELD2 = RECYCC ·YIELD2 .
M=9: RECYCC is a free model parameter, which can be used in the user supplied sputter model for any particular surface element.
Default: RECYCC = 1.

SPTPRM

free model parameter for user supplied sputtering models N=3, M=3.

Default: SPTPRM = 0.

Note:

All data in a block for "local reflection data" are, in general, independent of the type and species of the incident particle. Some are, however, copied identically NPHOTI+NATMI+NMOLI+NIONI+NPLSI times onto appropriately dimensioned arrays (with the same names). This is currently done for the parameters:

    ISRS  ISRC
    TRANSP(1,...)  TRANSP(2,...)
    RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
    RECYCS  RECYCC  SPTPRM

The first index is the species index, the second index is the surface number.

Overwriting local reflection data for some specific incident particle species can be done via a call to the entry RF0USR of the user supplied reflection routine REFUSR. This entry is called in the initialization phase of an EIRENE run, from subroutine REFLEC. Likewise, the entry SP0USR of the user supplied sputter routine SPTUSR is called in the initialization phase from subroutine SPUTER. For further details on user supplied surface interaction routines see section 3.3.

For example: RECYCT(3,5) is the recycling coefficient for incident species 3 onto surface no. 5.

If surface models are defined in input block 6 by the SURFMOD_modname label, (rather than individually for each surface in blocks 3a or 3b), then these variables may be made species dependent by adding an arbitrary number of lines to a SURFMOD-deck, each of which overwriting the specification for a particular particle species. The name of the species (blocks 4 and 5) must be uniquely determining one of the test particle or bulk species.

In the example of SURFMOD decks given above, e.g., the first such additional line overwrites the original value of TRANSP(1,...) (=0.0) for species D2 with the value 0.5. I.e., all surfaces, to which the reflection model modname=CARB_SPT_1153K is assigned, have a recycling coefficient =1.0 for all incident species, except for the D2 molecules, for which these surfaces are made transparent with probability 0.5, if these molecules are incident from the positive side. For T2 molecules incident from the negative side this transparency is set to 0.25. The other species dependent modifications in this SURFMOD deck are self-explaining.

2.7 Input data for Initial Distribution of Test Particles

General Remarks

The primary source (and: initial distribution of test particles, in time dependent mode) is given as a function Q(i, r, v ,t). Q is the density of the probability distribution from which the species index i , the starting point r , the velocity vector v and the starting time t are sampled in subroutine LOCATE. (i,r , v ,t) is the state of the starting particle, which then will be followed in subroutine FOLNEUT or FOLION. There are four types of sources (: "Strata") , namely Point sources, Line sources (to be written), Surface sources and Volume sources.

As the transport equation is linear, various sources of test particles can be treated subsequently and the responses can then be linearly superposed. EIRENE can handle up to NSTRA (see: PARMUSR, 3.1) different strata and prints output from each single stratum as well as from the sum over strata. The source strength is prescribed for each stratum separately (input card: FLUX(ISTRA),...). Subdividing the total source into such strata can be useful to increase the efficiency of the code, ("stratified source sampling", see any textbook on Monte Carlo integration), or if the contribution of such sub-sources is of interest by itself.

In order to facilitate sampling each stratum can be subdivided further into a number (NSRFSI) of "sub-strata". Random sampling is done by firstly sampling the sub-stratum ISRFSI, and then (conditional) the initial state of the test flight within this sub-stratum.

The total CPU resources (and storage for the census array, see section 2.13) are distributed over the strata. This distribution is controlled by some of the input flags in this input block, as described below. EIRENE stops working on a particular stratum if either the assigned CPU-time for this stratum has been reached (Message: "No further CPU time for this stratum"), or if all requested histories have been calculated (Message: "All histories for this stratum completed") or, in case of time dependent problems, if the storage on the census array of this stratum has been filled up (Message: "Census array filled for this stratum").

The Input Block

 *** 7. Data for primary sources
 NSTRAI
 INDSRC(ISTRA), ISTRA=1,NSTRAI
 ALLOC
 DO 70, ISTRAI=1,NSTRAI
   TXTSOU
   NLAVRP  NLAVRT  NLSYMP  NLSYMT
   NPTS  NINITL  NEMODS  NAMODS
   FLUX  SCALV   IVLSF  ISCLS  ISCLT  ISCL1  ISCL2  ISCL3  ISCLB  ISCLA
 *Species index distribution
   NLATM  NLMOL  NLION  NLPLS NLPHOT
   NSPEZ
 *Distribution in physical space
   NLPNT  NLLNE  NLSRF  NLVOL NLCNS
   NSRFSI
     DO 75, J=1,NSRFSI
       INUM    INDIM   INSOR   INGRDA(1)  INGRDE(1)
                               INGRDA(2)  INGRDE(2)
                               INGRDA(3)  INGRDE(3)
       SORWGT  SORLIM  SORIND  SOREXP  SORIFL
       NRSOR   NPSOR   NTSOR   NBSOR   NASOR   NISOR
       SORAD1  SORAD2  SORAD3  SORAD4  SORAD5  SORAD6
 75  CONTINUE
 Distribution in velocity space
   SORENI  SORENE  SORVDX  SORVDY  SORVDZ
   SORCOS  SORMAX  SORCTX  SORCTY  SORCTZ
 70 CONTINUE

Meaning of the Input Variables for primary sources

NSTRAI

Number of different sources (Strata), which are computed one after the other and are linearly superimposed at the end of the run.

(NSTRAI ≤ NSTRA, see "Parameter-Statements")

INDSRC

INDSRC(ISTRA)=0-5: the input data for stratum ISTRA are read here, but may be modified in some user routine (SAMUSR) or interface routine (INFCOP, at entry IF2COP(ISTRA))
INDSRC(ISTRA)=6: no input data for stratum ISTRA are read here. The definition of this stratum must be entirely in some problem specific routine (IF2COP, etc.). See section 3.4 for one such example, namely the default surface recycling source model as specified in coupled B2-EIRENE runs.
INDSRC(ISTRA)=-1: the input data for stratum ISTRA are read here, and no attempt is made to modify these. I.e. IF2COP(ISTRA) is not called.

ALLOC

Allocation of CPU-time to stratum weighted as

(1-ALLOC)*NPTS+ALLOC*FLUX

TXTSOU

Text to characterize the stratum (name of the source) on the printout file.

NLAVRP

= .TRUE. not in use

NLAVRT

= .TRUE. not in use

NLSYMP

= .TRUE.

Symmetrize profiles with respect to poloidal (y-) co-ordinate x², i.e., with respect to the poloidal surface x²= PSURF((NP2ND+1)/2) in case NP2ND is an odd integer, or with respect to the cell center x²= PZONE(NP2ND/2) in case NP2ND is an even integer.

NLSYMT

= .TRUE.

ditto, but for toroidal (z-) co-ordinate, i.e., for toroidal surface

TSURF((NT3RD+1)/2) or TZONE(NT3RD/2) respectively.

NPTS

NPTS > 0: Maximum number of test particle histories.

If there is more than one stratum (NSTRAI > 1) then the total CPU-time NTIME (input block 1) will be distributed proportional to NPTS to the single strata. NPTS is the maximum number of test-particles only if sufficient CPU time is available. Otherwise a message "NO FURTHER COMPUTATION TIME FOR THIS STRATUM" is printed and the particle loop for the respective stratum is stopped.

NPTS = 0: this stratum is "turned off".

NPTS < 0: no limitation in the number of particles. The entire CPU time assigned to this stratum will be used up. (NPTS is reset to the largest integer on the machine. Hence: some care is needed here in case of multiple strata, in combination with the ALLOC-options to assign CPU time to individual strata).

NINITL

If NINITL > 0 : seed for initialization of random number generator. The results for all those individual strata can be reproduced exactly for which the same number of test-flights is computed as in a previous run.

If NINITL = 0 : no initialization of random numbers for the particular stratum. In case of the first stratum, the default initialization is used. Runs can only be reproduced, if the same number of test-flights is computed for each stratum. Somewhat weakened correlation between subsequent runs as compared to the NINITL > 0 option.

If NINITL < 0: truly random initialization (determined by machine clock). These runs cannot be reproduced exactly. Subsequent runs are uncorrelated. (E.g.: recommended for stochastic approximation procedures in nonlinear applications.

NEMODS

Flag to select one of the preprogrammed source energy conditional distributions given the source particle's position and species. (See: "distribution in velocity space", below)

NAMODS

Flag to select one of the preprogrammed conditional source angular distributions given the position, species and energy of the source particle. (See: "distribution in velocity space", below)

FLUX, SCALV

SCALV=0: (default) FLUX = Source strength in Ampere.
FLUX is the scaling factor for all surface- or volume averaged tallies.
FLUX is an "atomic flux" (or: an "atomic ion flux"). Each source particle may carry a different flux NPRT(ISPZ) (initial weight) depending on the species ISPZ (see below: distribution for the species index). NPRT is specified in the blocks 4 and 5. The total "atomic" source particle flux for each stratum is scaled to be FLUX. For example, a H₂ molecule source, with NPRT_H₂ = 2, is treated as if a flux of FLUX/1.602E-19/2 H₂-molecules per second is emitted, resulting in an equivalent "atomic flux" FLUX/1.602E-19 per second.
SCALV ≠ 0: The default scaling of tallies with FLUX can be overruled by this flag. The common scaling factor for all surface- and volume averaged tallies is determined such that one particular tally (selected by the ISCL..-flags described below) has the prescribed value SCALV. This determines the scaling of all other volume averaged and surface averaged tallies. By this option, for example, one can set the neutral particle density to a prescribed value in one particular cell. Hence, one can prescribe the local Knudsen number for nonlinear applications including neutral-neutral interactions.

IVLSF

=1: The following ISCL..-flags select one particular volume averaged tally
=2: The following ISCL..-flags select one particular surface averaged tally

ISCLS

species index of selected tally

ISCLT

tally number of selected tally (refer to tables 5.2,5.3)

ISCL1, ISCL2, ISCL3, ISCLB, ISCLA

IVLSF=1: cell numbers NRCELL, NPCELL, NTCELL, NBLOCK, NACELL, respectively.
if (NPCELL = 0) or (NTCELL = 0), then ISCL1 = NCELL, the cell number in the 1-dimensional arrays (see end of section 2.2.1).
The additional cell region is specified by NRCELL=0, NPCELL=1, NTCELL=1, NBLOCK=NBMLT+1 (see section /refsec2.2) and the proper value of NACELL.
IVLSF=2: to be written

Distribution for the species index

NLATM: = .TRUE.
Atomic source. History starts in subroutine FOLNEUT with type index ITYP=1, species index ISPZ = IATM and initial weight NPRTA(IATM) (see block 4A)
NLMOL: = .TRUE.
Molecule source. History starts in subroutine FOLNEUT with type index ITYP=2, species index ISPZ = IMOL and initial weight NPRTM(IMOL) (see block 4B)
NLION: = .TRUE.
Test ion source. History starts in subroutine FOLION with type index ITYP=3, species index ISPZ = IION and initial weight NPRTI(IION) (see block 4C)
NLPHOT: = .TRUE.
Photon source. History starts in subroutine FOLNEUT with type index ITYP=0, species index ISPZ = IPHOT and initial weight NPRTPH(IPHOT) (see block 4D)
Not all options for direct photon sources are fully programmed. Currently we mostly use the bulk particle volume recombination source (see next) to simulated radiative decay from excited states as birth profile for (bound-bound) line-photons.
NLPLS: = .TRUE.
Bulk ion source. Initial co-ordinates of a bulk ion with species index ISPZ = IPLS and initial weight NPRTP(IPLS) (see block 5) are generated, then a surface reflection model or a volume re-combination model is called and atoms, molecules or test ions with species index either IATM, IMOL or IION are created.

One and only one of these five variables must be .TRUE. .

NSPEZ

Species index of the source particle

1 ≤ NSPEZ ≤ NATMI, NMOLI, NIONI, NPLSI: NSPEZ is the (fixed) species index of the source particle. No random sampling for the species index is done.
NSPEZ > NATMI, NMOLI, NIONI, NPLSI: (depending upon the type of the source particle) the species index is sampled from the distribution DATM, DMOL, DION, DPLS, respectively. The distributions DATM, DMOL, DION and DPLS are read in the block "Data for General Reflection Models".
NSPEZ = 0: the species index of the particle is directly sampled from the "analog distribution" WEISPZ, i.e., no biased source species sampling. WEISPZ is defined internally by the code.
NSPEZ < 0: The distributions DATM, DMOL, DION and DPLS are considered as biased source species distributions, whereas the analog (physical) distribution is provided automatically by the array WEISPZ from the source sampling routines SAMPNT, SAMLNE, SAMSRF or SAMVOL respectively. An appropriate weight correction is carried out after sampling from DATM, DMOL, DION or DPLS, respectively, in subroutine LOCATE.

The distribution WEISPZ is defined only for NLPLS=TRUE sources, and here only for volume re-combination sources, or for surface re-combination sources (see below, this section). It is set according to the local bulk ion composition.

(This, currently, does not work for volume recombination sources.)

In all other cases NSPEZ must be positive.

Distribution in physical space

NLPNT: = .TRUE. Point Source
NLLNE: = .TRUE. Line Source
NLSRF: = .TRUE. Surface Source
NLVOL: = .TRUE. Volume Source

One and only one of this 4 variables must be .TRUE.

Point source

Flags for the "distribution in physical space" not mentioned here are irrelevant for point sources.

NPNTSI

Total number of different points, over which the starting points for this stratum are distributed.

INUM

irrelevant; labelling index

SORWGT

Relative frequency for starting point labelled INUM The sum of SORWGT for all NPNTSI points is normalized to one internally.

NRSOR

> 0: x- or radial cell number NRCELL of the zone containing the point source.
= 0: NRCELL is found automatically from the "standard mesh" zoning.

NPSOR

ditto, for y- or poloidal cell number NPCELL

NTSOR

ditto, for z- or toroidal cell number NTCELL

NBSOR

standard mesh block number NBLOCK. Defaulted to NBLOCK = 1, if NBSOR ≤ 0

NASOR

additional cell number NACELL, if point source is located outside the standard mesh. Defaulted to NACELL = 0, if at least one of the variables NRCELL, NPCELL or NTCELL are larger than zero.

NISOR

polygon index IPOLG. Meaningless if NLPLG = .FALSE.

SORAD1

x-co-ordinate of source point X0

SORAD2

y-co-ordinate of source point Y0

SORAD3

z-co-ordinate of source point Z0

SORAD4, SORAD5, SORAD6

are the x,y and z coordinates of a vector C = (CRTX,CRTY,CRTZ) which may be used to distinguish one particular direction for the distribution in velocity space (see below). Internally this vector is normalized to length 1. Irrelevant for an isotropic velocity distribution.

Line source

to be written

Surface source

Flags for the "distribution in physical space" not mentioned here are irrelevant for surface sources.

NSRFSI

Total number of different surfaces, over which the starting points for this stratum are distributed.

INUM

irrelevant; labelling index

INDIM

= 0: source on "additional surface" ASURF (see block 3B)
= 1: source on "standard surface" RSURF, x- (or radial) mesh
(see block 2A and 3A)
= 2: source on "standard surface" PSURF, y- (or poloidal) mesh
(see block 2B and 3A)
= 3: source on "standard surface" TSURF, z- (or toroidal) mesh
(see block 2C and 3A)
= 4: source on a surface composed of one or more segments of radial and/or poloidal polygons. The further details of the spatial distribution are defined in code coupling routines, i.e., the code coupling routine IF1COP must be called. Special versions of IF1COP are available, e.g., in the code segments COUPLE_B2, COUPLE_B2.5 (coupling to BRAAMS multi-fluid plasma code) , COUPLE_DIVIMP (coupling to DIVIMP impurity ion kinetic transport code) or COUPLE_Ufile (TRANSP-code format).
The position on the surface is sampled from a (piecewise constant) step function defined from plasma fluxes onto that surface vs. arclength. The flags INSOR, INGRDA and INGRDE described below are set automatically in this option and hence need not be specified.

INSOR

number of the surface in the mesh RSURF, PSURF, TSURF or ASURF respectively. (Redundant in case INDIM=4)

INGRDA, INGRDE

same as IRPTA, IRPTE flags in input block 3a. Defines subrange on standard surfaces, on which the source is distributed. Irrelevant for sources on additional surfaces.

SORWGT

Relative frequency for starting points on surface labelled INUM. The sum of SORWGT for all NSRFSI surfaces is normalized to one internally.

SORLIM

= KLMN

if SORLIM ≤ 0, the user supplied Subroutine SAMUSR is called to sample all 3 initial co-ordinates (X0,Y0,Z0), see section 3.4.

if SORLIM > 0, then one of the preprogrammed options is used (the digits L,M and N are relevant only for surface sources). In this case:

N: Index to select one of the preprogrammed distributions in radial or x-direction on the surface.
M: Index to select one of the preprogrammed distributions in poloidal or y-direction on the surface.
L: Index to select one of the preprogrammed distributions in toroidal or z-direction on the surface.
K: Index to select one of the preprogrammed distributions for the starting time.
M,N,L = 0: The respective co-ordinate is computed from the 2 others and from the equation for surface number INSOR.

Thus, one and only one of these 3 digits must be equal to 0, because the birth-point for a surface source is determined already by two co-ordinates and the labelling o index of the surface.

If INDIM=0, any one of the 3 digits can be the 0, depending upon the particular equation for the surface ASURF(INSOR).

In case INDIM=1, one has to set N=0, and the poloidal (or y) and toroidal (or z) co-ordinate is sampled according to the flags M and L.

In case INDIM=2 one must specify M=0, and L=0 has to be set in case INDIM=3.

L,M,N = 1: δ-distribution at (a+b)/2
L,M,N = 2: Uniform distribution on the interval [a,b]
L,M,N = 3: Truncated exponential decay with decay length λ on the interval [a,b]. I.e. the sampling distribution reads:
f(x) = c ·exp( − x/ λ) if x ∈ [a,b] and f(x) = 0 elsewhere,
with normalizing constant
c = { λ(exp[ −a/ λ] − exp[−b/λ])}⁻¹
L,M,N = 4: Step-function (see below: Function STEP, subsection 2.7.1) (only one of either L or M or N should be 4)
K = 1: δ-distribution at TIME0 for time of particle birth. (A delta function source in time for the kinetic equation in integral form corresponds to an initial condition for time-dependent linear kinetic integro-differential equation).
K = 2: Uniform distribution in [TIME0,TIME0+DTIMV] for time of particle birth.
Default: K=2 in time-dependent mode (NTIME > 0) and K=1, TIME0=0 in time-independent mode (NTIME = 0), see section 2.1.

SORIND

Flag to choose one from the various step functions (option 4), which have been defined in the initialization phase. Up to NSTEP (PARMUSR, see section 3.1) step functions can be described there. SORIND is the labelling index of the selected step function.

Each step function ISTEP can consist of step functions for fluxes of up to NSPZ species, see 2.7.1. NSPZ depends upon the initialization of this function. By default the source species index NSPEZ is used when sampling from step functions.

New option (Aug. 2006): If SORIND ≥ 100, then the 3rd digit is used to choose the species index from step function ISTEP. I.e.: Let SORIND = LMN, then NM is used to sample from flux distribution for species ISPZ = L from step function ISTEP = NM.

SOREXP

Decay length λ in the exponential distribution (option 3)

SORIFL

The first of the 4 digits can be used to overrule the default orientation of the surface normal at the birth point, or if ILSIDE = 0 for this particular surface. If this digit is nonzero, a value 1 would lead to a test flight originating from the surface as if a particle has been incident onto this surface in the positive direction, and the value 2 means that this imaginary particle has been striking in the negative direction.

The last 3 digits of SORIFL act as LMN of the ILSWCH flag described in section 2.3B assuming incidence in the positive direction.

If any of this 3 digits equals zero, than the ILSWCH flag for this particular surface is activated. (See also: section 2.3B, input flag ILSWCH)

NRSOR