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
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
|
*DECK KINDRV
SUBROUTINE KINDRV(NEN,KEN,CMOD,NGR,NBM,NBFIS,NDG,NLF,ITY,NEL,
1 LL4,NUN,NUP,TTF,TTP,DT,IMPH,ICL1,ICL2,NADI,ADJ,MAXOUT,EPSOUT,
2 MAXINR,EPSINR,IFL,IPR,IEXP,INORM,IMPX,POWTOT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Driver to perform the space-time kinetics calculations.
*
*Copyright:
* Copyright (C) 2008 Ecole Polytechnique de Montreal.
*
*Author(s): D. Sekki
*
*Parameters: input
* NEN number of LCM objects used in the module.
* KEN addresses of LCM objects: (1) L_KINET; (2) L_MACROLIB;
* (3) L_TRACK; (4) L_SYSTEM; (5) L_MACROLIB.
* CMOD name of the assembly door (BIVAC or TRIVAC).
* NGR number of energy groups.
* NBM number of material mixtures.
* NBFIS number of fissile isotopes.
* NDG number of delayed-neutron groups.
* NLF number of Legendre orders for fluxes.
* ITY type of finite elements and tracking.
* NEL total number of finite elements.
* LL4 number of flux unknowns per energy group.
* NUN total number of unknowns per energy group.
* NUP number of precursor unknowns per delayed group.
* TTF value of theta-parameter for fluxes.
* TTP value of theta-parameter for precursors.
* DT current time increment.
* IMPH management of convergence histogram.
* ICL1 number of free iterations in one cycle of the inverse power
* method
* ICL2 number of accelerated iterations in one cycle
* NADI number of inner adi iterations per outer iteration
* ADJ flag for adjoint space-time kinetics calculation
* MAXOUT maximum number of outer iterations
* EPSOUT convergence criteria for the flux
* MAXINR maximum number of thermal iterations.
* EPSINR thermal iteration epsilon.
* IFL temporal integration scheme for fluxes.
* IPR temporal integration scheme for precursors.
* IEXP exponential transformation flag (=1 to activate).
* INORM type of flux normalization (=0: no normalization; =1: imposed
* factor; =2: maximum flux; =3 initial power).
* IMPX printing parameter (=0 for no print).
*
*Parameter: output
* POWTOT power.
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NEN,NGR,NBM,NBFIS,NDG,NLF,ITY,NEL,LL4,NUN,NUP,IMPH,ICL1,
1 ICL2,NADI,MAXOUT,MAXINR,IFL,IPR,IEXP,INORM,IMPX
TYPE(C_PTR) KEN(NEN)
REAL TTF,TTP,DT,EPSOUT,EPSINR,POWTOT
CHARACTER CMOD*12
LOGICAL ADJ
*----
* LOCAL VARIABLES
*----
PARAMETER(IOS=6)
INTEGER MAT(NEL),IDL(NEL),IDLPC(NEL)
REAL VOL(NEL),PDC(NDG),PMAX(NDG,NBFIS)
LOGICAL LNUD,LCHD
TYPE(C_PTR) IPMAC,IPSYS
REAL, DIMENSION(:), ALLOCATABLE :: DNF,AVG1,AVG2,WORK1,RM
REAL, DIMENSION(:,:), ALLOCATABLE :: EVECT,DNS,PHO,OVR,OMEGA
REAL, DIMENSION(:,:,:), ALLOCATABLE :: PC,CHI,SGF
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: SGD,CHD,SGO
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: SRC
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(EVECT(NUN,NGR),PC(NUP,NDG,NBFIS),SGD(NBM,NBFIS,NGR,NDG),
1 OMEGA(NBM,NGR))
*----
* RECOVER INFORMATION
*----
CALL KDRCPU(TK1)
TA1=TK1
IF(IMPX.GT.0) WRITE(IOS,1001)
CALL LCMGET(KEN(1),'E-KEFF',EVL)
CALL LCMGET(KEN(1),'LAMBDA-D',PDC)
CALL LCMGET(KEN(1),'E-IDLPC',IDLPC)
CALL LCMLEN(KEN(1),'OMEGA',ILONG,ITYLCM)
IF((IEXP.EQ.0).OR.(ILONG.EQ.0)) THEN
OMEGA(:NBM,:NGR)=0.0
ELSE
CALL LCMGET(KEN(1),'OMEGA',OMEGA)
ENDIF
CALL LCMGET(KEN(3),'VOLUME',VOL)
CALL LCMGET(KEN(3),'MATCOD',MAT)
CALL LCMGET(KEN(3),'KEYFLX',IDL)
*----
* RECOVER CROSS SECTIONS (BEGINNING-OF-STEP)
*----
ALLOCATE(DNF(NDG),DNS(NGR,NDG))
CALL LCMLEN(KEN(1),'BETA-D',LEN,ITYL)
LNUD=(LEN.EQ.NDG)
IF(LNUD) CALL LCMGET(KEN(1),'BETA-D',DNF)
CALL LCMLEN(KEN(1),'CHI-D',LEN,ITYL)
LCHD=(LEN.EQ.NGR*NDG)
IF(LCHD) CALL LCMGET(KEN(1),'CHI-D',DNS)
ALLOCATE(OVR(NBM,NGR),CHI(NBM,NBFIS,NGR),CHD(NBM,NBFIS,NGR,NDG),
1 SGF(NBM,NBFIS,NGR),SGO(NBM,NBFIS,NGR,NDG))
IF(NEN.EQ.4) THEN
IPMAC=KEN(2)
IPSYS=KEN(4)
ELSE IF(NEN.EQ.6) THEN
IPMAC=KEN(5)
IPSYS=KEN(6)
ENDIF
CALL KINXSD(IPMAC,NGR,NBM,NBFIS,NDG,EVL,DT,DNF,DNS,LNUD,LCHD,OVR,
1 CHI,CHD,SGF,SGO)
*----
* COMPUTE THE SOURCE TERM
*----
LL4B=LL4
IF((ITY.EQ.11).OR.(ITY.EQ.13)) LL4B=LL4*NLF/2
ALLOCATE(PHO(NUN,NGR),SRC(NUN,NGR))
CALL LCMGET(KEN(1),'E-PREC',PC)
CALL LCMGET(KEN(1),'E-VECTOR',PHO)
CALL KINSRC(KEN(3),IPSYS,CMOD,IMPX,IFL,IPR,IEXP,NGR,NBM,NBFIS,NDG,
1 ITY,LL4B,NUN,NUP,PDC,TTF,TTP,DT,ADJ,OVR,CHI,CHD,SGF,SGO,OMEGA,
2 PHO,PC,SRC)
*----
* RECOVER CROSS SECTIONS (END-OF-STEP)
*----
CALL KINXSD(KEN(2),NGR,NBM,NBFIS,NDG,EVL,DT,DNF,DNS,LNUD,LCHD,
1 OVR,CHI,CHD,SGF,SGD)
DEALLOCATE(DNS,DNF)
*----
* RECOVER THE BEGINNING-OF-STEP FLUX
*----
IF(IMPX.GT.0)THEN
CALL KDRCPU(TA2)
WRITE(IOS,1002) TA2-TA1
WRITE(IOS,1003)
ENDIF
DO 15 IGR=1,NGR
DO 10 IND=1,NUN
EVECT(IND,IGR)=PHO(IND,IGR)
10 CONTINUE
15 CONTINUE
*----
* COMPUTE THE FLUX SOLUTION
*----
IF(CMOD.EQ.'BIVAC')THEN
IF(ADJ) CALL XABORT('KINDRV: ADJOINT CALCULATION NOT IMPLEMENT'
1 //'ED WITH BIVAC.')
CALL KINSLB(KEN(3),KEN(4),KEN(1),LL4B,ITY,NUN,NGR,IFL,IPR,
1 IEXP,NBM,NBFIS,NDG,ICL1,ICL2,IMPX,IMPH,TITR,EPSOUT,MAXINR,
2 EPSINR,MAXOUT,PDC,TTF,TTP,DT,OVR,CHI,CHD,SGF,SGD,OMEGA,EVECT,
3 SRC)
ELSEIF(CMOD.EQ.'TRIVAC')THEN
CALL KINSLT(KEN(3),KEN(4),KEN(1),LL4B,ITY,NUN,NGR,IFL,IPR,
1 IEXP,NBM,NBFIS,NDG,ICL1,ICL2,IMPX,IMPH,TITR,EPSOUT,MAXINR,
2 EPSINR,NADI,ADJ,MAXOUT,PDC,TTF,TTP,DT,OVR,CHI,CHD,SGF,SGD,
3 OMEGA,EVECT,SRC)
ENDIF
DEALLOCATE(SRC)
*----
* COMPUTE THE PRECURSOR SOLUTION
*----
CALL KINPRC(KEN(3),KEN(4),CMOD,NGR,NBM,NBFIS,NDG,NEL,LL4,NUN,NUP,
1 MAT,VOL,IDLPC,EVECT,PHO,CHD,CHO,SGD,SGO,PDC,DT,ADJ,TTP,PC,IPR,
2 IEXP,OMEGA,IMPX)
CALL LCMPUT(KEN(1),'E-PREC',NDG*NUP*NBFIS,2,PC)
*----
* COMPUTE THE EXPONENTIAL TRANSFORMATION FACTORS
*----
IF(IEXP.EQ.1) THEN
ALLOCATE(WORK1(LL4),RM(LL4),AVG1(NBM),AVG2(NBM))
DO 35 IGR=1,NGR
DO 20 IBM=1,NBM
AVG1(IBM)=EXP(OMEGA(IBM,IGR)*DT)
20 CONTINUE
IF(CMOD.EQ.'BIVAC')THEN
CALL KINBLM(KEN(3),NBM,LL4,AVG1,EVECT(1,IGR),WORK1)
CALL MTLDLS('RM',KEN(3),KEN(4),LL4,1,WORK1)
ELSEIF(CMOD.EQ.'TRIVAC')THEN
CALL KINTLM(KEN(3),NBM,LL4,AVG1,EVECT(1,IGR),WORK1)
CALL LCMLEN(KEN(4),'RM',ILONG,ITYLCM)
CALL LCMGET(KEN(4),'RM',RM)
DO 25 IND=1,ILONG
FACT=RM(IND)
IF(FACT.EQ.0.0) CALL XABORT('KINDRV: SINGULAR RM.')
WORK1(IND)=WORK1(IND)/FACT
25 CONTINUE
ENDIF
DO 30 IND=1,LL4
EVECT(IND,IGR)=WORK1(IND)
30 CONTINUE
35 CONTINUE
CALL LCMPUT(KEN(1),'E-VECTOR',NGR*NUN,2,EVECT)
*
DO 60 IGR=1,NGR
AVG1(:NBM)=0.0
AVG2(:NBM)=0.0
DO 40 IEL=1,NEL
IBM=MAT(IEL)
IF(IBM.GT.0) THEN
AVG1(IBM)=AVG1(IBM)+VOL(IEL)*PHO(IDL(IEL),IGR)
AVG2(IBM)=AVG2(IBM)+VOL(IEL)*EVECT(IDL(IEL),IGR)
ENDIF
40 CONTINUE
DO 50 IBM=1,NBM
RATIO=MIN(10.0,ABS(AVG2(IBM)/AVG1(IBM)))
OMEGA(IBM,IGR)=LOG(RATIO)/DT
50 CONTINUE
IF(IMPX.GT.1) THEN
WRITE(IOS,1006) (OMEGA(IBM,IGR),IBM=1,NBM)
ENDIF
60 CONTINUE
CALL LCMPUT(KEN(1),'OMEGA',NBM*NGR,2,OMEGA)
DEALLOCATE(AVG2,AVG1,RM,WORK1)
ENDIF
*----
* COMPUTE AVERAGED FLUX VALUES.
*----
DO 70 IGR=1,NGR
IF(CMOD.EQ.'BIVAC')THEN
CALL FLDBIV(KEN(3),NEL,NUP,EVECT(1,IGR),MAT,VOL,IDL)
ELSEIF(CMOD.EQ.'TRIVAC')THEN
CALL FLDTRI(KEN(3),NEL,NUP,EVECT(1,IGR),MAT,VOL,IDL)
ENDIF
70 CONTINUE
CALL LCMPUT(KEN(1),'E-VECTOR',NGR*NUN,2,EVECT)
*----
* FIND THE MAXIMUM FLUX VALUE
*----
FMAX=0.0
IDMX=0
DO 85 IGR=1,NGR
DO 80 IEL=1,NEL
IND=IDL(IEL)
IF(IND.EQ.0) GO TO 80
IF(ABS(EVECT(IND,IGR)).GT.FMAX) THEN
FMAX=EVECT(IND,IGR)
IDMX=IEL
IGMX=IGR
ENDIF
80 CONTINUE
85 CONTINUE
IF(IDMX.EQ.0) CALL XABORT('KINDRV: UNABLE TO SET FMAX.')
IND=IDLPC(IDMX)
IF(IND.EQ.0) CALL XABORT('KINDRV: UNABLE TO SET PMAX.')
DO 95 IFIS=1,NBFIS
DO 90 IDG=1,NDG
PMAX(IDG,IFIS)=PC(IND,IDG,IFIS)
90 CONTINUE
95 CONTINUE
IF(IMPX.GT.0) THEN
WRITE(IOS,1004) FMAX,IDMX,IGMX
CALL KDRCPU(TK2)
WRITE(IOS,1005)TK2-TK1
ENDIF
CALL LCMPUT(KEN(1),'CTRL-FLUX',1,2,FMAX)
CALL LCMPUT(KEN(1),'CTRL-PREC',NDG*NBFIS,2,PMAX)
CALL LCMPUT(KEN(1),'CTRL-IDL',1,1,IDMX)
CALL LCMPUT(KEN(1),'CTRL-IGR',1,1,IGMX)
*----
* COMPUTE REACTOR POWER
*----
IF(INORM.EQ.3) THEN
CALL KINPOW(KEN(2),NGR,NBM,NUN,NEL,MAT,VOL,IDL,EVECT,POWTOT)
CALL LCMPUT(KEN(1),'E-POW',1,2,POWTOT)
IF(IMPX.GT.0) WRITE(6,*) 'REACTOR POWER (MW) =',POWTOT
ENDIF
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(PHO,SGO,SGF,CHD,CHI,OVR)
DEALLOCATE(OMEGA,SGD,PC,EVECT)
RETURN
*
1001 FORMAT(/1X,'=> ASSEMBLY OF THE SYSTEM MATRICES'/)
1002 FORMAT(/1X,'TOTAL CPU TIME USED FOR THE ASSEMBLING',
1 1X,'OF ALL SYSTEM MATRICES =',F6.3/)
1003 FORMAT(/1X,'=> COMPUTING THE KINETICS SOLUTION'/)
1004 FORMAT(/1X,'CONTROLLING PARAMETERS:',2X,'MAX-VA',
1 'L',1X,1PE12.5,3X,'IDL #',I5.5,3X,'IGR #',I2.2/)
1005 FORMAT(/1X,'TOTAL CPU TIME USED FOR KINETICS CALC',
1 'ULATIONS =',F10.3//1X,'=> SPACE-TIME',1X,
2 'KINETICS CALCULATION IS DONE.')
1006 FORMAT(39H KINDRV: MIXTURE-ORDERED OMEGA FACTORS:/(1P,10E14.6))
END
|
