1 |
C $Header: $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "SEAICE_OPTIONS.h" |
5 |
#ifdef ALLOW_AUTODIFF |
6 |
# include "AUTODIFF_OPTIONS.h" |
7 |
#endif |
8 |
|
9 |
C-- File seaice_krylov.F: seaice krylov dynamical solver S/R: |
10 |
|
11 |
CBOP |
12 |
C !ROUTINE: SEAICE_KRYLOV |
13 |
C !INTERFACE: |
14 |
SUBROUTINE SEAICE_KRYLOV( myTime, myIter, myThid ) |
15 |
|
16 |
C !DESCRIPTION: \bv |
17 |
C *==========================================================* |
18 |
C | SUBROUTINE SEAICE_KRYLOV |
19 |
C | o Picard solver for ice dynamics using a preconditioned |
20 |
C | KRYLOV (Generalized Minimum RESidual=GMRES) method for |
21 |
C | solving the linearised system following J.-F. Lemieux |
22 |
C | et al., JGR 113, doi:10.1029/2007JC004680, 2008. |
23 |
C *==========================================================* |
24 |
C | written by Martin Losch, Jan 2016 |
25 |
C *==========================================================* |
26 |
C \ev |
27 |
|
28 |
C !USES: |
29 |
IMPLICIT NONE |
30 |
|
31 |
C === Global variables === |
32 |
#include "SIZE.h" |
33 |
#include "EEPARAMS.h" |
34 |
#include "PARAMS.h" |
35 |
#include "DYNVARS.h" |
36 |
#include "GRID.h" |
37 |
#include "SEAICE_SIZE.h" |
38 |
#include "SEAICE_PARAMS.h" |
39 |
#include "SEAICE.h" |
40 |
|
41 |
C !INPUT/OUTPUT PARAMETERS: |
42 |
C === Routine arguments === |
43 |
C myTime :: Simulation time |
44 |
C myIter :: Simulation timestep number |
45 |
C myThid :: my Thread Id. number |
46 |
_RL myTime |
47 |
INTEGER myIter |
48 |
INTEGER myThid |
49 |
|
50 |
#ifdef SEAICE_ALLOW_KRYLOV |
51 |
C !FUNCTIONS: |
52 |
LOGICAL DIFFERENT_MULTIPLE |
53 |
EXTERNAL DIFFERENT_MULTIPLE |
54 |
|
55 |
C !LOCAL VARIABLES: |
56 |
C === Local variables === |
57 |
C i,j,bi,bj :: loop indices |
58 |
INTEGER i,j,bi,bj |
59 |
C loop indices |
60 |
INTEGER picardIter |
61 |
INTEGER krylovIter, krylovFails |
62 |
INTEGER krylovIterMax, picardIterMax |
63 |
INTEGER totalKrylovItersLoc, totalPicardItersLoc |
64 |
C FGMRES parameters |
65 |
C im :: size of Krylov space |
66 |
C ifgmres :: interation counter |
67 |
INTEGER im |
68 |
PARAMETER ( im = 50 ) |
69 |
INTEGER ifgmres |
70 |
C FGMRES flag that determines amount of output messages of fgmres |
71 |
INTEGER iOutFGMRES |
72 |
C FGMRES flag that indicates what fgmres wants us to do next |
73 |
INTEGER iCode |
74 |
_RL picardResidual |
75 |
_RL picardResidualKm1 |
76 |
C parameters to compute convergence criterion |
77 |
_RL krylovLinTol |
78 |
_RL FGMRESeps |
79 |
_RL picardTol |
80 |
C backward differences extrapolation factors |
81 |
_RL bdfFac, bdfAlpha |
82 |
C |
83 |
_RL recip_deltaT |
84 |
LOGICAL picardConverged, krylovConverged |
85 |
LOGICAL writeNow |
86 |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
87 |
|
88 |
_RL uIceLin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
89 |
_RL vIceLin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
90 |
C extra time level required for backward difference time stepping |
91 |
_RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
92 |
_RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
93 |
C u/vWork :: work arrays |
94 |
_RL uWork (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
95 |
_RL vWork (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
96 |
C u/vIceLHS :: left hand side of momentum equation (A*x) |
97 |
_RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
98 |
_RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
99 |
C u/vIceRHS :: right hand side of momentum equation (b) |
100 |
_RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
101 |
_RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
102 |
C helper array |
103 |
_RL resTmp (nVec,1,nSx,nSy) |
104 |
C work arrays |
105 |
_RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy) |
106 |
_RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy) |
107 |
_RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy) |
108 |
CEOP |
109 |
|
110 |
C Initialise |
111 |
picardIter = 0 |
112 |
krylovFails = 0 |
113 |
totalKrylovItersLoc = 0 |
114 |
picardConverged = .FALSE. |
115 |
picardTol = 0. _d 0 |
116 |
picardResidual = 0. _d 0 |
117 |
picardResidualKm1 = 0. _d 0 |
118 |
FGMRESeps = 0. _d 0 |
119 |
recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn |
120 |
|
121 |
krylovIterMax = SEAICElinearIterMax |
122 |
picardIterMax = SEAICEnonLinIterMax |
123 |
IF ( SEAICEusePicardAsPrecon ) THEN |
124 |
krylovIterMax = SEAICEpreconlinIter |
125 |
picardIterMax = SEAICEpreconNL_Iter |
126 |
ENDIF |
127 |
|
128 |
iOutFGMRES=0 |
129 |
C with iOutFgmres=1, seaice_fgmres prints the residual at each iteration |
130 |
IF ( debugLevel.GE.debLevC .AND. |
131 |
& .NOT.SEAICEusePicardAsPrecon .AND. |
132 |
& DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) ) |
133 |
& iOutFGMRES=1 |
134 |
|
135 |
C backward difference extrapolation factors |
136 |
bdfFac = 0. _d 0 |
137 |
IF ( SEAICEuseBDF2 ) THEN |
138 |
IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN |
139 |
bdfFac = 0. _d 0 |
140 |
ELSE |
141 |
bdfFac = 0.5 _d 0 |
142 |
ENDIF |
143 |
ENDIF |
144 |
bdfAlpha = 1. _d 0 + bdfFac |
145 |
|
146 |
DO bj=myByLo(myThid),myByHi(myThid) |
147 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
148 |
DO J=1-OLy,sNy+OLy |
149 |
DO I=1-OLx,sNx+OLx |
150 |
uIceLHS(I,J,bi,bj) = 0. _d 0 |
151 |
vIceLHS(I,J,bi,bj) = 0. _d 0 |
152 |
uIceRHS(I,J,bi,bj) = 0. _d 0 |
153 |
vIceRHS(I,J,bi,bj) = 0. _d 0 |
154 |
ENDDO |
155 |
ENDDO |
156 |
C cycle ice velocities |
157 |
DO J=1-OLy,sNy+OLy |
158 |
DO I=1-OLx,sNx+OLx |
159 |
duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha |
160 |
& + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac |
161 |
dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha |
162 |
& + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac |
163 |
uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj) |
164 |
vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj) |
165 |
uIceLin(I,J,bi,bj) = uIce(I,J,bi,bj) |
166 |
vIceLin(I,J,bi,bj) = vIce(I,J,bi,bj) |
167 |
ENDDO |
168 |
ENDDO |
169 |
C Compute things that do no change during the OL iteration: |
170 |
C sea-surface tilt and wind stress: |
171 |
C FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT |
172 |
DO J=1-OLy,sNy+OLy |
173 |
DO I=1-OLx,sNx+OLx |
174 |
FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj) |
175 |
& + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT |
176 |
FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj) |
177 |
& + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT |
178 |
ENDDO |
179 |
ENDDO |
180 |
CML ENDIF |
181 |
ENDDO |
182 |
ENDDO |
183 |
C Start nonlinear Picard iteration: outer loop iteration |
184 |
DO WHILE ( picardIter.LT.picardIterMax .AND. |
185 |
& .NOT.picardConverged ) |
186 |
picardIter = picardIter + 1 |
187 |
C smooth ice velocities in time for computation of |
188 |
C the non-linear drag coefficents |
189 |
DO bj=myByLo(myThid),myByHi(myThid) |
190 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
191 |
DO j=1-OLy,sNy+OLy |
192 |
DO i=1-OLx,sNx+OLx |
193 |
uIceLin(I,J,bi,bj) = 0.5 _d 0 * |
194 |
& (uIce(I,J,bi,bj) + uIceLin(I,J,bi,bj)) |
195 |
vIceLin(I,J,bi,bj) = 0.5 _d 0 * |
196 |
& (vIce(I,J,bi,bj) + vIceLin(I,J,bi,bj)) |
197 |
ENDDO |
198 |
ENDDO |
199 |
ENDDO |
200 |
ENDDO |
201 |
C u/vIce have changed in Picard iteration so that new drag |
202 |
C coefficients and viscosities are required (that will not change in |
203 |
C the Krylov iteration) |
204 |
CALL SEAICE_OCEANDRAG_COEFFS( |
205 |
I uIceLin, vIceLin, |
206 |
O DWATN, |
207 |
I 0, myTime, myIter, myThid ) |
208 |
CALL SEAICE_CALC_STRAINRATES( |
209 |
I uIceLin, vIceLin, |
210 |
O e11, e22, e12, |
211 |
I 0, myTime, myIter, myThid ) |
212 |
CALL SEAICE_CALC_VISCOSITIES( |
213 |
I e11, e22, e12, zMin, zMax, hEffM, press0, tensileStrength, |
214 |
O eta, etaZ, zeta, zetaZ, press, deltaC, |
215 |
I 0, myTime, myIter, myThid ) |
216 |
C compute rhs that does not change during Krylov iteration |
217 |
CALL SEAICE_CALC_RHS( |
218 |
O uIceRHS, vIceRHS, |
219 |
I picardIter, 0, myTime, myIter, myThid ) |
220 |
C compute rhs for initial residual |
221 |
CALL SEAICE_CALC_LHS( |
222 |
I uIce, vIce, |
223 |
O uIceLHS, vIceLHS, |
224 |
I picardIter, myTime, myIter, myThid ) |
225 |
C Calculate the residual |
226 |
DO bj=myByLo(myThid),myByHi(myThid) |
227 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
228 |
DO J=1,sNy |
229 |
DO I=1,sNx |
230 |
uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - uIceRHS(I,J,bi,bj) |
231 |
vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - vIceRHS(I,J,bi,bj) |
232 |
C save u/vIceLin as k-2nd step for linearization (does not work properly) |
233 |
CML uIceLin(I,J,bi,bj) = uIce(I,J,bi,bj) |
234 |
CML vIceLin(I,J,bi,bj) = vIce(I,J,bi,bj) |
235 |
ENDDO |
236 |
ENDDO |
237 |
ENDDO |
238 |
ENDDO |
239 |
C Important: Compute the norm of the residual using the same scalar |
240 |
C product that SEAICE_FGMRES does |
241 |
CALL SEAICE_MAP2VEC(nVec,uIceLHS,vIceLHS,resTmp,.TRUE.,myThid) |
242 |
CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp, |
243 |
& picardResidual,myThid) |
244 |
picardResidual = SQRT(picardResidual) |
245 |
|
246 |
C compute convergence criterion for linear preconditioned FGMRES |
247 |
krylovLinTol = JFNKgamma_lin_max |
248 |
C the best method is still not clear to me |
249 |
C this is described in Lemieux et al 2008 |
250 |
CML IF ( picardIter .EQ. 1 ) krylovLinTol = 1./10. |
251 |
CML IF ( picardIter .EQ. 2 ) krylovLinTol = 1./20. |
252 |
CML IF ( picardIter .EQ. 3 ) krylovLinTol = 1./20. |
253 |
CML IF ( picardIter .EQ. 4 ) krylovLinTol = 1./30. |
254 |
CML IF ( picardIter .EQ. 5 ) krylovLinTol = 1./30. |
255 |
CML IF ( picardIter .EQ. 6 ) krylovLinTol = 1./30. |
256 |
CML IF ( picardIter .EQ. 7 ) krylovLinTol = 1./30. |
257 |
CML IF ( picardIter .EQ. 8 ) krylovLinTol = 1./40. |
258 |
CML IF ( picardIter .EQ. 9 ) krylovLinTol = 1./50. |
259 |
CML IF ( picardIter .GT. 9 ) krylovLinTol = 1./80. |
260 |
C this is used with the JFNK solver, but the Picard-Krylov solver |
261 |
C converges too slowly for this scheme |
262 |
CML IF ( picardIter.GT.1.AND.picardResidual.LT.JFNKres_t ) THEN |
263 |
CMLC Eisenstat and Walker (1996), eq.(2.6) |
264 |
CML krylovLinTol = SEAICE_JFNKphi |
265 |
CML & *( picardResidual/picardResidualKm1 )**SEAICE_JFNKalpha |
266 |
CML krylovLinTol = min(JFNKgamma_lin_max, krylovLinTol) |
267 |
CML krylovLinTol = max(JFNKgamma_lin_min, krylovLinTol) |
268 |
CML ENDIF |
269 |
CML krylovLinTol = 1. _d -1 |
270 |
C save the residual for the next iteration |
271 |
picardResidualKm1 = picardResidual |
272 |
|
273 |
C The Krylov iteration uses FGMRES, the preconditioner is LSOR |
274 |
C for now. The code is adapted from SEAICE_LSR, but heavily stripped |
275 |
C down. |
276 |
C krylovIter is mapped into "its" in seaice_fgmres and is incremented |
277 |
C in that routine |
278 |
krylovIter = 0 |
279 |
iCode = 0 |
280 |
|
281 |
picardConverged = picardResidual.LT.picardTol |
282 |
|
283 |
C do Krylov loop only if convergence is not reached |
284 |
|
285 |
IF ( .NOT.picardConverged ) THEN |
286 |
|
287 |
C start Krylov iteration (FGMRES) |
288 |
|
289 |
krylovConverged = .FALSE. |
290 |
FGMRESeps = krylovLinTol * picardResidual |
291 |
CALL SEAICE_MAP2VEC(nVec,uIce,vIce,sol,.TRUE.,myThid) |
292 |
CALL SEAICE_MAP2VEC(nVec,uIceRHS,vIceRHS,rhs,.TRUE.,myThid) |
293 |
DO bj=myByLo(myThid),myByHi(myThid) |
294 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
295 |
DO j=1-OLy,sNy+OLy |
296 |
DO i=1-OLx,sNx+OLx |
297 |
uWork(i,j,bi,bj) = 0. _d 0 |
298 |
vWork(i,j,bi,bj) = 0. _d 0 |
299 |
ENDDO |
300 |
ENDDO |
301 |
ENDDO |
302 |
ENDDO |
303 |
DO WHILE ( .NOT.krylovConverged ) |
304 |
C solution vector sol = u/vIce |
305 |
C residual vector (rhs) Fu = u/vIceRHS |
306 |
C output work vectors wk1, -> input work vector wk2 |
307 |
|
308 |
C map results to wk2 |
309 |
CALL SEAICE_MAP2VEC(nVec,uWork,vWork,wk2,.TRUE.,myThid) |
310 |
|
311 |
CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter, |
312 |
U vv,w,wk1,wk2, |
313 |
I FGMRESeps,krylovIterMax,iOutFGMRES, |
314 |
U iCode, |
315 |
I myThid) |
316 |
C |
317 |
IF ( iCode .EQ. 0 ) THEN |
318 |
C map sol(ution) vector to u/vIce |
319 |
CALL SEAICE_MAP2VEC(nVec,uIce,vIce,sol,.FALSE.,myThid) |
320 |
CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid) |
321 |
ELSE |
322 |
C map work vector to du/vIce to either compute a preconditioner |
323 |
C solution (wk1=rhs) or a matrix times wk1 |
324 |
CALL SEAICE_MAP2VEC(nVec,uWork,vWork,wk1,.FALSE.,myThid) |
325 |
CALL EXCH_UV_XY_RL( uWork, vWork,.TRUE.,myThid) |
326 |
ENDIF |
327 |
|
328 |
C FGMRES returns iCode either asking for an new preconditioned vector |
329 |
C or product of matrix times vector. For iCode = 0, terminate |
330 |
C iteration |
331 |
IF (iCode.EQ.1) THEN |
332 |
C Call preconditioner |
333 |
IF ( SEAICEpreconLinIter .GT. 0 ) |
334 |
& CALL SEAICE_PRECONDITIONER( |
335 |
U uWork, vWork, |
336 |
I zeta, eta, etaZ, zetaZ, dwatn, |
337 |
I picardIter, krylovIter, myTime, myIter, myThid ) |
338 |
ELSEIF (iCode.GE.2) THEN |
339 |
C Compute lhs of equations (A*x) |
340 |
CALL SEAICE_CALC_STRAINRATES( |
341 |
I uWork, vWork, |
342 |
O e11, e22, e12, |
343 |
I krylovIter, myTime, myIter, myThid ) |
344 |
CALL SEAICE_CALC_LHS( |
345 |
I uWork, vWork, |
346 |
O uIceLHS, vIceLHS, |
347 |
I picardIter, myTime, myIter, myThid ) |
348 |
DO bj=myByLo(myThid),myByHi(myThid) |
349 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
350 |
DO j=1-OLy,sNy+OLy |
351 |
DO i=1-OLx,sNx+OLx |
352 |
uWork(i,j,bi,bj) = uIceLHS(i,j,bi,bj) |
353 |
vWork(i,j,bi,bj) = vIceLHS(i,j,bi,bj) |
354 |
ENDDO |
355 |
ENDDO |
356 |
ENDDO |
357 |
ENDDO |
358 |
ENDIF |
359 |
krylovConverged = iCode.EQ.0 |
360 |
C End of Krylov iterate |
361 |
ENDDO |
362 |
totalKrylovItersLoc = totalKrylovItersLoc + krylovIter |
363 |
C some output diagnostics |
364 |
IF ( debugLevel.GE.debLevA |
365 |
& .AND. .NOT.SEAICEusePicardAsPrecon ) THEN |
366 |
_BEGIN_MASTER( myThid ) |
367 |
totalPicardItersLoc = |
368 |
& picardIterMax*(myIter-nIter0)+picardIter |
369 |
WRITE(msgBuf,'(2A,2(1XI6),2E12.5)') |
370 |
& ' S/R SEAICE_KRYLOV: Picard iterate / total, ', |
371 |
& 'KRYLOVgamma_lin, initial norm = ', |
372 |
& picardIter, totalPicardItersLoc, |
373 |
& krylovLinTol,picardResidual |
374 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
375 |
& SQUEEZE_RIGHT, myThid ) |
376 |
WRITE(msgBuf,'(3(A,I6))') |
377 |
& ' S/R SEAICE_KRYLOV: Picard iterate / total = ', |
378 |
& picardIter, ' / ', totalPicardItersLoc, |
379 |
& ', Nb. of FGMRES iterations = ', krylovIter |
380 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
381 |
& SQUEEZE_RIGHT, myThid ) |
382 |
_END_MASTER( myThid ) |
383 |
ENDIF |
384 |
IF ( krylovIter.EQ.krylovIterMax ) THEN |
385 |
krylovFails = krylovFails + 1 |
386 |
ENDIF |
387 |
C Set the stopping criterion for the Picard iteration and the |
388 |
C criterion for the transition from accurate to approximate FGMRES |
389 |
IF ( picardIter .EQ. 1 ) THEN |
390 |
picardTol=SEAICEnonLinTol*picardResidual |
391 |
IF ( JFNKres_tFac .NE. UNSET_RL ) |
392 |
& JFNKres_t = picardResidual * JFNKres_tFac |
393 |
ENDIF |
394 |
ENDIF |
395 |
C end of Picard iterate |
396 |
ENDDO |
397 |
|
398 |
C-- Output diagnostics |
399 |
|
400 |
IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN |
401 |
C Count iterations |
402 |
totalJFNKtimeSteps = totalJFNKtimeSteps + 1 |
403 |
totalNewtonIters = totalNewtonIters + picardIter |
404 |
totalKrylovIters = totalKrylovIters + totalKrylovItersLoc |
405 |
C Record failure |
406 |
totalKrylovFails = totalKrylovFails + krylovFails |
407 |
IF ( picardIter .EQ. picardIterMax ) THEN |
408 |
totalNewtonfails = totalNewtonfails + 1 |
409 |
ENDIF |
410 |
ENDIF |
411 |
C Decide whether it is time to dump and reset the counter |
412 |
writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq, |
413 |
& myTime+deltaTClock, deltaTClock) |
414 |
#ifdef ALLOW_CAL |
415 |
IF ( useCAL ) THEN |
416 |
CALL CAL_TIME2DUMP( |
417 |
I zeroRL, SEAICE_monFreq, deltaTClock, |
418 |
U writeNow, |
419 |
I myTime+deltaTclock, myIter+1, myThid ) |
420 |
ENDIF |
421 |
#endif |
422 |
IF ( writeNow ) THEN |
423 |
_BEGIN_MASTER( myThid ) |
424 |
WRITE(msgBuf,'(A)') |
425 |
&' // =======================================================' |
426 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
427 |
& SQUEEZE_RIGHT, myThid ) |
428 |
WRITE(msgBuf,'(A)') ' // Begin KRYLOV statistics' |
429 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
430 |
& SQUEEZE_RIGHT, myThid ) |
431 |
WRITE(msgBuf,'(A)') |
432 |
&' // =======================================================' |
433 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
434 |
& SQUEEZE_RIGHT, myThid ) |
435 |
WRITE(msgBuf,'(A,I10)') |
436 |
& ' %KRYLOV_MON: time step = ', myIter+1 |
437 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
438 |
& SQUEEZE_RIGHT, myThid ) |
439 |
WRITE(msgBuf,'(A,I10)') |
440 |
& ' %KRYLOV_MON: Nb. of time steps = ', |
441 |
& totalJFNKtimeSteps |
442 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
443 |
& SQUEEZE_RIGHT, myThid ) |
444 |
WRITE(msgBuf,'(A,I10)') |
445 |
& ' %KRYLOV_MON: Nb. of Picard steps = ', totalNewtonIters |
446 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
447 |
& SQUEEZE_RIGHT, myThid ) |
448 |
WRITE(msgBuf,'(A,I10)') |
449 |
& ' %KRYLOV_MON: Nb. of Krylov steps = ', totalKrylovIters |
450 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
451 |
& SQUEEZE_RIGHT, myThid ) |
452 |
WRITE(msgBuf,'(A,I10)') |
453 |
& ' %KRYLOV_MON: Nb. of Picard failures = ', totalNewtonfails |
454 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
455 |
& SQUEEZE_RIGHT, myThid ) |
456 |
WRITE(msgBuf,'(A,I10)') |
457 |
& ' %KRYLOV_MON: Nb. of Krylov failures = ', totalKrylovFails |
458 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
459 |
& SQUEEZE_RIGHT, myThid ) |
460 |
WRITE(msgBuf,'(A)') |
461 |
&' // =======================================================' |
462 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
463 |
& SQUEEZE_RIGHT, myThid ) |
464 |
WRITE(msgBuf,'(A)') ' // End KRYLOV statistics' |
465 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
466 |
& SQUEEZE_RIGHT, myThid ) |
467 |
WRITE(msgBuf,'(A)') |
468 |
&' // =======================================================' |
469 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
470 |
& SQUEEZE_RIGHT, myThid ) |
471 |
_END_MASTER( myThid ) |
472 |
C reset and start again |
473 |
totalJFNKtimeSteps = 0 |
474 |
totalNewtonIters = 0 |
475 |
totalKrylovIters = 0 |
476 |
totalKrylovFails = 0 |
477 |
totalNewtonfails = 0 |
478 |
ENDIF |
479 |
|
480 |
C Print more debugging information |
481 |
IF ( debugLevel.GE.debLevA |
482 |
& .AND. .NOT.SEAICEusePicardAsPrecon ) THEN |
483 |
IF ( picardIter .EQ. picardIterMax ) THEN |
484 |
_BEGIN_MASTER( myThid ) |
485 |
WRITE(msgBuf,'(A,I10)') |
486 |
& ' S/R SEAICE_KRYLOV: Solver did not converge in timestep ', |
487 |
& myIter+1 |
488 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
489 |
& SQUEEZE_RIGHT, myThid ) |
490 |
_END_MASTER( myThid ) |
491 |
ENDIF |
492 |
IF ( krylovFails .GT. 0 ) THEN |
493 |
_BEGIN_MASTER( myThid ) |
494 |
WRITE(msgBuf,'(A,I4,A,I10)') |
495 |
& ' S/R SEAICE_KRYLOV: FGMRES did not converge ', |
496 |
& krylovFails, ' times in timestep ', myIter+1 |
497 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
498 |
& SQUEEZE_RIGHT, myThid ) |
499 |
_END_MASTER( myThid ) |
500 |
ENDIF |
501 |
_BEGIN_MASTER( myThid ) |
502 |
WRITE(msgBuf,'(A,I6,A,I10)') |
503 |
& ' S/R SEAICE_KRYLOV: Total number FGMRES iterations = ', |
504 |
& totalKrylovItersLoc, ' in timestep ', myIter+1 |
505 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
506 |
& SQUEEZE_RIGHT, myThid ) |
507 |
_END_MASTER( myThid ) |
508 |
ENDIF |
509 |
|
510 |
#endif /* SEAICE_ALLOW_KRYLOV */ |
511 |
|
512 |
RETURN |
513 |
END |