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
|
*DECK NSSFL2
SUBROUTINE NSSFL2(IPFLX,NUN,NG,NEL,NMIX,NALB,EPSOUT,MAXOUT,MAT,
1 XX,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,FD,IPRINT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Flux calculation for the coarse mesh finite differences 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.
* 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.
*
*Parameters: output
* KEFF effective multiplication factor
* EVECT neutron flux
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPFLX
INTEGER NUN,NG,NEL,NMIX,NALB,MAXOUT,IPRINT,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(:,:) :: 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=3*NEL
ALLOCATE(FUNKN(DIM,NG),A(DIM*NG,DIM*NG),B(DIM*NG,DIM*NG))
ALLOCATE(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)*3
FUNKN(IOF+1,:)=EVECT(IEL,:)
ENDDO
*----
* COMPUTE NODAL SOLUTION
*----
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('NSSFL2: 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
IOF2=(I-1)*DIM
IF(I == J) THEN
CALL NSS4TR(NEL,NMIX,MAT,XX,IQFR,QFR2,DIFF(:,I),SIGR(:,I),
1 FD(:,:,I,J),A11)
A(IOF1+1:IOF1+DIM,IOF1+1:IOF1+DIM)=A11(:,:)
ELSE
CALL NSS5TR(NEL,NMIX,MAT,IQFR,QFR2,SCAT(:,I,J),FD(:,:,I,J),
1 A11)
A(IOF2+1:IOF2+DIM,IOF1+1:IOF1+DIM)=-A11(:,:)
ENDIF
B(IOF2+1:IOF2+DIM,IOF1+1:IOF1+DIM)=0.0
NUM1=0
DO IEL=1,NEL
IBM=MAT(IEL)
B(IOF2+NUM1+1,IOF1+NUM1+1)=CHI(IBM,I)*SIGF(IBM,J)
NUM1=NUM1+3
ENDDO
ENDDO
ENDDO
DEALLOCATE(QFR2,A11)
*----
* SOLVE EIGENVALUE MATRIX SYSTEM
*----
CALL ALINV(DIM*NG,A,DIM*NG,IER)
IF(IER.NE.0) CALL XABORT('NSSFL2: 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+3
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('NSSFL2: NUN ERROR.')
DO IG=1,NG
DO KEL=1,NEL
IBM=MAT(KEL)
IOF=(KEL-1)*3
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))
EVECT(IOF3+KEL,IG)=-(DIFF(IBM,IG)/XX(KEL))*(FUNKN(IOF+2,IG)-
1 3.0*FUNKN(IOF+3,IG))
ENDDO
IBM=MAT(NEL)
IOF=(NEL-1)*3
EVECT(IOF3+NEL+1,IG)=-(DIFF(IBM,IG)/XX(NEL))*(FUNKN(IOF+2,IG)+
1 3.0*FUNKN(IOF+3,IG))
IF(IPRINT.GT.2) THEN
WRITE(6,'(/33H NSSFL2: AVERAGED FLUXES IN GROUP,I5)') IG
WRITE(6,'(1P,10e12.4)') (EVECT(I,IG),I=1,NEL)
WRITE(6,'(/39H NSSFL2: 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 NSSFL2: KEFF=,F11.8,12H OBTAINED IN,I5,11H ITERATIONS)
END
|
