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
|
*DECK NSSFL1
SUBROUTINE NSSFL1(IPFLX,NUN,NG,NEL,NMIX,NALB,ITRIAL,EPSOUT,MAXOUT,
1 MAT,XX,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,FD,IPRINT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Flux calculation for the nodal expansion method in Cartesian 1D
* geometry.
*
*Copyright:
* Copyright (C) 2021 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
* IPFLX nodal flux.
* NUN number of unknowns (=4*NEL+1).
* NG number of energy groups.
* NEL number of nodes in the nodal calculation.
* NMIX number of mixtures in the nodal calculation.
* NALB number of physical albedos.
* ITRIAL type of expansion functions in the nodal calculation
* (=1: polynomial; =2: hyperbolic).
* EPSOUT convergence epsilon for the power method.
* MAXOUT maximum number of iterations for the power method.
* MAT material mixtures.
* XX mesh spacings.
* IQFR boundary condition information.
* QFR albedo function information.
* DIFF diffusion coefficients
* SIGR removal cross sections.
* CHI fission spectra.
* SIGF nu times fission cross section.
* SCAT scattering cross section.
* BETA albedos.
* FD discontinuity factors.
* IPRINT edition flag.
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPFLX
INTEGER NUN,NG,NEL,NMIX,NALB,ITRIAL(NMIX,NG),MAXOUT,IPRINT,
1 MAT(NEL),IQFR(6,NEL)
REAL EPSOUT,XX(NEL),QFR(6,NEL),DIFF(NMIX,NG),SIGR(NMIX,NG),
1 CHI(NMIX,NG),SIGF(NMIX,NG),SCAT(NMIX,NG,NG),BETA(NALB,NG,NG),
2 FD(NMIX,2,NG,NG)
*----
* LOCAL VARIABLES
*----
TYPE(C_PTR) :: JPFLX
INTEGER :: DIM
REAL :: KEFF
REAL, ALLOCATABLE, DIMENSION(:) :: WORK
REAL, ALLOCATABLE, DIMENSION(:,:) :: EVECT,A,B,AI,A11,QFR2,FUNKN
*
ALB(X)=0.5*(1.0-X)/(1.0+X)
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(EVECT(NUN,NG))
DIM=5*NEL
ALLOCATE(FUNKN(DIM,NG),A(DIM*NG,DIM*NG),B(DIM*NG,DIM*NG))
ALLOCATE(WORK(NMIX),A11(DIM,DIM),QFR2(6,NEL))
*----
* INITIALIZATIONS
*----
CALL LCMLEN(IPFLX,'FLUX',ILONG,ITYLCM)
IF(ILONG == 0) THEN
JPFLX=LCMLID(IPFLX,'FLUX',NG)
EVECT(:NUN,:NG)=1.0
KEFF=1.0
ELSE
JPFLX=LCMGID(IPFLX,'FLUX')
DO IG=1,NG
CALL LCMLEL(JPFLX,IG,ILONG,ITYLCM)
IF(ILONG /= NUN) CALL XABORT('NSSFL3: INVALID FLUX.')
CALL LCMGDL(JPFLX,IG,EVECT(1,IG))
ENDDO
CALL LCMGET(IPFLX,'K-EFFECTIVE',KEFF)
ENDIF
FUNKN(:,:)=0.0
DO IEL=1,NEL
IOF=(IEL-1)*5
FUNKN(IOF+1,:)=EVECT(IEL,:)
ENDDO
*----
* COMPUTE NODAL SOLUTION
*----
DIM=5*NEL
A(:DIM*NG,:DIM*NG)=0.0
B(:DIM*NG,:DIM*NG)=0.0
QFR2(:6,:NEL)=0.0
DO J=1,NG
IOF1=(J-1)*DIM
DO I=1,NG
DO IQW=1,2
DO IEL=1,NEL
IALB=IQFR(IQW,IEL)
IF(IALB.GT.0) THEN
IF(IALB.GT.NALB) CALL XABORT('NSSFL1: BETA OVERFLOW.')
QFR2(IQW,IEL)=QFR(IQW,IEL)*ALB(BETA(IALB,I,J))
ELSE IF(I == J) THEN
QFR2(IQW,IEL)=QFR(IQW,IEL)
ELSE
QFR2(IQW,IEL)=0.0
ENDIF
ENDDO
ENDDO
DO IBM=1,NMIX
WORK(IBM)=CHI(IBM,I)*SIGF(IBM,J)
ENDDO
IOF2=(I-1)*DIM
IF(I == J) THEN
CALL NSS1TR(ITRIAL(1,J),NEL,NMIX,MAT,XX,IQFR,QFR2,DIFF(:,I),
1 SIGR(:,I),FD(:,:,I,J),A11)
A(IOF1+1:IOF1+DIM,IOF1+1:IOF1+DIM)=A11(:,:)
ELSE
CALL NSS2TR(ITRIAL(1,J),NEL,NMIX,MAT,XX,IQFR,QFR2,DIFF(:,J),
1 SIGR(:,J),SCAT(:,I,J),FD(:,:,I,J),A11)
A(IOF2+1:IOF2+DIM,IOF1+1:IOF1+DIM)=-A11(:,:)
ENDIF
CALL NSS3TR(ITRIAL(1,J),NEL,NMIX,MAT,XX,DIFF(:,J),SIGR(:,J),
1 WORK(:),A11)
B(IOF2+1:IOF2+DIM,IOF1+1:IOF1+DIM)=A11(:,:)
ENDDO
ENDDO
DEALLOCATE(QFR2,A11,WORK)
*----
* SOLVE EIGENVALUE MATRIX SYSTEM
*----
CALL ALINV(DIM*NG,A,DIM*NG,IER)
IF(IER.NE.0) CALL XABORT('NSSFL1: SINGULAR MATRIX')
ALLOCATE(AI(DIM*NG,DIM*NG))
AI(:DIM*NG,:DIM*NG)=MATMUL(A(:DIM*NG,:DIM*NG),B(:DIM*NG,:DIM*NG))
CALL AL1EIG(DIM*NG,AI,EPSOUT,MAXOUT,ITER,FUNKN,KEFF,IPRINT)
IF(IPRINT.GT.0) WRITE(6,10) KEFF,ITER
DEALLOCATE(AI,B,A)
*----
* NORMALIZE THE FLUX
*----
FLMAX=0.0
DO IG=1,NG
NUM1=0
DO IEL=1,NEL
IF(ABS(FUNKN(NUM1+1,IG)).GT.ABS(FLMAX)) FLMAX=FUNKN(NUM1+1,IG)
NUM1=NUM1+5
ENDDO
ENDDO
FUNKN(:,:)=FUNKN(:,:)/FLMAX
*----
* COMPUTE INTERFACE FLUXES AND CURRENTS
*----
IOF1=NEL
IOF2=2*NEL
IOF3=3*NEL
IF(IOF3+NEL+1.NE.NUN) CALL XABORT('NSSFL1: NUN ERROR.')
DO IG=1,NG
DO KEL=1,NEL
IBM=MAT(KEL)
IOF=(KEL-1)*5
EVECT(KEL,IG)=FUNKN(IOF+1,IG)
EVECT(IOF1+KEL,IG)=FUNKN(IOF+1,IG)+0.5*(-FUNKN(IOF+2,IG)+
1 FUNKN(IOF+3,IG))
EVECT(IOF2+KEL,IG)=FUNKN(IOF+1,IG)+0.5*(FUNKN(IOF+2,IG)+
1 FUNKN(IOF+3,IG))
IF(ITRIAL(IBM,IG).EQ.1) THEN
EVECT(IOF3+KEL,IG)=-(DIFF(IBM,IG)/XX(KEL))*(FUNKN(IOF+2,IG)-
1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)/2.0-FUNKN(IOF+5,IG)/5.0)
ELSE
ETA=XX(KEL)*SQRT(SIGR(IBM,IG)/DIFF(IBM,IG))
ALP1=ETA*COSH(ETA/2.0)-2.0*SINH(ETA/2.0)
EVECT(IOF1+KEL,IG)=EVECT(IOF1+KEL,IG)-FUNKN(IOF+4,IG)*
1 SINH(ETA/2.0)+FUNKN(IOF+5,IG)*ALP1/ETA
EVECT(IOF2+KEL,IG)=EVECT(IOF2+KEL,IG)+FUNKN(IOF+4,IG)*
1 SINH(ETA/2.0)+FUNKN(IOF+5,IG)*ALP1/ETA
EVECT(IOF3+KEL,IG)=-(DIFF(IBM,IG)/XX(KEL))*(FUNKN(IOF+2,IG)-
1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)*ETA*COSH(ETA/2.0)-
2 FUNKN(IOF+5,IG)*ETA*SINH(ETA/2.0))
ENDIF
ENDDO
IBM=MAT(NEL)
IOF=(NEL-1)*5
IF(ITRIAL(IBM,IG).EQ.1) THEN
EVECT(IOF3+NEL+1,IG)=-(DIFF(IBM,IG)/XX(NEL))*(FUNKN(IOF+2,IG)+
1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)/2.0+FUNKN(IOF+5,IG)/5.0)
ELSE
ETA=XX(NEL)*SQRT(SIGR(IBM,IG)/DIFF(IBM,IG))
EVECT(IOF3+NEL+1,IG)=-(DIFF(IBM,IG)/XX(NEL))*(FUNKN(IOF+2,IG)+
1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)*ETA*COSH(ETA/2.0)+
2 FUNKN(IOF+5,IG)*ETA*SINH(ETA/2.0))
ENDIF
IF(IPRINT.GT.2) THEN
WRITE(6,'(/33H NSSFL1: AVERAGED FLUXES IN GROUP,I5)') IG
WRITE(6,'(1P,10e12.4)') (EVECT(I,IG),I=1,NEL)
WRITE(6,'(/39H NSSFL1: SURFACIC NET CURRENTS IN GROUP,I5)') IG
WRITE(6,'(1P,10e12.4)') (EVECT(IOF3+I,IG),I=1,NEL+1)
ENDIF
ENDDO
*----
* SAVE SOLUTION
*----
JPFLX=LCMGID(IPFLX,'FLUX')
DO IG=1,NG
CALL LCMPDL(JPFLX,IG,NUN,2,EVECT(1,IG))
ENDDO
CALL LCMPUT(IPFLX,'K-EFFECTIVE',1,2,KEFF)
DEALLOCATE(FUNKN,EVECT)
RETURN
*
10 FORMAT(14H NSSFL1: KEFF=,F11.8,12H OBTAINED IN,I5,11H ITERATIONS)
END
|
