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
|
*DECK MCGDS2
SUBROUTINE MCGDS2(SUBDSC,LC,M,N,H,NOM,NZON,TR,SC,W,NFI,DIAGF,
1 DIAGQ,CA,CQ,PREV,NEXT,DINV2,A2,B2)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Calculation of contribution in second-order ACA coefficients on one
* track.
*
*Copyright:
* Copyright (C) 2002 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): I. Suslov and R. Le Tellier
*
*Parameters: input
* SUBDSC ACA coefficients calculation subroutine.
* LC dimension of vector MCU.
* M number of material mixtures.
* N number of elements for this track.
* H tracking widths.
* NOM integer tracking elements.
* NZON index-number of the mixture type assigned to each volume.
* TR macroscopic total cross section.
* SC macroscopic P0 scattering cross section.
* W weight associated with this track.
* NFI total number of volumes and surfaces for which specific values
* of the neutron flux and reactions rates are required.
*Parameters: input/output
* CA undefined.
* CQ undefined.
* DIAGQ undefined.
* DIAGF undefined.
*
*Parameters: scratch
* PREV undefined.
* NEXT undefined.
* DINV2 undefined.
* A2 undefined.
* B2 undefined.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*
*---
* SUBROUTINE ARGUMENTS
*---
INTEGER LC,M,N,NFI,NZON(NFI),NOM(N),PREV(N),NEXT(N)
DOUBLE PRECISION W,H(N),CA(LC),DIAGF(NFI),DINV2(N),A2(N),B2(N)
REAL TR(0:M),SC(0:M),DIAGQ(NFI),CQ(LC)
EXTERNAL SUBDSC
*---
* LOCAL VARIABLES
*---
INTEGER IBCV
DOUBLE PRECISION DMINV,DMAX
PARAMETER(DMINV=2.D-2,IBCV=-7)
DOUBLE PRECISION WW,AAW,CAW,AQW,CQW,CAWP,CQWP,
1 B,A,DINV,DH,DHP,BN,AN,DINVN,BP,AP,DINVP
INTEGER I,I1,NOMI,NZI,NOMIN,NZIN,NOMIP,NZIP,ICN,ICP
*
DMAX=1.D0/DMINV
WW=W
*---
* CALCULATE COEFFICIENTS OF THIS TRACK
*---
* MCGDS2A: Tabulated Exponentials
* MCGDS2E: Exact Exponentials
CALL SUBDSC(N,M,NFI,NOM,NZON,H,TR,SC,DINV2,B2,A2)
*----
* CONSTRUCTION OF ACA MATRICES
*---
* -------------------
* Left outer boundary
* -------------------
I1=2
ICN=NEXT(1)
ICP=PREV(1)
NZIP=IBCV
NOMI=NOM(1)
NZI=NZON(NOMI)
NOMIN=NOM(I1)
NZIN=NZON(NOMIN)
IF (NZI.NE.IBCV) THEN
* Other than void boundary condition (treated as white reflection)
DIAGF(NOMI)=DIAGF(NOMI)+WW
DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*(DINV2(I1)-1.D0))
IF(ICN.GT.0) THEN
CA(ICN)=CA(ICN)-WW*A2(I1)
CQ(ICN)=CQ(ICN)+REAL(W*B2(I1))
ENDIF
ENDIF
* ------------
* Volume Cells
* ------------
DHP=0.0
AP=0.0
BP=0.0
AN=0.0
BN=0.0
DO I=2,N-1
ICN=NEXT(I)
ICP=PREV(I)
NOMIP=NOMI
NZIP=NZI
NOMI=NOMIN
NZI=NZIN
NOMIN=NOM(I+1)
NZIN=NZON(NOMIN)
DINV=DINV2(I)
B=B2(I)
A=A2(I)
DH=0.0
IF (NZIN.EQ.IBCV) THEN
* next cell is a fixed boundary condition
DH=1.D0/(1.D0+DINV)
ELSEIF (NZIN.GE.0) THEN
* next cell is a volume
I1=I+1
DINVN=DINV2(I1)
BN=B2(I1)
AN=A2(I1)
IF (ABS(DINV+DINVN).LT.DMINV) THEN
DH=DMAX
ELSE
DH=1.D0/(DINV+DINVN)
ENDIF
ENDIF
IF (NZIP.EQ.IBCV) THEN
* previous cell is a fixed boundary condition
DHP=1.D0/(1.D0+DINV)
ELSEIF (NZIP.GE.0) THEN
* previous cell is a volume
I1=I-1
DINVP=DINV2(I1)
BP=B2(I1)
AP=A2(I1)
IF (ABS(DINV+DINVP).LT.DMINV) THEN
DHP=DMAX
ELSE
DHP=1.D0/(DINV+DINVP)
ENDIF
ENDIF
* assembling coefficients
IF ((NZIN.LT.0).AND.(NZIN.NE.IBCV)) THEN
* next cell is a surface with reflective boundary condition
AAW=0.D0
AQW=0.D0
CAW=0.D0
CQW=1.D0
ELSE
* next cell is a volume or a fixed boundary condition
AAW=DH*A
AQW=DH*B
IF (NZIN.GE.0) THEN
* next cell is a volume
CAW=DH*AN
CQW=DH*BN
ELSE
* next cell is a fixed boundary condition
CAW=0.D0
CQW=0.D0
ENDIF
ENDIF
IF ((NZIP.LT.0).AND.(NZIP.NE.IBCV)) THEN
* previous cell is a surface with reflective boundary condition
CAWP=0.D0
CQWP=1.D0
ELSE
* previous cell is a volume or a fixed boundary condition
AAW=AAW+DHP*A
AQW=AQW+DHP*B
IF (NZIP.GE.0) THEN
* previous cell is a volume
CAWP=DHP*AP
CQWP=DHP*BP
ELSE
* previous cell is a fixed boundary condition
CAWP=0.D0
CQWP=0.D0
ENDIF
ENDIF
* assembling matrices
DIAGF(NOMI)=DIAGF(NOMI)+AAW*WW
DIAGQ(NOMI)=DIAGQ(NOMI)-REAL(W*AQW)
IF(ICN.GT.0) THEN
* next cell is a volume different from this one
CA(ICN)=CA(ICN)-CAW*WW
CQ(ICN)=CQ(ICN)+REAL(W*CQW)
ELSE
* next cell is a voided boundary or a volume identical to this one
DIAGF(NOMI)=DIAGF(NOMI)-CAW*WW
DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*CQW)
ENDIF
IF(ICP.GT.0) THEN
* previous cell is a volume different from this one
CA(ICP)=CA(ICP)-CAWP*WW
CQ(ICP)=CQ(ICP)+REAL(W*CQWP)
ELSE
* previous cell is a voided boundary or a volume identical to this one
DIAGF(NOMI)=DIAGF(NOMI)-CAWP*WW
DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*CQWP)
ENDIF
ENDDO
* --------------------
* Right outer boundary
* --------------------
ICN=NEXT(N)
ICP=PREV(N)
NOMIP=NOMI
NZIP=NZI
NOMI=NOMIN
NZI=NZIN
NZIN=IBCV
IF (NZI.NE.IBCV) THEN
* Other than void boundary condition (treated as white reflection)
I1=N-1
DIAGF(NOMI)=DIAGF(NOMI)+WW
DIAGQ(NOMI)=DIAGQ(NOMI)+REAL(W*(DINV2(I1)-1.D0))
IF(ICP.GT.0) THEN
CA(ICP)=CA(ICP)-WW*A2(I1)
CQ(ICP)=CQ(ICP)+REAL(W*B2(I1))
ENDIF
ENDIF
*
RETURN
END
|
