Improvements and validation of the linear surface characteristics scheme

Annals of Nuclear Energy 36 (2009) 46–59

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Improvements and validation of the linear surface characteristics scheme S. Santandrea *, J.C. Jaboulay, P. Bellier, F. Fevotte, H. Golfier Direction de l’Energie Nucléaire, Département de Modélisation des Systèmes et Structures, Service des Etudes de Réacteurs et de Mathématiques Appliquées, CEA, DEN, DM2S-SERMA, F-91191 Gif Sur Yvette, France

a r t i c l e

i n f o

Article history: Received 20 June 2008 Received in revised form 24 October 2008 Accepted 27 October 2008 Available online 16 December 2008

a b s t r a c t In this paper we present the last improvements of the recently proposed linear surface (LS) characteristics scheme for unstructured meshes. First we introduce a new numerical tracking technique, specifically adapted to the LS method, which tailors transverse integration weights to take into account the geometrical discontinuities that appear along the pipe affected to every trajectory in classical characteristics schemes. Another development allows using the volumetric flux variation of the LS method to re-compute step-wise constant fluxes to be used in other parts of a computational scheme. This permits to take greater advantage of the higher precision of the LS method without necessarily conceiving specialized theories for all the modular functionalities of a spectral code such as APOLLO2. Moreover we present a multi-level domain decomposition method for solving the synthetic acceleration operator that is used to accelerate the free iterations for the LS method. We discuss all these new developments by illustrating some benchmarks results obtained with the LS method. This is done by detailed comparisons with Monte-Carlo calculations. In particular we show that the new method can be used not only as a reference tool, but also inside a suitable industrial calculation scheme. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Recently a new transport scheme for the characteristics method in unstructured meshes has been proposed (Santandrea and Mosca, 2007) and implemented in the current version of APOLLO2 code system. (Sanchez et al., 1988) The new method allows linear variations for fluxes and sources inside the computational meshes. This is based on a linear interpolation of fluxes over surface values. In Section 2 we briefly review the scheme. In this work we introduce many improvements to this method. First of all we adapt the tracking technique, which is at the base of every method of characteristics (MOC), to face the local geometrical discontinuities of the computational meshes. This problem was firstly considered in (Fevotte et al., 2007). The analysis carried out in that work showed that when tracking over an arbitrary geometry many deformations appear due to the different numerical tracking approximations. These deformations are of the order of the transverse integration step, and can only be reduced by a sufficiently fine tracking. In (Fevotte et al., 2007) a technique was proposed which consisted in local projections of all discontinuities, in a way that avoids propagating them onto all the geometry. Unfortunately this technique demands strong modifications of the stan-

* Corresponding author. Tel.: +33 (0) 1 69 08 81 78; fax: +33 (0) 1 69 08 94 90. E-mail addresses: [email protected] (S. Santandrea), [email protected] (J.C. Jaboulay), [email protected] (P. Bellier), [email protected] (F. Fevotte), herve.golfi[email protected] (H. Golfier). 0306-4549/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2008.10.007

dard sweeping algorithm, i.e. the algorithm that sweep tracking during transport. When comparing the computational overhead necessary to implement the method with the advantages it gives in terms of reduced transverse integration step, we find that, for the moment, the approach is not computationally efficient. In this work we propose a different approach. We adapt locally the transverse integration weights in order to take into account the possibility of only a limited number of geometrical discontinuities in the neighborhood of each intersection point of the trajectory with a computational region boundary. This simplified approach lacks the strong convergence property of the original algorithm but leaves unchanged the transport sweep. This proves to be more computationally efficient than previous strategies. In Section 3 we detail the problem and illustrate our improvement. Another important issue is how coupling the LS results with the rest of a code system such as APOLLO2 where the basic numerical representation for fluxes is the constant one. We introduce, in Section 4, an interpolation operator that allows to use the larger meshes of typical LS calculations and to couple them with the finer meshes typically used for nuclide evolution or self-shielding calculations. In Section 5 we show in realistic benchmark assembly calculations the results given by the above methods, and we compare them with reference Monte-Carlo calculations as well as with the standard step characteristics (SC) scheme. The aim of this paragraph is to show that reference benchmark calculations can be done with the LS method even in the delicate context of nuclide

47

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

evolution and evolutionary self-shielding calculations, where a finer mesh is necessary. In fact we show that a suitable operator can correctly interpolate the values of LS fluxes into the finer meshes without losing important physical information. A last important improvement is described into Sections 6 and 7, where adapted multi-level algorithms are developed to take full advantage of the cell-type configurations that are now routinely used for LS calculations. In these sections we also show on realistic calculations that their impact on computational time can be very important. Conclusions will follow in the last section. 2. The LS scheme In this paragraph we briefly recall the LS method, which has already been presented in full details in (Santandrea et al., 2008). This method uses numerical integrated surface quantities and a volumetric linear integration technique. In Sections 3 and 4 we will show how the computation of these quantities can be strongly improved with respect to the previous works. To introduce the LS scheme, we first discuss the interpolation technique at the center of our methodology. Before all we will introduce the quantities that our scheme is supposed to reproduce or conserve. We will give analytical or semi-analytical definitions without giving the explicit numerical counterparts. We will then show how the numerical interpolation works and discuss some of its difficulties. Finally we will substitute the interpolation scheme into the transport and balance equations to obtain a final numerical system that allows us to find the basic coefficients of the interpolation itself. We will suppose here that the system D we are considering is constituted by a set of homogeneous regions we denote by {Di}. We then want to interpolate a generic function / 2 C 0 ðDÞ, / 2 C 1 ðDi Þ 8i 2 f1; . . . ; ng. That is to say a continuous function on the entire domain and whose derivatives are also continuous inside each homogeneous region. The interpolation uses the following four integral quantities:

1 ~ /a ¼ 4 p Sa

Z

AðXÞ dX~

ð4pÞ

~ /surf ;i Þ ¼ /i;G ¼ Mi ð~

Z

dr s wðr s ; XÞ;

a2i

1 M a;i ~ /a ¼ 4pV

Z ð4pÞ

XZ a2i

dS? l Sa

a2i

where a is an index referring to surfaces while ‘‘i” is for regions, and A is the vector of real spherical harmonics (Sanchez and Chetaine, 2000). In Eqs. (1)–(3) ds? ¼ dr s jX  nj and ‘‘l” is the length of the chord intersected over region ‘‘i” by the trajectory passing on the point rs of the surface Sa and having direction X. Moreover the notation a e i (2) and (4) indicates that the sum has to be done over all surfaces that belong to the border of the region ‘‘i”. In Eq. (4), which is a simple analytical form of the balance equation, one can recognize the partial moment currents (~ J  ), the volume of the region (Vi), the emission source average (Qi) and the total cross-section (Ri). Finally t(a, rs, X) is the opposite surface (with respect to a) appearing on the trajectory passing through a in the point rs with direction X. The quantities in Eqs. (1)–(4) are, respectively, the flux surface average, the flux volume average in hypothesis of linear variation of the flux inside the region, the linear discontinuous moment fluxes and the flux volume average. In (2) we have used the notation

~ /a ; a 2 ig; /surf ;i ¼ f~

ð4pÞ

Z

ð4pÞ

Z

dr~ /ðrÞ

Vi

dX

XZ a2i

dS? Sa

Z

dx~ /ðr s þ xXÞ:

ð6Þ

½0;l

In the LS interpolation the flux moments inside a volume are linearly interpolated over each trajectory t crossing cell i in direction X making use of the quantities (3) so that we write

l  x ~ x ~ /ðxÞ ¼ ~ / ð0Þ þ / ðlÞ ; l l

