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
|
*DECK GPTMRA
SUBROUTINE GPTMRA(IPTRK,IPSYS,IPFLUP,LADJ,LL4,ITY,NUN,NGRP,ICL1,
1 ICL2,IMPX,NADI,MAXINR,NSTART,MAXX0,EPS2,EPSINR,EVAL,EVECT,ADECT,
2 EASS,SOUR,TKT,TKB,ZNORM,ITER)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Solution of a multigroup fixed source eigenvalue problem for the
* calculation of a gpt solution in Trivac. Use the preconditioned power
* method with GMRES(m) acceleration.
*
*Copyright:
* Copyright (C) 2019 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.
* IPFLUP L_FLUX pointer to the gpt solution
* LADJ flag set to .TRUE. for adjoint solution acceleration.
* LL4 order of the system matrices.
* ITY type of solution (2: classical Trivac; 3: Thomas-Raviart).
* NUN number of unknowns in each energy group.
* NGRP number of energy groups.
* ICL1 number of free up-scattering iterations in one cycle of the
* inverse power method.
* ICL2 number of accelerated up-scattering iterations in one cycle.
* IMPX print parameter (set to 0 for no printing).
* NADI initial number of inner ADI iterations per outer iteration.
* MAXINR maximum number of thermal iterations.
* NSTART restarts the GMRES method every NSTART iterations.
* MAXX0 maximum number of outer iterations
* EPS2 outer iteration convergence criterion
* EPSINR thermal iteration convergence criterion
* EVAL eigenvalue.
* EVECT unknown vector for the non perturbed direct flux
* ADECT unknown vector for the non perturbed adjoint flux
* SOUR fixed source
*
*Parameters: input/output
* EASS solution of the fixed source eigenvalue problem
* TKT CPU time spent to compute the solution of linear systems.
* TKB CPU time spent to compute the bilinear products.
* ZNORM Hotelling deflation accuracy.
* ITER number of iterations.
*
*-----------------------------------------------------------------------
*
USE GANLIB
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPTRK,IPSYS,IPFLUP
LOGICAL LADJ
INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,NADI,MAXINR,NSTART,MAXX0,
1 ITER
REAL EPS2,EPSINR,EVECT(NUN,NGRP),ADECT(NUN,NGRP),EASS(NUN*NGRP),
1 SOUR(NUN,NGRP),TKT,TKB,SDOT
DOUBLE PRECISION EVAL,ZNORM
*----
* LOCAL VARIABLES
*----
PARAMETER (IUNOUT=6)
*----
* ALLOCATABLE ARRAYS
*----
REAL, ALLOCATABLE, DIMENSION(:) :: RR,QQ,VV,GAR1
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: V,H
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: G,C,S,X
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(V(NUN*NGRP,NSTART+1),G(NSTART+1),H(NSTART+1,NSTART+1),
1 C(NSTART+1),S(NSTART+1),X(NUN*NGRP),GAR1(NUN*NGRP))
*----
* GLOBAL GMRES ITERATION.
*----
ALLOCATE(RR(NUN*NGRP),QQ(NUN*NGRP),VV(NUN*NGRP))
EPS1=EPS2*SQRT(SDOT(NUN*NGRP,SOUR,1,SOUR,1))
RHO=1.0E10
ITER=0
NITER=1
NNADI=NADI
DO WHILE((RHO.GT.EPS1).AND.(ITER.LT.MAXX0))
CALL GPTGRA(IPTRK,IPSYS,IPFLUP,LADJ,.TRUE.,LL4,ITY,NUN,NGRP,
1 ICL1,ICL2,IMPX,NNADI,MAXINR,EPSINR,EVAL,EVECT,ADECT,EASS(1),
2 SOUR(1,1),GAR1,JTER0,TKT,TKB,ZNORM,RR)
NITER=NITER+1
DO I=1,NUN*NGRP
X(I)=RR(I)
ENDDO
RHO=SQRT(DDOT(NUN*NGRP,X(1),1,X(1),1))
*----
* TEST FOR TERMINATION ON ENTRY
*----
IF(RHO.LT.EPS1) THEN
DEALLOCATE(VV,QQ,RR)
GO TO 100
ENDIF
*
V(:NUN*NGRP,:NSTART+1)=0.0D0
G(:NSTART+1)=0.0D0
H(:NSTART+1,:NSTART+1)=0.0D0
C(:NSTART+1)=0.0D0
S(:NSTART+1)=0.0D0
G(1)=RHO
DO I=1,NUN*NGRP
V(I,1)=X(I)/RHO
ENDDO
*----
* GMRES(1) ITERATION
*----
K=0
DO WHILE((RHO.GT.EPS1).AND.(K.LT.NSTART).AND.(ITER.LT.MAXX0))
K=K+1
ITER=ITER+1
IF(IMPX.GT.1) WRITE(IUNOUT,300) ITER,RHO,JTER0
DO I=1,NUN*NGRP
VV(I)=REAL(V(I,K))
QQ(I)=0.0
ENDDO
CALL GPTGRA(IPTRK,IPSYS,IPFLUP,LADJ,.TRUE.,LL4,ITY,NUN,NGRP,
1 ICL1,ICL2,IMPX,NNADI,MAXINR,EPSINR,EVAL,EVECT,ADECT,VV(1),
2 QQ(1),GAR1,JTER,TKT,TKB,ZNORM,RR)
IF(JTER.NE.JTER0) CALL XABORT('GPTMRA: INCONSISTENT PRECONDIT'
1 //'IONING IN GMRES.')
NITER=NITER+1
DO I=1,NUN*NGRP
V(I,K+1)=-RR(I)
ENDDO
*----
* MODIFIED GRAM-SCHMIDT
*----
DO J=1,K
HR=DDOT(NUN*NGRP,V(1,J),1,V(1,K+1),1)
H(J,K)=HR
DO I=1,NUN*NGRP
V(I,K+1)=V(I,K+1)-HR*V(I,J)
ENDDO
ENDDO
H(K+1,K)=SQRT(DDOT(NUN*NGRP,V(1,K+1),1,V(1,K+1),1))
*----
* REORTHOGONALIZE
*----
DO J=1,K
HR=DDOT(NUN*NGRP,V(1,J),1,V(1,K+1),1)
H(J,K)=H(J,K)+HR
DO I=1,NUN*NGRP
V(I,K+1)=V(I,K+1)-HR*V(I,J)
ENDDO
ENDDO
H(K+1,K)=SQRT(DDOT(NUN*NGRP,V(1,K+1),1,V(1,K+1),1))
*----
* WATCH OUT FOR HAPPY BREAKDOWN
*----
IF(H(K+1,K).NE.0.0) THEN
DO I=1,NUN*NGRP
V(I,K+1)=V(I,K+1)/H(K+1,K)
ENDDO
ENDIF
*----
* FORM AND STORE THE INFORMATION FOR THE NEW GIVENS ROTATION
*----
DO I=1,K-1
W1=C(I)*H(I,K)-S(I)*H(I+1,K)
W2=S(I)*H(I,K)+C(I)*H(I+1,K)
H(I,K)=W1
H(I+1,K)=W2
ENDDO
ZNU=SQRT(H(K,K)**2+H(K+1,K)**2)
IF(ZNU.NE.0.0) THEN
C(K)=H(K,K)/ZNU
S(K)=-H(K+1,K)/ZNU
H(K,K)=C(K)*H(K,K)-S(K)*H(K+1,K)
H(K+1,K)=0.0D0
W1=C(K)*G(K)-S(K)*G(K+1)
W2=S(K)*G(K)+C(K)*G(K+1)
G(K)=W1
G(K+1)=W2
ENDIF
*----
* UPDATE THE RESIDUAL NORM
*----
RHO=ABS(G(K+1))
ENDDO
*----
* AT THIS POINT EITHER K > NSTART OR RHO < EPS1.
* IT'S TIME TO COMPUTE X AND CYCLE.
*----
DO J=1,K
H(J,K+1)=G(J)
ENDDO
CALL ALSBD(K,1,H,IER,NSTART+1)
IF(IER.NE.0) CALL XABORT('GPTMRA: SINGULAR MATRIX.')
DO I=1,NUN*NGRP
EASS(I)=EASS(I)+REAL(DDOT(K,V(I,1),NUN*NGRP,H(1,K+1),1))
ENDDO
IF(K.EQ.NSTART) THEN
NNADI=NNADI+1
IF(IMPX.NE.0) WRITE (6,310) NNADI
ENDIF
ENDDO
DEALLOCATE(VV,QQ,RR)
*----
* SCRATCH STORAGE DEALLOCATION
*----
100 DEALLOCATE(GAR1,X,S,C,H,G,V)
RETURN
*
300 FORMAT(24H GPTMRA: OUTER ITERATION,I4,10H L2 NORM=,1P,E11.4,
1 28H (NB. OF THERMAL ITERATIONS=,I4,1H))
310 FORMAT(/53H GPTMRA: INCREASING THE NUMBER OF INNER ITERATIONS TO,
1 I3,36H ADI ITERATIONS PER OUTER ITERATION./)
END
|
