# Consistent Algorithms Marching Along Characteristics for the Numerical Solution of the Boltzmann Equation

Authors: Nilson C. Roberty; Rogerio C. Nunes
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

We introduce algorithms marching over a polygonal mesh with elements consistent with the propagation directions of the particle (radiation) flux. The decision for adopting this kind of mesh to solve the one-speed Boltzmann transport equation is due to characteristics of the domain of the transport operator which controls derivatives only in the direction of propagation of the particles (radiation) flux in the absorbing and scattering media. This a priori adaptivity has the advantages that it formulates a consistent scheme which makes appropriate the application of the Lax equivalence theorem framework to the problem. In this work, we present the main functional spaces involved in the formalism and a description of the algorithms for the mesh generation and the transport equation solution. Some numerical examples related to the solution of a transmission problem in a high-contrast model with absorption and scattering are presented. Also, a comparison with benchmarks problems for source and reactor criticality simulations shows the compatibility between calculations with the algorithms proposed here and theoretical results.


---

## Body

## 1. Introduction

A line propagation of radiation and particles is observed for low-density noncharged particles as neutron in nuclear reactors and X-rays orγ-rays in nuclear medicine and engineering nondestructive test applications. For charged particles, this collimation is found in accelerator and spallation target technologies for ADS nuclear plants. This natural behavior of collimated radiation and particles introduces in the numerical computation a new paradigm of operator, the one-speed Boltzmann transport operator, which has been extensively studied in the nuclear community Vladimirov [1], Duderstadt and Martin [2], Kaper et al. [3], Lewis and Miller Jr. [4], Dautray and Lions [5], Mokhtar-Kharroubi [6] and astrophysical community Chandrasekhar [7]. Outside the nuclear community, the analysis of radiation transfer and particle transport problems has several relevant applications in different areas such as population dynamics, heat transfer, remote sensing, global warming models, natural water radiative properties estimation and transmission, scattering and optical tomography.The huge number of degree of freedom associated with the seven-dimensional phase-space-time domain introduces the necessity of numerically preprocessing the equation in order to produce a solution computationally feasibly. This may introduce conceptual errors in the proposed approximated models. In this work, we are trying to avoid one of the most common mistakes in these models, that is, the disrespect with the natural characteristic of the solution, and also, develop a procedure for discriminating radiation flux of different orders of magnitudes. In image processing, the lower-order magnitude radiation flux is treated as noise, and if we have collimated radiation flux, angular averaged numerical methods will introduce errors. For planes and axis symmetric problems involving the linear Boltzmann transport equation, the adoption of a polygonal mesh with elements consistent with the directions of particles and radiation flux gives an a-priori adaptability that is necessary in the study of the propagation of absorbing and scattering media. This kind of mesh has been called, by us, a natural pixel partition of the domain of definition and has been used successfully in problems of image reconstruction from projection Reis and Roberty [8] and Montero et al. [9]. The use for the solution of the Boltzmann equation is now introduced by me and my students and collaborators Chaffin and Roberty [10], Montero [11], Montero et al. [12]. In the present work, we introduce a new formalism for the radiative transfer and linear particle transport problem. The procedure has been studied in two dimensional models, but extension for 3-dimensional problems is straightforward. Recent works explore the respect with the characteristic of the linear Boltzmann equation, but, contrary to the presented work, they are based on the ray tracing technique and not in the domain partition, as this work is. We cite Liu et al. [13], Chetaine et al. [14], Jevremovic et al. [15].This paper will be structured with the following sequence. The physical formulation for the stationary model is introduced in Section2. Some mathematical results that are important to the understanding of the solution given to this numerically implemented problem are presented in Section 3. The consequences of the necessity of a special basis that rotates to accomplish the regularities of the angular flux are discussed in Section 4. We introduce the implementation of the rotation pixel and the natural partitions elements in Section 5. The basis generated by this partition is the most relevant and original aspect found in the present work. Questions related to the algorithm implementation are introduced. The two-dimensional discrete equation used for algorithm implementation is shown in Section 7. The algorithms for the transillumination, source, and the eigenvalue problem are detailed in Section 7.1. The numerical results for the three simulated problems are presented and discussed in Section 8. Finally, we conclude and point out the advances introduced by the present work.


## 2. Physical Formulation for the Stationary Problem

The particle Boltzmann transport equation is a balance of particles (or photons) inside an absorbing, scattering, or multiplicative participating medium such as a reactor, a tomography, or the natural environment. To state the stationary equation, consider ann-dimensional convex region Ω,n=2,3, delimiting an region with internal radiation sources S(ω,x),(ω,x)∈Sn-1×Ω, subjected to externally generated parallel beams of particles. Here, Sn-1 is the set of directions in Rn. For the steady-state situation, with one-velocity particles or no spectral dependency in the case of radiation, the following formulation is obtained from the linear Boltzmann equation, Vladimirov [1], Kaper et al. [3], Dautray and Lions [5], which is usual for the mathematical modeling of the interaction of the radiation with the participating medium. Based on these initial assumptions, we consider the following boundary problem: for a given function g-, find the constant velocity radiation intensity ϕ∈L2(Sn-1×Ω) such that(2.1)ω⋅∇ϕ(ω,x)+σt(x)ϕ(ω,x)=∫Sn-1f(x,ω′⋅ω)ϕ(ω′,x)dω′+S(ω,x),(ω,x)∈Sn-1×Ω,(2.2)ϕ(ω,σ)=g-(ω,σ)for(ω,σ)∈Σ-,
where Σ±={(ω,σ)∈SN-1×Γ;±ν(σ)·ω>0} are, respectively, the phase space surfaces of the incident and emergent radiation on the physical surface Γ=∂Ω, which is the smooth boundary of the bounded, possibly convex region Ω⊂ℝn.The coefficientsσt and f in (2.1) are, respectively, the total extinction coefficient (absorption + outscattering) and the scattering and fission coefficient. We assume that σt∈L∞(Ω) and, since ω′·ω=cosθ0, with θ0 being the angle between the direction of incident radiation ω′ and the emergent scattered radiation ω, we also assume that(2.3)f(⋅,s)∈L∞(Ω),a.e.s∈[-1,1],f(x,⋅)∈L1([-1,1]),a.e.x∈Ω.
Under these assumptions, L2(Σ±;|ν(σ)·ω|dσdω) become the appropriated spaces for traces of the functions in the boundary value problem (2.1) and (2.2). For a more detailed discussion, see  Cipolatti [16] and the references therein.


## 3. Mathematical Preliminaries

