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
|
*DECK FLDBHR
SUBROUTINE FLDBHR(IPTRK,IPSYS,LADJ,LL4,ITY,NUN,NGRP,ICL1,ICL2,
1 IMPX,NADI,MAXINR,EPSINR,ITER,TKT,TKB,GRAD1)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Perform thermal (up-scattering) iterations in Bivac.
*
*Copyright:
* Copyright (C) 2002 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
* IPTRK L_TRACK pointer to the tracking information.
* IPSYS L_SYSTEM pointer to system matrices.
* LADJ flag set to .TRUE. for adjoint solution acceleration.
* LL4 order of the system matrices.
* ITY type of solution (2: classical Bivac; 3: Thomas-Raviart).
* NUN number of unknowns in each energy group.
* NGRP number of energy groups.
* ICL1 number of free iretations in one cycle of the up-scattering
* iterations.
* ICL2 number of accelerated up-scattering iterations in one cycle.
* IMPX print parameter (set to 0 for no printing).
* NADI number of inner ADI iterations per outer iteration (used with
* SPN approximations).
* MAXINR maximum number of thermal iterations.
* EPSINR thermal iteration epsilon.
*
*Parameters: input/output
* ITER actual number of thermal iterations.
* TKT CPU time spent to compute the solution of linear systems.
* TKB CPU time spent to compute the bilinear products.
* GRAD1 delta flux for this outer iteration.
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPTRK,IPSYS
INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,NADI,MAXINR,ITER
REAL EPSINR,TKT,TKB,GRAD1(NUN,NGRP)
LOGICAL LADJ
*----
* LOCAL VARIABLES
*----
PARAMETER(NSTATE=40)
INTEGER ISTATE(NSTATE)
CHARACTER TEXT12*12,TEXT3*3
*----
* ALLOCATABLE ARRAYS
*----
REAL, DIMENSION(:), ALLOCATABLE :: WORK2
REAL, DIMENSION(:,:), ALLOCATABLE :: GAR2
REAL, DIMENSION(:,:,:), ALLOCATABLE :: WORK
*----
* SCRATCH STORAGE ALLOCATION
*----
IF(MAXINR.EQ.0) RETURN
ALLOCATE(GAR2(NUN,NGRP),WORK(LL4,NGRP,3),WORK2(LL4))
*
CALL LCMGET(IPSYS,'STATE-VECTOR',ISTATE)
NAN=ISTATE(8)
NCTOT=ICL1+ICL2
IF(ICL2.EQ.0) THEN
NCPTM=NCTOT+1
ELSE
NCPTM=ICL1
ENDIF
DO 15 IGR=1,NGRP
DO 10 I=1,LL4
WORK(I,IGR,1)=0.0
WORK(I,IGR,2)=0.0
WORK(I,IGR,3)=GRAD1(I,IGR)
10 CONTINUE
15 CONTINUE
IGDEB=1
*----
* PERFORM THERMAL (UP-SCATTERING) ITERATIONS
*----
TEXT3='NO '
ITER=2
DO
CALL KDRCPU(TK1)
IF(LADJ) THEN
* ADJOINT SOLUTION
DO 35 IGR=IGDEB,NGRP
WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,IGR,3),
1 GAR2(1,IGR))
DO 30 JGR=1,NGRP
IF(JGR.EQ.IGR) GO TO 30
WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR
CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
IF(ILONG.EQ.0) GO TO 30
CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,JGR,3),WORK2(1))
DO 20 I=1,LL4
GAR2(I,IGR)=GAR2(I,IGR)-WORK2(I)
20 CONTINUE
30 CONTINUE
35 CONTINUE
DO 65 IGR=NGRP,IGDEB,-1
DO 50 JGR=NGRP,IGR+1,-1
WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR
CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
IF(ILONG.EQ.0) GO TO 50
CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,JGR),WORK2(1))
DO 40 I=1,LL4
GAR2(I,IGR)=GAR2(I,IGR)+WORK2(I)
40 CONTINUE
50 CONTINUE
CALL KDRCPU(TK2)
TKB=TKB+(TK2-TK1)
*
CALL KDRCPU(TK1)
WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
IF(ITY.EQ.11) THEN
* SIMPLIFIED PN BIVAC TRACKING.
IF(NAN.EQ.0) CALL XABORT('FLDBHR: SPN-ONLY ALGORITHM.')
CALL FLDBSS(TEXT12,IPTRK,IPSYS,LL4,NBMIX,NAN,GAR2(1,IGR),
1 NADI)
ELSE
CALL MTLDLS(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,IGR))
ENDIF
DO 60 I=1,LL4
WORK(I,IGR,1)=WORK(I,IGR,2)
WORK(I,IGR,2)=WORK(I,IGR,3)
WORK(I,IGR,3)=GRAD1(I,IGR)+(WORK(I,IGR,2)-GAR2(I,IGR))
60 CONTINUE
65 CONTINUE
ELSE
* DIRECT SOLUTION
DO 85 IGR=IGDEB,NGRP
WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,IGR,3),
1 GAR2(1,IGR))
DO 80 JGR=1,NGRP
IF(JGR.EQ.IGR) GO TO 80
WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
IF(ILONG.EQ.0) GO TO 80
CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK(1,JGR,3),WORK2(1))
DO 70 I=1,LL4
GAR2(I,IGR)=GAR2(I,IGR)-WORK2(I)
70 CONTINUE
80 CONTINUE
85 CONTINUE
DO 115 IGR=IGDEB,NGRP
DO 100 JGR=1,IGR-1
WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
IF(ILONG.EQ.0) GO TO 100
CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,JGR),WORK2(1))
DO 90 I=1,LL4
GAR2(I,IGR)=GAR2(I,IGR)+WORK2(I)
90 CONTINUE
100 CONTINUE
CALL KDRCPU(TK2)
TKB=TKB+(TK2-TK1)
*
CALL KDRCPU(TK1)
WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
IF(ITY.EQ.11) THEN
* SIMPLIFIED PN BIVAC TRACKING.
IF(NAN.EQ.0) CALL XABORT('FLDBHR: SPN-ONLY ALGORITHM.')
CALL FLDBSS(TEXT12,IPTRK,IPSYS,LL4,NBMIX,NAN,GAR2(1,IGR),
1 NADI)
ELSE
CALL MTLDLS(TEXT12,IPTRK,IPSYS,LL4,ITY,GAR2(1,IGR))
ENDIF
DO 110 I=1,LL4
WORK(I,IGR,1)=WORK(I,IGR,2)
WORK(I,IGR,2)=WORK(I,IGR,3)
WORK(I,IGR,3)=GRAD1(I,IGR)+(WORK(I,IGR,2)-GAR2(I,IGR))
110 CONTINUE
115 CONTINUE
ENDIF
IF(MOD(ITER-2,NCTOT).GE.NCPTM) THEN
CALL FLD2AC(NGRP,LL4,IGDEB,WORK,ZMU)
ELSE
ZMU=1.0
ENDIF
IGDEBO=IGDEB
DO 130 IGR=IGDEBO,NGRP
GINN=0.0
FINN=0.0
DO 120 I=1,LL4
GINN=MAX(GINN,ABS(WORK(I,IGR,2)-WORK(I,IGR,3)))
FINN=MAX(FINN,ABS(WORK(I,IGR,3)))
120 CONTINUE
GINN=GINN/FINN
IF((GINN.LT.EPSINR).AND.(IGDEB.EQ.IGR)) IGDEB=IGDEB+1
130 CONTINUE
CALL KDRCPU(TK2)
TKT=TKT+(TK2-TK1)
IF(GINN.LT.EPSINR) TEXT3='YES'
IF(IMPX.GT.2) WRITE(6,1000) ITER,GINN,EPSINR,IGDEB,ZMU,TEXT3
IF((GINN.LT.EPSINR).OR.(ITER.EQ.MAXINR)) GO TO 160
ITER=ITER+1
ENDDO
*----
* END OF THERMAL ITERATIONS
*----
160 DO 175 I=1,LL4
DO 170 IGR=1,NGRP
GRAD1(I,IGR)=WORK(I,IGR,3)
170 CONTINUE
175 CONTINUE
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(GAR2,WORK,WORK2)
RETURN
*
1000 FORMAT (10X,3HIN(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2,
1 7H IGDEB=, I13,6H ACCE=,0P,F12.5,12H CONVERGED=,A3)
END
|
