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
|
*DECK SNFT1C
SUBROUTINE SNFT1C(NREG,NMAT,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA,
1 PLZ,PL,MAT,VOL,SURF,TOTAL,IGAV,QEXT,LFIXUP,CURR,FLUX)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Perform one inner iteration for solving SN equations in 1D cylindrical
* geometry.
*
*Copyright:
* Copyright (C) 2005 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
* NREG number of regions.
* NMAT number of material mixtures.
* M2 number of axial $\\xi$ levels.
* NPQ number of SN directions in one octant.
* ISCAT anisotropy of one-speed sources:
* =1 isotropic sources;
* =2 linearly anisotropic sources.
* NSCT number of spherical harmonics components in the flux.
* JOP number of base points per axial level in one octant.
* U base points in $\\mu$ of the 1D quadrature. Used with
* zero-weight points.
* UPQ base points in $\\mu$ of the 2D SN quadrature.
* WPQ weights of the 2D SN quadrature.
* ALPHA angular redistribution terms.
* PLZ discrete values of the spherical harmonics corresponding
* to the 1D quadrature. Used with zero-weight points.
* MN moment-to-discrete matrix.
* DN discrete-to-moment matrix.
* MAT material mixture index in each region.
* VOL volumes of each region.
* SURF surfaces surrounding each region.
* TOTAL macroscopic total cross sections.
* IGAV type of condition at axial axis (=1 specular reflection;
* =2 zero-weight reflection; =3 averaged reflection).
* QEXT spherical harmonics components of the fixed source.
* LFIXUP flag to enable negative flux fixup.
*
*Parameters: input/output
* CURR entering current at input and leaving current at output.
*
*Parameters: output
* FLUX spherical harmonics components of the flux.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NREG,NMAT,M2,NPQ,ISCAT,NSCT,JOP(M2),MAT(NREG),IGAV
REAL U(M2),UPQ(NPQ),WPQ(NPQ),ALPHA(NPQ),PLZ(NSCT,M2),PL(NSCT,NPQ),
1 VOL(NREG),SURF(NREG+1),TOTAL(0:NMAT),QEXT(NSCT,NREG),CURR,
2 FLUX(NSCT,NREG)
LOGICAL LFIXUP
*----
* LOCAL VARIABLES
*----
PARAMETER(RLOG=1.0E-8,PI=3.141592654)
DOUBLE PRECISION Q,E2,AFB,Q1,Q2,PSIA,WSIA,CURSUM
REAL, ALLOCATABLE, DIMENSION(:) :: FLXB,AFGL
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(FLXB(M2),AFGL(NREG))
*----
* COMPUTE A NORMALIZATION CONSTANT.
*----
DENOM=0.0
DO 10 I=1,NPQ
IF(UPQ(I).GT.0.0) DENOM=DENOM+4.0*WPQ(I)*UPQ(I)
10 CONTINUE
*----
* OUTER LOOP OVER AXIAL LEVELS.
*----
FLUX(:NSCT,:NREG)=0.0
CURSUM=0.0D0
IPQ=0
DO 200 IP=1,M2
*----
* INITIALIZATION SWEEP (USING ZERO-WEIGHT POINTS). AFB IS THE EDGE
* FLUX VALUE AND AFGL(I) IS THE CENTERED FLUX VALUE.
*----
AFB=CURR/DENOM
DO 30 I=NREG,1,-1 ! Spatial sweep
Q=0.0D0
IBM=MAT(I)
IOF=0
DO 25 IL=0,ISCAT-1
DO 20 IM=0,IL
IF(MOD(IL+IM,2).EQ.1) GO TO 20
IOF=IOF+1
Q=Q+QEXT(IOF,I)*PLZ(IOF,IP)*(2.0*IL+1.0)/(4.0*PI)
20 CONTINUE
25 CONTINUE
Q1=-4.0D0*PI*SQRT(1.0-U(IP)*U(IP))/(SURF(I+1)-SURF(I))
E2=(Q-Q1*AFB)/(TOTAL(IBM)-Q1)
IF(LFIXUP.AND.(E2.LE.RLOG)) E2=0.0
AFB=2.0*E2-AFB
IF(LFIXUP.AND.(AFB.LE.RLOG)) AFB=0.0
AFGL(I)=REAL(E2)
IF(LFIXUP.AND.(AFGL(I).LE.RLOG)) AFGL(I)=0.0
30 CONTINUE
WSIA=0.0D0
IF(IGAV.EQ.2) THEN
PSIA=0.0D0
ELSE
PSIA=AFB
ENDIF
*----
* BACKWARD SWEEP (FROM EXTERNAL SURFACE TOWARD CENTRAL AXIS).
*----
ALPMAX=0.0
DO 80 M=2*JOP(IP),JOP(IP)+1,-1 ! Angular sweep
AFB=CURR/DENOM
ALPMIN=ALPHA(IPQ+M)
DO 70 I=NREG,1,-1 ! Spatial sweep
Q=0.0
IBM=MAT(I)
IOF=0
DO 45 IL=0,ISCAT-1
DO 40 IM=0,IL
IF(MOD(IL+IM,2).EQ.1) GO TO 40
IOF=IOF+1
Q=Q+QEXT(IOF,I)*PL(IOF,IPQ+M)*(2.0*IL+1.0)/(4.0*PI)
40 CONTINUE
45 CONTINUE
Q1=Q*VOL(I)-UPQ(IPQ+M)*(SURF(I)+SURF(I+1))*AFB+(SURF(I+1)-
1 SURF(I))*(ALPMIN+ALPMAX)*AFGL(I)/WPQ(IPQ+M)
Q2=TOTAL(IBM)*VOL(I)-2.0D0*UPQ(IPQ+M)*SURF(I)+2.0D0*
1 (SURF(I+1)-SURF(I))*ALPMIN/WPQ(IPQ+M)
E2=Q1/Q2
IF(LFIXUP.AND.(E2.LE.RLOG)) E2=0.0
AFB=2.0*E2-AFB ! IN SPACE
IF(LFIXUP.AND.(AFB.LE.RLOG)) AFB=0.0
AFGL(I)=2.0*REAL(E2)-AFGL(I) ! IN ANGLE
IF(LFIXUP.AND.(AFGL(I).LE.RLOG)) AFGL(I)=0.0
DO 60 K=1,NSCT
FLUX(K,I)=FLUX(K,I)+4.0*WPQ(IPQ+M)*REAL(E2)*PL(K,IPQ+M)
60 CONTINUE
70 CONTINUE
IF(IGAV.EQ.1) THEN
FLXB(2*JOP(IP)-M+1)=REAL(AFB)
ELSE IF(IGAV.EQ.2) THEN
PSIA=PSIA+WPQ(IPQ+M)*REAL(AFB)
WSIA=WSIA+WPQ(IPQ+M)
ENDIF
ALPMAX=ALPMIN
80 CONTINUE
IF(IGAV.EQ.2) PSIA=PSIA/WSIA
*----
* FORWARD SWEEP (FROM CENTRAL AXIS TOWARD EXTERNAL SURFACE).
*----
DO 130 M=JOP(IP),1,-1 ! Angular sweep
ALPMIN=ALPHA(IPQ+M)
IF(IGAV.EQ.1) THEN
AFB=FLXB(M)
ELSE IF(IGAV.EQ.2) THEN
AFB=PSIA
ELSE IF(IGAV.EQ.3) THEN
AFB=PSIA
ENDIF
DO 120 I=1,NREG ! Spatial sweep
Q=0.0
IBM=MAT(I)
IOF=0
DO 100 IL=0,ISCAT-1
DO 90 IM=0,IL
IF(MOD(IL+IM,2).EQ.1) GO TO 90
IOF=IOF+1
Q=Q+QEXT(IOF,I)*PL(IOF,IPQ+M)*(2.0*IL+1.0)/(4.0*PI)
90 CONTINUE
100 CONTINUE
Q1=Q*VOL(I)+UPQ(IPQ+M)*(SURF(I)+SURF(I+1))*AFB+(SURF(I+1)-
1 SURF(I))*(ALPMIN+ALPMAX)*AFGL(I)/WPQ(IPQ+M)
Q2=TOTAL(IBM)*VOL(I)+2.0D0*UPQ(IPQ+M)*SURF(I+1)+2.0D0*
1 (SURF(I+1)-SURF(I))*ALPMIN/WPQ(IPQ+M)
E2=Q1/Q2
IF(LFIXUP.AND.(E2.LE.RLOG)) E2=0.0
AFB=2.0*E2-AFB ! IN SPACE
IF(LFIXUP.AND.(AFB.LE.RLOG)) AFB=0.0
AFGL(I)=2.0*REAL(E2)-AFGL(I) ! IN ANGLE
IF(LFIXUP.AND.(AFGL(I).LE.RLOG)) AFGL(I)=0.0
DO 110 K=1,NSCT
FLUX(K,I)=FLUX(K,I)+4.0*WPQ(IPQ+M)*REAL(E2)*PL(K,IPQ+M)
110 CONTINUE
120 CONTINUE
CURSUM=CURSUM+4.0*WPQ(IPQ+M)*UPQ(IPQ+M)*AFB
ALPMAX=ALPMIN
130 CONTINUE
*----
* END OF OUTER LOOP OVER AXIAL LEVELS
*----
IPQ=IPQ+2*JOP(IP)
200 CONTINUE
IF(IPQ.NE.NPQ) CALL XABORT('SN1T1C: BUG.')
CURR=REAL(CURSUM)
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(AFGL,FLXB)
RETURN
END
|