ð7Þ

Mð~ /surf ;i Þ ¼ ~ /i;G ¼ ~ /i;C

;

ð2Þ ~ /a þ ~ /i;C  M i ð~ /a Þ; /a ¼ ~ ( )  X  ~þ ~ ~ ~ Ja  J a þ V i Q i =ðRi V i Þ; /i;C ¼

1 4pV

dX

ð4Þ

~ /tða;rs ;XÞ /a þ ~ 2

¼

Z

ð3Þ

ð1Þ dX

1 ~ /i ¼ 4pV

where x is the distance along the trajectory and l the chord length in region i. In formula (7) u(0), u(l) are the local surface values in the step-wise surface representation given by (3). The interpolation (7) is done via the ‘‘parametric” presence of X that is not a variable the flux moments depend on. This presence will turn out to be useful when a flux moment appears inside the angular integral in characteristics neutron equation. Even if expression (7) does not appear to be completely orthodox, and in fact flux values do not necessarily coincide in an arbitrary point for different angles X, we think it is an acceptable numerical approximation. Substituting (7) into (6) gives (2). Note also that (2) implicitly defines the average operator M. The suffixes ‘‘G” and ‘‘C” used into (2) and (4) denote two different volume averages. The first is the geometrical average of (2), the second is the conservative value in (4). The principal difficulty of the LS is to maintain the coherence between the geometrical average value (2) and the conservative value (4). It is important to note that Eqs. (2) and (3) are simply interpolative equations, that is to say they are completely uncorrelated to the physical meaning of the neutron flux w. The physical quantities that comes from the numerical model are given by Eqs. (1) and (4). Now consider the quantity Mð~ /a Þ, that is the geometrical average of the flux computed with the continuous flux moments and not with the discontinuous moments as it is done in (2). This quantity should be a good candidate for the geometrical average, but in fact it is not equal to the ~ /i;C except for very refined meshes. (Santandrea and Mosca, 2007). To impose this constraint we need a further degree of freedom that is in fact the discontinuity introduced in (3). In fact one can easily verify that whatever the value of ~ /C;i is

Sa

X

to indicate the discontinuous moment fluxes defined in (3) that belong to the internal part of region i. A similar notation is used in (3) for the continuous moments of (1). Formula (2) asks for more comment. In fact an analytical definition of the average flux could be written as

ð5Þ

ð8Þ

Formula (8) assures us that the conservative and the geometrical averages give the same value, in fact the one computed by the conservation Eq. (4). In conclusion our LS scheme interpolates flux moments, via the parametric presence of the angular variable by using the linear Lagrange-type formula (7). The interpolation is realized over the integral quantities (1)–(4). Since the two volume averages (2) and (4) are coherent between themselves, we will call our method conservative. In the following we will show how to compute the quantities w(rs, X) and ~ /Ci appearing in Eqs. (1) and (4). This will allow us to compute all the other quantities of system (1)–(3). The LS scheme is based on a linear interpolation between cell boundary values of the collision source in the transmission equation. To perform this integration one needs the value of the sources within each region that is expressed by

~s; Þ~ ~ qa ¼ diagðR /a þ ~ Sa ; i

ð9Þ

where this source is computed for each cell side because both crosssections and fluxes have been chosen to be discontinuous over cell

48

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

sides and are given by (3). For a trajectory t crossing cell i in direction X the emission density is calculated, coherently to (7), as

related to trajectory t and Sa is the area of surface a that is defined by

 x x qt ðx; XÞ ¼ qt ðl; XÞ þ qt ð0; XÞ 1  : l l

Sa ¼

ð10Þ

To derive from this the numerical transmission equation for the LS we first recall the integral transport equation, which we write as sðxin ;xout Þ

wout ðxout ; XÞ ¼ e

win ðxin ; XÞ þ

Z

xout

sðx0 ;xout Þ

e

0

0

qðx ; XÞdx ;

ð11Þ

xin

where x indicates the position along the trajectory, xin (xout) is the entering (exiting) point of the trajectory in the domain and s(x0 , xout) is the optical distance along the trajectory between point xout and x0 . Replacing (10) into (11) gives ðnÞ

ðnÞ

ðn1Þ

^t;i ð0; XÞ wout;i ðtÞ ¼ eRi Ri ðtÞ win;i ðtÞ þ bq   ðn1Þ ^t;i ðl; XÞ  q ^ðn1Þ ð0; XÞ ; þ b1 q t;i

ð12Þ

where the following notation is used:

^ ¼ q=R; q

b1 ¼ 1 

1  eRi Ri ; Ri R i

b ¼ 1  eRi Ri ¼ ð1  b1 ÞRi Ri ð13Þ

where Ri = xout  xin. It is important to note that (13) permits to tabulate only the function b1 and to deduce b and the exponential term appearing in (12) with simple operations. This is what is done in the transport sweep. Moreover we have written (12), and we will do the same in the following, with the iteration index ‘‘n”. We couple Eq. (12) with a numerical version of the balance Eq. (4) that we write as ðnÞ Ri V i ~ /i;C ¼

    X 1 ðn1Þ ~s;i ~ ð~ Ja  ~ J a ÞðnÞ þ V i diag R qext /i;G þ ~ i;G ; 2k þ 1 a2i

where ‘‘k” stands for the anisotropy order. Eq. (14) is an analytical balance equation where we approximate numerically some of the quantities involved into. The first of these numerical approximations is to use the following definitions of partial currents

Z

dSjX  n jwðnÞ ðr; XÞ 

a

X

x? ðtÞwðnÞ ;i ðtÞ;

ð15Þ

tkX;t\a

where the ‘‘+” and ‘‘” stand for outgoing or ingoing angular directions. The second approximation is to substitute the expression (2) for geometrical average fluxes, while surface moments are obtained by using (3). Eq. (14) is a truly balance equation because ~ /i;G is equal to the conservative value ~ /i;C as we have shown in (8). Eq. (14), even if not numerically coherent with the propagation equation, gives a conservative average value coherent with the computed currents. This matter is further discussed in (Santandrea et al., 2008). Ultimately the value of the emission density q within a region depends on the inhomogeneous sources, if any, and on the multigroup fluxes via fission and scattering. In the LS scheme we have used a cell linear source approximation obtained by linear interpolation from cell surface values (a surface being the common boundary between two adjacent regions; but it can also be defined by further subdivisions to improve convergence). In practice we compute and store the surface average flux moments,

Z Z 1 dS dXAk ðXÞwðr; XÞ 4pSa a ð4pÞ 1 X x ðtÞ  xn ? Ak ðXn Þwta ðXn Þ; Sa x 2S jX  nj n N |fflffl{zfflffl}

Z

dX

Z

ð4pÞ

a

dS 

X xn 2SN

xn

X x? ðtÞ : jX  nj t\Sa |fflffl{zfflffl} tkX

ð17Þ

dSt

Moreover, in (16) and (17) we have made use of an SN = {xn, Xn}n=1,N angular quadrature formula to replace exact integration with finite sums. Formula (16) shows that all integrals are numerically evaluated using the angular fluxes along numerical trajectories. These fluxes are obtained by applying formula (16) during the transport sweep, and they are used to update the surface source (9). In conclusion, the LS algorithm proceeds as follows: 1. Start with a guess surface source. 2. Make a transport sweep, calculate for each surface the moments of the flux with (16). Compute the current term in balance Eq. (14), and update volume fluxes. 3. Test convergence on volume and surface fluxes and eventually return to 1. The important case of cyclic trajectories can also be treated with specialized formulas (Santandrea et al., 2008). We end the section with an important observation. Forcing conservation makes the scheme not positive, since the correction in (3) can lead to negative surface and volume average values. This last eventuality is only limited to cases where the spatial mesh discretization is quite crude (Santandrea et al., 2008). 3. An improved tracking technique for the LS scheme

