IBD-2L : two line overrelaxation code for diffusion theory calculations

TÜRKİYE ATOM ENERJİSİ KURUMU

ÇEKMECE NÜKLEER ARAŞTIRMA VE EĞİTİM MERKEZİ

Ç .N .A .E .M . A .R - 284

IBD-2L : TWO LINE OVERRELAXATION CODE FOR DIFFUSION THEORY CALCULATIONS

Ulvi ADALIOĞLU

Nuclear Engineering Department

October 1990

SUMMARY

IBD-2L : TWO LINE OVERRELA£ATION CODE FOR DIFFUSION THEORY CALCULATIONS

One line block o v e r r el ax at i on method has been succ e ss fu ll y appl ie d to IBD code to accelerate the s olution of inner i te r at i ­ on p r ob l em of two d i m en s i o na l ,f e w group diff u si on equations r e ­ cently. By using same alg or i th m given for one line technique, a two line block o v e r r el a xa ti o n scheme was used for the inner ite- rat i o n .

The new code, IBD-2L, tested on VAX/11-750 machine proved to be faster than IBJ5-1L. The gain in CPU time is about 12 to 40 %, d e pe nding the c o m plexity of problem. It seems that the larger the size of mesh structure the larger the gain in c om pu ta ti on time for the same convergence criterion.

Ö Z E T

JBD-2L : DtFÜZYON HESAPLARI İÇİN İKİ SATIR OVERRELAKSASYON KODU

Bir kaç yıl önce iki boyutlu, birkaç gruplu difüzyon denklem lerinin çözümünü hız l a n dı rm a k için bir satır blok o v e r r e 1aksas- yon metodu IBD koduna başarıyla tatbik edilmişdi. Aynı a l g o r i t ­ mayı k u l la n ar a k bu sefer iç iterasyon için iki satır overrelak- zasyonu uygulanmıştır.

V A X / 1 1-750 makinesinde test edilen yeni kod, IBD-2L bir s a ­ tır k o dundan daha hızlıdır. P ro blemin k ar ma şı k lı ğı n a göre CPU zamanındaki kazanç % 12 ilâ 40 arasındadır. Ayni y akınsama k r i ­ teri için kafes yapısı büyüdükçe zaman kazancı da artıyor g ö r ü n ­ me k. t e d i r .

TABİLE OF CONTENTS

P a g e s

1. INTRODUCTION 1

2 * THEORY 1

2.1- Matrix formulation of Diffusion Equation 1 S

2.2- Partitioning of matrices A s 3

2.3- Algorithm and program 7

3. APPLICATIONS 8

4. RESULTS AND CONCLUSIONS 8

REFERENCES 10

Appendix 1- Input Data Deck

TABLES

Table 1- Some Characteristic Results obtained

by Two Codes 11

FIGURES

Figure 1- Mesh Structure for Difference Equation 2 Figure 2- Grouping and Ordering of Mesh Points 3

for Two Line Partitioning

1. INTRODUCTION

The block successive overrelaxation techniques have been w i ­ dely used to accelerate the solution of inner iteration problem of two dimensional neutron diffusion equations in past 30 years in nuclear field (1).

The simplest case of these techniques is one line o ve r r e ­ laxation and it was used to solve and accelerate the solution of a 2D diffusion code developed in ÇNAEM N u c l . Engineering Dep.

(2,3) .

The next case is the two line block overrelaxation to be used for the inner iteration acceleration. This report describes theory and the application of this technique, and the results o b t a i n e d .

2. THEORY

2.1- Matrix Formulation of Diffusion Equation

