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
|
*DECK SNT1DS
SUBROUTINE SNT1DS (IMPX,LX,NCODE,ZCODE,XXX,NLF,NSCT,U,W,ALPHA,
1 PLZ,PL,VOL,IDL,SURF,IQUAD,IL,IM)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Numbering corresponding to a 1-D spherical 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).
* NSCT maximum number of spherical harmonics moments of the flux.
* IQUAD quadrature type:
* =1 Gauss-Lobatto;
* =2 Gauss-Legendre.
*
*Parameters: output
* U base points in $\\xi$ of the axial quadrature.
* W weights for the quadrature in $\\mu$.
* 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,NSCT,IDL(LX),IQUAD,IL(NSCT),IM(NSCT)
REAL ZCODE(2),XXX(LX+1),U(NLF),W(NLF),PLZ(NSCT),PL(NSCT,NLF),
1 ALPHA(NLF),VOL(LX),SURF(LX+1)
*----
* LOCAL VARIABLES
*----
PARAMETER(PI=3.141592654)
*----
* GENERATE QUADRATURE BASE POINTS AND CORRESPONDING WEIGHTS.
*----
IF(MOD(NLF,2).EQ.1) CALL XABORT('SNT1DS: EVEN NLF EXPECTED.')
IF(IQUAD.EQ.1) THEN
CALL SNQU07(NLF,U,W) ! GAUSS-LOBATTO
ELSEIF(IQUAD.EQ.2) THEN
CALL ALGPT(NLF,-1.0,1.0,U,W) ! GAUSS-LEGENDRE
ELSE
CALL XABORT('SNT1DS: UNKNOWN QUADRATURE TYPE.')
ENDIF
*----
* COMPUTE ALPHA.
*----
SUMETA=0.0
DO IP=1,NLF
SUMETA=SUMETA-2.0*W(IP)*U(IP)
ALPHA(IP)=SUMETA
ENDDO
*----
* PRINT COSINES AND WEIGHTS.
*----
IF(IMPX.GT.1) THEN
WRITE(6,60) (U(M),M=1,NLF)
60 FORMAT(//,1X,'THE POSITIVE QUADRATURE COSINES FOLLOW'//
1 (1X,5E14.6))
WRITE(6,70) (W(M),M=1,NLF)
70 FORMAT(//,1X,'THE CORRESPONDING QUADRATURE WEIGHTS FOLLOW'//
1 (1X,5E14.6))
WRITE(6,76) (ALPHA(M),M=1,NLF)
76 FORMAT(//,1X,'THE ALPHAS FOLLOW'//(1X,5E14.6))
ENDIF
*----
* GENERATE LEGENDRE POLYNOMIALS FOR SCATTERING SOURCE.
*----
IL(:NSCT)=0
IM(:NSCT)=0
DO 150 M=1,NLF
PL(1,M)=1.0
IL(1)=0
IF(NSCT.GT.1) THEN
PL(2,M)=U(M)
IL(2)=1
DO 140 L=2,NSCT-1
PL(L+1,M)=((2.0*L-1.0)*U(M)*PL(L,M)-(L-1)*PL(L-1,M))/REAL(L)
IL(L+1)=L
140 CONTINUE
ENDIF
150 CONTINUE
PLZ(1)=1.0
IF(NSCT.GT.1) THEN
PLZ(2)=-1.0
DO 160 L=2,NSCT-1
PLZ(L+1)=(-(2.0*L-1.0)*PLZ(L)-(L-1)*PLZ(L-1))/REAL(L)
160 CONTINUE
ENDIF
*----
* COMPUTE SURFACES AND VOLUMES.
*----
DO 200 I=1,LX+1
SURF(I)=4.0*PI*XXX(I)*XXX(I)
200 CONTINUE
DO 210 I=1,LX
IDL(I)=(I-1)*NSCT+1
VOL(I)=4.0*PI*(XXX(I+1)**3-XXX(I)**3)/3.0
210 CONTINUE
*----
* SET BOUNDARY CONDITIONS.
*----
IF(NCODE(2).EQ.2) ZCODE(2)=1.0
*
RETURN
END
|