ð14Þ

J a ðXÞðnÞ ¼

1 4p

In this paragraph we discuss an important improvement of the tracking technique used for the LS method. This improvement is directly inspired by the analysis developed in (Fevotte et al., 2007). In fact one of the basic findings of that work was that much of the error due to transverse integration in characteristics methods is due to geometrical discontinuities. By ‘‘geometrical discontinuities” we mean a point which is either the extremity of a surface element, or the point of tangency between a surface element and a characteristic line in the direction of interest.

Ω

ω⊥

/ka ¼

t\Sa

ð16Þ

dSt

for each surface a and angular direction X. Here Ak is the spherical harmonic of order ‘‘k” dSt is the measure of the cross-sectional area

x⊥ Fig. 1. Sketch of the classical tracking scheme. We represent here a set of trajectories parallel to direction X intersecting an arbitrary region. We affect to each trajectory a pipe that is supposed to be centered on it and which has a x\ width. In the figure the pipes are represented by dotted lines.

49

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

The situation can be better discussed looking at Figs. 1 and 2. Fig. 1 represents an arbitrary region intersected by a uniformly spaced set of trajectories parallel to direction X. In Fig. 2 we can look at the local weights dSt of formula (16). In the classical tracking scheme we have dSt ¼ x? ðtÞ= cosðaÞ where x\(t) is the transversal integration step and a is the angle between the trajectory and the outgoing normal to the surface. The situation is described in Fig. 2. These weights should characterize the local surfaces intersected by each imaginary pipe related to every trajectory. Unfortunately this representation can be quite approximated as illustrated in Fig. 2. There are in particular two points where problems are clearly visible. Next to point 1 the end point of the region is not respected, while next to point 2 the circular surface is roughly approximated in the neighbor of the tangent limit. In order to strongly improve the tracking technique we propose to consider a different definition of the weight dSt. The strategy is to adapt this weight to the local form of the intersected equation. Let us suppose to describe the position vector ~ rðsÞ over the trajectory as

~ rðsÞ ¼ ~ a þ s~ e;

ð18Þ

where s is the local coordinate over the line, while ~ a and ~ e are, respectively, a point on the trajectory and the unit direction vector of the line. Note that the width between trajectories is given by x\. We can now distinguish two cases, according to the type of the intersected equation.

where D = eyfx + exfy is the determinant of the linear system constituted by Eqs. (18) and (19), we affect the following value to the intersection weight

dSt ¼ fg½ð1  t in ÞjDj þ g½t in jDjg= cosðaÞ;

ð20Þ

where a is the angle between the normal vector to the segment and ~ e and the function g is defined as

gðyÞ ¼



x? =2 if y > x? y otherwise

ð21Þ

:

The normal isffi taken on the side such that a < p/2, therefore qffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðaÞ ¼ fx2 þ fy2 =jDj. Formula (20) is obtained by the following reasoning. When a trajectory intersects an element we separately consider the two parts of the element that are divided by that intersection. The situation is depicted in Fig. 3, where only one trajectory intersected the element a. This situation is related to the second branch of function g in formula (21). When another previous or successive trajectory intersects the element, the local weight is chosen to be x\/(2*cos(a)), i.e. the local surface related to the half pipe. This situation corresponds to the first branch of (21). 3.2. Circles or arc of circles To treat the case of circles we take into consideration the local change of frame depicted in Fig. 4. We treat before the case of full circles. We write the equation of the circle as

3.1. Segments The parametric equation of a segment can be written as

~ rðtÞ ¼ ~ c þ t~ f;

tin ¼ ½ex ð~ c ~ aÞy  ey ð~ c ~ aÞx =D;

ð19Þ

where ~ c is the position vector of the first of the two limits of the segment, while ~ f is the displacement vector going from the first to the second limit, and t is supposed to be a local coordinate belonging to [0,1]. Given the intersection coordinate between the segment and the trajectory, whose value is

~ rðtÞ ¼ ~ c þ R~ eðhÞ;

ð22Þ

where ~ eðhÞ ¼ ðcos h; sin hÞ, R is the radius of the circle and ~ c is the center. We note that the reference frame used here, and described by Fig. 4, is locally defined in the sense it has the center into the ~ c point and that its axis are given by x\ and X. Two intersections are possible in this case, that are given by the following angles:

#1 ¼ 2p  cos1

x  loc ; R

#2 ¼ 2 p  # 1 ;

ð23Þ

where xloc is the local coordinate of the impact point of the trajectory on the x\ axis as described by Fig. 4, and is given xloc ¼ c ~ aÞx ey . Moreover all angles are measured with reð~ c ~ aÞy ex  ð~ spect to the x\ axis in the classical counter clock wise sense. We then define the following angles:

Ω

1

t-1

α

Ω

t+1

t

n 2

δSt

ω⊥

α

δSt=ω⊥/cos(α)

x⊥ Fig. 2. Drawbacks of the classical tracking scheme. This figure represents the transverse surfaces weight dSt of formula (16) in the case of the classical tracking scheme. This figure highlights many of the major drawbacks of the classical tracking technique. For example near point 1 one clearly see a geometrical paradox: the original end point of the region is not respected by the numerically computed surfaces, which intersect one each other and end in the following pipe limit. Another important deformation is next to point 2.

tin|Δ|

1-tin|]Δ|

Fig. 3. Illustrative sketch for formula (20). The element a is intersected only by trajectory t. The sum of two half side weights g(tinD) and g([1  tin] D) equals the total surface length.

50

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

θ2

Ω

t

ϑ2

x⊥ z2

θ1 z1 xloc

ϑ1 δSt

Fig. 4. Basic geometrical quantities involved into the integration of an arc of a circle. The arc of circle is the part of the circle comprised between the angles h1 and h2. The local reference frame is defined having the center into the circle center, and the two axis as being X and x\. An arbitrary trajectory t intersects orthogonally x\ at the local coordinate xloc, while the intersection points are determined by angles #1 and #2 . The imaginary pipe related to the trajectory t also intersects somehow the arc of circle. For example at the entering point the intersection of the pipe determines the two angles z1 and z2, that delimits the dSt surface over the circle.

 hx þ x i ? loc z1 ¼ cos1 min ;1 ; R  hx  x i ? loc z2 ¼ cos1 max ; 1 ; R

ð24Þ

and take for the tracking transverse weight in the case of circles the following expression:

dSt ¼ Rðz2  z1 Þ:

ð25Þ

As shown in Fig. 4 this weight defines the local arc over the circle that is intersected by the pipe related to the trajectory. As for arcs the situation is complicated by the presence of the extremes. All formulas before remains valid, but now h 2 ½h1 ; h2  where h1 ; h2 are the limits angles of the arc of circle and h1 2 ½0; 2p, which we suppose with no generality loss. Let first take into account the transverse weight related to the first intersection (23), once we suppose that it belongs to the arc. We define before all r1 ¼ maxðh1 ; pÞ; r2 ¼ minðh2 ; 2pÞ and then fix

z1 ¼ z2 ¼



¼ r1

if xloc  x? < R cosðr 1 Þ

ðAÞ;

¼ 2p  cos1 ½ðxloc  x? =2Þ=R otherwise ðBÞ; ðCÞ; ¼ r 2 if xloc þ x? > R cosðr 2 Þ

ð26Þ

¼ 2p  cos1 ½ðxloc þ x? =2Þ=R otherwise ðDÞ:

