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
|
*DECK SNT1DC
SUBROUTINE SNT1DC (IMPX,LX,NCODE,ZCODE,XXX,NLF,NPQ,NSCT,IQUAD,
1 JOP,U,W,TPQ,UPQ,VPQ,WPQ,ALPHA,PLZ,PL,VOL,IDL,SURF,IL,IM)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Numbering corresponding to a 1-D cylindrical geometry with discrete
* ordinates approximation of the flux.
*
*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
* IMPX print parameter.
* LX number of elements along the R axis.
* NCODE type of boundary condition applied on each side
* (i=1 R-; i=2 R+):
* =1 VOID; =2 REFL; =7 ZERO.
* ZCODE ZCODE(I) is the albedo corresponding to boundary condition
* 'VOID' on each side (ZCODE(I)=0.0 by default).
* XXX coordinates along the R axis.
* NLF order of the SN approximation (even number).
* NPQ number of SN directions in two octants.
* NSCT maximum number of spherical harmonics moments of the flux.
* IQUAD type of SN quadrature (1 Level symmetric, type IQUAD;
* 4 Gauss-Legendre and Gauss-Chebyshev; 10 product).
*
*Parameters: output
* JOP number of base points per axial level in one octant.
* U base points in $\\xi$ of the axial quadrature. Used with
* zero-weight points.
* W weights for the quadrature in $\\mu$.
* TPQ base points in $\\xi$ of the 2D SN quadrature.
* UPQ base points in $\\mu$ of the 2D SN quadrature.
* VPQ base points in $\\eta$ 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.
* PL discrete values of the spherical harmonics corresponding
* to the 2D SN quadrature.
* VOL volume of each element.
* IDL position of integrated fluxes into unknown vector.
* SURF surfaces.
* IL indexes (l) of each spherical harmonics in the
* interpolation basis.
* IM indexes (m) of each spherical harmonics in the
* interpolation basis.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER IMPX,LX,NCODE(2),NLF,NPQ,NSCT,IQUAD,JOP(NLF/2),IDL(LX),
1 IL(NSCT),IM(NSCT)
REAL ZCODE(2),XXX(LX+1),U(NLF/2),W(NLF/2),PLZ(NSCT,NLF/2),
1 PL(NSCT,NPQ),TPQ(NPQ),UPQ(NPQ),VPQ(NPQ),WPQ(NPQ),ALPHA(NPQ),
2 VOL(LX),SURF(LX+1)
*----
* LOCAL VARIABLES
*----
PARAMETER(PI=3.141592654)
REAL, ALLOCATABLE, DIMENSION(:) :: TPQ2,UPQ2,VPQ2,WPQ2
*----
* GENERATE QUADRATURE BASE POINTS AND CORRESPONDING WEIGHTS.
*----
IF(MOD(NLF,2).EQ.1) CALL XABORT('SNT1DC: EVEN NLF EXPECTED.')
M2=NLF/2
ALLOCATE(TPQ2(NPQ/2),UPQ2(NPQ/2),VPQ2(NPQ/2),WPQ2(NPQ/2))
IF(IQUAD.EQ.1) THEN
CALL SNQU01(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2)
ELSE IF(IQUAD.EQ.2) THEN
CALL SNQU02(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2)
ELSE IF(IQUAD.EQ.3) THEN
CALL SNQU03(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2)
ELSE IF(IQUAD.EQ.4) THEN
CALL SNQU04(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2)
ELSE IF(IQUAD.EQ.10) THEN
CALL SNQU10(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2)
ELSE
CALL XABORT('SNT1DC: UNKNOWN QUADRATURE TYPE.')
ENDIF
IPQ=0
IPR=0
WSUM=0.0
DO IP=1,M2
DO IQ=1,JOP(IP)
IPQ=IPQ+1
IPR=IPR+1
TPQ(IPQ)=TPQ2(IPR)
UPQ(IPQ)=UPQ2(IPR)
VPQ(IPQ)=VPQ2(IPR)
WPQ(IPQ)=WPQ2(IPR)
WSUM=WSUM+WPQ(IPQ)
ENDDO
DO IQ=JOP(IP)+1,2*JOP(IP)
IPQ=IPQ+1
JQ=IQ-(JOP(IP)+1)+1
JPQ=IPQ-2*JQ+1
TPQ(IPQ)=TPQ(JPQ)
UPQ(IPQ)=-UPQ(JPQ)
VPQ(IPQ)=VPQ(JPQ)
WPQ(IPQ)=WPQ(JPQ)
WSUM=WSUM+WPQ(IPQ)
ENDDO
ENDDO
DEALLOCATE(WPQ2,VPQ2,UPQ2,TPQ2)
IF(IPQ.NE.NPQ) CALL XABORT('SNT1DC: BAD VALUE ON NPQ.')
*----
* COMPUTE ALPHA.
*----
IPQ=0
DO IP=1,M2
SUMETA=0.0
DO IQ=1,2*JOP(IP)
IPQ=IPQ+1
ALPHA(IPQ)=SUMETA
SUMETA=SUMETA+WPQ(IPQ)*UPQ(IPQ)
ENDDO
ENDDO
*----
* PRINT COSINES AND WEIGHTS.
*----
IF(IMPX.GT.1) THEN
WRITE(6,'(/20H SNT1DC: WEIGHT SUM=,1P,E11.4)') WSUM
WRITE(6,60) (U(N),N=1,M2)
60 FORMAT(//,1X,'THE POSITIVE QUADRATURE COSINES FOLLOW'//
1 (1X,5E14.6))
WRITE(6,70) (W(N),N=1,M2)
70 FORMAT(//,1X,'THE CORRESPONDING QUADRATURE WEIGHTS FOLLOW'//
1 (1X,5E14.6))
WRITE(6,74) (UPQ(N),N=1,NPQ)
74 FORMAT(//,1X,'THE BASE POINTS (MU) FOLLOW'//(1X,5E14.6))
WRITE(6,75) (WPQ(N),N=1,NPQ)
75 FORMAT(//,1X,'THE WEIGHTS FOLLOW'//(1X,5E14.6))
WRITE(6,76) (ALPHA(N),N=1,NPQ)
76 FORMAT(//,1X,'THE ALPHAS FOLLOW'//(1X,5E14.6))
ENDIF
*----
* GENERATE SPHERICAL HARMONICS FOR SCATTERING SOURCE.
*----
IOF=0
IL(:NSCT)=0
IM(:NSCT)=0
DO 130 L=0,NLF-1
DO 120 M=0,L
IF(MOD(L+M,2).EQ.1) GO TO 120
IOF=IOF+1
IL(IOF)=L
IM(IOF)=M
IF(IOF.GT.NSCT) GO TO 140
DO 100 N=1,M2
PLZ(IOF,N)=PNSH(L,M,U(N),-SQRT(1.0-U(N)*U(N)),0.0)
100 CONTINUE
DO 110 N=1,NPQ
PL(IOF,N)=PNSH(L,M,TPQ(N),UPQ(N),VPQ(N))
110 CONTINUE
120 CONTINUE
130 CONTINUE
*----
* COMPUTE SURFACES AND VOLUMES.
*----
140 DO 150 I=1,LX+1
SURF(I)=2.0*PI*XXX(I)
150 CONTINUE
DO 160 I=1,LX
IDL(I)=(I-1)*NSCT+1
VOL(I)=PI*(XXX(I+1)*XXX(I+1)-XXX(I)*XXX(I))
160 CONTINUE
*----
* SET BOUNDARY CONDITIONS.
*----
IF(NCODE(2).EQ.2) ZCODE(2)=1.0
*
RETURN
END
|
