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
|
*DECK KINPRC
SUBROUTINE KINPRC(IPTRK,IPSYS,CMOD,NGR,NBM,NBFIS,NDG,NEL,LL4,NUN,
1 NUP,MAT,VOL,IDLPC,FLN,FLO,CHD,CHO,SGD,SGO,PDC,DT,ADJ,TTP,PC,IPR,
2 IEXP,OMEGA,IMPX)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the precursors unknowns for the current time step according
* to the pre-defined temporal integration scheme.
*
*Copyright:
* Copyright (C) 2008 Ecole Polytechnique de Montreal.
*
*Author(s): D. Sekki
*
*Parameters: input/output
* IPTRK pointer to L_TRACK object.
* IPSYS pointer to L_SYSTEM object.
* 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.
* 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.
* MAT mixture index assigned to each volume.
* VOL volume of each element.
* IDLPC position of averaged precursor values in unknown vector.
* FLN unknown flux vector at current time step.
* FLO unknown flux vector at previous time step.
* CHD current delayed fission spectrum.
* CHO previous delayed fission spectrum.
* SGD current delayed nu*fission macroscopic x-sections/keff.
* SGO previous delayed nu*fission macroscopic x-sections/keff.
* PDC precursor decay constants.
* DT current time increment.
* ADJ flag for adjoint space-time kinetics calculation.
* TTP value of theta-parameter for precursors.
* PC unknown vector for precursors.
* IPR integration scheme for precursors: =1 implicit;
* =2 Crank-Nicholson; =3 theta; =4 exponential.
* IEXP exponential transformation flag (=1 to activate).
* OMEGA exponential transformation parameter.
* IMPX printing parameter (=0 for no print).
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPTRK,IPSYS
INTEGER NGR,NBM,NBFIS,NDG,NEL,LL4,NUN,NUP,MAT(NEL),IDLPC(NEL),
1 IPR,IEXP,IMPX
REAL VOL(NEL),PDC(NDG),DT,TTP,PC(NUP,NDG,NBFIS),FLN(NUN,NGR),
1 FLO(NUN,NGR),CHD(NBM,NBFIS,NGR,NDG),CHO(NBM,NBFIS,NGR,NDG),
2 SGD(NBM,NBFIS,NGR,NDG),SGO(NBM,NBFIS,NGR,NDG),OMEGA(NBM,NGR)
CHARACTER CMOD*12
LOGICAL ADJ
*----
* LOCAL VARIABLES
*----
PARAMETER(IOS=6)
DOUBLE PRECISION DK,DTP,TK1(NDG),TK2(NDG),TK3(NDG)
REAL, DIMENSION(:), ALLOCATABLE :: GAR1,GAR2
REAL, DIMENSION(:,:), ALLOCATABLE :: XSEXP
REAL, DIMENSION(:), POINTER :: RM
TYPE(C_PTR) RM_PTR
*----
* COMPUTE THE KINETICS FACTORS
*----
TK1(:NDG)=0.0D0
TK2(:NDG)=0.0D0
TK3(:NDG)=0.0D0
DTP=9999.0D0
IF(IPR.EQ.2)THEN
* CRANK-NICHOLSON
DTP=0.5D0
ELSEIF(IPR.EQ.3)THEN
* THETA
DTP=DBLE(TTP)
ENDIF
DO 10 L=1,NDG
DK=PDC(L)*DT
IF(IPR.EQ.1)THEN
* IMPLICIT
TK1(L)=1.0D0/(1.0D0+DK)
TK2(L)=DT/(1.0D0+DK)
ELSEIF(IPR.EQ.4)THEN
* EXPONENTIAL
TK1(L)=DEXP(-DK)
TK2(L)=(1.0D0-(1.0D0-TK1(L))/DK)/PDC(L)
TK3(L)=((1.0D0-TK1(L))/DK-TK1(L))/PDC(L)
ELSE
* GENERAL
TK1(L)=(1.0D0-(1.0D0-DTP)*DK)/(1.0D0+DTP*DK)
TK2(L)=DTP*DT/(1.0D0+DTP*DK)
TK3(L)=(1.0D0-DTP)*DT/(1.0D0+DTP*DK)
ENDIF
10 CONTINUE
*----
* COMPUTE THE PRECURSOR UNKNOWN VECTOR
*----
IF(IMPX.GT.0)WRITE(IOS,1001)CMOD
ALLOCATE(GAR1(NUP),GAR2(NUP),XSEXP(NBM,NGR))
DO 220 IFIS=1,NBFIS
DO 210 IDG=1,NDG
DO 20 I=1,NUP
PC(I,IDG,IFIS)=REAL(TK1(IDG))*PC(I,IDG,IFIS)
20 CONTINUE
GAR2(:NUP)=0.0
IF(.NOT.ADJ) THEN
DO 40 IGR=1,NGR
DO 30 IBM=1,NBM
IF(IEXP.EQ.0) THEN
XSEXP(IBM,IGR)=SGD(IBM,IFIS,IGR,IDG)
ELSE
* exponential transformation
XSEXP(IBM,IGR)=SGD(IBM,IFIS,IGR,IDG)*EXP(OMEGA(IBM,IGR)*DT)
ENDIF
30 CONTINUE
40 CONTINUE
ELSE
DO 60 IGR=1,NGR
DO 50 IBM=1,NBM
IF(IEXP.EQ.0) THEN
XSEXP(IBM,IGR)=PDC(IDG)*CHD(IBM,IFIS,IGR,IDG)
ELSE
* exponential transformation
XSEXP(IBM,IGR)=PDC(IDG)*CHD(IBM,IFIS,IGR,IDG)*
1 EXP(OMEGA(IBM,IGR)*DT)
ENDIF
50 CONTINUE
60 CONTINUE
ENDIF
IF(CMOD.EQ.'BIVAC')THEN
ITY=1
DO 80 IGR=1,NGR
CALL KINBLM(IPTRK,NBM,NUP,XSEXP(1,IGR),FLN(1,IGR),GAR1)
DO 70 IND=1,NUP
GAR2(IND)=GAR2(IND)+GAR1(IND)
70 CONTINUE
80 CONTINUE
CALL MTLDLS('RM',IPTRK,IPSYS,LL4,1,GAR2)
ELSEIF(CMOD.EQ.'TRIVAC')THEN
DO 100 IGR=1,NGR
CALL KINTLM(IPTRK,NBM,NUP,XSEXP(1,IGR),FLN(1,IGR),GAR1)
DO 90 IND=1,NUP
GAR2(IND)=GAR2(IND)+GAR1(IND)
90 CONTINUE
100 CONTINUE
CALL LCMLEN(IPSYS,'RM',ILONG,ITYLCM)
CALL LCMGPD(IPSYS,'RM',RM_PTR)
CALL C_F_POINTER(RM_PTR,RM,(/ ILONG /))
DO 110 IND=1,ILONG
GAR2(IND)=GAR2(IND)/RM(IND)
110 CONTINUE
ENDIF
DO 120 IND=1,NUP
PC(IND,IDG,IFIS)=PC(IND,IDG,IFIS)+REAL(TK2(IDG))*GAR2(IND)
120 CONTINUE
IF(IPR.GT.1) THEN
GAR2(:NUP)=0.0
IF(CMOD.EQ.'BIVAC')THEN
IF(ADJ) CALL XABORT('KINPRC: ADJOINT CALCULATION NOT IMPLEME'
1 //'NTED WITH BIVAC.')
ITY=1
DO 140 IGR=1,NGR
CALL KINBLM(IPTRK,NBM,NUP,SGO(1,IFIS,IGR,IDG),FLO(1,IGR),
1 GAR1)
DO 130 IND=1,NUP
GAR2(IND)=GAR2(IND)+GAR1(IND)
130 CONTINUE
140 CONTINUE
CALL MTLDLS('RM',IPTRK,IPSYS,LL4,1,GAR2)
CALL FLDBIV(IPTRK,NEL,NUP,GAR2(1),MAT,VOL,IDLPC)
ELSEIF(CMOD.EQ.'TRIVAC')THEN
IF(.NOT.ADJ) THEN
DO 160 IGR=1,NGR
CALL KINTLM(IPTRK,NBM,NUP,SGO(1,IFIS,IGR,IDG),FLO(1,IGR),
1 GAR1)
DO 150 IND=1,NUP
GAR2(IND)=GAR2(IND)+GAR1(IND)
150 CONTINUE
160 CONTINUE
ELSE
DO 180 IGR=1,NGR
CALL KINTLM(IPTRK,NBM,NUP,CHO(1,IFIS,IGR,IDG),FLO(1,IGR),
1 GAR1)
DO 170 IND=1,NUP
GAR2(IND)=GAR2(IND)+PDC(IDG)*GAR1(IND)
170 CONTINUE
180 CONTINUE
ENDIF
CALL LCMLEN(IPSYS,'RM',ILONG,ITYLCM)
CALL LCMGPD(IPSYS,'RM',RM_PTR)
CALL C_F_POINTER(RM_PTR,RM,(/ ILONG /))
DO 190 IND=1,ILONG
GAR2(IND)=GAR2(IND)/RM(IND)
190 CONTINUE
ENDIF
DO 200 IND=1,NUP
PC(IND,IDG,IFIS)=PC(IND,IDG,IFIS)+REAL(TK3(IDG))*GAR2(IND)
200 CONTINUE
ENDIF
IF(CMOD.EQ.'BIVAC')THEN
CALL FLDBIV(IPTRK,NEL,NUP,PC(1,IDG,IFIS),MAT,VOL,IDLPC)
ELSEIF(CMOD.EQ.'TRIVAC')THEN
CALL FLDTRI(IPTRK,NEL,NUP,PC(1,IDG,IFIS),MAT,VOL,IDLPC)
ENDIF
210 CONTINUE
220 CONTINUE
DEALLOCATE(XSEXP,GAR1,GAR2)
*----
* EDITION
*----
IF(IMPX.GT.5) THEN
WRITE(IOS,1002)
DO 240 IFIS=1,NBFIS
DO 230 IDG=1,NDG
WRITE(IOS,1003) IDG,IFIS,(PC(IND,IDG,IFIS),IND=1,LL4)
230 CONTINUE
240 CONTINUE
ENDIF
IF(IMPX.GT.2) THEN
DO 260 IFIS=1,NBFIS
WRITE(IOS,1004) IFIS,(IDG,IDG=1,NDG)
DO 250 IEL=1,NEL
IND=IDLPC(IEL)
IF(IND.EQ.0) GO TO 250
WRITE(IOS,1005) IEL,(PC(IND,IDG,IFIS),IDG=1,NDG)
250 CONTINUE
WRITE(IOS,'(/)')
260 CONTINUE
ENDIF
RETURN
*
1001 FORMAT(/1X,'COMPUTING THE PRECURSOR UNKNOWN VECTOR',
1 1X,'ACCORDING TO THE TRACKING TYPE: ',A6/)
1002 FORMAT(/1X,'=> COMPUTED PRECURSOR UNKNOWN VECTOR')
1003 FORMAT(/17H PRECURSOR GROUP=,I5,18H FISSILE ISOTOPE=,I5/
1 (1P,8E14.5))
1004 FORMAT(/51H KINPRC: PRECURSOR UNKNOWN VECTOR (FISSILE ISOTOPE=,
1 I5,1H)/(9X,6I13,:))
1005 FORMAT(1X,I6,2X,1P,6E13.5,:/(9X,6E13.5,:))
END
|