The logic behind (26) is simple and can be retrieved by the reader by considering the signification of the different cases. Case (A) is when the previous trajectory t0 (such that x\ (t0 ) < x\ (t)) does not intersect the arc, while in (B) this event occurs. Similarly in case (C) the following trajectory does not intersect the arc, while this intersection is realized in (D). Similar formulas can be written for the second eventual intersection of (23). Two final comments are worthwhile. The first one is that with the definitions (20) and (25) the second sum in formula (17) directly normalizes as

Snum ¼ a

X X xn 2SN t\Sa tkXn

dSt ¼ Sa ;

ð27Þ

that is to say the integration is exact. Formula (27) is subject to only one condition, that at least one trajectory in tracking intersects surface a for the all the angles Xn. Since the fulfillment of this simple condition cannot always be granted we will have to consider that Snum a ð–Sa Þ is a numerical value, and the normalization (16) still remains necessary. Nevertheless, the capacity of satisfying (27) in almost all cases is an important property that gives our scheme a greater stability. Unfortunately it must be noticed that our locally adapted tracking technique palliates but does not solve the problem of taking into account geometrical discontinuities as it was done in (Fevotte et al., 2007). However, in that work this improvement was done at the price of strong modifications of the standard sweeping algorithm, so that the global computational overhead was too high to justify its use with respect to the standard free iterations. Our approach, on the contrary, does not change free iterations. All eventual precision improvement is obtained without any computational penalization. Secondly, we must point out that our treatment of discontinuities is not only approximated, in the sense we have specified before, but also incomplete in the sense that the transverse integration appears not only in formula (16) for the surface flux, but also in expression (15) for the current. In our treatment we only deal with the first but avoid to take into account the second. The careful reader will observe that applying the same approach to the currents, that is to say adapting the integration weight x\ of formula (15), is problematic. Actually we may induce ‘‘numerical” balance errors, because the weight for an entering current can be different to the weight of the exiting one for any intersection. Since we want to preserve conservation, we prefer to postpone to future works the possibility to take into consideration geometrical discontinuities also for the current term. Nevertheless it must be pointed out that the higher convergence property of the LS are related to the linear interpolation over surface fluxes, while the volume averages (to which currents are related) are only used to make a convenient per region linear rescaling. In fact, as it is showed in (Santandrea and Mosca, 2007) the balance rescaling reduce the observed numerical error, without changing the numerical observed rate of convergence. Therefore we can suppose that the integration of the surface fluxes hides the major part of error.

4. Coupling geometries with the linear surface interpolation In an integrated environment such as APOLLO2 (Sanchez et al., 1988), each solver has to be able to interact with the data coming from other modules. In the LS scheme we deal with a linear surface flux representation and a step-wise representation for cross-sections. In this section we show how we couple geometries where quantities are supposed to be step-wise constant, with geometries that adopt a linear surface representation for fluxes. Two steps are to be considered. The first one is the projection of the step-wise cross-sections of the fine mesh geometry to give the step-wise cross-sections in the coarse mesh geometry. The second step is to retrieve the constant fluxes in the fine mesh geometry from the coarse one. As for the first step, we have adopted in this work a quite crude assumption, by assuming that cross-sections still remain constant over the coarse mesh geometry, so that they are simply obtained by the classical formula

P

i2I RI ¼ P

Ri V i

iV i

;

ð28Þ

where the two sums are done over all ‘‘micro” regions belonging to ‘‘macro” region I. Formula (28) completely defines the projected cross-sections. This allows performing a linear surface calculation with the method described in Section 2. It is important to understand that assumption (28) does not prevent a detailed nuclide

51

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

evolution in each region of the fuel pins, coherently with the geometry described at the self-shielding level. But it asks however that the errors done by computing fluxes with the homogenised crosssection (28) can be neglected. This seemingly crude hypothesis reveals to be quite an acceptable one in many cases, even if we have found one severe case where it is in pain. We will comment more this point in the following section. The second step consists in retrieving fluxes and reaction rates over the original mesh. This can be done by computing the average flux in the fine mesh by supposing to have a linear variation over the coarse mesh. This last interpolation is given by the computed flux. This assumption leads us to write fine fluxes as

/i ¼

X

Nia; /a; ;

ð29Þ

/i ¼

1 4p V i

     1 la  lb 1 lb  la dX dSjX  njli /a   þ /b : 2 2 2lI 2lI ð4pÞ @i

Z

Z

ð31Þ After rearranging, expression (31) gives exactly (29), where coefficients N are given by

Nia; ¼

1 4pV i

Z

dS

a

  2la þ li : dXjX  njli 1  2lI ð2pþ Þ

Z

ð32Þ

A direct translation of (32) in a numerical form of the type of those of (15) and (16) is straightforward and we leave it to the reader. The microscopic fluxes are the final result that can be directly coupled with other methods.

a2I

which shows that we obtain the average flux of each ‘‘micro” region contained in I as a linear combination of inward surface fluxes of region I. To compute the coefficients N we need to define some important geometrical quantities that are represented in Fig. 5. In this Fig. 1 can see an arbitrary trajectory intersecting a macro region I and a micro region i (contained into I). This trajectory enters I through surface a and exits through surface b and determines three chords, la li and lb that are, respectively, the distance between the entrance of macro region I and the intersection of region i, the chord intersected over region i, and the distance between the exit point of region i and the exit point of region I. We will call also lI = la + li + lb. þ Moreover / i and /i indicate the entering and exiting fluxes of region i over trajectory t. Similarly /a and /b are the surface fluxes at the entrance and at the exit of macro region I over trajectory t. In a representation where fluxes have a linear surface representation with respect to macro surface values we can write

  la la /i ¼ /a 1  þ /b ; lI lI   lb lb þ /i ¼ /a þ /b 1  : lI lI

ð30Þ

We can then compute the average flux in region i coherently with (2), where we replace expressions (30). This in details leads to

5. Results and comparisons on physical test cases In this section we validate the LS characteristics scheme improvements with the help of numerical PWR assembly benchmarks. Two standard PWR lattice configurations (UO2 and MOX) are considered, whose rough descriptions are given in Fig. 6. The physical conditions of these two configurations are absolutely representative of realistic assemblies. The first UO2 lattice represents a classical PWR assembly with typical functioning temperatures. Also the display of water holes, indicated in blue in Fig. 6, is standard and is designed to produce a flat power distribution. The second assembly is a more difficult one, since the different plutonium enrichments create the conditions for much higher flux heterogeneities. For both assemblies the last cell column has a slightly enlarged pitch. The computational models presented in this paper are twofold: on one hand we propose the CEA point-wise Monte-Carlo code TRIPOLI-4 (Diop et al., 2007) and, on the other hand, the CEA reactor physics code APOLLO2. The validation process is based on the comparison of APOLLO2 calculations versus continuous-energy Monte-Carlo TRIPOLI-4 reference calculations. All the geometrical models used to make TRIPOLI-4 or APOLLO2 calculations have been produced by using the versatile tools offered by SILENE (Stankovski, 2008). All calculations use the same nuclear data

lβ li Ω

I

t

φβ

β

lα

φi+

φi−

i

φα α Fig. 5. Interpolation of a ‘‘microscopic” region ‘i’ starting from the values of the ‘‘macroscopic” surface fluxes of region I.

52

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

UO2 lattice

__

235

MOX lattice

U 4,9%

__ Pu 9,8% __ Pu 6,5% __ Pu 3,6%

Water hole

Cell pitch : 1,26 cm

Cell pitch : 1,26 cm

Fuel radius : 0,407 cm

Fuel radius : 0,415 cm

clad radius : 0,476 cm

