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
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
|
*DECK SNT1DP
SUBROUTINE SNT1DP (IMPX,LX,IELEM,NCODE,ZCODE,XXX,NLF,NSCT,U,W,PL,
1 VOL,IDL,LL4,NUN,LSHOOT,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM,IGLK,MN,
2 DN,IL,IM,IQUAD)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Numbering corresponding to a 1-D slab 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 and C. Bienvenue
*
*Parameters: input
* IMPX print parameter.
* LX number of elements along the X axis.
* IELEM measure of order of the spatial approximation polynomial:
* =1 constant - default for HODD;
* =2 linear - default for DG;
* >3 higher orders.
* NCODE type of boundary condition applied on each side
* (i=1 X-; i=2 X+):
* =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 number of $\\mu$ levels.
* NSCT maximum number of Legendre components in the flux.
* LSHOOT enablig flag for the shooting method.
* EELEM measure of order of the energy approximation polynomial:
* =1 constant - default for HODD;
* =2 linear - default for DG;
* >3 higher orders.
* IBFP type of energy proparation relation:
* =0 no Fokker-Planck term;
* =1 Galerkin type;
* =2 heuristic Przybylski and Ligou type.
* ISCHM method of spatial discretisation:
* =1 High-Order Diamond Differencing (HODD) - default;
* =2 Discontinuous Galerkin finite element method (DG);
* =3 Adaptive weighted method (AWD).
* ESCHM method of energy discretisation:
* =1 High-Order Diamond Differencing (HODD) - default;
* =2 Discontinuous Galerkin finite element method (DG);
* =3 Adaptive weighted method (AWD).
* IGLK angular interpolation type:
* =0 classical SN method.
* =1 Galerkin quadrature method (M = inv(D))
* =2 Galerkin quadrature method (D = inv(M))
* IQUAD quadrature type:
* =1 Gauss-Lobatto;
* =2 Gauss-Legendre.
*
*Parameters: output
* U base points in $\\mu$ of the 1D quadrature.
* W weights for the quadrature in $\\mu$.
* PL discrete values of the Legendre polynomials corresponding
* to the SN quadrature.
* VOL volume of each element.
* IDL position of integrated fluxes into unknown vector.
* LL4 number of unknowns being solved for, over the domain. This
* includes the various moments of the isotropic (and if present,
* anisotropic) flux.
* NUN total number of unknowns stored in the FLUX vector per group.
* This includes LL4 (see above) as well as any surface boundary
* fluxes, if present.
* WX spatial closure relation weighting factors.
* WE energy closure relation weighting factors.
* CST constants for the polynomial approximations.
* MN moment-to-discrete matrix.
* DN discrete-to-moment matrix.
* 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,IELEM,NCODE(2),NLF,NSCT,IDL(LX),LL4,NUN,EELEM,
1 IBFP,ISCHM,ESCHM,IL(NSCT),IM(NSCT),IQUAD,IGLK
REAL ZCODE(2),XXX(LX+1),U(NLF),W(NLF),VOL(LX),WX(IELEM+1),
1 WE(EELEM+1),CST(MAX(IELEM,EELEM)),PX,PE,MN(NLF,NSCT),
2 DN(NSCT,NLF),PL(NSCT,NLF)
LOGICAL LSHOOT
*----
* LOCAL VARIABLES
*----
DOUBLE PRECISION MND(NLF,NSCT)
*----
* GENERATE QUADRATURE BASE POINTS AND CORRESPONDING WEIGHTS.
*----
IF(MOD(NLF,2).EQ.1) CALL XABORT('SNT1DP: 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('SNT1DP: UNKNOWN QUADRATURE TYPE.')
ENDIF
*----
* PRINT COSINES AND WEIGHTS.
*----
IF(IMPX.GT.1) THEN
WRITE(6,60) (U(M),M=1,NLF)
60 FORMAT(//,1X,'THE 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))
ENDIF
*----
* GENERATE LEGENDRE POLYNOMIALS FOR SCATTERING SOURCE.
*----
DO 145 M=1,NLF
PL(1,M)=1.0
IF(NSCT.GT.1) THEN
PL(2,M)=U(M)
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))/L
140 CONTINUE
ENDIF
145 CONTINUE
*----
* GENERATE MAPPING MATRIX FOR GALERKIN QUADRATURE METHOD
*----
MN(:NLF,:NSCT)=0.0D0
DN(:NSCT,:NLF)=0.0D0
IL(:NSCT)=0
IM(:NSCT)=0
IF(IGLK.EQ.1) THEN
DO L=0,NSCT-1
IL(L+1)=L
DO N=1,NLF
DN(L+1,N)=W(N)*PL(L+1,N)
ENDDO
ENDDO
MND=DN
CALL ALINVD(NLF,MND,NLF,IER)
IF(IER.NE.0) CALL XABORT('SNT1DP: SINGULAR MATRIX.')
MN = REAL(MND)
ELSEIF(IGLK.EQ.2) THEN
DO L=0,NSCT-1
IL(L+1)=L
DO N=1,NLF
MN(N,L+1)=(2.0*L+1.0)/2.0*PL(L+1,N)
ENDDO
ENDDO
MND=MN
CALL ALINVD(NLF,MND,NLF,IER)
IF(IER.NE.0) CALL XABORT('SNT1DP: SINGULAR MATRIX.')
DN = REAL(MND)
ELSE
DO L=0,NSCT-1
IL(L+1)=L
DO N=1,NLF
DN(L+1,N)=W(N)*PL(L+1,N)
MN(N,L+1)=(2.0*L+1.0)/2.0*PL(L+1,N)
ENDDO
ENDDO
ENDIF
*----
* GENERATE THE WEIGHTING PARAMETERS OF THE CLOSURE RELATION.
*----
PX=1
PE=1
IF(ISCHM.EQ.1.OR.ISCHM.EQ.3) THEN
PX=1
ELSEIF(ISCHM.EQ.2) THEN
PX=0
ELSE
CALL XABORT('SNT1DP: UNKNOWN TYPE OF SPATIAL CLOSURE RELATION.')
ENDIF
IF(MOD(IELEM,2).EQ.1) THEN
WX(1)=-PX
WX(2:IELEM+1:2)=1+PX
IF(IELEM.GE.2) WX(3:IELEM+1:2)=1-PX
ELSE
WX(1)=PX
WX(2:IELEM+1:2)=1-PX
IF(IELEM.GE.2) WX(3:IELEM+1:2)=1+PX
ENDIF
IF(IBFP.NE.0) THEN
IF(ESCHM.EQ.1.OR.ESCHM.EQ.3) THEN
PE=1
ELSEIF(ESCHM.EQ.2) THEN
PE=0
ELSE
CALL XABORT('SNT1DP: UNKNOWN TYPE OF ENERGY CLOSURE RELATION.')
ENDIF
IF(MOD(EELEM,2).EQ.1) THEN
WE(1)=-PE
WE(2:EELEM+1:2)=1+PE
IF(EELEM.GE.2) WE(3:EELEM+1:2)=1-PE
ELSE
WE(1)=PE
WE(2:EELEM+1:2)=1-PE
IF(EELEM.GE.2) WE(3:EELEM+1:2)=1+PE
ENDIF
ENDIF
! NORMALIZED LEGENDRE POLYNOMIAL CONSTANTS
DO IEL=1,MAX(IELEM,EELEM)
CST(IEL)=SQRT(2.0*IEL-1.0)
ENDDO
*----
* COMPUTE VOLUMES, ISOTROPIC FLUX INDICES and UNKNOWNS.
*----
NM=IELEM*EELEM
NMX=EELEM
LL4=LX*NSCT*NM
NUN=LL4
DO I=1,LX
IDL(I)=(I-1)*NSCT*NM+1
VOL(I)=XXX(I+1)-XXX(I)
ENDDO
IF(.NOT.LSHOOT) NUN=NUN+NMX*NLF
*----
* SET BOUNDARY CONDITIONS.
*----
IF(NCODE(1).EQ.2) ZCODE(1)=1.0
IF(NCODE(2).EQ.2) ZCODE(2)=1.0
IF((NCODE(1).EQ.4).OR.(NCODE(2).EQ.4)) THEN
IF((NCODE(1).NE.4).OR.(NCODE(2).NE.4)) CALL XABORT('SNT1DP: '
1 //'INVALID TRANSLATION BOUNDARY CONDITIONS.')
ZCODE(1)=1.0
ZCODE(2)=1.0
ENDIF
IF((ZCODE(1).EQ.0.0).AND.(ZCODE(2).NE.0.0)) CALL XABORT('SNT1DP:'
1 //' CANNOT SUPPORT LEFT VACUUM BC IF RIGHT BC IS NOT VACUUM.')
*
RETURN
END
|
