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
|
*DECK SYBRHL
SUBROUTINE SYBRHL(IPRT,NSOUT,NREG,G,PROB)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Stamm'ler normalisation of collision, escape and transmission
* probabilities.
*
*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): R. Roy
*
*Parameters: input
* IPRT print parameter (equal to zero for no print).
* NSOUT number of surfaces.
* NREG number of regions.
*
*Parameters: input/output
* G surface and volume array on input and
* renormalization factors at output.
* PROB collision, escape and transmission probabilities on input
* and normalized collision, escape and transmission
* probabilities at output.
*
*Reference:
* R. Roy et.al.,
* Normalization techniques for CP matrices, Physor-90,
* Marseille/France, v 2, p ix-40 (1990).
* A. Villarino, R.J.J.Stamm'ler and A.A.Ferri and J.J.Casal
* Helios: Angularly dependent collision probabilities E.,
* Nucl.Sci.Eng. 112,16-31, 1992.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER IPRT,NSOUT,NREG
REAL G(NSOUT+NREG),PROB((NSOUT+NREG)*(NSOUT*NREG+1)/2)
*----
* LOCAL VARIABLES
*----
INTEGER CPTLB,CPTAC,CTOT
PARAMETER (CPTLB=3, CPTAC=3, CTOT=CPTAC+CPTLB, IUNOUT=6,
1 EPSCON=1.0E-6, NITMAX=20)
DOUBLE PRECISION WFSPAD,WFSP,NOM,DENOM,DMU
LOGICAL NOTCON
REAL, ALLOCATABLE, DIMENSION(:,:) :: WEIG
*
INDPOS(I,J)=MAX(I,J)*(MAX(I,J)-1)/2+MIN(I,J)
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(WEIG(NSOUT+NREG,3))
*----
* INITIALISATION OF WEIGHTS
*----
NOTCON=.FALSE.
DO 10 IR=1,NSOUT+NREG
WEIG(IR,1)=0.0
WEIG(IR,2)=0.5
WEIG(IR,3)=0.5
10 CONTINUE
*----
* MAIN ITERATION LOOP
*----
IF(IPRT.GT.8) THEN
WRITE(IUNOUT,'(/30H SYBRHL: NORMALIZATION FACTORS)')
WRITE(IUNOUT,'(1X,A24)') 'ITER. MU ERROR '
ENDIF
NIT=0
20 NIT=NIT+1
IF(NIT.GT.NITMAX) THEN
NOTCON=.TRUE.
WRITE(IUNOUT,'(31H SYBRHL: WEIGHTS NOT CONVERGED.)')
GO TO 80
ENDIF
DO 40 IR=1,NSOUT+NREG
WFSPAD=G(IR)+PROB(INDPOS(IR,IR))*WEIG(IR,3)
WFSP=PROB(INDPOS(IR,IR))
DO 30 JR=1,NSOUT+NREG
WFSPAD=WFSPAD-WEIG(JR,3)*PROB(INDPOS(IR,JR))
WFSP=WFSP+PROB(INDPOS(IR,JR))
30 CONTINUE
WEIG(IR,3)=REAL(WFSPAD/WFSP)
40 CONTINUE
*----
* ACCELERATION BY RESIDUAL MINIMIZATION
*----
IF(MOD(NIT-1,CTOT).GE.CPTAC) THEN
NOM = 0.0D0
DENOM = 0.0D0
DO 50 IR=1,NSOUT+NREG
R1= WEIG(IR,2) - WEIG(IR,1)
R2= WEIG(IR,3) - WEIG(IR,2)
NOM = NOM + R1*(R2-R1)
DENOM = DENOM + (R2-R1)*(R2-R1)
50 CONTINUE
IF(DENOM.EQ.0.0D0) THEN
DMU=1.0D0
ELSE
DMU=-NOM/DENOM
ENDIF
ZMU=REAL(DMU)
IF(ZMU.GT.10.0 .OR. ZMU.LT.0.0) THEN
IF( IPRT.GT.2 ) WRITE(IUNOUT,'(I3,1P,G12.4,A)') NIT,ZMU,
> ' =MU / SYBRHL: NON ACCELERATION'
ZMU=1.0
ENDIF
DO 60 IR=1,NSOUT+NREG
WEIG(IR,3)=WEIG(IR,2)+ZMU*(WEIG(IR,3)-WEIG(IR,2))
WEIG(IR,2)=WEIG(IR,1)+ZMU*(WEIG(IR,2)-WEIG(IR,1))
60 CONTINUE
ELSE
ZMU = 1.0
ENDIF
*----
* CALCULATIONS OF SQUARE DISTANCE BETWEEN 2 ITERATIONS AND UPDATING
* OF THE SOLUTION
*----
TOTCON = 0.0
DO 70 IR=1,NSOUT+NREG
TMPCON=ABS(WEIG(IR,3)-WEIG(IR,2))/WEIG(IR,3)
TOTCON=MAX(TMPCON,TOTCON)
WEIG(IR,1)=WEIG(IR,2)
WEIG(IR,2)=WEIG(IR,3)
70 CONTINUE
IF(IPRT.GT.8) WRITE(IUNOUT,'(I3,F9.5,E15.7)') NIT,ZMU,TOTCON
*----
* CONVERGENCE TEST
*----
IF(TOTCON.LT.EPSCON) GO TO 80
GO TO 20
*----
* RENORMALIZE "PIJ" SYMMETRIC MATRIX
*----
80 IPRB=0
DO 95 IR=1,NSOUT+NREG
G(IR)=WEIG(IR,1)
DO 90 JR=1,IR
IPRB=IPRB+1
PROB(IPRB)=PROB(IPRB)*(WEIG(IR,1)+WEIG(JR,1))
90 CONTINUE
95 CONTINUE
*----
* PRINT WEIGHT FACTORS IF THERE IS A PROBLEM
*----
IF(NOTCON .OR. (IPRT.GE.15)) THEN
WRITE(IUNOUT,'(24H SURFACE WEIGHTS FACTORS/)')
NSURC=1
DO 100 IP =1,(9+NSOUT)/10
NSURM=MIN(NSOUT,NSURC+9)
WRITE(IUNOUT,'(10X,10(A5,I6)/)')
> (' SUR ',IR,IR= NSURC, NSURM)
WRITE(IUNOUT,'(10H WEIGHT ,10F11.5)')
> (WEIG(IR,1),IR=NSURC,NSURM)
NSURC=NSURC+10
100 CONTINUE
WRITE(IUNOUT,'(24H VOLUME WEIGHTS FACTORS/)')
NVOLC=NSOUT+1
DO 110 IP=1,(9+NREG)/10
NVOLM=MIN(NSOUT+NREG,NVOLC+9)
WRITE(IUNOUT,'(10X,10(A5,I6)/)')
> (' VOL ',IR,IR=NVOLC-NSOUT,NVOLM-NSOUT)
WRITE(IUNOUT,'(10H WEIGHT ,10F11.5)')
> (WEIG(IR,1),IR=NVOLC,NVOLM)
NVOLC=NVOLC+10
110 CONTINUE
ENDIF
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(WEIG)
RETURN
END
|