clad radius : 0,475 cm

TFuel : 552˚C

TFuel : 650˚C

Tclad : 400˚C

Tclad : 335˚C

Tmod : 311˚C

Tmod : 306˚C

Fig. 6. PWR assemblies used in benchmarks.

JEFF3.1 (Koning et al., 2006). APOLLO2 calculations are based on the SHEM 281 group structure (Hfaiedh and Santamarina, 2005) and the resonant mixture self-shielding treatment (Coste et al., 2005). The resonances overlapping phenomena are taken into account from 2.2 Kev up to 22.5 eV for the UO2 lattice and from 7.5 KeV up to 22.5 eV for the MOX lattice. This treatment proves to be essential to obtain a cross-section library able to grant a benchmark quality result. Another important remark concerns the collision anisotropy. In all the following calculations we have used a P1 cross-section library. This anisotropy gives the same overall precision on the pin-by-pin absorption or production reaction rates than a P0 transport corrected cross-section library. But the behavior of the total reactivity is different. In fact P0 cross-section shows much compensation effects between different isotopes. These compensations are much more reduced for the P1 cross-sections. Paradoxically the error on reactivity increases with P1 cross-sections. We consider that for a benchmark purpose it is better to present a P1 calculation even if the total discrepancies on reactivity are higher. The 2D neutron flux calculation is performed using the LS and SC characteristics schemes to evaluate their absolute and relative accuracy and performance. Using the SC scheme, two different spatial discretizations are defined for the UO2 and MOX lattice as illustrated in Fig. 7. The MOX lattice spatial discretization employs

sectors in fuel to represent a nonazimuthally symmetrical flux (this does not matter significantly in UO2). All calculations showed here have been done on Intel Xeon 2.6 GHz PCs. As far as LS characteristics scheme is concerned, UO2 and MOX lattice are described with the same discretization but with fewer meshes (in fact some meshes in the moderator have been removed). Two descriptions are considered: the first with four rings in fuel and the second with two rings in the fuel. Since the LS scheme has a linear flux variation inside rings we can reduce their number in fuel rods. We will keep the original four ring geometry for the self-shielding and depletion models which use flat fluxes. The coupling tool presented in Section 4 will permit to link selfshielding and depletion (four rings model), with the LS calculation (two rings model). The four rings model permits to validate the coupling tool. The main TRIPOLI-4 results are given in Table 1; the first two lines give the Kinf obtained, with their co respective statistical uncertainties. We note that these uncertainties are very low (less than 20 pcm). The two following lines give the root mean square of the standard deviations in pin-by-pin absorption and fission reaction rates distributions for both lattice types. The last two lines gives the maximum standard deviation of these same pin-by-pin reaction distributions. We can note that we dispose here of very well converged Monte-Carlo results, and that the following

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

53

spatial discretization with step-wise characteristics scheme MOX lattice (1578 regions)

UO2 lattice (801 regions)

spatial discretization with Linear-Surface characteristics scheme UO2 and MOX lattice, 2 rings (180 regions)

UO2 and MOX lattice, 4 rings (258 regions)

Fig. 7. UO2 and MOX lattice spatial discretization.

comparisons with deterministic calculations could be pushed well beyond the percent. We provide in Table 2 a first batch of comparisons. Here we evaluate the impact of the new tracking technique presented in Section 3 as compared to the old tracking technique in the LS scheme, these errors are computed by taking as reference the TRIPOLI-4 results.

Table 1 TRIPOLI-4 reference results. For the two benchmark PWR assembly cases the computed reference Kinf are affected by very low statistical uncertainties (less than 20 pcm). Also the error on pin-by-pin integrated fission and absorption reaction rates are very low, by far less than a percent. The statistical uncertainties are shown for the reactivity, for the root mean square error (RMS) and for the absolute maximal errors in pin-by-pin integrated reaction rates.

kinf (T4) r (pcm)

UO2

MOX

1.31247 17

1.12257 16

RMS

r(uRa) r(umRf)

0.07% 0.08%

0.09% 0.10%

MAX

r(uRa) r(umRf)

0.09% 0.11%

0.11% 0.14%

Both the two rings and the four rings models are compared. All calculations are made with a quadrature angular formula consisting of a product of 24 uniform azimutal angles and two polar Bickley– Naylor angles. (Sanchez et al., 2002) The transverse integration is made with a uniform step of 0.1 cm. Infinite multiplication factor, pin-by-pin absorption and production rates distributions are compared. We recall that since there are no differences in the sweeping algorithm, the computational times are the same for the calculation made with the old and the new tracking. In fact the only differences in computational times are related to the tracking phase itself, which is increased of 20–30%, but still remains negligible with respect to the rest (less than 0.3 s in any case). A last remark is necessary which concerns LS results. To improve results we can adopt an automatic surface division process that allows having a better representation of surface fluxes. In fact, since the surface flux representation is constant per surface, we have to reduce surface lengths in order to have a better spatial flux variation. We do not detail here the full study of this effect but we point out that all the LS calculations we show in the following have been done by re-sizing all surfaces in such a way that their length is next to 0.6 cm for the UO2 lattice and 0.4 cm for the MOX lattice. On the base of our experience this prescription gives good results at a reasonable cost.

54

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

Table 2 Effectiveness of the new tracking technique and impact of a reduced mesh. In this table we compare APOLLO2 results obtained with the old and the new tracking and with the two or four rings geometry with respect to reference TRIPOLI-4 calculations. The comparisons are done for the reactivity, for the root mean square error (RMS) and for the absolute maximal errors in pin-by-pin integrated reaction rates. The reduction of average and absolute maximal error on pin-wise integrated absorption and fissions reaction rates is particularly apparent for the MOX assembly. Thanks to the higher order convergence property of the LS method, we can reduce the number of rings in every pin without introducing any perceptible error. Geometry

Two rings

Four rings

UO2

MOX

UO2

MOX

Tracking

Old

New

Old

New

Old

New

Old

New

Kinf (AP2) Dqpcm

1.30993 148

1.30980 155

1.12067 150

1.12046 167

1.30965 164

1.30955 170

1.12040 172

1.12018 190

RMS

DuRa DumRf

0.23% 0.29%

0.18% 0.24%

0.32% 0.42%

0.23% 0.27%

0.23% 0.29%

0.18% 0.24%

0.31% 0.41%

0.22% 0.27%

MAX

DuRa DumRf

0.54% 0.63%

0.48% 0.57%

1.24% 1.62%

0.73% 0.94%

0.54% 0.62%

0.48% 0.56%

1.19% 1.61%

0.71% 0.97%

Solver CPU time (s)

36

54

We note that the new tracking technique improves pin-by-pin reaction rates distribution and reduces maximal discrepancies for rates distributions. This is particularly apparent for the absolute maximal errors made in the MOX assembly: here the new tracking reduces the errors by almost a factor two. We can also evaluate the impact of reducing in the LS scheme the mesh size in fuel pin from four rings to two rings. In Table 2 we show that there are slight differences in reaction rates distributions, no significant impact in reactivity and CPU time is reduced by 24% in average. Comparisons between the LS and SC characteristics schemes are presented in Table 3. LS calculations are performed with the new tracking technique and with only two rings in the fuel rod (corresponding to 180 regions) while SC calculations use the geometries shown in Fig. 7. The impact of the solver on reactivity is limited to 30 pcm and reaction rates distributions are very similar. LS characteristics scheme reduces CPU time from 33% up to 54%. It is important to note that very important gains in computational times are obtained with the new LS scheme. We now look at depletion calculation. In this context we fully apply the interpolation technique of the Section 4. Actually in what follows we will apply for both the SC and LS calculation the same depletion and evolutionary self-shielding meshes. But for the LS scheme we adopt, for the MOC calculation, the two rings geometries described above, while the connection with the rest of the computational scheme is done by interpolation. The SC scheme is considered as the reference depletion calculation. The Fig. 8 shows the discrepancies in reactivity between LS and SC characteristics schemes. LS calculations use the four and two rings geometry. Concerning the UO2 lattice, there are slight differ-