Two dimensional, multigroup neutron diffusion equations to be solved are G g g g \ 7 g h g - V . (D V<i> ) + [E + / E 3 <t> = a fa-jy s i (2.1) S = 1,2, >G

where G is the maximum number of energy groups, and the other symbols have meanings wellknown.

The in Fig. finite d

F q . (2.1) is integrated over the mesh structure shown 2.1 and simplified in order to obtain the corresponding

~

2

-Figure 1- Mesh Structure for difference Equation

The five point finite difference equations obtained are

B

B

e

e

s

B

~ Ö

Ş

e

e

,s

s

s

ç

= --- —

S'

+, s

£ S

t

e - 1,2,

*G

s

where S s and S s are the total fisyon and group f i ,j si , j

scattering sources at point (i , j )

,

- 3

-g S g

A 9 = S (2.3)

g

if the mesh points are numbered succesively. A s are MNxMN d i ­ mensional three diagonal matrices with two non-zero off-diagonal

g g

entry lines. (j) and S s are MNxl dimensional solution and s ou r ­ ce vectors, respectively.

2 .2- Partitioning of matrices A g s

The mesh points of successive two lines can be grouped, as shown in Fig. 2*2, in order to obtain a particular partitioning

g

of matrices A s. The ordering of mesh points is shown in the fi­ gure . 1

1

1 i

11

----Figure 2- Grouping and ordering of mesh points for two

lire partitionin g

Notice that for a two line mesh point grouping the number of rows, M should be an even number.

Using this grouping and ordering one can easily obtain a new form of E q . (2.3), i.e. ,dropping group index g

-4 r 1 ~ B E S 1 1 1 1 0 D B E

_♦

s

s

i

( 2 . 4 )

where the submatrices B s, D s, and E s are of the form of

i i i B = j a 1 j -e 1 j — c 1 j -d 2 j a 2j 0 -c 2 j -b 3 j 0 a 3 j -e 3 j -b 2N-2 j -b G 2N-1 j 2Î- -2 j j a -e j 2N-1 j 2‘- - l j j -b -d a 2N j

-

5

We see that matrix group equations (2.4) have dimensions which are half of the corresponding equations used in one line o v e r r e l a x a t i o n . On the contrary, submatrices E s or D s have

i i

dimensions twice of those of corresponding matrices (see R e f . 3). All couplings within one grouping of mesh points, that is

-6-for a particular j of row groups are expressed with the o ff -d i ­ agonal entries of B s. The coupling with other rows, namely j-1

j

or j+1 are given by the matrices E s and D s.

j J

Different ordering» other than the one chosen here, of mesh points produces matrix equations exactly same form of E q .(2.4). But submatrices becomes quite different and give rise to an inc­ rease of number of arithmetic operations, and hence CPU time. Therefore one must be very careful for the choice of ordering in order not to get that kind of submatrices.

The successive overrelaxation technique applied to the E q . (2,,4) gives the following set of equations:

g , s(m) g ,g(m) • g g ( m - l ) g B _♦ - D <b E © + S i i i i-1 i i + 1 i ,g(m) r g(m) g(m-l) , g(m <|> _{U) <P} + 4> i b L i İ i i = l , 2 , ...- , M / 2

where m shows the iteration n u m b e r .

E q s . (2.5) can be converted into the familiar form of the one line overrelaxation technique :

(m) Y (m ) = - N Y N (m - 1 ) Y

i , i i i,i-1 i-1 _ij, i + 1 i + 1

(m ) _{p (m) (m - 1 ) - (m~ 1)} Y = i U) j Y - Y b L i i + Y i

for each energy group. Here

= N

1,0 M/2,M/2+l

- 7 -B is factored as such i T B = S S i i,i i,i

and R is the diagonal matrix whose diagonal elements are i , i

the diagorial entries of S . And also, as definitions i > i -1 -1 N = ■ R B R i , i i , i i i , i N = -1 R D -1 R i , i-1 _i i i - 1 , i-1 N = -1 R E -1 R i , i + 1 i , i i i + 1 , i+1

The factorization of matrices N is very easy and given as i,i

a product of

? two known matrices, that is

N

-1 = ( R

T -1

S ) . ( S R ) (2.8)

i i >i i,i i,i i,i

The new solut ion vector (j) is again:

, (m) -1

i (m)

<!> = R Y (2.9)

i _i i

2.3- Algori thm and program

-8-re with some modifications.

The subroutines COZ2 and KATSA are modified for the new m a ­ trix formulation. The prog r am itself is also mod if i ed a cc or d in g to new dimensional requirements.

Input is changed in order to give axial b uc kl in g besides the reactor physical height. (See Appendix 1.)

The inner iteration process utilizes same c on vergence c r i t e ­ ria and steps as defined in Ref. 3 for IBD-1L code, namely flux c on v er g en ce c r iterion is £ . It becomes smaller than £ after

4 3

few iterations. Its initial values are taken as 0.02 or 0.01 and final values are 0.0002 or 0,0001 for sample problems in this r e p o r t .

3. A P PL I CA T I O N S

This new code IBD-2L has been tested with several problems on V A X /1 1- 7 50 machine. The results were compared with those o b ­ tained by IBD-1L code.

i- The first problem is the TR-2 small core conf i gu ra t io n shown in Figure 2.3. The core is reflected from one side by B e ­ ri Ilium blocks, and the other side has two A lu minium blocks lo- catod at each corner. The control rod BC-2 is fulljr inserted while three other rods are out of the core for this test run.

A four group, 24 x 30 mesh point structure in X-Y g eo metry m o d e l l i n g has been used for c riticality calculations. 53 regions are defined for 10 m a terials used.

i i - The second problem is the first problem but with smaller mesh intervals or larger mesh size.. The new mesh point numbers are 48x58, R e g i o n and material numbers are again 53 and 1.0, r e s ­ pect. i v e l y .

4. R E S UL T S AND C O N C L US I ON S

Table 1 summarizes the results obtained. It seems that the a greement is quite s a t is f a c t or y with respect to m ul ti pl i ca ti o n factors, They diverges each other at fifth digit. On the c o nt r ar y the average group fluxes change at third digit.

-9“ Y Water reflector Be ||Be !Be 1 Be 1!B S l ; jj BS2 B C 1 I 1 1 I lB C 2 l A1 _{i !} ! A 1 Be : Be blocks A1 : A1 blocks BC : Control rod BS : Shim rod

Blanks are standart fuel elements

X

