1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
|
*DECK SNF
SUBROUTINE SNF(KPSYS,IPTRK,IFTRAK,IMPX,NGEFF,NGIND,IDIR,NREG,
1 NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,TITR,NBS,KPSOU1,KPSOU2,
2 FLUXC,EVALRHO)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Solve N-group transport equation for fluxes using the discrete
* ordinates (SN) method.
*
*Copyright:
* Copyright (C) 2007 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* KPSYS pointer to the assembly LCM object (L_PIJ signature). KPSYS is
* an array of directories.
* IPTRK pointer to the tracking (L_TRACK signature).
* IFTRAK not used.
* IMPX print flag (equal to zero for no print).
* NGEFF number of energy groups processed in parallel.
* NGIND energy group indices assign to the NGEFF set.
* IDIR not used.
* NREG total number of regions for which specific values of the
* neutron flux and reactions rates are required.
* NBMIX number of mixtures.
* NUN total number of unknowns in vectors SUNKNO and FUNKNO.
* MAT index-number of the mixture type assigned to each volume.
* VOL volumes.
* KEYFLX position of averaged flux elements in FUNKNO vector.
* SUNKNO input source vector.
* TITR title.
* NBS
* KPSOU1
* KPSOU2
*
*Parameters: input/output
* FUNKNO unknown vector.
* FLUXC flux at the cutoff energy.
* EVALRHO dominance ratio.
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
CHARACTER TITR*72
TYPE(C_PTR) KPSYS(NGEFF),IPTRK,KPSOU1(NGEFF),KPSOU2(NGEFF)
INTEGER NGEFF,NGIND(NGEFF),IFTRAK,IMPX,IDIR,NREG,NBMIX,NUN,
1 MAT(NREG),KEYFLX(NREG),NBS(NGEFF)
REAL VOL(NREG),FUNKNO(NUN,NGEFF),SUNKNO(NUN,NGEFF)
REAL,OPTIONAL :: FLUXC(NREG)
REAL,OPTIONAL :: EVALRHO
*----
* LOCAL VARIABLES
*----
PARAMETER (IUNOUT=6,NSTATE=40)
INTEGER IPAR(NSTATE)
LOGICAL LIVO
DOUBLE PRECISION F1,F2,R1,R2,DMU
*----
* ALLOCATABLE ARRAYS
*----
REAL, ALLOCATABLE, DIMENSION(:) :: FGAR,TEST
REAL, ALLOCATABLE, DIMENSION(:,:) :: OLD1,OLD2,OLD3
LOGICAL, ALLOCATABLE, DIMENSION(:) :: INCONV
*----
* RECOVER SN SPECIFIC PARAMETERS
*----
IF(IMPX.GT.2) THEN
WRITE(IUNOUT,'(//6H SNF: ,A72)') TITR
CALL KDRCPU(TK1)
ENDIF
IF(IDIR.NE.0) CALL XABORT('SNF: EXPECTING IDIR=0')
IF(IFTRAK.NE.0) CALL XABORT('SNF: EXPECTING IFTRAK=0')
CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR)
NSTART=IPAR(20)
MAXIT=IPAR(22)
LIVO=(IPAR(23).EQ.1)
ICL1=IPAR(24)
ICL2=IPAR(25)
IBFP=IPAR(31)
NFOU=IPAR(34)
CALL LCMGET(IPTRK,'EPSI',EPSINR)
IF(IMPX.GT.3) THEN
ALLOCATE(FGAR(NREG))
DO II=1,NGEFF
FGAR(:NREG)=0.0
DO I=1,NREG
IF(KEYFLX(I).NE.0) FGAR(I)=SUNKNO(KEYFLX(I),II)
ENDDO
WRITE(IUNOUT,'(/33H N E U T R O N S O U R C E S (,I5,
1 3H ):)') NGIND(II)
WRITE(IUNOUT,'(1P,6(5X,E15.7))') (FGAR(I),I=1,NREG)
ENDDO
DEALLOCATE(FGAR)
ENDIF
*
IF(NSTART.GT.0) THEN
*----
* GMRES(M) INNER ITERATION LOOP FOR ONE-SPEED TRANSPORT EQUATION
*----
CALL SNGMRE (KPSYS,NGIND,IPTRK,IMPX,NGEFF,NREG,NBMIX,NUN,
1 NSTART,MAXIT,EPSINR,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,NBS,KPSOU1,
2 KPSOU2,FLUXC)
ELSE
*----
* LIVOLANT INNER ITERATION LOOP FOR ONE-SPEED TRANSPORT EQUATION
*----
LNCONV=NGEFF
ALLOCATE(INCONV(NGEFF),OLD1(NUN,NGEFF),OLD2(NUN,NGEFF),
1 TEST(NGEFF),OLD3(NUN,NGEFF))
*
INCONV(:NGEFF)=.TRUE.
TEST(:NGEFF)=0.0
*
OLD1(:NUN,:NGEFF)=0.0
OLD2(:NUN,:NGEFF)=0.0
IF(NFOU.GT.0) OLD3(:NUN,:NGEFF)=0.0
*
ITER=0
10 ITER=ITER+1
IF(ITER.GT.MAXIT) THEN
WRITE(IUNOUT,'(40H SNF: MAXIMUM NUMBER OF ONE-SPEED ITERAT,
1 12HION REACHED.)')
IF(NFOU.GT.0) WRITE(6,310) TRHO1
GO TO 70
ENDIF
*
IF(NFOU.GT.0)THEN
OLD1(:NUN,:NGEFF)= OLD2(:NUN,:NGEFF)
OLD2(:NUN,:NGEFF)= OLD3(:NUN,:NGEFF)
OLD3(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF)
ELSE
OLD1(:NUN,:NGEFF)= OLD2(:NUN,:NGEFF)
OLD2(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF)
ENDIF
*----
* UPDATE THE FIXED SOURCE AND COMPUTE THE FLUX
*----
CALL SNFLUX(KPSYS,INCONV,NGIND,IPTRK,IMPX,NGEFF,NREG,
1 NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,ITER,NBS,KPSOU1,
2 KPSOU2,FLUXC)
*----
* LOOP OVER ENERGY GROUPS
*----
DO 60 II=1,NGEFF
IF(INCONV(II)) THEN
*----
* FOURIER ANALYSIS NUMERICAL EIGENVALUE CALCULATION
*----
IF(NFOU.GT.0)THEN
TRHO1=0.0
IF(ITER.GT.3)THEN
TOT1=0.0
TOT2=0.0
DO IR=1,NREG
IND=KEYFLX(IR)
IF(IND.GT.0)THEN
TOT1=TOT1+(FUNKNO(IND,II)-OLD3(IND,II))**2
TOT2=TOT2+( OLD2(IND,II)-OLD1(IND,II))**2
ENDIF
ENDDO
TOT1 = SQRT(TOT1)
TOT2 = SQRT(TOT2)
TRHO1 = SQRT(TOT1/TOT2)
EVALRHO=TRHO1
ENDIF
ENDIF
*----
* VARIATIONAL ACCELERATION. LIVOLANT INNER ITERATION LOOP FOR ONE-GROUP
* TRANSPORT EQUATION.
*----
DMU=1.0D0
IF(LIVO.AND.(MOD(ITER-1,ICL1+ICL2).GE.ICL1)) THEN
F1=0.0
F2=0.0
DO 30 I=1,NUN
R1=OLD2(I,II)-OLD1(I,II)
R2=FUNKNO(I,II)-OLD2(I,II)
F1=F1+R1*(R2-R1)
F2=F2+(R2-R1)*(R2-R1)
30 CONTINUE
DMU=-F1/F2
IF(DMU.GT.0.0) THEN
RDMU=REAL(DMU)
DO 40 I=1,NUN
FUNKNO(I,II)=OLD2(I,II)+RDMU*(FUNKNO(I,II)-OLD2(I,II))
OLD2(I,II)=OLD1(I,II)+RDMU*(OLD2(I,II)-OLD1(I,II))
40 CONTINUE
ENDIF
ENDIF
*----
* CALCULATE ERROR AND TEST FOR CONVERGENCE
*----
AAA=0.0
BBB=0.0
DO 50 I=1,NREG
IF(KEYFLX(I).EQ.0) GO TO 50
AAA=MAX(AAA,ABS(FUNKNO(KEYFLX(I),II)-OLD2(KEYFLX(I),II)))
BBB=MAX(BBB,ABS(FUNKNO(KEYFLX(I),II)))
50 CONTINUE
IF(IMPX.GT.2) WRITE(IUNOUT,300) NGIND(II),ITER,AAA,BBB,
1 AAA/BBB,DMU
IF(IMPX.GT.5) THEN
ALLOCATE(FGAR(NREG))
FGAR(:NREG)=0.0
DO I=1,NREG
IF(KEYFLX(I).NE.0) FGAR(I)=FUNKNO(KEYFLX(I),II)
ENDDO
WRITE(IUNOUT,'(//33H N E U T R O N F L U X E S :)')
WRITE(IUNOUT,'(1P,6(5X,E15.7))') (FGAR(I),I=1,NREG)
DEALLOCATE(FGAR)
ENDIF
IF(AAA.LE.0.1*EPSINR*BBB) THEN
LNCONV=LNCONV-1
INCONV(II)=.FALSE.
ENDIF
IF(ITER.EQ.1) TEST(II)=AAA
! Be careful if the value of ITER is changed below.
IF((ITER.GE.10).AND.(AAA.GT.TEST(II))) THEN
WRITE(IUNOUT,'(39H SNF: UNABLE TO CONVERGE ONE-SPEED ITER,
1 15HATIONS IN GROUP,I5,1H.)') NGIND(II)
LNCONV=LNCONV-1
INCONV(II)=.FALSE.
ENDIF
ENDIF
60 CONTINUE
IF((NFOU.GT.0).AND.(LNCONV.EQ.0)) WRITE(6,310) TRHO1
IF(LNCONV.EQ.0) GO TO 70
GO TO 10
*----
* CONVERGENCE OF ONE-SPEED ITERATIONS IN ALL NGEFF GROUPS
*----
70 IF(IMPX.GT.1) WRITE(IUNOUT,'(29H SNF: NUMBER OF ONE-SPEED ITE,
1 8HRATIONS=,I5,1H.)') ITER
DEALLOCATE(OLD3,TEST,OLD2,OLD1,INCONV)
ENDIF
IF(IMPX.GT.2) THEN
CALL KDRCPU(TK2)
WRITE(IUNOUT,'(15H SNF: CPU TIME=,1P,E11.3,8H SECOND./)')
1 TK2-TK1
ENDIF
RETURN
*
300 FORMAT(11H SNF: GROUP,I5,20H ONE-SPEED ITERATION,I4,8H ERROR=,
1 1P,E11.4,5H OVER,E11.4,5H PREC,E12.4,22H ACCELERATION FACTOR=,
2 0P,F7.3)
310 FORMAT (44H SNF: EIGENVALUE FOR FOURIER ANALYSIS, RHO= ,E13.6)
END
|