47

71

ences (max: 70 pcm for the four rings geometry and 50 pcm for the two rings geometry) and CPU time is reduced by a factor 2. As for the MOX lattice, no significant discrepancies are noticed for the four rings geometry (around 50 pcm between SC and LS) but an important deviation is shown for the LS two rings geometry (around 190 pcm at 70 GWd/t). To understand these differences we have given a detailed look at reaction rates in fuel rods. They have shown important variations in internal rings (theses variations are less important with UO2). The constant cross-section variation we impose in our scheme for the internal region (see formula (28)) is not sufficient to describe these variations. To reduce this deviation, the internal ring must be described with more regions (four rings geometry). Finally in Table 4 we show how minimal, maximal and root mean square errors change during evolution for the four and the two rings geometry at 0, 40 and 80 GWd/t that are, respectively, the beginning, the half and the end of fuel life. This table shows that the two rings geometry errors degrade very strongly at the end of depletion for the MOX case. This is of course a result of error accumulations starting from the very beginning. As for the UO2 assembly our interpolation technique gives very good results all along the evolution. To improve our approach we envisage in future works to use a polynomial function to describe the neutron cross-section in internal ring. In conclusion, the LS scheme and the new developments presented in previous sections allow a reduction of CPU time by a factor 2 with the same accuracy obtained with SC characteristics. New developments (polynomial flux in internal rings) are underway to reduce the discrepancies observed for MOX lattice with the LS two rings geometry. 6. Advances in the DPN synthetic acceleration

Table 3 Comparisons between the two ring geometry LS calculation with the four ring geometry SC calculation. TRIPOLI-4 Monte-Carlo calculations are taken as reference. The comparisons are done for the reactivity, for the root mean square error (RMS) and for the absolute maximal errors in pin-by-pin integrated reaction rates. Results show that the two calculations have similar precisions, but the LS is notably faster. Characteristics

UO2

MOX

SC

LS

SC

LS

kinf (AP2) Dqpcm

1.31013 136

1.30980 155

1.12085 136

1.12046 167

RMS

DuR a Dum Rf

0.16% 0.22%

0.18% 0.24%

0.20% 0.25%

0.23% 0.27%

MAX

DuR a Dum Rf

0.39% 0.46%

0.48% 0.57%

0.71% 0.90%

0.73% 0.94%

54

36

118

54

Solver CPU time (s)

Every numerical transport method that is based on free iterations must necessarily accelerate these iterations to be able to solve realistic problems.(Adams and Larsen, 2002) For the MOC in APOLLO2 a synthetic method has been proposed and used since many years (Sanchez and Chetaine, 2000; Santandrea and Sanchez, 2002), we call it the DPN method. The basic assumption of the DPN method is to suppose that the angular flux in a given surface b can be approximated as

wb ðr; XÞ  ~ AðXÞdðX; bþ Þ~ wb þ þ ~ AðXÞdðX; b Þ~ wb ;

ð33Þ

where the following notation has been used

dðX; a Þ ¼



¼ 1 if current a has the same direction as X; ¼ 0 otherwise:

55

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

UO2 lattice

MOX lattice

0

0 -20

-10

disc repanc ies in reac tivity (pc m)

discrepancies in reactivity (pcm)

-40 -20 -30

-40

-50

LS (2 rings) / SC -60

LS (4 rings) / SC

-70

-80

0

10000

20000

30000

40000 50000 BU (MWd/t)

60000

70000

80000

Solver CPU time: SC: 1H 26’ LS: 40’ (4 rings) LS: 33’ (2 rings)

-60 -80 -100 -120 -140 -160

LS (2 rings) / SC

-180

LS (4 rings) / SC

-200

0

10000

20000

30000

40000 50000 BU (MWd/t)

60000

70000

80000

Solver CPU time: SC: 1H 56’ LS: 55’ (4 rings) LS: 41’ (2 rings)

Fig. 8. Comparison between SC and LS during depletion calculation.

Table 4 Comparisons during depletion between LS calculation with two ring and four ring geometries and SC four ring geometry. SC calculations are taken as reference. The comparisons are done for the root mean square error (rms) and for the maximal and minimal errors in pin-by-pin integrated production rates. Results show that the two calculations have similar precisions with the UO2 assembly, but discrepancies are higher with the MOX lattice and the LS two ring geometry. UO2

0 GWd/t

DmRfu

MOX DmRfu

40 GWd/t Two rings

Four rings

Two rings

Four rings

Two rings

rms max min

0.08% 0.18% 0.19%

0.10% 0.22% 0.20%

0.05% 0.12% 0.11%

0.06% 0.15% 0.13%

0.04% 0.09% 0.08%

0.04% 0.10% 0.08%

rms max min

0.13% 0.25% 0.31%

0.16% 0.36% 0.39%

0.07% 0.14% 0.17%

0.12% 0.20% 0.35%

0.04% 0.07% 0.10%

0.40% 0.36% 0.80%

The assumption (33) permits to express the average flux moments on the surface a as

Z Z Z Z 1 1 A ~ A~ waþ þ A~ A~ wa dS dX~ dS dX~ 4pSa a 4pSa a ð2pþÞ ð2pÞ ¼ Aa;þ ~ ð34Þ wa;þ þ Aa; ~ wa; ¼ Aa ~ wa :

~ /a ¼

Formula (34) also defines a DPN geometrical averaging operator A. The DPN problem conserves the current moments and couples them with the integral transport Eq. (11) to give

~ J þa ¼ ¼

Z

dS Sa

X

Z

^ aþ ~ dXjX  nj~ AðXÞwðr; XÞ ¼ A waþ

ð2pþÞ

^b Þ ðT ab ~ wb þ Eab q

ð35Þ

b2i

where the escape (E) and transmission (T) coefficients are factors that are defined coherently with the transport transmission (12) and q is the source vector (9). Eq. (35) is finally coupled with (14) and this closes also the final synthetic problem. Many algebraic manipulations are necessary to give the DPN method a form adapted to a fast solution (Santandrea et al., 2008). The final form of the DPN equations can be cast as

~ w þ Sext : wþ ¼ T ~

80 GWd/t

Four rings

ð36Þ

Eq. (36) is written for each computational region and simply gives the outgoing DPN components in function of ingoing ones and of the external sources contributions (collision sources have already been eliminated). Before looking into the detail of what we propose here, it is essential to stress an important point. The linear problem (36) has a very specific form and in TDT much work has been done to optimize its solution. Many techniques have been conceived that take full advantage of the response matrix form (Santandrea and Sanchez, 2005). In our following improvement we will conserve all these properties for the final problem. The main of these properties is that the profile of the response matrix problem is only region dependent. We take advantage of this by a suitable region-wise ordering that minimizes memory swap during the fundamental matrix by vector operation. Another fundamental property relies on the fact that an adapted ILU0 decomposition exists that demands almost no supplementary storage, and that has a very reduced computational cost. We shall show that the improvement we propose in this paper preserves the response matrix form, while inducing further optimization. The first result we want to show in this paragraph is that it is possible to take advantage of the very low region coupling of cell-type geometries to reduce the computational effort involved in the domain decomposition technique for this case. It is worth noting that even if what we propose is only related to the

