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
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
|
*DECK MCTRK
SUBROUTINE MCTRK(IPTRK,IPRINT,NFREG,NFSUR,NDIM,NMIX,ANGBC,ITYPBC,
1 MAXMSH,NUCELL,MXGSUR,MXGREG,MAXPIN,BCRT,ICODE,
2 ALBEDO,IUNFLD,DGMESH,XYZL,INDEX,IDREG,DCMESH,
3 ITPIN,DRAPIN,ISEED,NGRP,NL,NFM,NDEL,NED,MATCOD,
4 XST,XSS,XSN2N,XSN3N,XSSNN,XSNUSI,XSCHI,XSEDI,
5 MIX,ISONBR,NU,POS,ITALLY,NBSCO,NMERGE,NGCOND,
6 IMERGE,INDGRP,SCORE1,SCORE2)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Simulation of a single particle from source to death.
*
*Copyright:
* Copyright (C) 2008 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): B. Arsenault
*
*Parameters: input
* IPTRK pointer to the TRACKING data structure.
* IPRINT print level.
* NFREG number of regions.
* NFSUR number of surfaces.
* NDIM problem dimensions.
* NMIX number of mixtures in the geometry.
* ANGBC angular treatment for boundary conditions (=0 isotropic;
* =1 specular).
* ITYPBC type of boundary.
* MAXMSH maximum number of elements in MESH array.
* NUCELL number of cell after unfolding in $X$, $Y$ and $Z$ directions.
* MXGSUR maximum number of surfaces for any geometry.
* MXGREG maximum number of region for any geometry.
* MAXPIN maximum number of pins in a cell.
* BCRT reflection/translation array.
* ICODE albedo index array.
* ALBEDO albedo array.
* IUNFLD description of unfolded geometry.
* DGMESH meshing vector for global geometry.
* XYZL Cartesian boundary coordinates.
* ISEED the seed for the generation of random numbers.
* NGRP number of energy groups.
* NL number of Legendre orders required in the estimations
* (NL=1 or higher).
* NFM number of fissile isotopes.
* NDEL number of delayed precursor groups.
* NED number of extra edit vectors.
* XST total macroscopic cross sections for each mixture and energy
* group.
* XSS total scattering cross sections for each mixture and energy
* group.
* XSN2N N2N macroscopic cross sections for each mixture and energy
* group.
* XSN3N N3N macroscopic cross sections for each mixture and energy
* group.
* MATCOD region material.
* XSSNN in-group and out-of-group macroscopic transfert cross sections
* for each mixture.
* XSNUSI the values of Nu time the fission cross sections for each
* isotope per mixture and energy group.
* XSCHI the values of fission spectrum per isotope per mixture for
* each energy group.
* XSEDI extra edit cross sections for each mixture and energy group.
* MIX the mixture number where the fission occurs.
* ISONBR the isotopic number where the fission occurs.
* NU the value of the particle weight.
* POS location of the particle in the x, y, and z directions.
* ITALLY type of tally (=0 no tally; =1 score effective
* multiplication factor; =2 also score macrolib information).
* NBSCO number of macrolib related scores.
* NMERGE number of homogenized regions.
* NGCOND number of condensed energy groups.
* IMERGE homogenized regions indices.
* INDGRP condensed groups indices.
* INDEX undefined.
* IDREG undefined.
* DCMESH undefined.
* ITPIN undefined.
* DRAPIN undefined.
*
*Parameters: input/output
* SCORE1 score for total flux and effective multiplication factor.
* SCORE2 macrolib score matrix.
*
*-----------------------------------------------------------------------
*
USE GANLIB
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPTRK
INTEGER IPRINT,NFREG,NFSUR,NDIM,NMIX,ANGBC,ITYPBC,MAXMSH,
1 NUCELL(3),MXGSUR,MXGREG,MAXPIN,MIX,ISONBR,
2 BCRT(NFSUR),ICODE(6),IUNFLD(2,NUCELL(1),NUCELL(2),NUCELL(3)),
3 INDEX(5,-MXGSUR:MXGREG,2),IDREG(MXGREG,2),ITPIN(3,MAXPIN),
4 ISEED,NGRP,NL,NFM,NDEL,NED,MATCOD(NFREG),ITALLY,NBSCO,NMERGE,
5 NGCOND,IMERGE(NFREG),INDGRP(NGRP)
REAL ALBEDO(6),XST(NMIX,NGRP),XSS(NMIX,NGRP,NL),XSN2N(NMIX,NGRP),
1 XSN3N(NMIX,NGRP),XSSNN(NMIX,NGRP,NGRP,NL),
2 XSNUSI(NMIX,NFM,NGRP,1+NDEL),XSCHI(NMIX,NFM,NGRP,1+NDEL),
3 XSEDI(NMIX,NGRP,NED),NU,SCORE1(3),SCORE2(NBSCO,NMERGE,NGCOND)
DOUBLE PRECISION DGMESH(-1:MAXMSH,4),XYZL(2,NDIM),
1 DCMESH(-1:MAXMSH,4,2),DRAPIN(-1:4,MAXPIN),POS(3)
*----
* LOCAL VARIABLES
*----
INTEGER NBIND
PARAMETER (NBIND=7)
DOUBLE PRECISION XDRCST,PI,TWOPI
INTEGER INDX(NBIND,0:2),ODIR(3),IREG,IDIRC(2),MESHC(4,2),NSURC(2),
1 NREGC(2),NTPIN,IDIR,ILEV,IDS,IDSO,IPSP,IEV,II,IGROUP,IGR,
2 ISUR,ISUR2,IBM,IDIM,IFIRST
DOUBLE PRECISION CELLPO(3,2),PINCEN(3),LENGTH,MU,MUS,VIS,PHI,
2 VDIR(3),OMEGA(3),ALB,NUSIGF0
LOGICAL TRACK,ACTIVE,KILLED,VIRTUAL
PARAMETER (ACTIVE=.TRUE.,KILLED=.FALSE.)
REAL RAND,XSM
DOUBLE PRECISION PSCAT,PN2N,PN3N,XPROB,XX
*----
* DATA VALUES
*----
INTEGER INDOS(2,3),INDIC(-3:3)
DATA INDOS / 2,3,
1 3,1,
2 1,2 /
DATA INDIC / 5,3,1,0,2,4,6/
*----
* LOAD CONSTANTS
*----
PI=XDRCST('Pi',' ')
TWOPI=2.D0*PI
IFIRST=1
*----
* TRACK THE NEUTRON PATH
*----
DO ILEV=0,2
DO IDIR=1,3
INDX(IDIR,ILEV)=1
ENDDO
DO IDIR=4,NBIND
INDX(IDIR,ILEV)=0
ENDDO
ENDDO
ODIR(:3)=1
TRACK=ACTIVE
LENGTH=0.D0
IDSO=0
VIRTUAL=.FALSE.
DO WHILE(TRACK.EQV.ACTIVE)
IF(IDSO.EQ.0) THEN
IF(LENGTH.EQ.0.D0) THEN
IF(MIX.EQ.-1) THEN
IGR=1
ELSE
CALL RANDF(ISEED,IFIRST,RAND)
XPROB=0.D0
IGR=0
DO WHILE((RAND.GE.XPROB).AND.(IGR.LT.NGRP))
IGR=IGR+1
XPROB=XPROB+XSCHI(MIX,ISONBR,IGR,1)
ENDDO
ENDIF
*----
* SET THE TRACK DIRECTION FOR THE FIRST FLIGHT
*----
XSM=0.0
DO II=1,NMIX
XSM=MAX(XSM,XST(II,IGR))
ENDDO
CALL RANDF(ISEED,IFIRST,RAND)
MU=2.D0*RAND-1.D0
MUS=DSQRT(1.D0-MU*MU)
CALL RANDF(ISEED,IFIRST,RAND)
PHI=TWOPI*RAND
VDIR(1)=MUS*COS(PHI)
VDIR(2)=MUS*SIN(PHI)
VDIR(3)=MU
ELSE
*----
* IF IT IS A VIRTUAL COLLISION, KEEP THE SAME DIRECTION,
* IF IT IS NOT A VIRTUAL COLLISION IT CONSISTS OF AN ISOTROPIC
* SCATTERING REACTION.
*----
IF(.NOT.VIRTUAL) THEN
CALL RANDF(ISEED,IFIRST,RAND)
PHI=TWOPI*RAND
CALL RANDF(ISEED,IFIRST,RAND)
MU=COS(PI*(2.D0*RAND-1.D0))
MUS=DSQRT(1.D0-MU*MU)
VIS=DSQRT(1.D0-VDIR(1)*VDIR(1))
OMEGA(1)=VDIR(1)*MU-VIS*MUS*COS(PHI)
OMEGA(2)=VDIR(2)*MU+MUS*(VDIR(1)*VDIR(2)*
1 COS(PHI)-VDIR(3)*SIN(PHI))/VIS
OMEGA(3)=VDIR(3)*MU+MUS*(VDIR(1)*VDIR(3)*
1 COS(PHI)+VDIR(2)*SIN(PHI))/VIS
DO II=1,3
VDIR(II)=OMEGA(II)
ENDDO
ENDIF
ENDIF
*----
* SAMPLE THE FREE PATH DISTANCE
*----
CALL RANDF(ISEED,IFIRST,RAND)
LENGTH=-LOG(RAND)/XSM
ENDIF
*----
* LOCATE THE NEUTRON IN THE GEOMETRY
*----
CALL MCTPTR(IPTRK,IPRINT,NDIM,MAXMSH,ITYPBC,NUCELL,MXGSUR,
1 MXGREG,MAXPIN,IUNFLD,DGMESH,XYZL,NBIND,POS,LENGTH,VDIR,
2 ODIR,IDS,IDSO,IREG,INDX,IDIRC,MESHC,NSURC,NREGC,NTPIN,
3 CELLPO,PINCEN,INDEX,IDREG,DCMESH,ITPIN,DRAPIN)
*---
* TALLY PROCESSING
*---
IF(ITALLY.GT.0) THEN
CALL MCTALLY(ITALLY,NFREG,NMIX,NGRP,NL,NFM,NDEL,NED,NBSCO,
1 NMERGE,NGCOND,IREG,IGR,NU,MATCOD,IMERGE,INDGRP,XSM,XST,XSS,
2 XSN2N,XSN3N,XSSNN,XSNUSI,XSCHI,XSEDI,SCORE1,SCORE2)
ENDIF
*----
* AN INTERACTION HAS BEEN DETECTED, DETECT IF IT IS A VIRTUAL OR A
* REAL COLLISION
*----
IF(IREG.GT.0) THEN
IBM=MATCOD(IREG)
CALL RANDF(ISEED,IFIRST,RAND)
VIRTUAL=(RAND.LE.((XSM-XST(IBM,IGR))/XSM))
*---
* DETERMINE THE TYPE OF REACTION
*---
IF(.NOT.VIRTUAL) THEN
CALL RANDF(ISEED,IFIRST,RAND)
PSCAT=XSS(IBM,IGR,1)/XST(IBM,IGR)
PN2N=XSN2N(IBM,IGR)/XST(IBM,IGR)
PN3N=XSN3N(IBM,IGR)/XST(IBM,IGR)
IF(RAND.LE.PSCAT+PN2N+PN3N) THEN
*---
* ISOTROPIC SCATTERING OR NxN EVENT
*---
IF(RAND.GT.PSCAT+PN2N) THEN
NU=3.0*NU
ELSE IF(RAND.GT.PSCAT) THEN
NU=2.0*NU
ENDIF
XX = 0.D0
CALL RANDF(ISEED,IFIRST,RAND)
IGROUP=NGRP
DO WHILE((RAND.GT.XX).AND.(IGROUP.GE.1))
XX=XX+XSSNN(IBM,IGROUP,IGR,1)/XSS(IBM,IGR,1)
IGROUP=IGROUP-1
ENDDO
IGR=IGROUP+1
TRACK=ACTIVE
XSM=0.0
DO II=1,NMIX
XSM=MAX(XSM,XST(II,IGR))
ENDDO
ELSE
*---
* CAPTURE OR FISSION EVENT
*---
TRACK = KILLED
MIX = IBM
NUSIGF0=0.0D0
DO II=1,NFM
NUSIGF0=NUSIGF0+XSNUSI(IBM,II,IGR,1)
ENDDO
IF((NUSIGF0.GT.0.0).AND.(NFM.EQ.1)) THEN
NU=NU*XSNUSI(IBM,1,IGR,1)/(XST(IBM,IGR)-XSS(IBM,IGR,1)-
1 XSN2N(IBM,IGR)-XSN3N(IBM,IGR))
ISONBR=1
ELSE IF(NUSIGF0.GT.0.0) THEN
CALL RANDF(ISEED,IFIRST,RAND)
XPROB=0.D0
ISONBR=0
DO WHILE((RAND.GE.XPROB).AND.(ISONBR.LT.NFM))
ISONBR=ISONBR+1
XPROB=XPROB+XSNUSI(IBM,ISONBR,IGR,1)/NUSIGF0
ENDDO
NU=NU*XSNUSI(IBM,ISONBR,IGR,1)/(XST(IBM,IGR)-
1 XSS(IBM,IGR,1)-2.0*XSN2N(IBM,IGR)-3.0*XSN3N(IBM,IGR))
ELSE
NU=0.0
ISONBR=0
ENDIF
ENDIF
ENDIF
ELSE
*----
* A BOUNDARY CONDITION HAS BEEN ENCOUNTERED
*----
IF(ITYPBC.EQ.0) THEN
*----
* CARTESIAN BOUNDARY
*----
ISUR=-IREG
ISUR2=BCRT(ISUR)
IREG=-ISUR2
IF((ISUR2.EQ.ISUR).AND.(ANGBC.EQ.1)) THEN
* SPECULAR REFLECTIVE BOUNDARY CONDITION
IF(IPRINT.GT.4) WRITE(6,*) 'SPECULAR REFLECTION ON ',ISUR
ALB=ALBEDO(-ICODE(INDIC(IDSO)))
IF(ALB.EQ.0.0) THEN
* no leakage
ISONBR = 0
TRACK = KILLED
ELSE IF(ALB.EQ.1.0) THEN
VDIR(IDS)=-VDIR(IDS)
ELSE
CALL RANDF(ISEED,IFIRST,RAND)
TRACK=(RAND.LE.ALB)
VDIR(IDS)=-VDIR(IDS)
ENDIF
ELSE IF(ISUR2.NE.ISUR) THEN
* PERIODIC BOUNDARY CONDITION
IF(IPRINT.GT.4) WRITE(6,*) 'BC TRANSLATION FROM ',ISUR,
1 ' TO',ISUR2
IPSP=-ISUR
IEV=-INDIC(IDSO)
IF(IPRINT.GT.99) CALL MCTPSP(IPTRK,POS,IPSP,IEV)
IF(IDSO.GT.0) THEN
POS(IDS)=POS(IDS)-(XYZL(2,IDS)-XYZL(1,IDS))
ELSE
POS(IDS)=POS(IDS)+(XYZL(2,IDS)-XYZL(1,IDS))
ENDIF
IDSO=-IDSO
ELSE IF((ISUR2.EQ.ISUR).AND.(ANGBC.EQ.0)) THEN
* WHITE BOUNDARY CONDITION
ALB=ALBEDO(-ICODE(INDIC(IDSO)))
IF(ALB.EQ.0.0) THEN
* no leakage
ISONBR = 0
TRACK = KILLED
ELSE
* otherwise, choose randomly according to albedo
CALL RANDF(ISEED,IFIRST,RAND)
TRACK=(RAND.LE.ALB)
IF(TRACK.EQV.ACTIVE) THEN
* FOR ISOTROPIC BC, CHOOSE RANDOMLY THE REENTERING
* DIRECTION
CALL RANDF(ISEED,IFIRST,RAND)
MU=2.D0*RAND-1.D0
MUS=DSQRT(1.D0-MU*MU)
CALL RANDF(ISEED,IFIRST,RAND)
PHI=PI*RAND
VDIR(IDS)=-SIGN(1,IDSO)*MUS*SIN(PHI)
VDIR(INDOS(1,IDS))=MUS*COS(PHI)
VDIR(INDOS(2,IDS))=MU
DO IDIM=1,NDIM
IF(IDIM.NE.IDS) THEN
CALL RANDF(ISEED,IFIRST,RAND)
POS(IDIM)=XYZL(1,IDIM)+RAND*(XYZL(2,IDIM)-
1 XYZL(1,IDIM))
ENDIF
ENDDO
CALL RANDF(ISEED,IFIRST,RAND)
LENGTH=-LOG(RAND)/XSM
ELSE
ISONBR = 0
ENDIF
ENDIF
ELSE
CALL XABORT('MCTRK: INVALID TYPE OF BOUNDARY CONDITION.')
ENDIF
ELSE
*----
* CYLINDRICAL BOUNDARY
*----
CALL XABORT('MCTRK: CYLINDRICAL BOUNDARY NOT IMPLEMENTED.')
ENDIF
IF(IPRINT.GT.99) THEN
IPSP=-IREG
IEV=-INDIC(IDSO)
CALL MCTPSP(IPTRK,POS,IPSP,IEV)
ENDIF
ENDIF
*----
* SAVE NEUTRON PATHS IN NXT TABLE FOR PSP DISPLAY
*----
IF(IPRINT.GT.99) THEN
IF(TRACK.EQV.KILLED) THEN
IPSP=-ABS(IREG)
ELSE
IPSP=ABS(IREG)
ENDIF
CALL MCTPSP(IPTRK,POS,IPSP,1)
ENDIF
ENDDO
RETURN
END
|
