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
|
!
!-----------------------------------------------------------------------
!
!Purpose:
! Calculation of cellwise scattering-reduced collision, escape and
! transmission probabilities for the current iteration method in the
! multicell surfacic approximation.
!
!Copyright:
! Copyright (C) 2025 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
! IPSYS pointer to the system matrices.
! IPTRK pointer to the tracking (L_TRACK signature).
! IFTRAK tracking file unit.
! IMPX print flag (equal to zero for no print).
! NREG total number of merged regions for which specific values
! of the neutron flux and reactions rates are required.
! NBMIX number of mixtures.
! SIGT0 total macroscopic cross sections ordered by mixture.
! SIGW0 within-group scattering macroscopic cross section ordered
! by mixture.
! NBATCH number of tracks dispached in eack OpenMP core.
! TITREC title.
! NALBP number of multigroup physical albedos.
! ALBP multigroup physical albedos.
!
!-----------------------------------------------------------------------
!
SUBROUTINE MUSA(IPSYS,IPTRK,IFTRAK,IMPX,NREG,NBMIX,SIGT0,SIGW0,NBATCH, &
& TITREC,NALBP,ALBP)
USE GANLIB
!----
! SUBROUTINE ARGUMENTS
!----
TYPE(C_PTR) IPSYS,IPTRK
INTEGER IFTRAK,IMPX,NREG,NBMIX,NBATCH,NALBP
REAL SIGT0(0:NBMIX),SIGW0(0:NBMIX),ALBP(NALBP)
CHARACTER TITREC*72
!----
! LOCAL VARIABLES
!----
PARAMETER (EPS1=1.0E-4,NSTATE=40)
TYPE(C_PTR) JPTRK,KPTRK
INTEGER ISTATT(NSTATE),NNPSYS(1)
CHARACTER TITRE2*72
logical LSKIP
!----
! ALLOCATABLE ARRAYS
!----
INTEGER, ALLOCATABLE, DIMENSION(:) :: MATALB,NMC_NODE,NMC_SURF,MAT2,IGEN,INUM
REAL, ALLOCATABLE, DIMENSION(:) :: SIGT2,SIGW2,WORK
REAL, ALLOCATABLE, DIMENSION(:,:) :: VOLSUR
REAL, POINTER, DIMENSION(:) :: PSSW,PSJW,PISW,PIJW
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DPROB,DPROBX
TYPE(C_PTR) :: PSSW_PTR,PSJW_PTR,PISW_PTR,PIJW_PTR
!
IND(I,J) = MAX(I+J3+1,J+J3+1)*(MAX(I+J3+1,J+J3+1)-1)/2 &
& + MIN(I+J3+1,J+J3+1)
!
WPR(I,J)= REAL(DPROB( IND(I,J),1 ) / DPROB( IND(I,0),1 ))
!----
! BICKLEY FLAG
!----
SAVE IBICKL
DATA IBICKL/0/
!----
! RECOVER BICKLEY TABLES
!----
IF(IBICKL.EQ.0) THEN
CALL XDRTA2
IBICKL=1
ENDIF
!----
! RECOVER SALT SPECIFIC PARAMETERS
!----
REWIND IFTRAK
CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATT)
IF(NREG.NE.ISTATT(1)) THEN
CALL XABORT('MUSA: STATE VECTOR HAS INVALID # OF ZONES.')
ENDIF
NMACRO=ISTATT(24) ! NGEN
NMCEL=NMACRO
NMERGE=NMACRO
NGEN=NMACRO
ALLOCATE(IGEN(NMERGE),INUM(NMCEL))
DO IK=1,NMERGE
IGEN(IK)=IK
ENDDO
DO IK=1,NMCEL
INUM(IK)=IK
ENDDO
IF(NMACRO.EQ.0) CALL XABORT('MUSA: MUST MODULE TRACKING IS MANDATORY.')
ALLOCATE(NMC_NODE(NMACRO+1),NMC_SURF(NMACRO+1))
JPTRK=LCMGID(IPTRK,'MACRO-TRACK')
CALL LCMGET(IPTRK,'NMC_NODE',NMC_NODE)
CALL LCMGET(IPTRK,'NMC_SURF',NMC_SURF)
NMIX=NMC_SURF(NMACRO+1)
NIFR=NMC_SURF(NMACRO+1)
!----
! LOOP OVER MACRO GEOMETRIES AND COMPUTE PIJ MATRICES USING EXCELP
!----
J1=0
NMIX=0
NPIJ=0
NPIS=0
NPSS=0
DO IMACRO=1,NMACRO
J2=NMC_NODE(IMACRO+1)-NMC_NODE(IMACRO)
J3=NMC_SURF(IMACRO+1)-NMC_SURF(IMACRO)
J1=J1+J2
NMIX=NMIX+J3
NPIJ=NPIJ+J2*J2
NPIS=NPIS+J2*J3
NPSS=NPSS+J3*J3
ENDDO
IF(J1.NE.NREG) CALL XABORT('MUSA: INVALID NREG.')
IF(NMIX.NE.NMC_SURF(NMACRO+1)) CALL XABORT('MUSA: INVALID NMIX.')
PIJW_PTR=LCMARA(NPIJ)
PISW_PTR=LCMARA(NPIS)
PSJW_PTR=LCMARA(NPIS)
PSSW_PTR=LCMARA(NPSS)
CALL C_F_POINTER(PIJW_PTR,PIJW,(/ NPIJ /))
CALL C_F_POINTER(PISW_PTR,PISW,(/ NPIS /))
CALL C_F_POINTER(PSJW_PTR,PSJW,(/ NPIS /))
CALL C_F_POINTER(PSSW_PTR,PSSW,(/ NPSS /))
J1=0
IPIJ=0
IPIS=0
IPSS=0
DO IMACRO=1,NMACRO
J2=NMC_NODE(IMACRO+1)-NMC_NODE(IMACRO)
J3=NMC_SURF(IMACRO+1)-NMC_SURF(IMACRO)
N2PRO=(J2+J3+1)**2
WRITE(TITRE2,'(A,9H -- MACRO,I5.5)') TRIM(TITREC),IMACRO
KPTRK=LCMGIL(JPTRK,IMACRO)
KNORM=4 ! use HELIOS-type normalization
NNPSYS(1)=1
ALLOCATE(MAT2(J2),SIGT2(J2),SIGW2(J2))
ALLOCATE(MATALB(-J3:J2),VOLSUR(-J3:J2,1),DPROB(N2PRO,1),DPROBX(N2PRO,1))
CALL LCMGET(KPTRK,'MATCOD',MAT2)
CALL EXCELP(KPTRK,IFTRAK,IMPX,J3,J2,NBMIX,MAT2,KNORM,SIGT0,1,N2PRO, &
& 1,NNPSYS(1),NBATCH,TITRE2,NALBP,ALBP,MATALB,VOLSUR,DPROB,DPROBX)
!----
! CHECK IF SCATTERING REDUCTION IS REQUIRED
!----
DO I=1,J2
SIGT2(I)=SIGT0(MAT2(I)) ! sigt by node
SIGW2(I)=SIGW0(MAT2(I)) ! sigw by node
ENDDO
LSKIP=.TRUE.
DO I=1,J2
LSKIP=LSKIP.AND.(SIGW2(I).EQ.0.0)
ENDDO
!----
! SCATTERING REDUCTION IF LSKIP=.FALSE.
!----
IF(LSKIP) THEN
! DO NOT PERFORM SCATTERING REDUCTION.
DO I=1,J2
DO J=1,J2
IF(SIGT2(J).EQ.0.0) THEN
PIJW(IPIJ+(J-1)*J2+I)=WPR(I,J)
ELSE
PIJW(IPIJ+(J-1)*J2+I)=WPR(I,J)/SIGT2(J)
ENDIF
ENDDO
ENDDO
DO I=1,J2
DO JC=1,J3
PISW(IPIS+(JC-1)*J2+I)=WPR(I,-JC)
IF(SIGT2(I).EQ.0.0) THEN
PSJW(IPIS+(I-1)*J3+JC)=WPR(-JC,I)
ELSE
PSJW(IPIS+(I-1)*J3+JC)=WPR(-JC,I)/SIGT2(I)
ENDIF
ENDDO
ENDDO
DO IC=1,J3
DO JC=1,J3
PSSW(IPSS+(JC-1)*J3+IC)=WPR(-IC,-JC)
ENDDO
ENDDO
ELSE
! COMPUTE THE SCATTERING-REDUCED COLLISION AND ESCAPE MATRICES.
DO I=1,J2
DO J=1,J2
IF(SIGT2(J).EQ.0.0) THEN
PIJW(IPIJ+(J-1)*J2+I)=0.0
ELSE
PIJW(IPIJ+(J-1)*J2+I)=-WPR(I,J)*SIGW2(J)/SIGT2(J)
ENDIF
ENDDO
PIJW(IPIJ+(I-1)*J2+I)=1.0+PIJW(IPIJ+(I-1)*J2+I)
ENDDO
CALL ALINV(J2,PIJW(IPIJ+1),J2,IER)
IF(IER.NE.0) CALL XABORT('MUSA: SINGULAR MATRIX.')
ALLOCATE(WORK(J2))
DO I=1,J2
DO K=1,J2
WORK(K)=PIJW(IPIJ+(K-1)*J2+I)
ENDDO
DO J=1,J2
WGAR=0.0
DO K=1,J2
IF(SIGT2(J).EQ.0.0) THEN
WGAR=WGAR+WORK(K)*WPR(K,J)
ELSE
WGAR=WGAR+WORK(K)*WPR(K,J)/SIGT2(J)
ENDIF
ENDDO
PIJW(IPIJ+(J-1)*J2+I)=WGAR
ENDDO
DO JC=1,J3
WGAR=0.0
DO K=1,J2
WGAR=WGAR+WORK(K)*WPR(K,-JC)
ENDDO
PISW(IPIS+(JC-1)*J2+I)=WGAR
ENDDO
ENDDO
DEALLOCATE(WORK)
!
! COMPUTE THE SCATTERING-REDUCED COLLISION PROBABILITY MATRIX
! FOR INCOMING NEUTRONS.
DO IC=1,J3
DO J=1,J2
IF(SIGT2(J).EQ.0.0) THEN
WGAR=WPR(-IC,J)
ELSE
WGAR=WPR(-IC,J)/SIGT2(J)
ENDIF
DO K=1,J2
IF(SIGT2(K).NE.0.0) THEN
WGAR=WGAR+WPR(-IC,K)*PIJW(IPIJ+(J-1)*J2+K)*SIGW2(K)/SIGT2(K)
ENDIF
ENDDO
PSJW(IPIS+(J-1)*J3+IC)=WGAR
ENDDO
ENDDO
!
! COMPUTE THE SCATTERING-REDUCED TRANSMISSION PROBABILITY MATRIX.
DO IC=1,J3
DO JC=1,J3
WGAR=WPR(-IC,-JC)
DO K=1,J2
IF(SIGT2(K).NE.0.0) THEN
WGAR=WGAR+WPR(-IC,K)*PISW(IPIS+(JC-1)*J2+K)*SIGW2(K)/SIGT2(K)
ENDIF
ENDDO
PSSW(IPSS+(JC-1)*J3+IC)=WGAR
ENDDO
ENDDO
ENDIF
DEALLOCATE(DPROBX,DPROB,VOLSUR,MATALB)
IF(IMPX.GE.8) THEN
IF(LSKIP) THEN
IN=1
ELSE
IN=2
ENDIF
CALL SYBPRX(IN,J3,J2,IMACRO,SIGT2,SIGW2,PIJW(IPIJ+1),PISW(IPIS+1), &
& PSJW(IPIS+1),PSSW(IPSS+1))
ENDIF
DEALLOCATE(SIGW2,SIGT2,MAT2)
J1=J1+J2
IPIJ=IPIJ+J2*J2
IPIS=IPIS+J2*J3
IPSS=IPSS+J3*J3
ENDDO
! end of SYB004 equivalent
CALL LCMPPD(IPSYS,'PSSW$SYBIL',NPSS,2,PSSW_PTR)
CALL LCMPPD(IPSYS,'PSJW$SYBIL',NPIS,2,PSJW_PTR)
CALL LCMPPD(IPSYS,'PISW$SYBIL',NPIS,2,PISW_PTR)
CALL LCMPPD(IPSYS,'PIJW$SYBIL',NPIJ,2,PIJW_PTR)
IF(IMPX.GT.1) THEN
WRITE(6,'(/31H MUSA: PIJ INFORMATION IN GROUP)')
CALL LCMLIB(IPSYS)
ENDIF
RETURN
END SUBROUTINE MUSA
|