56

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

response-matrix form, and in this sense it is applicable to both LS and SC DPN operators, in practice it is of true interest only for the first one. In fact only for the LS scheme, a cell-type mesh is precise enough to be of practical use. We will call this technique a cell adapted domain decomposition technique (CADDT). For a generic layer region (let say the jth one) in a cell geometry we will write (36) as

wþ;int j

¼

T iij w;int j

þ

T iej w;ext j

þ

Sþ;int j

ð37Þ

;ext wþ;ext ¼ T eij w;int þ T ee þ Sjþ;ext j wj j j

where the label ‘‘ext” or ‘‘int” denotes the external or internal frontier of the region with respect to the centre. The situation is described in Fig. 9. Of course the global system consists of the union of all layers of each cell and of all the cells of the domain. The coupling between each cell being done by the currents located at the interfaces between cells. Moreover one can recognize that the transmission matrix for each layer ‘‘j” has been partitioned according to the distinction between the internal and external currents, so that the final transmission matrix can be written as

Tj ¼

T iij

T iej

T eij

T ee j

! ð38Þ

;

We will try now to apply the principle of domain decomposition to the DPN but in a different way with respect to the general approach of (Santandrea and Sanchez (2002)). The logic of our procedure is to show that the dependence on internal layers can be eliminated in a way that is computational efficient and by introducing a minimum fill in into the original system. Keeping in mind that objective, let us suppose that we have already done the work for the (j  1)th layer so that we can write þ;ext ^þ ~ ¼ T^ j1 ~ wj1 w;ext j1 þ Sj1 ;

ð39Þ

where only the DPN components over the external boundary of layer ‘‘j-1” remain and the hat sign over the T transmission matrix and the source vector S is to denote that they are different from the original quantities in (36). The system of Eq. (39) can be directly coupled with Eq. (37). In fact, paying attention to our notation, we see that ~;int of (37), and analogously the vector ~ wþ;ext j1 of (39) coincides with wj ;ext þ;int ~ ~ wj1 of (39) is the same vector as wj of (37). This means that, with local algebraic operations we can easily eliminated the more internal currents. In fact by defining

Aþjint ¼ ðI  T iij T^ j1 Þ1 ðT iij ^Sþj1 þ Sþjint Þ; |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}

Njint

wþe

we can write the problem for the jth layer in the same form as (39) where the quantities T and S now reads

T^ j ¼ T eij T^ j1 T^ iej þ T ee j ;

ð41Þ

^Sþ ¼ T ei N þ Sþ : jext j jint j1

ψj+,ext -,int

!

wþi

 ¼

T ee

ðn1Þ 

T ei

T ie

ψj-,ext

ψj+,int j-th layer Fig. 9. Description of the variables of the system (37) related to the jth layer of a cell geometry.

we wi

T ii



 þ

Se Si

ðn1Þ :

ð43Þ

and we can compute the transmission matrix and the source for the following nth level as

S ð40Þ

ð42Þ

Eqs. (42) and (39) allow to recover the internal unknowns and so to complete the solution. It is easy to check that the only supplementary storage needed per layer with respect to the original mav trix is that due to the quantity T in j;j1 . This quantity can be as great as half of the original transmission matrix in the limit case when each cell geometry has a very high number of layers. It is important to note that the set of matrix operations involved in (41) are not computationally expensive since we deal here with locally defined quantities. In contrast with a standard elimination procedure, where the amount of computational work grows quadratically with the number of unknowns that are eliminated, here this computational effort grows only linearly. This is due to the fact that the v basic operation we do in our algorithm is to compute T in j;j1 whose 2 cost is proportional to (Slayer) where Slayer is the number of surfaces per layer. Typically this number is limited to a maximum of four, so that the total cost of the algorithm is proportional to 16*Nlayer where Nlayer is the total number of layers in the cell. This is the main advantage of the technique we propose here. The other advantage is due to the fact that the storage needed for the standard algorithm grows quadratically with the number of unknowns, while it grows only linearly for the algorithm we propose. When no cell structure exists in the system, as it is the case for many experimental reactors, the principle of domain decomposition can of course be applied in a multi-level framework. Let us suppose we have computed a transmission matrix at (n  1)th level over N total levels. We next denote as internal those currents of this level that do not belong to the grid of the n-th level, the others being denoted as external. We can then write the transmission problem of this level as

ðn1Þ

T^ iej ¼ ðI  T iij T^ j1 Þ1 T iej

ψj

~ ¼ T^ iej ~ w;int wj;ext þ Aþjint : j

ðn1Þ T ðnÞ ¼ T ee þ T ei

v T in j;j1

¼ T^ j1 Aþjint þ ^Sþj1 ;

The system (41) allows us to eliminate the dependence on the internal DPN components of the ring j. This procedure can be applied iteratively until we reach the sides of the square cell. This gives us a system where only the unknowns related to these surfaces are present. Once this final reduced problem is solved we can come back by using Eq. (39) and the following equation:

ðnÞ

¼

Sðn1Þ e

þ

ðnÞ T ei ðI

ðn1Þ 1 ðn1Þ Þ T ie ; ðn1Þ 1 ðn1Þ T ii Þ Si :

ðI  T ii



ð44Þ

This procedure can be carried on until the last level. Formula (44) gives us the way to compute the transmission matrices of all levels (which is done only once for all), and the sources that are computed iteratively during the calculation. This allows us to arrive to the last level, where the problem is completely solved. Once the DPN components of the last level are known, we have to retrieve back all components of previous levels until we reach the basic one and compute fluxes. For each level this is done through the following equation:

wþi ¼ ðI  T ii Þ1 ½T ie we þ Si ;

ð45Þ

where the distinction between internal and external components is to be evaluated at each level following the criterion we explain for matrix (43) and where the internal component of the source is that appearing in (44). Once all DPN components are available, the flux moments can be computed from the basic level equations. The situation is schematically represented in Fig. 10. We conclude with a final remark. Domain decomposition techniques are general techniques (Saad, 1996) and many numerical

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

57

Fig. 10. Multilevel DPN algorithm. The starting point in the multilevel domain decomposition method is to suppose to have a response matrix problem at, say, the (N  1)th level. By using the unknowns partition that relates the (N  1)th with the Nth level we decompose the transmission matrix as in (43). Thanks to (44) we can define quantities for the N-th level, and go on until we reach the last level where the solution is computed on the unknowns that belong to it. Then thanks to (45) we can come back until the first level is reached, where a balance equation can be used to recover the final result.

libraries can be found that solve linear matrix problems on the base of them. The reason why we do not use them here is that specialized storage and preconditioning techniques exist that are typical of a response matrix problem. To our knowledge it is not possible to find them in standard libraries. We refer the interest reader to (Santandrea and Sanchez, 2005). Here we only briefly re-call some of these aspects. The first of them is what we have called before the region-wise reordering, which we impose to all

the matrices appearing in the domain decomposition algorithm. Another important optimization consists in the building phase of the multi-level method. In our implementation we compute all the matrices necessary to solve a multi-group DPN in a preliminary phase where many energy groups (eventually all) matrices are computed at once. This means that in the initial phase the algorithm is organized in such a way that the multi-group calculation is performed with a full scalar operations, i.e. loops over groups

Fig. 11. Computational mesh for the Osiris reactor.

58

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59

are the most internal ones. At the same time during the building phase also the ordering of the DPN unknowns is not the same of that used in the final group by group iteration phase. It is clearly impossible to ask an independent library this level of optimizations. This can also be said with the respect to the very specific situation and treatment appearing in the CADDT technique.