The stationary transport (2.1) is composed of the following operators. The operator(3.1)A:W(Sn-1×Ω)⟶L2(Sn-1×Ω)ϕ(ω,x)⟼A[ϕ](ω,x)=ω⋅∇ϕ(ω,x),
where derivatives are in the sense of distributions. The space H=L2(Sn-1×Ω) is the closure of C(Sn-1×Ω¯) with respect to the norm (3.2)‖ϕ‖H=(∫Sn-1∫Ω|ϕ(ω,x)|2dxdω)1/2,
where dω denotes the measure on Sn-1 associated with the Lebesgue measure ℝn. The space W is defined from H as (3.3)W={ϕ∈H;Aϕ∈H}.
W is a Hilbert space for the norm (3.4)‖ϕ‖W=(∫Sn-1∫Ω[|ϕ(t,ω,x)|2+|Aϕ(t,ω,x)|2]dωdx)1/2.The second operator(3.5)B:H⟶Hϕ(ω,x)⟼[Bϕ](ω,x)=σt(x)ϕ(ω,x)
is a bounded operator which has continuous inverse for σt(x)>0,x∈Ω. The total macroscopic cross-section is, in applications, piecewise smooth (piecewise and infinitely differentiable), bounded, and strictly positive on Ω with (3.6)0<inf⁡{σt(x);x∈Ω}≤σt(x)≤c0=sup⁡⁡{σt(x);x∈Ω}.
The streaming operator is the sum(3.7)Lσ=A+B:W(Sn-1×Ω)⟶L2(Sn-1×Ω)ϕ(ω,x)⟼Lσ[ϕ](ω,x)=ω⋅∇ϕ(ω,x)+σt(x)ϕ(ω,x).The third operator(3.8)C:H⟶Hϕ(ω,x)⟼[Cϕ](ω,x)=-∫Sn-1f(x,ω′⋅ω)ϕ(ω′,x)dω′
describes the anisotropic scattering and fission process. Kernel f can be expanded in an absolutely and uniformly convergent series (3.9)f(x,ω⋅ω′)=∑k=1∞2k-14πσk(x)Pk-1(ω⋅ω′)+1λσf(x),
where Pk-1 is the surface harmonic of degree k-1 (Legendre polynomials). The coefficients σk, k≥1 are nonnegative, and the first coefficient verifies the inequality (3.10)0≤inf⁡⁡{(σt(x)-σ1(x);x∈Ω)}≤c0-c1,wherec1=inf⁡⁡{(σ1(x);x∈Ω)},
and (1/λ)σf is the nonnegative fission cross section divided by the criticality parameter. With this, the collision operator takes its usual form(3.11)C[ϕ](ω,x)=∫Sn-1(σs(ω⋅ω′,x)+1λσf(x))ϕ(ω′,x)dω′.The operatorC with this kernel is compact from L2(Sn-1)↦L2(Sn-1), but not from L2(Sn-1×Ω)↦L2(Sn-1×Ω), and this is the main mathematical problem in this engineering problem.The trace operator(3.12)γ-:W̃={ϕ∈W:∫Σ|ϕ(ω,σ)|2|ω⋅ν(σ)|dσdω<+∞}⟶L2(Σ-,|ω⋅ν(σ)|dσdω)
is continuous and surjective, where W̃ is a Hilbert space with the norm ∥ϕ∥W̃2=∥ϕ∥W2+∫Σ|ϕ(ω,σ)|2|ω·ν(σ)|dσdω.Given the incident flux by the boundary condition (2.2),  g-∈L2(Σ-,|ω·ν(σ)|dσdω), the following prolongation of g- inside the phase space is defined:
(3.13)Jσ[g-](ω,x)=exp⁡⁡(-∫0τ-(ω,x)σt(x-sω)ds)g-(ω,x-ωτ-(ω,x)),
where τ±(ω,x)=sup{t≥0:x±tω∈Ω}are the distance to the emergent and incident boundary, respectively. Note that this definition implies that Jσ is the solution operator of the problem Lσ[ϕ]=0,ϕ|Σ-=g-.These operators can be used to formulate the first-order boundary non-homogeneous stationary direct problem given (2.1) and (2.2) as the following equivalent homogeneous problem: to find φ=ϕ-J[g-]∈W∩dom(A) such that(3.14)(Lσ-C)[φ]=S-C[Jσ[g-]]=h∈H,
where dom(A)={φ∈W:φ=0onΣ-}.Theorem 3.1.
For eachh∈H, the boundary value problem (3.14) has a unique solution witch is given by
(3.15)ϕ=(I-Lσ-1C)-1Lσ-1[h],
and verifies the following estimate:
(3.16)‖ϕ‖≤(1-exp⁡(-c0d))‖(I-(Lσ)-1S)-1‖‖h‖,
where d is the diameter of Ω.Proof.
See  Chapter (III) and Lemma 3.1 of Vladimirov [1].In Sections IV thought VII of Vladimirov [1], operators Lσ:D(Lσ)↦L2(Sn-1×Ω) and C:L2(Sn-1×Ω)↦L2(Sn-1×Ω) are investigated for more general Sobolev spaces. We can easily adapt them to our Hilbert space framework.Lemma 3.2.
The operatorLσ-1∘C:L2(Sn-1×Ω)↦L2(Sn-1×Ω) is compact.Proof.
See Lemma 5.3 of Vladimirov [1].This result is fundamental for the reformulation of boundary value problem given by (2.1) and (2.2) in terms of an integral equation based on a compact perturbation of the identity operator used in this work: to find ϕ(ω,x) such that(3.17)ϕ=Lσ-1∘C[ϕ]+Jσ[g-],(ω,x)∈Sn-1×Ω.Remark 3.3.
The study of the numerical behavior of compact operators is an very well known subject. The spectral properties of compact operators are very similar to that of finite-dimensional operators, fact that make analysis very easy. For this kind of operator, the subspaces associated with a eigenvalue are finite dimensional. It is based on this behavior our interest in model this problem using integral equation (3.17). A situation similar also occurs in models based on the Laplacian operator paradigm, which fortunately, for bounded domains, is an operator with compact resolvent. Since the one-velocity transport equation has a more complex structure than the Laplacian, it is necessary to do the composition of the inverse of the streaming operator, Lσ-1, and the collision operator, C, to obtain a compact operator. Note that Lσ-1 presupposes and line integrals through the ray paths of particles inside the domain, and is this fact that will have a remarkable influence on the numerical behavior of the problem. We will see that this connects compactness and characteristic basis.


## 4. Relation between the Integral Operator and the Characteristic Basis