Figure 3- TR-2 Small Core Configuration

For large size problems average group fluxes differ at e a r­ lier than third digit.

Comparing CPU times we see that the two line overrelaxation has an ddvantage over the one line scheme. The gains in CPU times are about 12 to 40 % for these examples.

The iterates of overrelaxation factor, M is obtained by b

solving homogeneous equation which is E q . (2.4) without sourc term. Cone calculate* spectral radius of Gauss-Seidel matrix o homogeneous equation in each iteration. After certain number o iterations the true «pectral radius can be obtained and optimum overrelaxation facto) are calculated. Hence, one does not need to solve homogeneous equation after that point of iteration process. Therefore one reduces number of floating point calculations. This f* ;»ture of IBD-2L code saved about o m minute of CPU time for the urge size problem n o .2.

Hcnc IBD-2L is quite satisfactory aru faster code IBD-1L i >r diffusion theory calculations with reasonable integral and point r< suits.

than good

-10-REFERENCES

1. L. A. Hageman, "Numerical Methods and Techniques used in the Two-Dimensional Neutron Diffusion Program, P D Q - 5 " , WAPD-TM-364, Feb. 1963.

2. U. Adalıoğlu, R. T u n ç e l , "IBD-1: Two Dimensional, M u lt i ­ group Neutron Diffusion C o d e ” , ÇNAEM A.R. No. 233, Oct. 1985.

3. U. Adalıoğlu, R. Tunçel, "IBD-1L : One Line O v er re la xa ­ tion Code for Diffusion Theory Calculations (An improved version of IBD-1 C o d e ) ” , ÇNAEM A .R .-263, July 1989.

1 a bl e 1-So me Ch ar ac te ri st ic Re su lt s ob ta in ed by Tw o Co de s 11

APPENDIX 1- INPUT DATA DECK 1 . card : I R R X 1,I R R X 2,I R H Z 1,IRHZ2 415 2 . card : E P S 2,E P S 3,E P S 4,NRR 3 F 1 0 .7,12 3 . card : I B , İ C A P , IGEO,IAD,IGUC,IBAK,ICEY, K M A X,I M A X,J M A X,K OMSA.IBIS 1213 4 . card : ( I S 1 ( I ) , I S 2 ( I ) , J S 1 ( I ) , J S 2 ( I ), 1 = 1 , IB 1814 ( 4 +I B ). card : KK(K) i I M A T (K ) , K = 1,IB 1814 ( 4 + 2 J B ). c a r d : DELRI(I), 1 = 1 , IMAX 8 F 1 0.5 < . . .> c a r d : DELZJ(I), 1=1,JMAX 8 F 1 0.5 card : X (K ), K =1,KMAX 8 F 1 0.5 ( • • •i c a r d : X N U (K.N ),D F (K,N ),S A (K,N ),S F (K,N ) : SS(K,M,N), M = 1,KMAX N=l,KOMSA, K = 1,KMAX 8 F 1 0.5 < . . .> card : P E l 4.7 ( . . . ) c a r d : HYUK E l 4.7 ( . . .* c r d : BKARE E l 4.7

Inf ( m a t i o n abou some of the input parameters :

IB : number f regions

I GEO : 1 for ; 1 indrical geometry 0 for cartesian geometry

 2

-ICÂP = 1 + ICEY = 1 half core (Y axis s y m .) in cartesian geometry, or whole core in cy l i n d r i ­ cal geometry

İCAP = 1 + ICEY = 1 half core (X axis sym.) in cartesian geometry İCAP = 1 IAD I GUO : I BAR : KMAX : I MAX : JMAX : ROMSA : IBIS : IS1(I) : IS 2 ( I ) : J S 1 (I ) : J S 2 ( I ) : TRPX1 : IRPX2 : IRHZ1 : IRHZ2 : KK(K) : I MAT('K) : D E L H I (I ): D E L Z J (T ): X(K) :

+ ICEY = 1 quarter core in cartesian geometry, or half core in cylindrical geometry 0 for flux and adjoint flux calculation

> 2 for flux calculation only

1 for power and source distribution calculation 1 B 2 z is used for X-Y geometry

0 for cylindrical geometry case number of groups

number of mesh points in X or r axis number of mesh points in Y or Z axis number of compositions

0 for relative flux calculation

1 for relative and normalised (power) flux left side boundary of region I

right side boundary of region I bottom side boundary of region I top side boundary of region I left boundary of core region right boundary of core region bottom boundary of core region top boundary of core region mesh number of region K

composition number in region K

mesh intervals in region I on X or r axis mesh intervals in region I on Y or Z axis f rac t ion of fission neutrons in group K

A 3

-XNU(K,N): for group K end composition N

D F (K ,N ) : diffusion coeffient for group K and composition N SA(K,N) : E « for group K and composition N

S F (K ,N ) : Ei for group K and composition N

S S (K ,M ,N ): E» from group

K

P : reactor power in watts for IGUC -1,or IBIS=1 cases , otherwise no need for P

HYUK : height of reactor core in Cm for X-Y geometry

-2