7. Results We present in this paragraph two realistic numerical examples, and we test separately the two domain decomposition techniques we describe above. Firstly we show the performance of the CADDT over the well known C5G7 benchmark (NEA, 2003). We propose here to look at the results obtained with the computational mesh, the spatial and angular quadratures used in (Santandrea and Mosca, 2007) to test the efficiency of the LS method. We only recall that the mesh consists of 2516 computational regions, that the reactivity error with respect to the reference Monte-Carlo calculation was 5 pcm, while the maximum pin-wise relative error was of 1.08%. It is important to note, as shown in (Santandrea and Mosca, 2007), that the mesh consists of many cells of the type described by Fig. 9, with no sectors. The good performance of the LS method even with no sectors is due its numerical representation capabilities, and has been discussed in (Santandrea et al., 2008). The impact of using the new technique is to reduce the computational time from 34.7 s to 19.5 s on the global time. In fact the impact of the technique is limited to the synthetic operator. When one looks directly at the incidence of this technique on the acceleration operator, one finds that the computational time is reduced by more than a factor 2. A second result, we briefly discuss, is the application of the complete multi-level domain decomposition technique to the case of the Osiris reactor (Lokhov, 2008). Osiris is an open-core light water material testing reactor with mixed H2O–Be reflector. The thermal power of this reactor is 70 MW. It runs with a 19.75% enriched U3Si2AL fuel and can produce up to 4.5  1014 n cm2 s1 fast neutrons (>0.1 MeV) and about 3  1014 n cm2 s1 thermal neutrons (<0.625 eV) in the total core volume. The core operation is performed with six control elements having an absorbing (up) and a fissile (down) parts. The neutron absorber in the control element is Hafnium. In Fig. 11 a sketch of the reactor is visible, where the six Hafnium control elements are clearly distinguishable. The fuel is formed by plates. We will not discuss here in detail the physics of the problem, we only need to know that results compare very well with MonteCarlo calculations and that the geometry does not have any cellular cluster since we deal here with plate fuel assemblies, while the core geometry is very heterogeneous. The calculation is done with a geometry constituted by 40.000 regions and with a 172 group cross-section library that is issued from the self-shielding procedures of APOLLO2. The calculation is done by using the external synthetic acceleration scheme and it converges in five external iterations. When we used no multi-level technique the total computational time is of 26,400 s, while it decreases to 19,500 s with our domain decomposition technique. We have used here a four level grid, which is obtained by the standard Spectral Recursive Bisection algorithm. (Saad, 1996) The coarser grid is an aggregation of about 16 elementary regions, while intermediary grids are chosen to be made by 2, 4 and 8 elementary regions. 8. Conclusions In this paper we have introduced new numerical techniques aimed at improving the precision and the solution speed of the

LS scheme. The first of them palliates the negative effects of geometrical discontinuities treatment of standard tracking approaches. To meet the terms of this task, the transverse integration weights that are usually affected to each trajectory are tailored to the local situation appearing at each intersection. In fact these weights are normally equal to the constant transverse cross-sectional area divided by the cosine of the angle between the normal to the intersected surface and the trajectory angle. This can give a very rough evaluation of the area of the surface intersected by the imaginary pipe related to the trajectory. In our work we give them the exact value of this area. A second development deals with the important issue of the coupling of any given method with any other one in a modular environment such as the one offered by the spectral code APOLLO2. In this framework the accepted convention is that both fluxes and cross-sections are step-wise constant. To fully take profit of the LS scheme we have developed an interpolation operator that can correctly project linear surface fluxes, computed over a coarse mesh geometry, into a fine mesh geometry, adapted to the step-wise assumption necessary to correctly carry out selfshielding or nuclide evolution calculations. These two first developments have then been deeply tested over PWR benchmark assembly calculations, and errors with respect to reference TRIPOLI-4 Monte-Carlo calculations have been shown. Comparisons to test the coupling with evolution and self-shielding calculations have also been carried out. Results show the overall good performances of the new techniques. Only for MOX assemblies we find that during nuclide evolution a not negligible reactivity derive appears in the LS calculation using our projection/reconstruction technique. We will try in the next future to improve this approach by implementing a polynomial (at a first attempt a linear one) interpolation for cross-sections inside the fuel pins. Finally an adapted version of a multi-level domain decomposition technique is presented. In fact we present a general algorithm that only relies on response matrix form of the DPN algorithm and a second version that is fully specialized to the case of cell-type geometries. Results over the C5G7 benchmark and over the Osiris experimental reactor show the very good performances of these techniques. Acknowledgement The authors thank Alexey Lokhov for having provided the Osiris calculations.

References Adams, M.L., Larsen, E., 2002. Fast iterative methods for discrete ordinates particle transport calculations. Prog. Nucl. Energy 40 (1), 3–159. Coste, M., Aggeri, A., Huot, N., 2005. New developments in resonant mixture selfshielding treatment with APOLLO2 code. M&C2005, Avignon (France), September 12–15. Diop, C.M., Petit, O., Dumonteil, E., Hugot, F.X., Lee, Y.K., Mazzolo, A., Trama, J.C., 2007. TRIPOLI-4: a 3D continuous-energy Monte Carlo transport code. PHYTRA1, Marrakech (Maroc), March 14–16. Fevotte, F., Santandrea, S., Sanchez, R., 2007. Advance transverse integration for the method of characteristics. M&C and SNA 2007, April 15–19, Monterey CA. Hfaiedh, N., Santamarina, A., 2005. Determination of the optimized SHEM mesh for neutron transport calculations. M&C2005, Avignon (France), September 12–15. Koning, A., Forrest, R., Kellett, M., Mills, R., Henriksson, H., Rugama, Y., 2006. The JEFF-3.1 Nuclear Data Library. JEFF Report 21, NEA No. 6190, Data Bank, November. Lokhov, A., 2008. Personal communication. NEA, 2003. Benchmark on Deterministic Transport Calculations Without Spatial Homogenization. A 2D/3D MOX Fuel Benchmark. NEA/NSC/DOC(2003)16. Saad, Y., 1996. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, Boston, Massachusetts. Sanchez, R., Mondot, J., Stankovski, Z., Cossic, A., Zmijarevic, I., 1988. APOLLO II: a user-oriented, portable, modular code for multigroup transport assembly calculations. Nucl. Sci. Eng. 100, 352–362.

S. Santandrea et al. / Annals of Nuclear Energy 36 (2009) 46–59 Sanchez, R., Chetaine, A., 2000. A synthetic acceleration for a two dimensional characteristics method in unstructured meshes. Nucl. Sci. Eng. 136, 122–139. Sanchez, R., Mao, L., Santandrea, S., 2002. Treatment of boundary conditions in trajectory based deterministic transport methods. Nucl. Sci. Eng. 140, 23–50. Santandrea, S., Sanchez, R., 2002. Acceleration techniques for the characteristics methods in unstructured meshes. Ann. Nucl. Energy 29, 323–352. Santandrea, S., Sanchez, R., 2005. Analysis and improvements of the DPN acceleration technique for the method of characteristics in unstructured meshes. Ann. Nucl. Energy 32, 163–193.

59

Santandrea, S., Mosca, P., 2007. Linear Surface Characteristics scheme for the neutron transport equation in unstructured geometries. M&C and SNA 2007, April 15–19, Monterey CA. Santandrea, S., Sanchez, R., Mosca, P., 2008. A linear surface characteristics scheme for neutron transport in unstructured meshes. Nucl. Sci. Eng. 160, 22–40. Stankovski, Z., 2008. Mutant Components and Scripting in Silene 2D/3D Pre&Post Processing GUI. Physor 2008, September 14–18, Interlaken (Switzerland).