Due to its very nature, the transport equation for the angular flux may be viewed as a system of ordinary differential equations coupled by the collision operator. Each angular field is a different field which has its own direction of propagation as its characteristic directions. This means what characteristic means: the discontinuities of each angular flux field will separately propagate in its own direction. This observation introduces a natural way to solve the transport equation, Askew [17], that will be computationally more expensive but will produce more accurate calculations. For this, consider the mapping in the angular phase space Sn-1×Rn(4.1)Q(ω,x)=(ω′,xω,xω⊥)
with(4.2)ω′=ω,xω=ω⋅x,xω⊥=x-(ω⋅x)ω.
This is an isometry that rotates the axis x to a new direction xω coincident with the propagation direction. For each fixed direction ω∈Sn-1, it introduces a new coordinate system (xω,xω⊥) with the second coordinate in the projection Πω[Ω¯] of the physical domain Ω in a subspace of Rn made with directions ω⊥ and evidences the ordinary character of each angular flux in the equation(4.3)∂ϕ∂xω+σtϕ=Lσ[ϕ]=C[ϕ]+S,(ω,xω,xω⊥)∈Q[Sn-1]×Q[Ω],
subject to the same incident boundary condition(4.4)ϕ(ω,xω-,xω⊥-)=g-(ω,xω-,xω⊥-),(ω,x)∈Σ-.
Note that this incident boundary condition can be interpreted as an initial value for this system of ordinary differential equations coupled through the collision operator. For simplicity, we have used the same notation ϕ(ω,xω,xω⊥) for the rotated angular flux Q[ϕ](ω,xω,xω⊥).As pointed out in Section3., the inversion of the operator Lσ and its composition with special class of collision operator found in nuclear applications gives a compact operator(4.5)Lσ-1C=Lσ-1Cscattering+1λLσ-1Cfission:L2(Sn-1×Ω)⟶L2(Sn-1×Ω).
Of course this important property will pose the direct problem involving recovery of the angular flux from appropriate parameter data σ, σs, and σf. Some more functionals defined in terms of the rotated coordinates systems need to be introduced in order to make comprehensive the approximations that will be done to solve the problem. First of all, the optical distance from interior point x to a point x0∈Rn(4.6)α(x,x0)=∫xω0x⋅ωσt(Xω,Xω⊥)dXω,
where ω=(x-x0)/∥x-x0∥ is the direction that aligns x and x0, and the extension for points in the possibly bare medium outside the physical domain Ω¯ may be done by assuming σ(x)=0,x∈Rn∖Ω¯. When x0 is a boundary point, it is given by x±(ω,x)=x±τ±(ω,x)ω which is the point in the boundary where the line passing x with direction ω crosses the emergent or the incident boundary, respectively. It may be used to define the attenuation from the incident boundary(4.7)a-(ω,xω,xω⊥)=exp⁡(-∫xω-xωσt(Xω,Xω⊥)dXω).
We also define the right inverse of the trace operator, the extension operator as(4.8)E[g-](ω,x)=g-(ω,x-ωτ-(ω,x));(ω,x)∈Sn-1×Ω,
and the extension with attenuation from the incident boundary(4.9)E[g-](ω,x)a-(ω,x)=g-(ω,x-ωτ-(ω,x))a-(ω,x)=J[g-](ω,x);(ω,x)∈Sn-1×Ω.
By noting that the attenuation is the inverse of the integrating factor for the system (4.3), we may integrate from the incident boundary to the interior point x to obtain an integral equation which is an explicit expression of (3.17): (4.10)ϕ(ω,xω,xω⊥)=∫xω-xωa-(ω,xω,xω⊥)a-(ω,Xω,Xω⊥)⋅S(ω,Xω,Xω⊥)dXω+g-(x-ωτ-(ω,x))⋅a-(ω,xω,xω⊥)⋅∫xω-xωa-(ω,xω,xω⊥)a-(ω,Xω,Xω⊥)⋅∫Sn-1(σs(ω⋅ω′,Xω,Xω⊥)+1λσf(Xω,Xω⊥))ϕ(ω′,Xω,Xω⊥)dω′dXω,
for all (xω,xω⊥)∈[xω-,xω+]×Πω[Ω¯],ω∈Sn-1.The integration may also be done between two arbitrary points, inside or outside the domainΩ¯, x, and x0. In this case we must extend the cross-sections outside the domain as zero. By defining ω=(x-x0)/|x-x0| and a0(ω,xω,xω⊥)=exp(-α(x,x0)), we obtain (4.11)ϕ(ω,xω,xω⊥)=∫xω0xωa0(ω,xω,xω⊥)a0(ω,Xω,Xω⊥)⋅S(ω,Xω,Xω⊥)dXω+g-(x-ωτ-(ω,x))⋅a0(ω,xω,xω⊥)⋅∫xω0xωa0(ω,xω,xω⊥)a0(ω,Xω,Xω⊥)⋅∫Sn-1(σs(ω⋅ω′,Xω,Xω⊥)+1λσf(Xω,Xω⊥))ϕ(ω′,Xω,Xω⊥)dω′dXω,
for all (xω,xω⊥)∈[xω0,xω+]×Πω[Ω¯],ω∈Sn-1.Remark 4.1.
Problems given by the equations(1)
(2.1) and (2.2),(2)
(3.17),(3)
(4.3) and (4.4), and(4)
(4.10) and (4.11) are alternative formulations to the one-velocity transport equation model and express different aspects of this equation. Note that the initial value character of the problem will be accomplished if we construct the solution by using a method that marches along each one of the characteristics.Remark 4.2.
The consistence that will be appropriated to the application of the Lax equivalence theorem framework that has been formulated for time transient problems will be accomplished by the adoption of a rotated pixel basis in the next section. The compact operatorLσ-1C is subcritical or critical, which means the inverse of its eigenvalue is not less than one. So this gives stability to the class of problems been investigated.As a consequence of consistency, the Lax equivalence theorem applies, Valtchev and Roberty [18], and may be stated as follows.Theorem 4.3.
For a well-posed initial value problem, a consistent numerical scheme is convergent if and only if it is stable.


## 5. The Natural Partition

