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
|
*DECK SNFT1S
SUBROUTINE SNFT1S(NREG,NMAT,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL,
1 SURF,XXX,TOTAL,IGAV,QEXT,LFIXUP,CURR,FLUX)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Perform one inner iteration for solving SN equations in 1D spherical
* 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.
* NLF number of SN directions.
* NSCT number of Legendre components in the flux.
* =1 isotropic sources;
* =2 linearly anisotropic sources.
* U base points in $\\mu$ of the 1D quadrature.
* W weights of the 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.
* XXX radii.
* 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,NLF,NSCT,MAT(NREG),IGAV
REAL U(NLF),W(NLF),ALPHA(NLF),PLZ(NSCT),PL(NSCT,NLF),VOL(NREG),
1 SURF(NREG+1),XXX(NREG+1),TOTAL(0:NMAT),QEXT(NSCT,NREG),CURR,
2 FLUX(NSCT,NREG)
LOGICAL LFIXUP
*----
* LOCAL VARIABLES
*----
DOUBLE PRECISION Q,E2,AFB,Q1,Q2,PSIA,WSIA,CURSUM
PARAMETER(RLOG=1.0E-8)
REAL, ALLOCATABLE, DIMENSION(:) :: AFGL,FLXB
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(AFGL(NREG),FLXB(NLF/2))
*----
* COMPUTE A NORMALIZATION CONSTANT.
*----
DENOM=0.0
DO 10 I=1,NLF
IF(U(I).GT.0.0) DENOM=DENOM+W(I)*U(I)
10 CONTINUE
*----
* 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 20 IL=0,NSCT-1
Q=Q+QEXT(IL+1,I)*PLZ(IL+1)*(2.0*IL+1.0)/2.0
20 CONTINUE
Q1=(XXX(I+1)-XXX(I))*Q+2.0*AFB
Q2=(XXX(I+1)-XXX(I))*TOTAL(IBM)+2.0
E2=Q1/Q2
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
IF(IGAV.EQ.2) THEN
PSIA=0.0D0
WSIA=0.0D0
ELSE
PSIA=AFB
ENDIF
*----
* BACKWARD SWEEP (FROM EXTERNAL SURFACE TOWARD CENTRAL AXIS).
*----
FLUX(:NSCT,:NREG)=0.0
CURSUM=0.0D0
ALPMIN=0.0
DO 80 IP=1,NLF/2
AFB=CURR/DENOM
ALPMAX=ALPHA(IP)
DO 70 I=NREG,1,-1 ! Spatial sweep
Q=0.0
IBM=MAT(I)
IOF=0
DO 40 IL=0,NSCT-1
Q=Q+QEXT(IL+1,I)*PL(IL+1,IP)*(2.0*IL+1.0)/2.0
40 CONTINUE
Q1=Q*VOL(I)-U(IP)*(SURF(I)+SURF(I+1))*AFB+0.5D0*(SURF(I+1)-
1 SURF(I))*(ALPMIN+ALPMAX)*AFGL(I)/W(IP)
Q2=TOTAL(IBM)*VOL(I)-2.0D0*U(IP)*SURF(I)+(SURF(I+1)-SURF(I))*
1 ALPMAX/W(IP)
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)+W(IP)*REAL(E2)*PL(K,IP)
60 CONTINUE
70 CONTINUE
IF(IGAV.EQ.1) THEN
FLXB(IP)=REAL(AFB)
ELSE IF(IGAV.EQ.2) THEN
PSIA=PSIA+W(IP)*REAL(AFB)
WSIA=WSIA+W(IP)
ENDIF
ALPMIN=ALPMAX
80 CONTINUE
IF(IGAV.EQ.2) PSIA=PSIA/WSIA
*----
* FORWARD SWEEP (FROM CENTRAL AXIS TOWARD EXTERNAL SURFACE).
*----
DO 130 IP=1+NLF/2,NLF
ALPMAX=ALPHA(IP)
IF(IGAV.EQ.1) THEN
AFB=FLXB(NLF-IP+1)
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 90 IL=0,NSCT-1
Q=Q+QEXT(IL+1,I)*PL(IL+1,IP)*(2.0*IL+1.0)/2.0
90 CONTINUE
Q1=Q*VOL(I)+U(IP)*(SURF(I)+SURF(I+1))*AFB+0.5D0*(SURF(I+1)-
1 SURF(I))*(ALPMIN+ALPMAX)*AFGL(I)/W(IP)
Q2=TOTAL(IBM)*VOL(I)+2.0D0*U(IP)*SURF(I+1)+(SURF(I+1)-SURF(I))*
1 ALPMAX/W(IP)
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)+W(IP)*REAL(E2)*PL(K,IP)
110 CONTINUE
120 CONTINUE
CURSUM=CURSUM+W(IP)*U(IP)*AFB
ALPMIN=ALPMAX
130 CONTINUE
CURR=REAL(CURSUM)
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(FLXB,AFGL)
RETURN
END
|
