- Email: [email protected]
Using the Schwinger inverse method for solutions of inverse transport problems in two-dimensional cylindrical geometries
Annals of Nuclear Energy 36 (2009) 966–973
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Using the Schwinger inverse method for solutions of inverse transport problems in two-dimensional cylindrical geometries Keith C. Bledsoe a,*, Jeffrey A. Favorite b, Tunc Aldemir a a b
Nuclear Engineering Program, The Ohio State University, Columbus, OH 43210, USA Applied Physics Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
a r t i c l e
i n f o
Article history: Received 27 October 2008 Received in revised form 26 February 2009 Accepted 27 February 2009 Available online 29 April 2009
a b s t r a c t The Schwinger method for solving inverse transport problems is applied to the problems of interface location identification, shield material identification, source isotope weight fraction identification, and material mass density identification (separately) in multilayered two-dimensional cylindrical gamma radiation source/shield systems. The method is iterative and estimates unknown interface locations, source isotope weight fractions, and material densities directly, while the unknown shield material is identified by estimating its total macroscopic gamma-ray cross sections. The energies of discrete gamma-ray lines emitted by the source are assumed to be known, while the unscattered flux of the lines is assumed to be measured at points external to the system. In numerical test cases, the Schwinger method correctly identifies the unknowns when the same deterministic ray-tracing code is used for both the parameter estimation process and simulation of the measured data. With realistic simulation of the measured data using a Monte Carlo code, the method produces more ambiguous results for interface location, shield material identification, and material density identification. The method works well for source weight fraction identification with measured data simulated by Monte Carlo. In addition to the application to more realistic (two-dimensional) problems, this paper extends previous work on the Schwinger inverse method by using surface formulas for unknown interface locations, automatic correction attempts for violated constraints, and ray-tracing instead of discrete-ordinates for transport calculations. Published by Elsevier Ltd.
1. Introduction The Schwinger inverse method was recently presented (Favorite, 2004, 2005; Favorite and Bledsoe, 2006; Bledsoe et al., 2007) for obtaining solutions to the problems of determining, in a gamma radiation source/shield system: (a) a set of unknown internal interface locations, (b) the isotopic content of an unknown source, (c) interface locations and source isotopic content simultaneously, (d) the composition of an unknown layer in the shield, and (e) unknown material mass densities. Favorite (2004) presented the mathematical background for problems (a–d) and documented numerical solutions for problems (a) and (b) in one-dimensional spherical geometries. Numerical results for the unknown shield layer problem were presented in Favorite and Bledsoe (2006), also for one-dimensional spherical geometries. In Favorite (2005) numerical results for the mass density problem were presented for one-dimensional spherical problems. The method was successful for identification of the unknowns in each of those papers. In Bledsoe et al. (2007) the method was used successfully for material identification in cylindrical systems, but the geometries were small
* Corresponding author. Address: X-1-TA, MS P365, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. E-mail address: [email protected] (K.C. Bledsoe). 0306-4549/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.anucene.2009.02.014
and the required leakage measurements would be extremely difficult to make. In this paper, the Schwinger inverse method is implemented and tested in realistic two-dimensional cylindrical geometries for the problems of identifying (separately) unknown interface locations, a shield material, source isotope weight fractions, and material mass densities. In addition to the application to multidimensional problems, this paper extends previous work on the Schwinger inverse method by using surface formulas for unknown interface locations, automatic correction attempts for violated constraints, and ray-tracing instead of discrete-ordinates for transport calculations. 2. Methodology of the Schwinger inverse method 2.1. General approach The Schwinger inverse method is derived from a perturbationtheory approach to the inverse transport problem (Favorite, 2004). Instead of calculating the effect of a system perturbation on a quantity of interest, the quantity of interest is assumed to be given (from a measurement) and the Schwinger functional is manipulated to produce an equation describing the system perturbation. This equation is applied iteratively to determine the unknown parameters in the system. The quantity of interest is the
967
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
flux of a set of discrete gamma-ray lines measured at a point (or the sum of fluxes measured at points) external to the geometry, which implies that the scattering term in the transport equation is assumed to be negligible. The mathematics of the Schwinger inverse method are presented in detail in Favorite (2004) and Favorite and Bledsoe (2006). Here a brief review of those results will be given. Begin by considering a radiation shield surrounding a gamma-ray source. Both the source and the shield may be multilayered but the layers are assumed to be homogeneous. The source emits gamma rays at discrete energies or lines, which can be resolved quite well using a high-purity germanium detector. It is assumed that there is no scattering into the energy lines. In Favorite (2004), the actual (unknown) configuration was considered the perturbed system, while an estimate of the configuration was the unperturbed system. The Schwinger functional was used to derive an expression for the quantity of interest in the system (in our case, unscattered scalar fluxes measured at points external to the system) that is second-order accurate with respect to the difference between the unperturbed and perturbed scalar fluxes. By manipulating this expression, equations were derived to solve for the unknown system parameter. ^ Þ be the forward angular flux For a discrete energy line, let wðr; X ^ Þ the adjoint angu^ , w ðr; X of gamma rays at position r and angle X ^ , Rt(r) the total lar flux of gamma rays at position r and angle X macroscopic cross section at position r, q(r) the source, at position r, of gamma rays (c/cm3 s), M the sum of the calculated flux of photons over all detector points (c/cm2 s), and M0 the sum of the measured flux of photons over all detector points (c/cm2 s). 2.2. Unknown interface locations If the system contains unknown radii and/or unknown heights, the equation for the unknowns Drn is (Favorite, 2004)
" # N N X 1 M X M M0 0 0 DRt;n In ðr n Þ Dq I ðr Þ ¼ ; R R ^ w q M0 n¼1 n s;n n M0 dV dX n¼1
ð1Þ
where N is the total number of unknown interfaces in the system and In ðr 0n Þ and Is;n ðr0n Þ are defined by
In ðr0n Þ
Z
r 0n
dV
Z
^ w w dX
ð2Þ
rn
DRt;n Rt;k Rt;kþ1 and
Dqn qk qkþ1 ;
Z
r0n
dV
Z
^ w : dX
ð3Þ
2.3. Unknown shield material For the unknown shield material problem, the following equation can be derived for the perturbed cross section for a discrete energy line (Favorite, 2004; Favorite and Bledsoe, 2006):
R0t ¼ R
In ðr0n Þ Dr n
Z
dS
Z
^ w w dX
ð4Þ
^ w ; dX
ð5Þ
and
Is;n ðr 0n Þ Dr n
Z
dS
Z
R where dS represents an integral over the unperturbed radial or axial surface r = rn. These approximations eliminate the inaccuracies introduced by truncation of a Taylor series. The terms DRt;n and Dqn are the differences in cross section and source density across interface n; they are defined by
R
R
^ w q M M 0 dX þ Rt : R ^ w w M0 dV dX DV dV
ð8Þ
In Eq. (8) it is assumed that the unknown cross section is constant over the unknown region. The integral in the numerator in Eq. (8) is over the entire problem, but the integral in the denominator is over the unknown material region DV only. After Eq. (8) converges, the calculated cross sections are used to identify the unknown material by comparison with known cross sections for a set of materials in a library. The method used is a root-mean-squared (RMS) difference between the calculated and known cross sections. This process is described in Favorite and Bledsoe (2006). The 40 materials included in the library are listed in Table 1 (this is the same set used in Favorite and Bledsoe (2006)). The cross sections used in the library were obtained from the MCPLIB04 continuous-energy photon cross section library, (White, 2003) used by the Monte Carlo code MCNP (X-5 Monte Carlo Team, 2003). 2.4. Unknown source isotope weight fractions For the unknown source isotope weight fractions problem, the following equation can be derived for a discrete energy line (Favorite, 2004):
R
qs NA
R ^ wq dV dX
rn
In Favorite (2004), Eq. (1) was written explicitly in terms of Drn r0n r n by using Taylor expansions of In ðr0n Þ and Is;n ðr 0n Þ about rn. It was noted in that paper that the truncation of these Taylor expansions may lead to inaccuracies. Here a different approach is used: Eqs. (2) and (3) are approximated using the surface-layer formulas (Rahnema, 1996)
ð7Þ
where region k is inside or below interface n and region k + 1 is outside or above interface n. If measurements are taken of a total of G discrete energy lines, the left side of Eq. (1) can be recast as a G N matrix and the right side can be recast as a G 1 vector. The system can then be solved for the update to the unknown interface location, Drn, n = 1, . . ., N (Favorite, 2004).
and
Is;n ðr 0n Þ
ð6Þ
¼
"Z
dV
Z
^ w w dX
Vs
J X rt;j j¼1
M qi Df j M 0 Aj Aj
Z
dV
Z
# ^ w Dfi dX
Vs
M M0 ; M0 ð9Þ
where fj is the weight fraction of isotope j in the source and Dfj fj0 fj . In addition, J is the number of isotopes with unknown abundances in the source material, Aj the gram atomic weight of isotope j, qs the mass density of the source material, NA Avogadro’s number, rt,j the microscopic total cross section of isotope j, and qi the source strength of isotope i [c/(atom of isotope i) s]. Eq. (9) can be recast as a G J matrix equation with a G 1 right-hand side and solved for Dfj. 2.5. Unknown mass densities Now suppose that the geometry of the system is known but the material mass densities are not. Favorite (2005) presented the equation that can be solved for Dqn, the update for the mass density of material n:
968
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
Table 1 Materials and densities included in the shield material library. Density (g/cm3)
Material Foam Polyurethane Wood 2 Wood 1 Wax Fiberglass Water Seawater Natural rubber Lucite C-phenolic Kevlar Beryllium Dry ground a
R
Material a
0.021 0.10 0.224 0.50 0.93 1.00 1.00 1.025 1.10 1.16 1.40 1.45 1.85 1.90
Generic shield Ordinary concrete Portland concrete Wet ground Aluminum Beryllium oxide Titanium Vanadium Type 304 stainless steel Iron Carbon steel Type 316 stainless steel Cadmium
Density (g/cm3)
Material
Density (g/cm3)
1.91 2.35 2.35 2.60 2.70 3.025 4.50 5.80 7.86 7.86 7.86 7.98 8.65
Nickel Cobalt Copper Bismuth Silver Niobium FS-85 Lead Thorium Tantalum Gold Tungsten Neptunium Platinum
8.902 8.90 8.96 9.80 10.5 10.7 11.4 11.7 16.6 19.3 19.3 20.4 21.4
5.1% hydrogen, 23.5% carbon, 5.1% nitrogen, 36.4% oxygen, 18.3% aluminum, and 15% lead.
R
N X Dq
1
n
^w q qn dV dX n¼1 M M0 ; ¼ M0
Rt;n
Z
dV
Z
^ w w dX
Vn
M q M0 n
Z
dV
Z
^ w dX
Vn
ð10Þ
where qn is the mass density of material n in the current iteration and N is the total number of materials with unknown densities. The volume integrals in Eq. (10) are over the entire region containing material n. If a total of G energy lines are measured, Eq. (10) can be recast as a G N matrix equation with a G 1 right-hand side. 3. Implementation for finite cylinders For the implementation of the Schwinger inverse method in two-dimensional cylindrical geometries, the adjoint flux integrals, forward-adjoint inner product integrals, and detector fluxes in Eqs. (1) and (8)–(10) were solved for using a deterministic ray-tracing code (Favorite et al., 2009). In the only previous paper in which we considered cylindrical geometries (Bledsoe et al., 2007), a discrete-ordinates code (Alcouffe et al., 2005) was used to perform these calculations. Due to the ray-effects problem inherent in discrete-ordinates calculations (Lewis and Miller, 1984), the scalar flux at external points could not be used as the quantity of interest in that paper. The total (4p) leakage from the system, an extremely difficult quantity to measure, was used instead. The forward-adjoint ray-tracing described in Favorite et al. (2009) allows the application of the Schwinger method to realistic two-dimensional problems. A Fortran code was written to read the geometry and materials from a user-generated input file, use this information in the raytrace routine to calculate fluxes and inner products, and then perform the necessary calculations to solve Eqs. (1) and (8)–(10) for the unknowns. In addition to geometry and materials, the usergenerated input file contains information on which type of inverse problem is to be solved, which radii, heights, shield material, weight fractions, or densities are unknown, and the measured fluxes for the discrete energy lines considered. This input file will be referred to as the Schwinger input file. The two-dimensional spatial discretization for the ray-trace calculations used mesh boundaries at each material interface. The ray-tracing routine evaluates uncollided point detector fluxes as well as uncollided adjoint and forward-adjoint product integrals in volumes and on surfaces using a deterministic quadrature with thousands of angles in a spherical coordinate system centered at the detector location (Favorite et al., 2009). For the problems in this paper, 1000 angular mesh points were used in both the polar and azimuthal directions.
Interface identification in two-dimensional geometries presents several interesting difficulties that don’t arise in one-dimensional spheres. Consider the geometry shown in Fig. 1. Suppose that the location of the source outer radius is unknown, but it is known to lie somewhere between the inner void cavity (2.0 cm) and the outer edge of the geometry (8.0 cm). Thus its location with respect to the radii of the top and bottom aluminum layers is unknown. In this case it is possible that a postulated source radius is greater than 7.0 cm, meaning it is on the wrong side of the aluminum radius (Fig. 2a). Another possibility is that the source radius is estimated to be equal to the radius of the aluminum regions at 7.0 cm (Fig. 2b). These cases illustrate that at times interfaces may cross and the mesh structure must be changed during the iterative process. Details on these issues are given in the appendix. Various constraints exist that the solutions of Eqs. (1) and (8)– (10) must satisfy. For interface location problems, interface values cannot be negative nor may two radii or heights cross except in those instances that the user allows for the crossing of a pair (discussion of the crossing of interfaces is given in the appendix). For source isotope identification, isotope weight fractions cannot be negative. Similarly for shield identification, calculated total cross sections for the unknown material cannot be negative, and for mass density problems, calculated densities cannot be negative. In the previous implementations of the Schwinger method in Favorite (2004), Favorite and Bledsoe (2006), and Bledsoe et al. (2007), a violation of any of these constraints led to a failure of the method. In the implementation of the method for this paper, however, a technique was used to try to ‘‘fix” the calculated results
Heights (cm) 10.0 9.0 7.0
Detector Location 1 (r,z) = (10.0, 7.0) cm
Source 3.0 Aluminum 1.0 Nickel 2.0
6.0 7.0 8.0
Detector Location 2 (r,z) = (0.0, -1.0) cm Fig. 1. Sample cylindrical geometry with cutaway to show internal structure. The central cavity is void.
969
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
Fig. 2. Two possible initial guesses for unknown source radius for the geometry given in Fig. 1.
so that the constraints were not violated. If in a given iteration the results violate the problem constraints, an analysis of the matrix equation [matrix forms of Eqs. (1) and (8)–(10)] used to solve for the unknowns is performed. In the Schwinger method the matrix equations are solved using singular value decomposition (SVD) (Press et al., 1994), a method that decomposes any matrix into a product of three matrices, one of which is a square diagonal matrix of positive or zero values called singular values. If an element of this matrix is zero, the original matrix (the matrix on the left side of Eqs. (1) and (8)–(10)) is singular. The original matrix can also be treated as singular if the ratio of the minimum to maximum singular value is below some preset threshold. If the matrix on the left side of Eqs. (1) and (8)–(10) is initially singular when the constraints are violated, the singular value threshold is adjusted so that the matrix is guaranteed to be nonsingular, and the SVD algorithm is carried out again. If, on the other hand, the matrix is initially nonsingular, the threshold is adjusted so that the matrix is guaranteed to be singular and the calculations are performed again. This method of automatically adjusting the singular value threshold was first used with the Schwinger method in Favorite (2005). If the fix of adjusting the singular value threshold is unsuccessful, an attempt is made at damping the update to the parameter values. Each of the calculated updates is multiplied by a damping factor, denoted a, that is equal to the minimum of unity or unity divided by the norm of the set of updates. This decreases the magnitude of the change of the parameter values from the previous generation and thus increases the likelihood that constraints are satisfied. Suppose, however, that after the damping factor is introduced the constraints are still violated. In this case, the damping factor is kept while the process of making the nonsingular matrix singular (or vice versa) is repeated. If once again the solutions violate the constraints (meaning a violation has occurred three times), the damping factor is adjusted to a/10 and the calculations are repeated. If the constraints are violated four times, the new damping factor is kept while the routine again switches the singular matrix to nonsingular (or vice versa) and recalculates the unknowns. If constraints are violated five times the method stops.
(HEU) source with 94.73 wt% 235U and 5.27 wt% 238U and a mass density of 18.74 g/cm3, which was similar to the composition of the Godiva assembly that was used for critical experiments at Los Alamos National Laboratory (Malenfant, 1981). Above and below the source is a layer of aluminum (density 2.70 g/cm3) and on the outside is a layer of nickel (density 8.902 g/cm3). The measured energy lines are the four strongest gamma emission lines for uranium (144, 186, 766, and 1001 keV). The quantity of interest for the test problems is the sum of the scalar fluxes for each of these lines over two detector points. The first detector point was at a radius of 10.0 cm and a height of 7.0 cm, and the second was at a radius of 0.0 cm and a height of 1.0 cm. The logarithmic derivatives of the fluxes with respect to source/ shield interface locations were analyzed to see if this is truly a twodimensional problem. The derivatives ranged between 4.7% and 15.7% for the four energy lines with respect to the outer source radius (6.0 cm). For the top of the source (9.0 cm), they ranged between 0.08% and 1.48%, and for the bottom of the source (1.0 cm) they ranged between 40.6% and 66.23%. Thus the results depend on both the axial and radial source locations, and this is a two-dimensional problem. The measured fluxes were simulated using two different methods. The first set was simulated by using the ray-trace code with the same angular quadrature as used by the inverse method transport calculations. Thus this set of measurements was exactly consistent with the iterative identification processes. These measurements will be referred to as ray-trace measurements. The second set of measurements was simulated stochastically using MCNP in order to represent more realistic measurements. These will be referred to as Monte Carlo measurements. The simulated total fluxes are shown in Table 2. For each test problem the convergence criterion was a relative difference of 0.01% between calculated and measured fluxes for
Table 2 Simulated flux measurements (a) for the test cases. Line energy (keV)
4. Test cases A set of numerical test cases will illustrate the Schwinger inverse method in a cylindrical geometry for each type of inverse problem. For the test cases, the geometry of Fig. 1 will again be considered. The geometry contains a 16-kg high-enriched uranium
144 186 766 1001 a b
Monte Carlob (%)
Ray-trace 2
5.5317 10 5.7794 103 5.2982 100 1.9521 101
Expressed in c/cm2 s. 1r relative standard deviations are given.
5.38 102 ± 8.51 5.78 103 ± 6.15 5.34 100 ± 5.34 1.99 101 ± 3.93
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
each energy line. In the event that this could not be met, four other stopping criteria were implemented: 1. A stop due to continued violation of constraints, as outlined in Section 3. 2. A test to see if the iterative results are diverging instead of converging toward the solution. To test this, the Fortran code checks the value of the root-mean-square difference between the calculated and measured fluxes over the previous ten iterations. If this difference increases every time the iteration number increases, the problem appears to be diverging and the calculations are stopped (this test is only administered after ten iterations have been completed). 3. The fluxes are changing very little with each iteration. To quantify this, the method checks the last ten values of the rootmean-square difference between calculated and measured fluxes and determines its slope over these iterations. The value of this slope is then divided by the average root-mean-square difference between calculated and measured fluxes over the ten iterations. If the resulting relative slope is less than a prescribed value (0.01), the method stops. 4. The same values for all the unknowns are calculated twice (i.e. in any two iterations). The script looks at all input files for repeats. In this case ‘‘same” means that the first seven digits are equal. Once the method stops, the set of calculated unknowns for the iteration with the lowest v2 difference between the calculated and measured fluxes is considered the solution. This v2 difference is given by
v2 ¼
G X
g
M
g¼1
rg0
2 Mg0
;
ð11Þ
where rg0 is the measurement uncertainty for line g. The set of calculated values for the unknowns corresponding to the iteration with the lowest v2 difference is called the best model. 4.1. Unknown interface locations (see Table 3) For Case 1 suppose the location of all the axial interfaces are known but the location of the radial interfaces at 2.0 cm and 6.0 cm are unknown. The outer source radius is known to be greater than 2.0 cm and is inside the outer edge of the geometry at 8.0 cm, but its location with respect to the aluminum radius at 7.0 cm is unknown. Here the outer source radius was postulated to be 7.2 cm (a situation similar to Fig. 2a). The inner source radius at 2.0 cm was guessed to be 3.0 cm. The Schwinger input file in this case indicated that the second and third radii were permitted to cross. The crossing did occur, and with ray-trace measurements convergence occurred in seven iterations. The outer source radius was calculated very accurately to be 6.0001 cm, but the inner radius was calculated to be a somewhat inaccurate 2.3691 cm. The reason for this is that the source is so thick that the measured flux is insensitive to the inner radius. With MCNP measurements, the outer source radius was calculated accurately at 6.054 cm, but the inner radius had a value of 5.0168 cm. Although this is far from the actual location of the inner radius, the agreement between calculated and measured fluxes was very good, with a v2 difference of 0.2371. Thus this model belongs to the family of solutions of the inverse problem. The method took 14 iterations to converge with MCNP measurements. In Case 2 suppose that a radial interface and axial interface are unknown. The outer source radius at 6.0 cm was unknown and guessed to be 6.3 cm, while the height at 1.0 cm was unknown and guessed to be 0.8 cm. Using ray-trace measurements, very
accurate values of 6.002 cm and 1.0023 cm were calculated in 16 iterations. With Monte Carlo measurements, the unknown height was calculated with high accuracy to be 1.044 cm. The unknown radius was calculated to be 6.218 cm. These values correspond to a v2 difference of 0.06075. Thus once again the method converged to a member of the family of solutions. 4.2. Unknown shield material In Case 3 the bottom aluminum portion of the shield was unknown and postulated to be nickel. The convergence of the Schwinger inverse method for this case can be seen in Fig. 3. When ray-trace measurements were used the method converged to nearly the exact cross sections of aluminum. Aluminum was clearly identified as the unknown material, as shown in Table 4, where the RMS difference for aluminum (0.00017%) was three orders of magnitude smaller than its nearest competitor, wet ground (0.501%). With Monte Carlo measurements, the 186-keV line converged to approximately 4% of the aluminum cross section, but the 1001-keV line converged to only within about 8.2%. The 144keV and 766-keV lines (not shown in Fig. 3) converged to within 5% and 3%, respectively, of the aluminum cross sections. As a result of these differences, aluminum was one of six materials identified as likely candidates for the unknown. Aluminum had the lowest RMS difference at 1.150%, just better than the nearest competitor, wet ground (1.228%). The other likely materials were beryllium oxide (2.785%), ordinary concrete (2.963%), Portland concrete (2.968%), and dry ground (7.511%). 4.3. Unknown source weight fractions Case 4 featured unknown weight fractions for the source material, with both guessed at 0.50000. With ray-trace measurements, the correct weight fractions of 0.9473 for 235U and 0.0527 for 238 U were calculated in a single iteration. Using Monte Carlo measurements, weight fractions of 0.946545 for 235U and 0.053455 for 238 U were calculated in three iterations. This model corresponds to a v2 difference of 0.1349. The iterative history for the case of Monte Carlo measurements is shown in Table 5. 4.4. Unknown material densities In Case 5 the densities of the source and aluminum layers were unknown. The source was postulated to have density 13.00 g/cm3
Percent Difference from True Cross Section
970
350% 186 keV; Ray-Trace
300%
1001 keV; Ray-Trace 186 keV; Monte Carlo
250%
1001 keV; Monte Carlo
200%
3.00% 0.00% -3.00%
150%
-6.00%
100%
-9.00% 3
4
5
50% 0% -50% 0
1
2
3
4
5
Iteration Fig. 3. Convergence of the Schwinger inverse method for Case 3. The unknown material was guessed to be nickel in iteration 0; it was actually aluminum. On the large graph the cross sections all appear to converge to 0%; differences can be seen in the inset.
971
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973 Table 3 Results of interface identification problems. Descriptor
Model Outer radii (cm)
Case 1
Actual Model Initial Converged – RT Converged – MC
Void
Source
Aluminum
Nickel
2.00 3.00 2.3691 5.0168
6.00 7.20 6.0001 6.054
7.00 7.00 7.00 7.00
8.00 8.00 8.00 8.00
Aluminum
Source
Void
Source
Aluminum
1.00 1.00 1.00 1.00
3.00 3.00 3.00 3.00
7.00 7.00 7.00 7.00
9.00 9.00 9.00 9.00
10.00 10.00 10.00 10.00
Void
Source
Aluminum
Nickel
2.00 2.00 2.00 2.00 Heights (cm)
6.00 6.30 6.002 6.218
7.00 7.00 7.00 7.00
8.00 8.00 8.00 8.00
Aluminum
Source
Void
Source
Aluminum
1.00 0.8 1.0023 1.044
3.00 3.00 3.00 3.00
7.00 7.00 7.00 7.00
9.00 9.00 9.00 9.00
10.00 10.00 10.00 10.00
Heights (cm)
Actual Model Initial Converged – RT Converged – MC
Outer Radii (cm)
Case 2
Actual Model Initial Converged – RT Converged – MC
Actual Model Initial Converged – RT Converged – MC Italics represent known (fixed) quantities.
Table 4 Materials with the lowest RMS difference for Case 3. Density (g/cm3)
RMS difference (%)
Ray-trace measurements 1 Aluminum 2 Wet ground 3 Portland concrete 4 Ordinary concrete 5 Beryllium oxide
2.70 2.60 2.35 2.35 3.025
0.00017 0.501 2.472 2.492 2.687
Monte Carlo measurements 1 Aluminum 2 Wet ground 3 Beryllium oxide 4 Ordinary concrete 5 Portland concrete 6 Dry ground
2.70 2.60 3.025 2.35 2.35 1.90
1.150 1.228 2.785 2.963 2.968 7.511
Order
Material
and aluminum was postulated to have density 5.00 g/cm3. With ray-trace measurements the source density was calculated to be 18.742 g/cm3 and the aluminum density was calculated to be 2.701 g/cm3. With Monte Carlo measurements the source density was calculated to be much higher than the actual source density, with a value of 25.469 g/cm3, but the aluminum density was calculated to be an accurate 2.758 g/cm3. These values have differences of 35.91% and 2.15%, respectively, to the actual densities. Though the source density is far from the actual density, these
numbers result in good agreement between calculated and measured fluxes, with a v2 difference of 0.1361. The iterative histories for this problem are presented in Table 6. 5. Summary and conclusions The Schwinger inverse method was originally applied to the inverse transport problems of identifying unknown interface locations and unknown source isotope weight fractions in spherical radiation source/shield systems in Favorite (2004). In Favorite and Bledsoe (2006), the method was implemented for shield identification on numerical test cases in one-dimensional spherical geometries and in an analytic monodirectional slab problem. Favorite (2005) presented numerical solutions for unknown mass density problems in one-dimensional spheres. The method was then applied for material identification in a two-dimensional cylinder in Bledsoe et al. (2007), but the geometry considered was unrealistically small and the quantity of interest was the total (4p) leakage from the system. In this paper, the method is applied to interface location identification, shield material identification, source weight fraction identification, and mass density identification in realistic finite twodimensional cylindrical geometries. In addition to the application to more realistic (two-dimensional) problems, this paper extends previous work on the Schwinger inverse method by using surface formulas for unknown interface locations, automatic correction at-
Table 5 Iteration history of source weight fraction identification for Case 4; Monte Carlo measurements. Iteration
235
Difference from actual weight fraction (%)
238
Difference from actual weight fraction (%)
0 1 2 3
0.5000 0.944535 0.946540 0.946545
47.22 0.29 0.08 0.08
0.5000 0.055465 0.053460 0.053455
848.77 5.25 1.44 1.33
U weight fraction
U weight fraction
972
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
Table 6 Iteration history of density identification for Case 5. Difference from actual source density (%)
Aluminum density (g/cm3)
Difference from actual aluminum density (%)
Ray-trace measurements 0 13.000 1 11.553 2 11.043 3 17.606 4 18.742
30.63 38.35 41.07 6.05 0.01
5.000 5.844 3.386 2.354 2.701
85.19 116.44 25.41 12.81 0.04
Monte Carlo measurements 0 13.000 1 31.722 2 31.459 3 26.020 4 25.604 5 25.456 6 25.469
30.63 69.27 67.87 38.85 36.63 35.84 35.91
5.000 0.421 3.269 2.863 2.761 2.758 2.758
85.19 84.41 21.07 6.04 2.26 2.15 2.15
Iteration
Source density (g/cm3)
tempts for violated constraints, and ray-tracing instead of discreteordinates for transport calculations, which enables us to use the scalar flux at a point as the quantity of interest. For the interface location problem, the issue addressed was whether the Schwinger inverse method could successfully identify unknown radii or heights of material interfaces in the interior of the source/shield system. In Favorite (2004), the answer to this question was yes for identification of unknown radii. Here the answer was yes for both radii and heights. When measurements were used that were perfectly consistent with the transport calculations of the iterative method, the Schwinger method converged nearly perfectly to the correct interface locations in both cases considered. When less consistent, more realistic, measurements were used, the Schwinger method calculated one interface very accurately and one inaccurately in both cases. However, in these cases the method converged to a model that belonged to a family of solutions that corresponds to small v2 differences between calculated and measured fluxes. For the shield identification problem, two issues were addressed in the numerical test cases considered here. The first was whether the method would converge to the correct cross section set. When perfectly consistent measurements were used, the Schwinger method converged nearly perfectly to the correct cross sections in all cases. When less consistent, more realistic simulated measurements were used, the method converged to within 15% for each energy line. The second issue considered was whether the calculated cross sections could be used to identify the unknown material. In the spherical problems of Favorite and Bledsoe (2006), the answer was yes. Here the answer was the same. With perfectly consistent measurements the actual material was correctly identified as the only suitable candidate in both test cases. When less consistent measurements were used, the difference between the converged cross sections and the actual cross sections of the unknown material led to some ambiguity. Six possible shield materials (including the correct one) were identified. For the source identification problem the issue addressed was whether the method could identify the correct weight fractions of the isotopes in the photon source. In the spherical test cases of Favorite (2004), the answer was yes. Here the answer is also yes. When perfectly consistent measurements were used the exact weight fractions were identified in one iteration for both test cases. Using less accurate measurements the method was still highly successful, with errors smaller than 2% for the calculated weight fractions. With these measurements, convergence required three iterations. In mass density identification problems the issue addressed was whether the method could identify the correct densities of the
materials in the source and shield. In the spherical test cases of Favorite (2005), the answer was yes. Here the answer was the same. When perfectly consistent measurements were used, the Schwinger method calculated densities to within 0.04% of the actual densities. With more realistic measurements the Schwinger method found a source density that differed from the actual source density by 35.91%. The shield density differed from the actual density by 2.15%. Although the source density differed greatly from the actual density, the final model output matched the measurements to within 1%. In this paper, the Schwinger inverse method was applied to realistic two-dimensional cylindrical inverse transport problems. By using a ray-tracing algorithm, the ray effects inherent in the discrete-ordinates code used in Favorite (2004), Favorite and Bledsoe (2006) and Bledsoe et al. (2007) have been eliminated. Since ray effects become worse as the distance between the source and detector increases, in our previous attempts at solving inverse problems we were forced to either use the total (4p) leakage from the system (which tends to wash out ray effects, but measurements of which are prohibitively difficult) as the quantity of interest or make the geometries unrealistically small. With ray effects eliminated in this paper, we were able to explore much larger geometries and use the sum of the scalar fluxes at a few points external to the geometry as the quantities of interest. There are several potential applications of the Schwinger inverse method. Identification of radioactive source constituents and dimensions is of enormous concern for national security. Identifying and determining the spatial locations of materials is also important in radioactive waste assay. The identification of shield material given a known radioactive source has application in industry, where, for instance, the composition of mined materials or manufactured alloys may be needed. Acknowledgements The authors wish to thank D.I. Ketcheson for contributions to the interface-location problem and the ray-tracing code and D.E. Vaughan for providing the damping algorithm. Appendix A. Identifying material interface locations in twodimensional cylinders and the method for the crossing of radii or heights This appendix describes how the Fortran code identifies the unknown material interface locations and how it can accommodate the crossing of two interface radii or heights. These considerations are necessary for the solution of two-dimensional problems. Consider the example shown in Fig. 1. Suppose first that in the initially
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
postulated system the source radius is equal to the radius of the aluminum layers. This example illustrates that with two-dimensional methods it is not sufficient to simply state in the Schwinger input file that an entire radius is unknown (as was the case with one-dimensional methods). In this case, the aluminum radii are known to be at 7.0 cm and should not change, but the source radius postulated to be at 7.0 cm will change. If the Schwinger input file simply stated that the 7.0 cm radius was unknown, the Fortran code would wrongly assume that all interfaces at 7.0 cm were unknown. To avoid this problem the Schwinger input file contains an array that has an element for each mesh. If the outer radius of the mesh is unknown (or, for height identification, if the top of the mesh is unknown) the element of the array corresponding to that mesh is assigned a value of 1. If the outer radius is known, the element is assigned a value of 0. In this way the Fortran code is able to discern exactly where the unknown boundary is. This array will be referred to as the array of unknown interfaces. For the unknown source radius problem considered here, the array of unknown interfaces is initially given by
0 0
0
0 1 0 0 1 0 0 1 0 0 0
0
Suppose the source radius postulated to be equal to the radius of the aluminum layers is calculated by the Schwinger method to have a different location. The result is a source radius moved to the inside or outside of the radius of the aluminum layers. A fourth mesh must now be created because the aluminum and source radii are no longer equal (if the source radius is moved inward, the situation is similar to Fig. 1, if it is moved outward, the situation is similar to Fig. 2a. With the additional mesh, the array of unknown interfaces must be updated to the form (assuming the source radius is moved inward)
0 0
0 0
0 1 0 0 0 1 0 0 0 1 0 0 0 0
0 0
The Fortran code was written to automatically handle this, and the iterative process continues until the unknown is identified. Suppose now that the source is postulated to have a radius greater than that of the aluminum layers. It is thus necessary for the radii of the aluminum layers and source layer to cross if the unknown is to be identified. In one dimension, the crossing of any two radii is not allowed in the Schwinger inverse method
973
(Favorite, 2004) but for interface identification in two dimensions the user may specify two interfaces that may cross (this specification occurs in the Schwinger input file). For the geometry shown in Fig. 1, the Schwinger input file would indicate that the second and third radii may cross. Note that when this crossing occurs, the array of unknown interfaces changes from
0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 to
0 0
0 0
0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 The Fortran code was written to handle this crossing and reassign the array of unknowns. References Alcouffe, R.E., Baker, R.S., Dahl, J.A., Turner, S.A., Ward, R.C., 2005. PARTISN: A TimeDependent, Parallel Neutral Particle Transport Code System. LA-UR-05-3925, Los Alamos National Laboratory. Bledsoe, K.C., Favorite, J.A., Aldemir, T., 2007. Material identification in finite cylindrical geometries using the Schwinger inverse method. Trans. Am. Nucl. Soc. 96, 545–547. Favorite, J.A., 2004. Using the Schwinger variational functional for the solution of inverse transport problems. Nucl. Sci. Eng. 146, 51–70. Favorite, J.A., 2005. Identification of unknown densities in a source/shield system using the Schwinger inverse method. Trans. Am. Nucl. Soc. 93, 588–589. Favorite, J.A., Bledsoe, K.C., 2006. Identification of an unknown material in a radiation shield using the Schwinger inverse method. Nucl. Sci. Eng. 152, 106– 117. Favorite, J.A., Bledsoe, K.C., Ketcheson, D.I., 2009. Surface and volume integrals of uncollided adjoint fluxes and forward-adjoint flux products. Nucl. Sci. Eng., in press. Lewis, E.E., Miller, W.F., 1984. Computational Methods of Neutron Transport. John Wiley and Sons (Chapter 4). Malenfant, R.E., 1981. Los Alamos Critical Assemblies Facility, LA-8762-MC(UC-46). Los Alamos Scientific Laboratory, Los Alamos, New Mexico. X-5 Monte Carlo Team, 2003. MCNP-A General Monte Carlo N-Particle Transport Code, Version 5. LA-UR-03-1987, vol. I. Los Alamos National Laboratory. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1994. Numerical Recipes in FORTRAN: The Art of Scientific Computing, second ed. (reprinted with corrections), (Chapter 2). Cambridge University Press. Rahnema, F., 1996. On the internal interface perturbation. Ann. Nucl. Energy 23, 1401. White, M.C., 2003. Photoatomic Data Library MCPLIB04: A New Photoatomic Library Based on Data from ENDF/B-VI Release 8. LA-UR-03-1019 (Revised), Los Alamos National Laboratory (February).
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Using the Schwinger inverse method for solutions of inverse transport problems in two-dimensional cylindrical geometries Keith C. Bledsoe a,*, Jeffrey A. Favorite b, Tunc Aldemir a a b
Nuclear Engineering Program, The Ohio State University, Columbus, OH 43210, USA Applied Physics Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
a r t i c l e
i n f o
Article history: Received 27 October 2008 Received in revised form 26 February 2009 Accepted 27 February 2009 Available online 29 April 2009
a b s t r a c t The Schwinger method for solving inverse transport problems is applied to the problems of interface location identification, shield material identification, source isotope weight fraction identification, and material mass density identification (separately) in multilayered two-dimensional cylindrical gamma radiation source/shield systems. The method is iterative and estimates unknown interface locations, source isotope weight fractions, and material densities directly, while the unknown shield material is identified by estimating its total macroscopic gamma-ray cross sections. The energies of discrete gamma-ray lines emitted by the source are assumed to be known, while the unscattered flux of the lines is assumed to be measured at points external to the system. In numerical test cases, the Schwinger method correctly identifies the unknowns when the same deterministic ray-tracing code is used for both the parameter estimation process and simulation of the measured data. With realistic simulation of the measured data using a Monte Carlo code, the method produces more ambiguous results for interface location, shield material identification, and material density identification. The method works well for source weight fraction identification with measured data simulated by Monte Carlo. In addition to the application to more realistic (two-dimensional) problems, this paper extends previous work on the Schwinger inverse method by using surface formulas for unknown interface locations, automatic correction attempts for violated constraints, and ray-tracing instead of discrete-ordinates for transport calculations. Published by Elsevier Ltd.
1. Introduction The Schwinger inverse method was recently presented (Favorite, 2004, 2005; Favorite and Bledsoe, 2006; Bledsoe et al., 2007) for obtaining solutions to the problems of determining, in a gamma radiation source/shield system: (a) a set of unknown internal interface locations, (b) the isotopic content of an unknown source, (c) interface locations and source isotopic content simultaneously, (d) the composition of an unknown layer in the shield, and (e) unknown material mass densities. Favorite (2004) presented the mathematical background for problems (a–d) and documented numerical solutions for problems (a) and (b) in one-dimensional spherical geometries. Numerical results for the unknown shield layer problem were presented in Favorite and Bledsoe (2006), also for one-dimensional spherical geometries. In Favorite (2005) numerical results for the mass density problem were presented for one-dimensional spherical problems. The method was successful for identification of the unknowns in each of those papers. In Bledsoe et al. (2007) the method was used successfully for material identification in cylindrical systems, but the geometries were small
* Corresponding author. Address: X-1-TA, MS P365, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. E-mail address: [email protected] (K.C. Bledsoe). 0306-4549/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.anucene.2009.02.014
and the required leakage measurements would be extremely difficult to make. In this paper, the Schwinger inverse method is implemented and tested in realistic two-dimensional cylindrical geometries for the problems of identifying (separately) unknown interface locations, a shield material, source isotope weight fractions, and material mass densities. In addition to the application to multidimensional problems, this paper extends previous work on the Schwinger inverse method by using surface formulas for unknown interface locations, automatic correction attempts for violated constraints, and ray-tracing instead of discrete-ordinates for transport calculations. 2. Methodology of the Schwinger inverse method 2.1. General approach The Schwinger inverse method is derived from a perturbationtheory approach to the inverse transport problem (Favorite, 2004). Instead of calculating the effect of a system perturbation on a quantity of interest, the quantity of interest is assumed to be given (from a measurement) and the Schwinger functional is manipulated to produce an equation describing the system perturbation. This equation is applied iteratively to determine the unknown parameters in the system. The quantity of interest is the
967
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
flux of a set of discrete gamma-ray lines measured at a point (or the sum of fluxes measured at points) external to the geometry, which implies that the scattering term in the transport equation is assumed to be negligible. The mathematics of the Schwinger inverse method are presented in detail in Favorite (2004) and Favorite and Bledsoe (2006). Here a brief review of those results will be given. Begin by considering a radiation shield surrounding a gamma-ray source. Both the source and the shield may be multilayered but the layers are assumed to be homogeneous. The source emits gamma rays at discrete energies or lines, which can be resolved quite well using a high-purity germanium detector. It is assumed that there is no scattering into the energy lines. In Favorite (2004), the actual (unknown) configuration was considered the perturbed system, while an estimate of the configuration was the unperturbed system. The Schwinger functional was used to derive an expression for the quantity of interest in the system (in our case, unscattered scalar fluxes measured at points external to the system) that is second-order accurate with respect to the difference between the unperturbed and perturbed scalar fluxes. By manipulating this expression, equations were derived to solve for the unknown system parameter. ^ Þ be the forward angular flux For a discrete energy line, let wðr; X ^ Þ the adjoint angu^ , w ðr; X of gamma rays at position r and angle X ^ , Rt(r) the total lar flux of gamma rays at position r and angle X macroscopic cross section at position r, q(r) the source, at position r, of gamma rays (c/cm3 s), M the sum of the calculated flux of photons over all detector points (c/cm2 s), and M0 the sum of the measured flux of photons over all detector points (c/cm2 s). 2.2. Unknown interface locations If the system contains unknown radii and/or unknown heights, the equation for the unknowns Drn is (Favorite, 2004)
" # N N X 1 M X M M0 0 0 DRt;n In ðr n Þ Dq I ðr Þ ¼ ; R R ^ w q M0 n¼1 n s;n n M0 dV dX n¼1
ð1Þ
where N is the total number of unknown interfaces in the system and In ðr 0n Þ and Is;n ðr0n Þ are defined by
In ðr0n Þ
Z
r 0n
dV
Z
^ w w dX
ð2Þ
rn
DRt;n Rt;k Rt;kþ1 and
Dqn qk qkþ1 ;
Z
r0n
dV
Z
^ w : dX
ð3Þ
2.3. Unknown shield material For the unknown shield material problem, the following equation can be derived for the perturbed cross section for a discrete energy line (Favorite, 2004; Favorite and Bledsoe, 2006):
R0t ¼ R
In ðr0n Þ Dr n
Z
dS
Z
^ w w dX
ð4Þ
^ w ; dX
ð5Þ
and
Is;n ðr 0n Þ Dr n
Z
dS
Z
R where dS represents an integral over the unperturbed radial or axial surface r = rn. These approximations eliminate the inaccuracies introduced by truncation of a Taylor series. The terms DRt;n and Dqn are the differences in cross section and source density across interface n; they are defined by
R
R
^ w q M M 0 dX þ Rt : R ^ w w M0 dV dX DV dV
ð8Þ
In Eq. (8) it is assumed that the unknown cross section is constant over the unknown region. The integral in the numerator in Eq. (8) is over the entire problem, but the integral in the denominator is over the unknown material region DV only. After Eq. (8) converges, the calculated cross sections are used to identify the unknown material by comparison with known cross sections for a set of materials in a library. The method used is a root-mean-squared (RMS) difference between the calculated and known cross sections. This process is described in Favorite and Bledsoe (2006). The 40 materials included in the library are listed in Table 1 (this is the same set used in Favorite and Bledsoe (2006)). The cross sections used in the library were obtained from the MCPLIB04 continuous-energy photon cross section library, (White, 2003) used by the Monte Carlo code MCNP (X-5 Monte Carlo Team, 2003). 2.4. Unknown source isotope weight fractions For the unknown source isotope weight fractions problem, the following equation can be derived for a discrete energy line (Favorite, 2004):
R
qs NA
R ^ wq dV dX
rn
In Favorite (2004), Eq. (1) was written explicitly in terms of Drn r0n r n by using Taylor expansions of In ðr0n Þ and Is;n ðr 0n Þ about rn. It was noted in that paper that the truncation of these Taylor expansions may lead to inaccuracies. Here a different approach is used: Eqs. (2) and (3) are approximated using the surface-layer formulas (Rahnema, 1996)
ð7Þ
where region k is inside or below interface n and region k + 1 is outside or above interface n. If measurements are taken of a total of G discrete energy lines, the left side of Eq. (1) can be recast as a G N matrix and the right side can be recast as a G 1 vector. The system can then be solved for the update to the unknown interface location, Drn, n = 1, . . ., N (Favorite, 2004).
and
Is;n ðr 0n Þ
ð6Þ
¼
"Z
dV
Z
^ w w dX
Vs
J X rt;j j¼1
M qi Df j M 0 Aj Aj
Z
dV
Z
# ^ w Dfi dX
Vs
M M0 ; M0 ð9Þ
where fj is the weight fraction of isotope j in the source and Dfj fj0 fj . In addition, J is the number of isotopes with unknown abundances in the source material, Aj the gram atomic weight of isotope j, qs the mass density of the source material, NA Avogadro’s number, rt,j the microscopic total cross section of isotope j, and qi the source strength of isotope i [c/(atom of isotope i) s]. Eq. (9) can be recast as a G J matrix equation with a G 1 right-hand side and solved for Dfj. 2.5. Unknown mass densities Now suppose that the geometry of the system is known but the material mass densities are not. Favorite (2005) presented the equation that can be solved for Dqn, the update for the mass density of material n:
968
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
Table 1 Materials and densities included in the shield material library. Density (g/cm3)
Material Foam Polyurethane Wood 2 Wood 1 Wax Fiberglass Water Seawater Natural rubber Lucite C-phenolic Kevlar Beryllium Dry ground a
R
Material a
0.021 0.10 0.224 0.50 0.93 1.00 1.00 1.025 1.10 1.16 1.40 1.45 1.85 1.90
Generic shield Ordinary concrete Portland concrete Wet ground Aluminum Beryllium oxide Titanium Vanadium Type 304 stainless steel Iron Carbon steel Type 316 stainless steel Cadmium
Density (g/cm3)
Material
Density (g/cm3)
1.91 2.35 2.35 2.60 2.70 3.025 4.50 5.80 7.86 7.86 7.86 7.98 8.65
Nickel Cobalt Copper Bismuth Silver Niobium FS-85 Lead Thorium Tantalum Gold Tungsten Neptunium Platinum
8.902 8.90 8.96 9.80 10.5 10.7 11.4 11.7 16.6 19.3 19.3 20.4 21.4
5.1% hydrogen, 23.5% carbon, 5.1% nitrogen, 36.4% oxygen, 18.3% aluminum, and 15% lead.
R
N X Dq
1
n
^w q qn dV dX n¼1 M M0 ; ¼ M0
Rt;n
Z
dV
Z
^ w w dX
Vn
M q M0 n
Z
dV
Z
^ w dX
Vn
ð10Þ
where qn is the mass density of material n in the current iteration and N is the total number of materials with unknown densities. The volume integrals in Eq. (10) are over the entire region containing material n. If a total of G energy lines are measured, Eq. (10) can be recast as a G N matrix equation with a G 1 right-hand side. 3. Implementation for finite cylinders For the implementation of the Schwinger inverse method in two-dimensional cylindrical geometries, the adjoint flux integrals, forward-adjoint inner product integrals, and detector fluxes in Eqs. (1) and (8)–(10) were solved for using a deterministic ray-tracing code (Favorite et al., 2009). In the only previous paper in which we considered cylindrical geometries (Bledsoe et al., 2007), a discrete-ordinates code (Alcouffe et al., 2005) was used to perform these calculations. Due to the ray-effects problem inherent in discrete-ordinates calculations (Lewis and Miller, 1984), the scalar flux at external points could not be used as the quantity of interest in that paper. The total (4p) leakage from the system, an extremely difficult quantity to measure, was used instead. The forward-adjoint ray-tracing described in Favorite et al. (2009) allows the application of the Schwinger method to realistic two-dimensional problems. A Fortran code was written to read the geometry and materials from a user-generated input file, use this information in the raytrace routine to calculate fluxes and inner products, and then perform the necessary calculations to solve Eqs. (1) and (8)–(10) for the unknowns. In addition to geometry and materials, the usergenerated input file contains information on which type of inverse problem is to be solved, which radii, heights, shield material, weight fractions, or densities are unknown, and the measured fluxes for the discrete energy lines considered. This input file will be referred to as the Schwinger input file. The two-dimensional spatial discretization for the ray-trace calculations used mesh boundaries at each material interface. The ray-tracing routine evaluates uncollided point detector fluxes as well as uncollided adjoint and forward-adjoint product integrals in volumes and on surfaces using a deterministic quadrature with thousands of angles in a spherical coordinate system centered at the detector location (Favorite et al., 2009). For the problems in this paper, 1000 angular mesh points were used in both the polar and azimuthal directions.
Interface identification in two-dimensional geometries presents several interesting difficulties that don’t arise in one-dimensional spheres. Consider the geometry shown in Fig. 1. Suppose that the location of the source outer radius is unknown, but it is known to lie somewhere between the inner void cavity (2.0 cm) and the outer edge of the geometry (8.0 cm). Thus its location with respect to the radii of the top and bottom aluminum layers is unknown. In this case it is possible that a postulated source radius is greater than 7.0 cm, meaning it is on the wrong side of the aluminum radius (Fig. 2a). Another possibility is that the source radius is estimated to be equal to the radius of the aluminum regions at 7.0 cm (Fig. 2b). These cases illustrate that at times interfaces may cross and the mesh structure must be changed during the iterative process. Details on these issues are given in the appendix. Various constraints exist that the solutions of Eqs. (1) and (8)– (10) must satisfy. For interface location problems, interface values cannot be negative nor may two radii or heights cross except in those instances that the user allows for the crossing of a pair (discussion of the crossing of interfaces is given in the appendix). For source isotope identification, isotope weight fractions cannot be negative. Similarly for shield identification, calculated total cross sections for the unknown material cannot be negative, and for mass density problems, calculated densities cannot be negative. In the previous implementations of the Schwinger method in Favorite (2004), Favorite and Bledsoe (2006), and Bledsoe et al. (2007), a violation of any of these constraints led to a failure of the method. In the implementation of the method for this paper, however, a technique was used to try to ‘‘fix” the calculated results
Heights (cm) 10.0 9.0 7.0
Detector Location 1 (r,z) = (10.0, 7.0) cm
Source 3.0 Aluminum 1.0 Nickel 2.0
6.0 7.0 8.0
Detector Location 2 (r,z) = (0.0, -1.0) cm Fig. 1. Sample cylindrical geometry with cutaway to show internal structure. The central cavity is void.
969
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
Fig. 2. Two possible initial guesses for unknown source radius for the geometry given in Fig. 1.
so that the constraints were not violated. If in a given iteration the results violate the problem constraints, an analysis of the matrix equation [matrix forms of Eqs. (1) and (8)–(10)] used to solve for the unknowns is performed. In the Schwinger method the matrix equations are solved using singular value decomposition (SVD) (Press et al., 1994), a method that decomposes any matrix into a product of three matrices, one of which is a square diagonal matrix of positive or zero values called singular values. If an element of this matrix is zero, the original matrix (the matrix on the left side of Eqs. (1) and (8)–(10)) is singular. The original matrix can also be treated as singular if the ratio of the minimum to maximum singular value is below some preset threshold. If the matrix on the left side of Eqs. (1) and (8)–(10) is initially singular when the constraints are violated, the singular value threshold is adjusted so that the matrix is guaranteed to be nonsingular, and the SVD algorithm is carried out again. If, on the other hand, the matrix is initially nonsingular, the threshold is adjusted so that the matrix is guaranteed to be singular and the calculations are performed again. This method of automatically adjusting the singular value threshold was first used with the Schwinger method in Favorite (2005). If the fix of adjusting the singular value threshold is unsuccessful, an attempt is made at damping the update to the parameter values. Each of the calculated updates is multiplied by a damping factor, denoted a, that is equal to the minimum of unity or unity divided by the norm of the set of updates. This decreases the magnitude of the change of the parameter values from the previous generation and thus increases the likelihood that constraints are satisfied. Suppose, however, that after the damping factor is introduced the constraints are still violated. In this case, the damping factor is kept while the process of making the nonsingular matrix singular (or vice versa) is repeated. If once again the solutions violate the constraints (meaning a violation has occurred three times), the damping factor is adjusted to a/10 and the calculations are repeated. If the constraints are violated four times, the new damping factor is kept while the routine again switches the singular matrix to nonsingular (or vice versa) and recalculates the unknowns. If constraints are violated five times the method stops.
(HEU) source with 94.73 wt% 235U and 5.27 wt% 238U and a mass density of 18.74 g/cm3, which was similar to the composition of the Godiva assembly that was used for critical experiments at Los Alamos National Laboratory (Malenfant, 1981). Above and below the source is a layer of aluminum (density 2.70 g/cm3) and on the outside is a layer of nickel (density 8.902 g/cm3). The measured energy lines are the four strongest gamma emission lines for uranium (144, 186, 766, and 1001 keV). The quantity of interest for the test problems is the sum of the scalar fluxes for each of these lines over two detector points. The first detector point was at a radius of 10.0 cm and a height of 7.0 cm, and the second was at a radius of 0.0 cm and a height of 1.0 cm. The logarithmic derivatives of the fluxes with respect to source/ shield interface locations were analyzed to see if this is truly a twodimensional problem. The derivatives ranged between 4.7% and 15.7% for the four energy lines with respect to the outer source radius (6.0 cm). For the top of the source (9.0 cm), they ranged between 0.08% and 1.48%, and for the bottom of the source (1.0 cm) they ranged between 40.6% and 66.23%. Thus the results depend on both the axial and radial source locations, and this is a two-dimensional problem. The measured fluxes were simulated using two different methods. The first set was simulated by using the ray-trace code with the same angular quadrature as used by the inverse method transport calculations. Thus this set of measurements was exactly consistent with the iterative identification processes. These measurements will be referred to as ray-trace measurements. The second set of measurements was simulated stochastically using MCNP in order to represent more realistic measurements. These will be referred to as Monte Carlo measurements. The simulated total fluxes are shown in Table 2. For each test problem the convergence criterion was a relative difference of 0.01% between calculated and measured fluxes for
Table 2 Simulated flux measurements (a) for the test cases. Line energy (keV)
4. Test cases A set of numerical test cases will illustrate the Schwinger inverse method in a cylindrical geometry for each type of inverse problem. For the test cases, the geometry of Fig. 1 will again be considered. The geometry contains a 16-kg high-enriched uranium
144 186 766 1001 a b
Monte Carlob (%)
Ray-trace 2
5.5317 10 5.7794 103 5.2982 100 1.9521 101
Expressed in c/cm2 s. 1r relative standard deviations are given.
5.38 102 ± 8.51 5.78 103 ± 6.15 5.34 100 ± 5.34 1.99 101 ± 3.93
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
each energy line. In the event that this could not be met, four other stopping criteria were implemented: 1. A stop due to continued violation of constraints, as outlined in Section 3. 2. A test to see if the iterative results are diverging instead of converging toward the solution. To test this, the Fortran code checks the value of the root-mean-square difference between the calculated and measured fluxes over the previous ten iterations. If this difference increases every time the iteration number increases, the problem appears to be diverging and the calculations are stopped (this test is only administered after ten iterations have been completed). 3. The fluxes are changing very little with each iteration. To quantify this, the method checks the last ten values of the rootmean-square difference between calculated and measured fluxes and determines its slope over these iterations. The value of this slope is then divided by the average root-mean-square difference between calculated and measured fluxes over the ten iterations. If the resulting relative slope is less than a prescribed value (0.01), the method stops. 4. The same values for all the unknowns are calculated twice (i.e. in any two iterations). The script looks at all input files for repeats. In this case ‘‘same” means that the first seven digits are equal. Once the method stops, the set of calculated unknowns for the iteration with the lowest v2 difference between the calculated and measured fluxes is considered the solution. This v2 difference is given by
v2 ¼
G X
g
M
g¼1
rg0
2 Mg0
;
ð11Þ
where rg0 is the measurement uncertainty for line g. The set of calculated values for the unknowns corresponding to the iteration with the lowest v2 difference is called the best model. 4.1. Unknown interface locations (see Table 3) For Case 1 suppose the location of all the axial interfaces are known but the location of the radial interfaces at 2.0 cm and 6.0 cm are unknown. The outer source radius is known to be greater than 2.0 cm and is inside the outer edge of the geometry at 8.0 cm, but its location with respect to the aluminum radius at 7.0 cm is unknown. Here the outer source radius was postulated to be 7.2 cm (a situation similar to Fig. 2a). The inner source radius at 2.0 cm was guessed to be 3.0 cm. The Schwinger input file in this case indicated that the second and third radii were permitted to cross. The crossing did occur, and with ray-trace measurements convergence occurred in seven iterations. The outer source radius was calculated very accurately to be 6.0001 cm, but the inner radius was calculated to be a somewhat inaccurate 2.3691 cm. The reason for this is that the source is so thick that the measured flux is insensitive to the inner radius. With MCNP measurements, the outer source radius was calculated accurately at 6.054 cm, but the inner radius had a value of 5.0168 cm. Although this is far from the actual location of the inner radius, the agreement between calculated and measured fluxes was very good, with a v2 difference of 0.2371. Thus this model belongs to the family of solutions of the inverse problem. The method took 14 iterations to converge with MCNP measurements. In Case 2 suppose that a radial interface and axial interface are unknown. The outer source radius at 6.0 cm was unknown and guessed to be 6.3 cm, while the height at 1.0 cm was unknown and guessed to be 0.8 cm. Using ray-trace measurements, very
accurate values of 6.002 cm and 1.0023 cm were calculated in 16 iterations. With Monte Carlo measurements, the unknown height was calculated with high accuracy to be 1.044 cm. The unknown radius was calculated to be 6.218 cm. These values correspond to a v2 difference of 0.06075. Thus once again the method converged to a member of the family of solutions. 4.2. Unknown shield material In Case 3 the bottom aluminum portion of the shield was unknown and postulated to be nickel. The convergence of the Schwinger inverse method for this case can be seen in Fig. 3. When ray-trace measurements were used the method converged to nearly the exact cross sections of aluminum. Aluminum was clearly identified as the unknown material, as shown in Table 4, where the RMS difference for aluminum (0.00017%) was three orders of magnitude smaller than its nearest competitor, wet ground (0.501%). With Monte Carlo measurements, the 186-keV line converged to approximately 4% of the aluminum cross section, but the 1001-keV line converged to only within about 8.2%. The 144keV and 766-keV lines (not shown in Fig. 3) converged to within 5% and 3%, respectively, of the aluminum cross sections. As a result of these differences, aluminum was one of six materials identified as likely candidates for the unknown. Aluminum had the lowest RMS difference at 1.150%, just better than the nearest competitor, wet ground (1.228%). The other likely materials were beryllium oxide (2.785%), ordinary concrete (2.963%), Portland concrete (2.968%), and dry ground (7.511%). 4.3. Unknown source weight fractions Case 4 featured unknown weight fractions for the source material, with both guessed at 0.50000. With ray-trace measurements, the correct weight fractions of 0.9473 for 235U and 0.0527 for 238 U were calculated in a single iteration. Using Monte Carlo measurements, weight fractions of 0.946545 for 235U and 0.053455 for 238 U were calculated in three iterations. This model corresponds to a v2 difference of 0.1349. The iterative history for the case of Monte Carlo measurements is shown in Table 5. 4.4. Unknown material densities In Case 5 the densities of the source and aluminum layers were unknown. The source was postulated to have density 13.00 g/cm3
Percent Difference from True Cross Section
970
350% 186 keV; Ray-Trace
300%
1001 keV; Ray-Trace 186 keV; Monte Carlo
250%
1001 keV; Monte Carlo
200%
3.00% 0.00% -3.00%
150%
-6.00%
100%
-9.00% 3
4
5
50% 0% -50% 0
1
2
3
4
5
Iteration Fig. 3. Convergence of the Schwinger inverse method for Case 3. The unknown material was guessed to be nickel in iteration 0; it was actually aluminum. On the large graph the cross sections all appear to converge to 0%; differences can be seen in the inset.
971
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973 Table 3 Results of interface identification problems. Descriptor
Model Outer radii (cm)
Case 1
Actual Model Initial Converged – RT Converged – MC
Void
Source
Aluminum
Nickel
2.00 3.00 2.3691 5.0168
6.00 7.20 6.0001 6.054
7.00 7.00 7.00 7.00
8.00 8.00 8.00 8.00
Aluminum
Source
Void
Source
Aluminum
1.00 1.00 1.00 1.00
3.00 3.00 3.00 3.00
7.00 7.00 7.00 7.00
9.00 9.00 9.00 9.00
10.00 10.00 10.00 10.00
Void
Source
Aluminum
Nickel
2.00 2.00 2.00 2.00 Heights (cm)
6.00 6.30 6.002 6.218
7.00 7.00 7.00 7.00
8.00 8.00 8.00 8.00
Aluminum
Source
Void
Source
Aluminum
1.00 0.8 1.0023 1.044
3.00 3.00 3.00 3.00
7.00 7.00 7.00 7.00
9.00 9.00 9.00 9.00
10.00 10.00 10.00 10.00
Heights (cm)
Actual Model Initial Converged – RT Converged – MC
Outer Radii (cm)
Case 2
Actual Model Initial Converged – RT Converged – MC
Actual Model Initial Converged – RT Converged – MC Italics represent known (fixed) quantities.
Table 4 Materials with the lowest RMS difference for Case 3. Density (g/cm3)
RMS difference (%)
Ray-trace measurements 1 Aluminum 2 Wet ground 3 Portland concrete 4 Ordinary concrete 5 Beryllium oxide
2.70 2.60 2.35 2.35 3.025
0.00017 0.501 2.472 2.492 2.687
Monte Carlo measurements 1 Aluminum 2 Wet ground 3 Beryllium oxide 4 Ordinary concrete 5 Portland concrete 6 Dry ground
2.70 2.60 3.025 2.35 2.35 1.90
1.150 1.228 2.785 2.963 2.968 7.511
Order
Material
and aluminum was postulated to have density 5.00 g/cm3. With ray-trace measurements the source density was calculated to be 18.742 g/cm3 and the aluminum density was calculated to be 2.701 g/cm3. With Monte Carlo measurements the source density was calculated to be much higher than the actual source density, with a value of 25.469 g/cm3, but the aluminum density was calculated to be an accurate 2.758 g/cm3. These values have differences of 35.91% and 2.15%, respectively, to the actual densities. Though the source density is far from the actual density, these
numbers result in good agreement between calculated and measured fluxes, with a v2 difference of 0.1361. The iterative histories for this problem are presented in Table 6. 5. Summary and conclusions The Schwinger inverse method was originally applied to the inverse transport problems of identifying unknown interface locations and unknown source isotope weight fractions in spherical radiation source/shield systems in Favorite (2004). In Favorite and Bledsoe (2006), the method was implemented for shield identification on numerical test cases in one-dimensional spherical geometries and in an analytic monodirectional slab problem. Favorite (2005) presented numerical solutions for unknown mass density problems in one-dimensional spheres. The method was then applied for material identification in a two-dimensional cylinder in Bledsoe et al. (2007), but the geometry considered was unrealistically small and the quantity of interest was the total (4p) leakage from the system. In this paper, the method is applied to interface location identification, shield material identification, source weight fraction identification, and mass density identification in realistic finite twodimensional cylindrical geometries. In addition to the application to more realistic (two-dimensional) problems, this paper extends previous work on the Schwinger inverse method by using surface formulas for unknown interface locations, automatic correction at-
Table 5 Iteration history of source weight fraction identification for Case 4; Monte Carlo measurements. Iteration
235
Difference from actual weight fraction (%)
238
Difference from actual weight fraction (%)
0 1 2 3
0.5000 0.944535 0.946540 0.946545
47.22 0.29 0.08 0.08
0.5000 0.055465 0.053460 0.053455
848.77 5.25 1.44 1.33
U weight fraction
U weight fraction
972
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
Table 6 Iteration history of density identification for Case 5. Difference from actual source density (%)
Aluminum density (g/cm3)
Difference from actual aluminum density (%)
Ray-trace measurements 0 13.000 1 11.553 2 11.043 3 17.606 4 18.742
30.63 38.35 41.07 6.05 0.01
5.000 5.844 3.386 2.354 2.701
85.19 116.44 25.41 12.81 0.04
Monte Carlo measurements 0 13.000 1 31.722 2 31.459 3 26.020 4 25.604 5 25.456 6 25.469
30.63 69.27 67.87 38.85 36.63 35.84 35.91
5.000 0.421 3.269 2.863 2.761 2.758 2.758
85.19 84.41 21.07 6.04 2.26 2.15 2.15
Iteration
Source density (g/cm3)
tempts for violated constraints, and ray-tracing instead of discreteordinates for transport calculations, which enables us to use the scalar flux at a point as the quantity of interest. For the interface location problem, the issue addressed was whether the Schwinger inverse method could successfully identify unknown radii or heights of material interfaces in the interior of the source/shield system. In Favorite (2004), the answer to this question was yes for identification of unknown radii. Here the answer was yes for both radii and heights. When measurements were used that were perfectly consistent with the transport calculations of the iterative method, the Schwinger method converged nearly perfectly to the correct interface locations in both cases considered. When less consistent, more realistic, measurements were used, the Schwinger method calculated one interface very accurately and one inaccurately in both cases. However, in these cases the method converged to a model that belonged to a family of solutions that corresponds to small v2 differences between calculated and measured fluxes. For the shield identification problem, two issues were addressed in the numerical test cases considered here. The first was whether the method would converge to the correct cross section set. When perfectly consistent measurements were used, the Schwinger method converged nearly perfectly to the correct cross sections in all cases. When less consistent, more realistic simulated measurements were used, the method converged to within 15% for each energy line. The second issue considered was whether the calculated cross sections could be used to identify the unknown material. In the spherical problems of Favorite and Bledsoe (2006), the answer was yes. Here the answer was the same. With perfectly consistent measurements the actual material was correctly identified as the only suitable candidate in both test cases. When less consistent measurements were used, the difference between the converged cross sections and the actual cross sections of the unknown material led to some ambiguity. Six possible shield materials (including the correct one) were identified. For the source identification problem the issue addressed was whether the method could identify the correct weight fractions of the isotopes in the photon source. In the spherical test cases of Favorite (2004), the answer was yes. Here the answer is also yes. When perfectly consistent measurements were used the exact weight fractions were identified in one iteration for both test cases. Using less accurate measurements the method was still highly successful, with errors smaller than 2% for the calculated weight fractions. With these measurements, convergence required three iterations. In mass density identification problems the issue addressed was whether the method could identify the correct densities of the
materials in the source and shield. In the spherical test cases of Favorite (2005), the answer was yes. Here the answer was the same. When perfectly consistent measurements were used, the Schwinger method calculated densities to within 0.04% of the actual densities. With more realistic measurements the Schwinger method found a source density that differed from the actual source density by 35.91%. The shield density differed from the actual density by 2.15%. Although the source density differed greatly from the actual density, the final model output matched the measurements to within 1%. In this paper, the Schwinger inverse method was applied to realistic two-dimensional cylindrical inverse transport problems. By using a ray-tracing algorithm, the ray effects inherent in the discrete-ordinates code used in Favorite (2004), Favorite and Bledsoe (2006) and Bledsoe et al. (2007) have been eliminated. Since ray effects become worse as the distance between the source and detector increases, in our previous attempts at solving inverse problems we were forced to either use the total (4p) leakage from the system (which tends to wash out ray effects, but measurements of which are prohibitively difficult) as the quantity of interest or make the geometries unrealistically small. With ray effects eliminated in this paper, we were able to explore much larger geometries and use the sum of the scalar fluxes at a few points external to the geometry as the quantities of interest. There are several potential applications of the Schwinger inverse method. Identification of radioactive source constituents and dimensions is of enormous concern for national security. Identifying and determining the spatial locations of materials is also important in radioactive waste assay. The identification of shield material given a known radioactive source has application in industry, where, for instance, the composition of mined materials or manufactured alloys may be needed. Acknowledgements The authors wish to thank D.I. Ketcheson for contributions to the interface-location problem and the ray-tracing code and D.E. Vaughan for providing the damping algorithm. Appendix A. Identifying material interface locations in twodimensional cylinders and the method for the crossing of radii or heights This appendix describes how the Fortran code identifies the unknown material interface locations and how it can accommodate the crossing of two interface radii or heights. These considerations are necessary for the solution of two-dimensional problems. Consider the example shown in Fig. 1. Suppose first that in the initially
K.C. Bledsoe et al. / Annals of Nuclear Energy 36 (2009) 966–973
postulated system the source radius is equal to the radius of the aluminum layers. This example illustrates that with two-dimensional methods it is not sufficient to simply state in the Schwinger input file that an entire radius is unknown (as was the case with one-dimensional methods). In this case, the aluminum radii are known to be at 7.0 cm and should not change, but the source radius postulated to be at 7.0 cm will change. If the Schwinger input file simply stated that the 7.0 cm radius was unknown, the Fortran code would wrongly assume that all interfaces at 7.0 cm were unknown. To avoid this problem the Schwinger input file contains an array that has an element for each mesh. If the outer radius of the mesh is unknown (or, for height identification, if the top of the mesh is unknown) the element of the array corresponding to that mesh is assigned a value of 1. If the outer radius is known, the element is assigned a value of 0. In this way the Fortran code is able to discern exactly where the unknown boundary is. This array will be referred to as the array of unknown interfaces. For the unknown source radius problem considered here, the array of unknown interfaces is initially given by
0 0
0
0 1 0 0 1 0 0 1 0 0 0
0
Suppose the source radius postulated to be equal to the radius of the aluminum layers is calculated by the Schwinger method to have a different location. The result is a source radius moved to the inside or outside of the radius of the aluminum layers. A fourth mesh must now be created because the aluminum and source radii are no longer equal (if the source radius is moved inward, the situation is similar to Fig. 1, if it is moved outward, the situation is similar to Fig. 2a. With the additional mesh, the array of unknown interfaces must be updated to the form (assuming the source radius is moved inward)
0 0
0 0
0 1 0 0 0 1 0 0 0 1 0 0 0 0
0 0
The Fortran code was written to automatically handle this, and the iterative process continues until the unknown is identified. Suppose now that the source is postulated to have a radius greater than that of the aluminum layers. It is thus necessary for the radii of the aluminum layers and source layer to cross if the unknown is to be identified. In one dimension, the crossing of any two radii is not allowed in the Schwinger inverse method
973
(Favorite, 2004) but for interface identification in two dimensions the user may specify two interfaces that may cross (this specification occurs in the Schwinger input file). For the geometry shown in Fig. 1, the Schwinger input file would indicate that the second and third radii may cross. Note that when this crossing occurs, the array of unknown interfaces changes from
0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 to
0 0
0 0
0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 The Fortran code was written to handle this crossing and reassign the array of unknowns. References Alcouffe, R.E., Baker, R.S., Dahl, J.A., Turner, S.A., Ward, R.C., 2005. PARTISN: A TimeDependent, Parallel Neutral Particle Transport Code System. LA-UR-05-3925, Los Alamos National Laboratory. Bledsoe, K.C., Favorite, J.A., Aldemir, T., 2007. Material identification in finite cylindrical geometries using the Schwinger inverse method. Trans. Am. Nucl. Soc. 96, 545–547. Favorite, J.A., 2004. Using the Schwinger variational functional for the solution of inverse transport problems. Nucl. Sci. Eng. 146, 51–70. Favorite, J.A., 2005. Identification of unknown densities in a source/shield system using the Schwinger inverse method. Trans. Am. Nucl. Soc. 93, 588–589. Favorite, J.A., Bledsoe, K.C., 2006. Identification of an unknown material in a radiation shield using the Schwinger inverse method. Nucl. Sci. Eng. 152, 106– 117. Favorite, J.A., Bledsoe, K.C., Ketcheson, D.I., 2009. Surface and volume integrals of uncollided adjoint fluxes and forward-adjoint flux products. Nucl. Sci. Eng., in press. Lewis, E.E., Miller, W.F., 1984. Computational Methods of Neutron Transport. John Wiley and Sons (Chapter 4). Malenfant, R.E., 1981. Los Alamos Critical Assemblies Facility, LA-8762-MC(UC-46). Los Alamos Scientific Laboratory, Los Alamos, New Mexico. X-5 Monte Carlo Team, 2003. MCNP-A General Monte Carlo N-Particle Transport Code, Version 5. LA-UR-03-1987, vol. I. Los Alamos National Laboratory. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1994. Numerical Recipes in FORTRAN: The Art of Scientific Computing, second ed. (reprinted with corrections), (Chapter 2). Cambridge University Press. Rahnema, F., 1996. On the internal interface perturbation. Ann. Nucl. Energy 23, 1401. White, M.C., 2003. Photoatomic Data Library MCPLIB04: A New Photoatomic Library Based on Data from ENDF/B-VI Release 8. LA-UR-03-1019 (Revised), Los Alamos National Laboratory (February).