The integral form of the transport equation (4.11) has been written in a different coordinate system for each direction of propagation ω. In order to attempt the regularities requirements of each angular flux ϕ, we must introduce a grid with one direction as ω. This is impossible unless we choose a set of discrete ordinates as representative of all angular directions. This must be done a priori. Let S be the surface of the unity sphere in R3 and (5.1)ω=(ω1,ω2,ω3)=(sin(θ)cos⁡(ψ),sin(θ)sin(ψ),cos⁡(θ)),ψ∈[0,2π),θ∈[0,π),
where θ and ψ are the usual angles of the spherical coordinates system. To be more specific, let us consider the two-dimensional equation embedded in R3. So, we choose a set of (5.2)Sji={ωji=(sin(θi)cos⁡(ψj),sin(θi)sin(ψj),cos⁡(θi)),j=1:2J,i=1:2I}
of equally spaced angular angles ψj=((j-1)/2J)2π and at least symmetrical with the angle θ=π/2 of 2I collocations directions out of the plane generated by the vectors (cos(ψ),sin(ψ)), that is, out of the plane xy. We now may use these directions to introduces grids, one for each direction, or in the case of a two-dimensional problem, only in the plane directions Πj={cos(ψj),sin(ψj)}. Note that in the two dimensional problem, since the problem is embedded in the three-dimensional space, angular flux is defined for all directions in the S2 sphere but has spatial variations only in the plane xy. This introduces, for each direction out of the xy plane, the necessity of a sin(θi) factor correction. This gives the possibility of extending this formulation to investigate axis symmetric problems, but, from now on, we will fix in the two-dimensional problem. We have introduced 2J coordinate systems. Without loss of generality, let us consider each one with a rectangular equally spaced grid. These coordinates systems are supposed to be equally rotated with respect to the others. We have chose an even number for J, so, directions of propagations appear in groups of four, to complain with the four quadrants.We know that it is necessary to combine the different directions in order to calculate the collision operator for groups of pair the directions. Since the cross sections are properties common for all angular flux, the complete intercessions of the 2J rotated grid is the appropriate partition consistent with the numerical solution of the transport equation. To fix ideas, let us introduce the appropriate notation in R2. The direction ωj contains one strip(5.3)χnj(x,y)={1,ifnj≤-xsin⁡(ψj)+ycos⁡⁡(ψj)<1+nj,0,else,
which may be intercept with a strip in one direction in ωj⊥(5.4)χmj(x,y)={1,ifmj≤xcos⁡⁡(ψj)+ysin⁡(ψj)<1+mj,0,else,
to form a rotated pixel (5.5)Pjmjnj(x,y)={1,(x,y)∈χnj∩χmj,0,else.
The regular spacing hxj,hyj=1 has been adopted without loss of generality for simplicity. Note that in the direction ωj=cos(ψj)ex+sin(ψj)ey,the integrated (4.11) must be collocated for points (ωj,x0ωj,xωj⊥) and (ωj,xωj,xωj⊥) that will be referenced by the indexes (j,mj,nj) and (j,mj+1,nj), respectively. This is a consequence of the fact that ωj aligns the points and x0ωj⊥=xωj⊥. As the marching direction is ωj, the equation index mj will grow from its value corresponding to a point xωj in the incident boundary, where the boundary condition is imposed, to the corresponding point in the emergent boundary. This happens for each transversal index nj in each direction ωj. We now face the problem of computing the collision integral in (4.11) with a consistent discretization. Note that the measure dω′dXω(5.6)∫xω0xωa0(ω,xω,xω⊥)a0(ω,Xω,Xω⊥)⋅∫Sn-1(σs(ω⋅ω′,Xω,Xω⊥)+1λσf(Xω,Xω⊥)dω′dXω)
will corresponds to a volumetric Lebesgue measure when the two integrals interchanges, that is, when the Fubini theorem is applicable. The natural way to obtain this in the discretized approximation is by considering the partition presented by Roberty [19] that makes a complete intersection of the pixels Pjmjnj. This may be done with the complete interception of the strips defining the pixels to form the natural pixels(5.7)Enj1,…,nj2J(x,y)=∏j=12Jχnj(x,y).
This interception produces a domain partition which is consistent in the sense that if we increase the number of directions and reduce the size of the spacing hj=1 by using a more refined mesh, the error in the approximation to the exact stationary transport equation will approache zero accordingly. A particle or photon collected in one interception can only come from pixels aligned with the current position thought the discrete directions considered in the model. The approximation made can answer the question “from where the particle came?”, in a way consistent with the exact equation.Calculations of the vertex, area, centroid, and all geometric information are a classical computational geometry problem of domain decomposition, O’Rourke [20], which may present a great complexity in dimensions greater than two.As an optimization problem, the mesh generation is the determination of the convex polygonal sets limiting the maximum area that follows linear restrictions of the type:Ax≤b, where x is a vector in ℝ2, b is a vector in some ℝNb related with the strips nj, and A is a 2×Nb matrix related with the direction angle θj. The search process of the elements in the mesh executes the following steps:(1)
decomposition of the domain in search subdomains(2)
determination of the minimum and maximum strips for the search subdomain(3)
search for the true intersections,(4)
calculus of elements vertices, edges, area, and pertinence index for all directions(5)
assemble subdomains to form the natural pixel partition.


### 5.1. Decomposition of the Domain in Search Subdomains

The greater majority of the combinations(n1,…,n2J) are false intersections and can be a prior avoided by restricting the search for parameters on regions for which we can expect true intersections. We call these regions search subdomains, and they can be a pixel for some prescribed direction j, an intersection of some given strip with a sector, or any other type of polygonal subdomain. With them, we avoid a search with a large number of false intersections in the step (3). As an example in Figure 1 we show a detail of the mesh J64M128 for huge numbers of direction (64) and pixels (128×128), which is typical of image processing and reconstruction problems. Since the rotational symmetry allows us to calculate the interception only in a J/π sector, we may have less computation than the case of doing only two interceptions, as illustrated by Figure 2. Finally, in Figure 3, we present the complete mesh used to simulate our solution with test and benchmark problems.Figure 1
A detail in the J32M64 mesh.Figure 2
Sector mesh generated strip by strip in thej=16 and M=15 problem.Figure 3
The J6M24 mesh.


### 5.2. Determination of the Minimum and Maximum Strips for the Search Subdomain

This is a complementary procedure for minimizing the determination of false intersection.


### 5.3. Search for the True Intersections

In each domain of search, every strip combinations must be tested until the area of the domain is totally filled with natural pixels.


### 5.4. Calculus of Element Vertexes, Edges, Area, and Pertinence Index for all Directions

This calculation is done simultaneously with the search for true intersections. Note that not all intersections of lines limiting the strips are vertexes of the polygons and not all lines across a vertex are a polygonal edge.


### 5.5. Assemble Subdomains to Form the Natural Pixel Partition

In each subdomain, the pertinence indexes are calculated independently and are assembled to form the matrix of pertinence index for the entire mesh.For more information, please see Roberty [19] and in the references there. In Figure 4, the dependence of the numbers of polygonal  with the numbers of views and slices is presented. Note that the almost exponential growth observed for growing numbers of slice and views introduces memories problems in the computational implementation of the present methodology. In Figure 5, we present a comparison of mesh J6M3 with the classical Delaunay triangulation made for the same vertex number. Note that it is inconsistent with the natural partition, that is, it is not a triangular decomposition of it.Figure 4
Numbers of polygonal dependence with the numbers of views and slices.A comparison between the natural and the Delaunay mesh.
(a)(b)


## 5.1. Decomposition of the Domain in Search Subdomains

The greater majority of the combinations(n1,…,n2J) are false intersections and can be a prior avoided by restricting the search for parameters on regions for which we can expect true intersections. We call these regions search subdomains, and they can be a pixel for some prescribed direction j, an intersection of some given strip with a sector, or any other type of polygonal subdomain. With them, we avoid a search with a large number of false intersections in the step (3). As an example in Figure 1 we show a detail of the mesh J64M128 for huge numbers of direction (64) and pixels (128×128), which is typical of image processing and reconstruction problems. Since the rotational symmetry allows us to calculate the interception only in a J/π sector, we may have less computation than the case of doing only two interceptions, as illustrated by Figure 2. Finally, in Figure 3, we present the complete mesh used to simulate our solution with test and benchmark problems.Figure 1
A detail in the J32M64 mesh.Figure 2
Sector mesh generated strip by strip in thej=16 and M=15 problem.Figure 3
The J6M24 mesh.


## 5.2. Determination of the Minimum and Maximum Strips for the Search Subdomain

This is a complementary procedure for minimizing the determination of false intersection.


## 5.3. Search for the True Intersections

In each domain of search, every strip combinations must be tested until the area of the domain is totally filled with natural pixels.


## 5.4. Calculus of Element Vertexes, Edges, Area, and Pertinence Index for all Directions

This calculation is done simultaneously with the search for true intersections. Note that not all intersections of lines limiting the strips are vertexes of the polygons and not all lines across a vertex are a polygonal edge.


## 5.5. Assemble Subdomains to Form the Natural Pixel Partition

In each subdomain, the pertinence indexes are calculated independently and are assembled to form the matrix of pertinence index for the entire mesh.For more information, please see Roberty [19] and in the references there. In Figure 4, the dependence of the numbers of polygonal  with the numbers of views and slices is presented. Note that the almost exponential growth observed for growing numbers of slice and views introduces memories problems in the computational implementation of the present methodology. In Figure 5, we present a comparison of mesh J6M3 with the classical Delaunay triangulation made for the same vertex number. Note that it is inconsistent with the natural partition, that is, it is not a triangular decomposition of it.Figure 4
Numbers of polygonal dependence with the numbers of views and slices.A comparison between the natural and the Delaunay mesh.
(a)(b)


## 6. The Discrete Equation for the Two-Dimensional Problem

As we mention before, we may consider the integral form of the transport (4.11) for direction ωj, j=1,2J, in the pixel Pjmjnj, nj=1,2M and mj=1,2M. Since we are focusing on the two dimensional problem, it necessary introduces the sin(θi) correction for each out of the xy plane direction. Note that the xy plane corresponds to the case sin(θi)=sin(π/2)=1. There is by hypotheses no variation of the angular flux in the z=xy⊥ direction, but it still remains dependent of azimuthal angle. Let us rebuild the notation to make this behavior explicit. Since ω3·∇xϕ=∂ϕ/∂z=0 and all other functional variations in this direction are also zero, we rewrite ωji=sin(ψi)·ωj and obtain by straightforward calculations, the optical distance for points with projection xy along the ω3 direction as(6.1)α(xy,xy0,ψi)=∫xy0⋅ωjxy⋅ωjσ(xyω′,xyω⊥(xy,xy0))dxyωj′sin(ψi).
Since all other parameters, as the cross sections, and sources are linear measures of the change in the particle population in the direction ωji, it must be corrected accordingly.As a consequence of the divergence theorem, the angular flux in directionωji has average value on all lateral segments (surfaces in three dimension), that is, the segments that are transversal to the strip χmj(x,y), parallel to ωj for all mj, are equal to zero. We complement this by introducing the following hypotheses:  (1)
that angular flux is well represented by its average value in each transversal section to the stripχnj(x,y),(2)
the cross sections are known by its average values in the natural pixel(6.2)χen1,…,n2J=∏j=12Jχnj,(3)
the average cross sections inside each rotated rectangular pixelPjmjnj are the area (volume in three dimensions) average weight of each natural pixel cross section value,(4)
to accomplish with high cross section values, we integrate the exponentials inside the rotated pixelPjmjnj,(5)
the mean value of angular flux inside the pixelPjmjnj is used in the calculation of scattering and fission reaction rate.


## 7. The Method of Successive Approximations

With respect to the intensity of the extinction cross sectionσt, we observe two important class of problems:  (1)
there exists a transillumination of the mediumΩ, which means that the extinction cross sections not sufficient to produce an optical path that shields the radiation flux. In this class of problems, a flux radiation in the incident boundary propagates through the medium and can be collected in the emergent boundary. Its modification will be used to interpret the optical properties of the medium inside Ω,(2)
there is no transillumination, which means that the extinction cross sections are sufficient to produce an optical path that completely shields the radiation flux. In this case, the shielded radiation is must be replaced by an internal source or the fission in the case of reactors.Since the extinction coefficients may present high values, we considered exponential flux attenuation inside the rotated pixel and calculated the reaction rates using the average flux expressed in terms of the incident and emergent flux in the rotated pixel. Equation (4.11) becomes (7.1)ϕmj+1,nj,j=ϕmj,nj,jE1,mj,nj,j+qmj,nj,jE2,mj,nj,jformj,nj=1,2M;j=1,2J.
Here, the source with contribution of the scattering and which forms fission inside the rotated pixel is (7.2)qmj,nj,j=∑k=12J∑e∈Pjmjnjae⋅[σs,e+1λσf,e]⋅[E3,m(e,k),n(e,k),k)ϕm(e,k),n(e,k),k+E4,m(e,k)+1,n(e,k),kϕm(e,k)+1,n(e,k),k];k=1,2J,E1,mj,nj,j=exp⁡⁡(-σmj,nj,j),E2,mj,nj,j=1-exp⁡⁡(-σmj,nj,j)σmj,nj,j,E3,mj,nj,j=1-(1+σmj,nj,j)exp⁡⁡(-σmj,nj,j)σmj,nj,j(1-exp⁡⁡(-σmj,nj,j)),E4,mj,nj,j,i=-1+σmj,nj,j+exp⁡⁡(-σmj,nj,j)σmj,nj,j(1-exp⁡⁡(-σmj,nj,j)),σmj,nj,j=∑e∈Pjmjnjaeσe.Remark 7.1.
Note that the correction for extending this methodology for axis symmetric problems is(7.3)σe,i=σesin(ψi)
as mentioned in the introduction of the natural partition in Section 5.The localization of the natural elementen1,…,n2J inside the rotated pixel Pmj,nj,j is made by the incident matrix. As we mentioned before, given one pair (e,j), one may calculate the rotate integer coordinates (m(e,j),n(e,j),j) by extracting the integer part of the natural element centroid in each rotate coordinate system. This is a procedure similar to the retro projection that is done in image reconstruction algorithms. The function relating the two integer coordinate systems is a discrete version of transformation (4.2). It is one-to-one between (e,j)→(m,n,j). It is not convenient to make it onto since the implementation of an procedure for the localization of the boundary of Ω will be computationally expensive. We will start the march in the algorithm in the position labeled with mj=1, which may correspond to a pixel outside the domain, and march in one of the directions j until the pixel labeled 2M, which also may be outside the domain, is arrived. This introduces no problem if we take care in extending outside the domain the cross sections value by zero, the flux by its boundary value and calculates the exponentials accordingly.As mentioned before, these equations are consistent and introduce no spatial error when the properties are constant; in that case only, the angular discretization error will be presented. The simulation problems that we will solve in Section8.3 have material parameters very close to this situation.As we pointed out before, the marching scheme equation (7.1) is a discretization of the first-order form of the first order system representing the stationary transport equation that has initial value in the incident boundary. The fact that the total cross-section is bigger than the scattering assures contractility of the compact operator related to the scattering kernel, and from this, we will obtain the stability and as a consequence a convergent scheme.


### 7.1. Algorithms Based on Successive Approximation Method

Properties of collision operator equation (3.11) and of the boundary conditions allow us to distinguish the basic class of problems for the stationary transport equation:(7.4)ϕ=(Ea-)[g-]+Lσ-1Cscattering[ϕ]+1λLσ-1Cfission[ϕ]+Lσ-1S;(ω,x)∈Sn-1×Ω.(1)
The collision operator and the external source localized insideΩ are negligible. In this case, the equation is integrable, and we have the direct problem for the X-ray transform.(2)
The fission part of the collision operatorLσ-1Cfission has the first positive eigenvalue greater than one, and the domain Ω is subcritical external source problem. The incident flux may be zero or not. There are two important cases of nonfission domain, that is, when there is no fissile material inside the domain and σf is properly zero. The external source inside the domain problem, with q≠0 and g-=0 and transillumination problem with S=0 and g-≠0. The transillumination problem has a greater interest for tomography and nondestructive testing.(3)
The fission part of the collision operatorLσ-1Cfission is defined with support inside Ω, and we want to find his first positive eigenvalue and its positive eigenfunction. Once we find the eigenvalue, we give to κ the same value and obtain a new fission operator with first eigenvalue equal one, which means that we have found the critical material composition for the domain Ω.


### 7.2. The External Source Inside the Domain Problem

The algorithm for the angular flux calculation in the external source inside the domain problem follows the following steps.(1)
Establish an error tolerance criteria.(2)
Input material properties and source:σt,σs,S.(3)
Input mesh parameters:M, J, I.(4)
Generate mesh: elements area, centroid, and vertex and incidence integers matrices for natural elements inside the rotated pixels.(5)
Adjust geometric scale between the regular mesh and the spatial distribution of reacting material.(6)
Calculate the attenuation, the exponentials, and external source for each rotated pixel.(7)
Initialize the flux and reaction rates.(8)
For iteration equal 1 to maximum number of iterations do.(9)
Calculate scattering reaction rate based with the previous iteration flux estimate.(10)
Calculate the new iteration flux with (7.1) with zero fission part and zero incident flux.(11)
Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the iteration before the  maximum number of iterations.(12)
End the iteration do.


### 7.3. First Positive Eigenvalue and Eigen Angular Flux Determination Problems

The algorithm for first positive eigenvalue and eigen angular flux in the criticality determination problem follows the following steps.(1)
Establish an error tolerance criteria.(2)
Input material properties and source:σt,σs,σf.(3)
Input mesh parameters:M, J, I.(4)
Generate mesh: elements area, centroid and vertex and incidence integers matrices for natural elements inside the rotated pixels.(5)
Adjust geometric scale between the regular mesh and the spatial distribution of reacting material.(6)
Calculate the attenuation and the exponentials for each rotated pixel.(7)
Set the first eigenvalue to one:λ0=1.(8)
Initializes with the fluxϕ0 normalized by the fission part of the collision operator.(9)
For outer iteration equal 1 to maximum number of outer iterations do.(10)
Since the reactions rate is calculated by collecting the contribution of each natural pixel to its averaged value inside the rotated pixels, take care to zero the rates before start.(11)
Calculate scattering outer reaction rates and the fission outer reaction rate based on the previous outer iteration flux.(12)
Calculate the new external outer iteration flux with (7.1) with zero external source and zero incident flux.(13)
Initialize the internal loop fluxϕ0 with the outer loop flux ϕiterouter and normalizes it with the fission operator.(14)
For inner iterations equal 1 to maximum number of inner iterations do.(15)
Take care to zero the scattering reaction rate, but not the reaction fission rate, since it will not change inside the internal loop.(16)
Calculate scattering reaction rate based on the previous inner iteration calculated flux.(17)
Calculate the new internal iteration flux with (7.1) with zero external source and zero incident flux and external loop value of the fission reaction rate.(18)
Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the inner iteration before the maximum number of inner iterations. Note that we may alternatively use the fission reaction rate error as criteria.(19)
End inner iteration do.(20)
Take the convergent inner loop value as the outer loop flux.(21)
Calculates the new eigenvalue value as a ratio between the fission rate based on the actual outer iteration flux and the fission rate based on the previous outer iteration flux.(22)
Verify if the difference between the old eigenvalue and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the outer iteration before the maximum number of outer iterations.(23)
End outer iteration do.Remark 7.2.
The first and the inner part of the second algorithm are based on the contracting behavior of the scattering collision operator, that is, we have a compact operator with spectral radius less than one and the rate of convergence will be proportional to this ratio. Perhaps other rates of convergence may be found with other algorithms. The outer loop of the second algorithm is also based on compactness. The problem is to determine the first eigenvalue in the pencil problem(7.5)(I-Lσ-1Cscattering)[ϕ]=λLσ-1Cfission[ϕ]
in L2(sn-1×Ω), that is, a perturbation of the identity by a compact operator in a compact weigh operator pencil problem. Its convergence with a decreasing sequence of eigenvalues may be proved by using the fact that the multiplicative region inside Ω has compact support and the composite operator is positive and compact. The fact that the first eigen angular flux is positive almost everywhere in Sn-1×Ω is a classical result based on Jentzsch’s theorem for compact operators  which leaves nonnegative those functions that are non negative. Once again the consistence has a positive consequence, since numerical artifacts have not been introduced in the approximation, and good properties of the continuous equation are preserved.


## 7.1. Algorithms Based on Successive Approximation Method

Properties of collision operator equation (3.11) and of the boundary conditions allow us to distinguish the basic class of problems for the stationary transport equation:(7.4)ϕ=(Ea-)[g-]+Lσ-1Cscattering[ϕ]+1λLσ-1Cfission[ϕ]+Lσ-1S;(ω,x)∈Sn-1×Ω.(1)
The collision operator and the external source localized insideΩ are negligible. In this case, the equation is integrable, and we have the direct problem for the X-ray transform.(2)
The fission part of the collision operatorLσ-1Cfission has the first positive eigenvalue greater than one, and the domain Ω is subcritical external source problem. The incident flux may be zero or not. There are two important cases of nonfission domain, that is, when there is no fissile material inside the domain and σf is properly zero. The external source inside the domain problem, with q≠0 and g-=0 and transillumination problem with S=0 and g-≠0. The transillumination problem has a greater interest for tomography and nondestructive testing.(3)
The fission part of the collision operatorLσ-1Cfission is defined with support inside Ω, and we want to find his first positive eigenvalue and its positive eigenfunction. Once we find the eigenvalue, we give to κ the same value and obtain a new fission operator with first eigenvalue equal one, which means that we have found the critical material composition for the domain Ω.


## 7.2. The External Source Inside the Domain Problem

The algorithm for the angular flux calculation in the external source inside the domain problem follows the following steps.(1)
Establish an error tolerance criteria.(2)
Input material properties and source:σt,σs,S.(3)
Input mesh parameters:M, J, I.(4)
Generate mesh: elements area, centroid, and vertex and incidence integers matrices for natural elements inside the rotated pixels.(5)
Adjust geometric scale between the regular mesh and the spatial distribution of reacting material.(6)
Calculate the attenuation, the exponentials, and external source for each rotated pixel.(7)
Initialize the flux and reaction rates.(8)
For iteration equal 1 to maximum number of iterations do.(9)
Calculate scattering reaction rate based with the previous iteration flux estimate.(10)
Calculate the new iteration flux with (7.1) with zero fission part and zero incident flux.(11)
Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the iteration before the  maximum number of iterations.(12)
End the iteration do.


## 7.3. First Positive Eigenvalue and Eigen Angular Flux Determination Problems

The algorithm for first positive eigenvalue and eigen angular flux in the criticality determination problem follows the following steps.(1)
Establish an error tolerance criteria.(2)
Input material properties and source:σt,σs,σf.(3)
Input mesh parameters:M, J, I.(4)
Generate mesh: elements area, centroid and vertex and incidence integers matrices for natural elements inside the rotated pixels.(5)
Adjust geometric scale between the regular mesh and the spatial distribution of reacting material.(6)
Calculate the attenuation and the exponentials for each rotated pixel.(7)
Set the first eigenvalue to one:λ0=1.(8)
Initializes with the fluxϕ0 normalized by the fission part of the collision operator.(9)
For outer iteration equal 1 to maximum number of outer iterations do.(10)
Since the reactions rate is calculated by collecting the contribution of each natural pixel to its averaged value inside the rotated pixels, take care to zero the rates before start.(11)
Calculate scattering outer reaction rates and the fission outer reaction rate based on the previous outer iteration flux.(12)
Calculate the new external outer iteration flux with (7.1) with zero external source and zero incident flux.(13)
Initialize the internal loop fluxϕ0 with the outer loop flux ϕiterouter and normalizes it with the fission operator.(14)
For inner iterations equal 1 to maximum number of inner iterations do.(15)
Take care to zero the scattering reaction rate, but not the reaction fission rate, since it will not change inside the internal loop.(16)
Calculate scattering reaction rate based on the previous inner iteration calculated flux.(17)
Calculate the new internal iteration flux with (7.1) with zero external source and zero incident flux and external loop value of the fission reaction rate.(18)
Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the inner iteration before the maximum number of inner iterations. Note that we may alternatively use the fission reaction rate error as criteria.(19)
End inner iteration do.(20)
Take the convergent inner loop value as the outer loop flux.(21)
Calculates the new eigenvalue value as a ratio between the fission rate based on the actual outer iteration flux and the fission rate based on the previous outer iteration flux.(22)
Verify if the difference between the old eigenvalue and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the outer iteration before the maximum number of outer iterations.(23)
End outer iteration do.Remark 7.2.
The first and the inner part of the second algorithm are based on the contracting behavior of the scattering collision operator, that is, we have a compact operator with spectral radius less than one and the rate of convergence will be proportional to this ratio. Perhaps other rates of convergence may be found with other algorithms. The outer loop of the second algorithm is also based on compactness. The problem is to determine the first eigenvalue in the pencil problem(7.5)(I-Lσ-1Cscattering)[ϕ]=λLσ-1Cfission[ϕ]
in L2(sn-1×Ω), that is, a perturbation of the identity by a compact operator in a compact weigh operator pencil problem. Its convergence with a decreasing sequence of eigenvalues may be proved by using the fact that the multiplicative region inside Ω has compact support and the composite operator is positive and compact. The fact that the first eigen angular flux is positive almost everywhere in Sn-1×Ω is a classical result based on Jentzsch’s theorem for compact operators  which leaves nonnegative those functions that are non negative. Once again the consistence has a positive consequence, since numerical artifacts have not been introduced in the approximation, and good properties of the continuous equation are preserved.


## 8. Numerical Simulations

### 8.1. The Transillumination Test Problem

In order to start the evaluation of the method proposed, we first study a subcritical transillumination test problem with no fission and no external sources. The flux for test is(8.1)ϕtest=exp⁡⁡(-0.0125⋅y2+0.05⋅(x-M))forj=24andzerosforotherviews.
The algorithm is used to solve problem (2.1) with source (8.2)S(ω,x)=ω⋅∇ϕtest(ω,x)+σt(x)ϕtest(ω,x)-∫Sn-1f(x,ω′⋅ω)ϕtest(ω′,x)dω′,(ω,x)∈Sn-1×Ω,
and incident boundary condition
(8.3)ϕ(ω,σ)=ϕtest(ω,σ)for(ω,σ)∈Σ-.
A high-contrast inclusion problem as shown in Figure 6 in a 24 views with 30×30 rectangular pixels/view model is solved with the algorithm proposed here which is based on (7.4) accordingly. Material properties are given in Table 1. As can be seen from Figure 7 comparing the given and the reconstructed solution showa that they are almost equal. From Figure 8 showing the absolute error for two views 9 and 24, the error is of order 10-8, bigger in the region with the high-contrast inclusion, but with the same magnitude. The absolute error variation with J varying between 6 to 16 and M from 10 to 24 is showed in Figure 9.Table 1
Absorption and scattering coefficient in the transillumination model.
PropertySymbolValueSection radius (cm)M15Circular inclusion radius (cm)M1.5Circular inclusion center (cm)(x1,x2)(3.0,1.5)Absorption cross-section blue (nonperturbed) (cm−1)σa0.0111Absorption cross-section green (cm−1)σa0.0555Scattering cross-section blue (non perturbed) (cm−1)σs0.0014Scattering cross-section green (cm−1)σs0.0167Figure 6
High-contrast nonhomogeneous model with 24 views and30×30 pixels view.Given and reconstructed view 24 angular test flux in the model.
(a)(b)Absolute error in view 24 and 9 angular test flux for the model.
(a)(b)Figure 9
Maximum error variation forJ and M.


### 8.2. A Transillumination Problem

The same transillumination problem with no fission and no external sources is now investigated. The incident radiation boundary condition isg(n2J,2J)=1, for n2J=1,…,2M and g(nj,j)=0 for all others views. Figure 10 shows selected particles (radiation) flux for j=9 and j=24. We can see the depleted effect in the main direction of the incident flux j=24 due to absorption and outscattering and the buildup effect in the other direction j=9 due to the scattered flux. Also, it is possible to note the differences in magnitudes in the flux for the incident direction, j=24, with magnitude order of 1 and the other directions where there is no incident radiation, with magnitude order of 0,02, one up to 50-times the other. Another important information that can be determinated by the solver presented here is emergent radiation (flux) map shown in Figure 11. This set of figures graphically defines the flux in the emergent boundary Σ+(8.4)ϕ=g+on⁡Σ+,
and the Cauchy data for the problem (8.5)Cqa,f={(g,g+)on⁡Σ-×Σ+},
characterizes the graph for the albedo operator (the incident to emergent flux mapping) (8.6)Aqa,f:L2(Σ-,|ω⋅ν(σ)|dσdω)⟼L2(Σ+,|ω⋅ν(σ)|dσdω).Selected radiation (particles) flux for illumination on view 24 in a high-contrast nonhomogeneous problem.
(a)(b)Figure 11
Emergent flux for illumination on view 24 in the high-contrast model problem.In the inverse problem, we ask if it is possible to determine the coefficients functionsσt and f from the a priori knowledge of the albedo operator, that is, data such as that found in Figure 11 for different incident radiation g-. Finally, the convergent sequence {ϕiter;iter=1,…,31,…} is tested as a Cauchy sequence in the complete Hilbert space H. Figure 12 shows that before the eleventh iteration, we achieved convergence of order 0,008. If we claim for 9,8e-008, it is achieved in the model problem in the 31 iteration.Figure 12
Cauchy error for 24 views and30×30 pixels view.


### 8.3. Reactor with the EIR-2A and EIR-2B Data

The EIR-2A and EIR-2B are two problems found in Mordant [21] with data about values on cross sections, source and geometric domain definition presented in Table 2 and represented in Figure 13. The total and scattering cross-section is the same for the two problems. The source values is for the EIR-2A problem, and in this case the fission cross-section is zero. The fission cross-section is for the EIR-2B problem, and in this case the source term is zero. Results obtained are in accordance with the theory and are showed in figures. In all cases, the number of axial directions is I=1. Figure 14 shows the angular flux for direction j=1 in the source problem. Figure 15 shows the maximum Cauchy error flux for the angular flux function, where a high convergence rate is observed. The Figure 16 shows the maximum Cauchy error in the first eigenvalue function, where a not-so-high convergence rate is observed. Finally, Figure 17 shows the first eigen angular flux for direction j=1. These results showing a monotonically decreasing sequence of eigenvalues with the smallest accumulation point at 1 and a nonnegative angular flux are in complete agreement with theoretical results pointed out in Remark 7.2. We must finally note that the continuity in the direction perpendicular to ω1 has been artificially introduced by the graphic algorithm used to present the results. In reality, the calculations does not impose this restriction on the function.Table 2
Cross sections and sources values for EIR problems.
Spatial regionTotalScatteringFissionSource|x|≤48,|x|≥30,|y|≤43.90.89.00.00|x|≤30,|y|≤43,|y|≥25.90.89.00.00x<30, x≥0, y≥0, y<25.70.66.431.00x>-30, x<0, y≥0, y<25.65.50.00.00x>-30, x≤0, y>-25, y<0.60.53.791.00x>0, x<30, y>-25, y<0.48.20.00.00EIR-2A and EIR-2B benchmark problems with data as in Table2.
(a)(b)Figure 14
Angular flux in J6M24 EIR 2A.Figure 15
Cauchy error in flux convergence for J6M24 EIR 2A.Figure 16
First eigenvalue in J6M24 EIR-2B.Figure 17
First eigenfunction in J6M24 EIR-2B.


## 8.1. The Transillumination Test Problem

In order to start the evaluation of the method proposed, we first study a subcritical transillumination test problem with no fission and no external sources. The flux for test is(8.1)ϕtest=exp⁡⁡(-0.0125⋅y2+0.05⋅(x-M))forj=24andzerosforotherviews.
The algorithm is used to solve problem (2.1) with source (8.2)S(ω,x)=ω⋅∇ϕtest(ω,x)+σt(x)ϕtest(ω,x)-∫Sn-1f(x,ω′⋅ω)ϕtest(ω′,x)dω′,(ω,x)∈Sn-1×Ω,
and incident boundary condition
(8.3)ϕ(ω,σ)=ϕtest(ω,σ)for(ω,σ)∈Σ-.
A high-contrast inclusion problem as shown in Figure 6 in a 24 views with 30×30 rectangular pixels/view model is solved with the algorithm proposed here which is based on (7.4) accordingly. Material properties are given in Table 1. As can be seen from Figure 7 comparing the given and the reconstructed solution showa that they are almost equal. From Figure 8 showing the absolute error for two views 9 and 24, the error is of order 10-8, bigger in the region with the high-contrast inclusion, but with the same magnitude. The absolute error variation with J varying between 6 to 16 and M from 10 to 24 is showed in Figure 9.Table 1
Absorption and scattering coefficient in the transillumination model.
PropertySymbolValueSection radius (cm)M15Circular inclusion radius (cm)M1.5Circular inclusion center (cm)(x1,x2)(3.0,1.5)Absorption cross-section blue (nonperturbed) (cm−1)σa0.0111Absorption cross-section green (cm−1)σa0.0555Scattering cross-section blue (non perturbed) (cm−1)σs0.0014Scattering cross-section green (cm−1)σs0.0167Figure 6
High-contrast nonhomogeneous model with 24 views and30×30 pixels view.Given and reconstructed view 24 angular test flux in the model.
(a)(b)Absolute error in view 24 and 9 angular test flux for the model.
(a)(b)Figure 9
Maximum error variation forJ and M.


## 8.2. A Transillumination Problem

The same transillumination problem with no fission and no external sources is now investigated. The incident radiation boundary condition isg(n2J,2J)=1, for n2J=1,…,2M and g(nj,j)=0 for all others views. Figure 10 shows selected particles (radiation) flux for j=9 and j=24. We can see the depleted effect in the main direction of the incident flux j=24 due to absorption and outscattering and the buildup effect in the other direction j=9 due to the scattered flux. Also, it is possible to note the differences in magnitudes in the flux for the incident direction, j=24, with magnitude order of 1 and the other directions where there is no incident radiation, with magnitude order of 0,02, one up to 50-times the other. Another important information that can be determinated by the solver presented here is emergent radiation (flux) map shown in Figure 11. This set of figures graphically defines the flux in the emergent boundary Σ+(8.4)ϕ=g+on⁡Σ+,
and the Cauchy data for the problem (8.5)Cqa,f={(g,g+)on⁡Σ-×Σ+},
characterizes the graph for the albedo operator (the incident to emergent flux mapping) (8.6)Aqa,f:L2(Σ-,|ω⋅ν(σ)|dσdω)⟼L2(Σ+,|ω⋅ν(σ)|dσdω).Selected radiation (particles) flux for illumination on view 24 in a high-contrast nonhomogeneous problem.
(a)(b)Figure 11
Emergent flux for illumination on view 24 in the high-contrast model problem.In the inverse problem, we ask if it is possible to determine the coefficients functionsσt and f from the a priori knowledge of the albedo operator, that is, data such as that found in Figure 11 for different incident radiation g-. Finally, the convergent sequence {ϕiter;iter=1,…,31,…} is tested as a Cauchy sequence in the complete Hilbert space H. Figure 12 shows that before the eleventh iteration, we achieved convergence of order 0,008. If we claim for 9,8e-008, it is achieved in the model problem in the 31 iteration.Figure 12
Cauchy error for 24 views and30×30 pixels view.


## 8.3. Reactor with the EIR-2A and EIR-2B Data

The EIR-2A and EIR-2B are two problems found in Mordant [21] with data about values on cross sections, source and geometric domain definition presented in Table 2 and represented in Figure 13. The total and scattering cross-section is the same for the two problems. The source values is for the EIR-2A problem, and in this case the fission cross-section is zero. The fission cross-section is for the EIR-2B problem, and in this case the source term is zero. Results obtained are in accordance with the theory and are showed in figures. In all cases, the number of axial directions is I=1. Figure 14 shows the angular flux for direction j=1 in the source problem. Figure 15 shows the maximum Cauchy error flux for the angular flux function, where a high convergence rate is observed. The Figure 16 shows the maximum Cauchy error in the first eigenvalue function, where a not-so-high convergence rate is observed. Finally, Figure 17 shows the first eigen angular flux for direction j=1. These results showing a monotonically decreasing sequence of eigenvalues with the smallest accumulation point at 1 and a nonnegative angular flux are in complete agreement with theoretical results pointed out in Remark 7.2. We must finally note that the continuity in the direction perpendicular to ω1 has been artificially introduced by the graphic algorithm used to present the results. In reality, the calculations does not impose this restriction on the function.Table 2
Cross sections and sources values for EIR problems.
Spatial regionTotalScatteringFissionSource|x|≤48,|x|≥30,|y|≤43.90.89.00.00|x|≤30,|y|≤43,|y|≥25.90.89.00.00x<30, x≥0, y≥0, y<25.70.66.431.00x>-30, x<0, y≥0, y<25.65.50.00.00x>-30, x≤0, y>-25, y<0.60.53.791.00x>0, x<30, y>-25, y<0.48.20.00.00EIR-2A and EIR-2B benchmark problems with data as in Table2.
(a)(b)Figure 14
Angular flux in J6M24 EIR 2A.Figure 15
Cauchy error in flux convergence for J6M24 EIR 2A.Figure 16
First eigenvalue in J6M24 EIR-2B.Figure 17
First eigenfunction in J6M24 EIR-2B.


## 9. Conclusions

The present methodology utilizes the natural basis concept, originally developed for image reconstruction, in the numerical solution of very classical and technological problems found in reactor physics calculations. The algorithms presented here may be considered as giving almost exact solutions in case where material properties and sources are constants in a small number of distinct regions inside the domain, as observed in the simulations EIR-2A and EIR-2B implemented. It does not array system matrices, and as a consequence few bigger problems signify more processing time. The algorithms are deduced by introducing, as less as possible, artifacts related with preprocessing the transport equation by doing inconsistent averages and simplifications hypotheses. It respects the basic characteristic of the transport equation, that is, no continuity requirements in the directions transverse to the propagation direction. It is also based on a scheme consistent with the exact equation. The algorithm implementation also shows compatibility with the present status of computer development.



---

*Source: 126309-2011-04-06.xml*