/[MITgcm]/MITgcm/pkg/seaice/seaice_jfnk.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_jfnk.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by mlosch, Tue Nov 6 13:18:14 2012 UTC revision 1.5 by mlosch, Wed Nov 7 09:56:23 2012 UTC
# Line 51  C     myThid :: my Thread Id. number Line 51  C     myThid :: my Thread Id. number
51  #if ( (defined SEAICE_CGRID) && \  #if ( (defined SEAICE_CGRID) && \
52        (defined SEAICE_ALLOW_JFNK) && \        (defined SEAICE_ALLOW_JFNK) && \
53        (defined SEAICE_ALLOW_DYNAMICS) )        (defined SEAICE_ALLOW_DYNAMICS) )
54    C     !FUNCTIONS:
55          LOGICAL  DIFFERENT_MULTIPLE
56          EXTERNAL DIFFERENT_MULTIPLE
57    
58  C     i,j,bi,bj :: loop indices  C     i,j,bi,bj :: loop indices
59        INTEGER i,j,bi,bj        INTEGER i,j,bi,bj
60  C     loop indices  C     loop indices
61        INTEGER newtonIter, newtonIterFail        INTEGER newtonIter
62        INTEGER krylovIter, krylovIterFail        INTEGER krylovIter, krylovFails
63        INTEGER totalKrylovIter        INTEGER totalKrylovItersLoc
64  C     FGMRES flag that indicates what to do next  C     FGMRES flag that determines amount of output messages of fgmres
65          INTEGER iOutFGMRES
66    C     FGMRES flag that indicates what fgmres wants us to do next
67        INTEGER iCode        INTEGER iCode
68        _RL     JFNKresidual, JFNKresidualTile(nSx,nSy)        _RL     JFNKresidual, JFNKresidualTile(nSx,nSy)
69        _RL     JFNKresidualKm1        _RL     JFNKresidualKm1
# Line 86  C     zeta, eta, and DWATN, press Line 91  C     zeta, eta, and DWATN, press
91  CEOP  CEOP
92    
93  C     Initialise  C     Initialise
94        newtonIter      = 0        newtonIter          = 0
95        newtonIterFail  = 0        krylovFails         = 0
96        krylovIterFail  = 0        totalKrylovItersLoc = 0
97        totalKrylovIter = 0        JFNKconverged       = .FALSE.
98        JFNKconverged   = .FALSE.        JFNKtol             = 0. _d 0
99        JFNKtol         = 0. _d 0        JFNKresidual        = 0. _d 0
100        JFNKresidual    = 0. _d 0        JFNKresidualKm1     = 0. _d 0
101        JFNKresidualKm1 = 0. _d 0        FGMRESeps           = 0. _d 0
102        FGMRESeps       = 0. _d 0        recip_deltaT        = 1. _d 0 / SEAICE_deltaTdyn
103        recip_deltaT    = 1. _d 0 / SEAICE_deltaTdyn  
104          iOutFGMRES=0
105    C     iOutFgmres=1 gives a little bit of output
106          IF ( debugLevel.GE.debLevA .AND.
107         &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
108         &     iOutFGMRES=1
109    
110  C      C    
111        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
112         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 194  C     in that routine Line 205  C     in that routine
205         krylovIter    = 0         krylovIter    = 0
206         iCode         = 0         iCode         = 0
207         IF ( debugLevel.GE.debLevA ) THEN           IF ( debugLevel.GE.debLevA ) THEN  
208            _BEGIN_MASTER( myThid )
209          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
210       &       ' S/R SEAICE_JFNK: newtonIter,',       &       ' S/R SEAICE_JFNK: newtonIter,',
211       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',
# Line 201  C     in that routine Line 213  C     in that routine
213       &       JFNKgamma_lin, JFNKresidual       &       JFNKgamma_lin, JFNKresidual
214          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
215       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
216            _END_MASTER( myThid )
217         ENDIF         ENDIF
218  C  C
219         JFNKconverged = JFNKresidual.LT.JFNKtol         JFNKconverged = JFNKresidual.LT.JFNKtol
# Line 221  C Line 234  C
234           CALL SEAICE_FGMRES_DRIVER(           CALL SEAICE_FGMRES_DRIVER(
235       I        uIceRes, vIceRes,       I        uIceRes, vIceRes,
236       U        duIce, dvIce, iCode,       U        duIce, dvIce, iCode,
237       I        FGMRESeps,         I        FGMRESeps, iOutFGMRES,
238       I        newtonIter, krylovIter, myTime, myIter, myThid )       I        newtonIter, krylovIter, myTime, myIter, myThid )
239  C     FGMRES returns iCode either asking for an new preconditioned vector  C     FGMRES returns iCode either asking for an new preconditioned vector
240  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate
# Line 242  C     Compute Jacobian times vector Line 255  C     Compute Jacobian times vector
255           krylovConverged = iCode.EQ.0           krylovConverged = iCode.EQ.0
256  C     End of Krylov iterate  C     End of Krylov iterate
257          ENDDO          ENDDO
258          totalKrylovIter = totalKrylovIter + krylovIter          totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
259  C     some output diagnostics  C     some output diagnostics
260          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
261             _BEGIN_MASTER( myThid )
262           WRITE(msgBuf,'(3(A,I6))')           WRITE(msgBuf,'(3(A,I6))')
263       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,
264       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
265       &        ', Nb. of FGMRES iterations = ', krylovIter       &        ', Nb. of FGMRES iterations = ', krylovIter
266           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
267       &        SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
268             _END_MASTER( myThid )
269          ENDIF          ENDIF
270          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
271           krylovIterFail = krylovIterFail + 1           krylovFails = krylovFails + 1
272          ENDIF          ENDIF
273  C     Update linear solution vector and return to Newton iteration  C     Update linear solution vector and return to Newton iteration
274          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 274  C     Set the stopping criterion for the Line 289  C     Set the stopping criterion for the
289         ENDIF         ENDIF
290  C     end of Newton iterate  C     end of Newton iterate
291        ENDDO        ENDDO
292  C     some output diagnostics  C
293        IF ( debugLevel.GE.debLevA ) THEN  C--   Output diagnostics
294    C
295    C     Count iterations
296          totalJFNKtimeSteps = totalJFNKtimeSteps + 1
297          totalNewtonIters   = totalNewtonIters + newtonIter
298          totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
299  C     Record failure  C     Record failure
300          totalKrylovFails   = totalKrylovFails + krylovFails
301          IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
302           totalNewtonFails = totalNewtonFails + 1
303          ENDIF
304    C     Decide whether it is time to dump and reset the counter
305          IF ( DIFFERENT_MULTIPLE(SEAICE_monFreq,myTime+deltaTClock,
306         &     deltaTClock) ) THEN
307           _BEGIN_MASTER( myThid )
308           WRITE(msgBuf,'(A)')
309         &' // ======================================================='
310           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
311         &      SQUEEZE_RIGHT, myThid )
312           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
313           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
314         &      SQUEEZE_RIGHT, myThid )
315           WRITE(msgBuf,'(A)')
316         &' // ======================================================='
317           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
318         &      SQUEEZE_RIGHT, myThid )
319           WRITE(msgBuf,'(A,I10)')
320         &      ' %JFNK_MON: time step              = ', myIter+1
321           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
322         &      SQUEEZE_RIGHT, myThid )
323           WRITE(msgBuf,'(A,I10)')
324         &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
325           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
326         &      SQUEEZE_RIGHT, myThid )
327           WRITE(msgBuf,'(A,I10)')
328         &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
329           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
330         &      SQUEEZE_RIGHT, myThid )
331           WRITE(msgBuf,'(A,I10)')
332         &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
333           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
334         &      SQUEEZE_RIGHT, myThid )
335           WRITE(msgBuf,'(A,I10)')
336         &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
337           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
338         &      SQUEEZE_RIGHT, myThid )
339           WRITE(msgBuf,'(A,I10)')
340         &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
341           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
342         &      SQUEEZE_RIGHT, myThid )
343           WRITE(msgBuf,'(A)')
344         &' // ======================================================='
345           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
346         &      SQUEEZE_RIGHT, myThid )
347           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
348           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
349         &      SQUEEZE_RIGHT, myThid )
350           WRITE(msgBuf,'(A)')
351         &' // ======================================================='
352           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
353         &      SQUEEZE_RIGHT, myThid )
354           _END_MASTER( myThid )
355    C     reset and start again
356           totalJFNKtimeSteps = 0
357           totalNewtonIters   = 0
358           totalKrylovIters   = 0
359           totalKrylovFails   = 0
360           totalNewtonFails   = 0
361          ENDIF
362    
363    C     Print more debugging information
364          IF ( debugLevel.GE.debLevA ) THEN
365         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
366          newtonIterFail = newtonIterFail + 1          _BEGIN_MASTER( myThid )
367          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
368       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
369       &       myIter       &       myIter+1
370          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
371       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
372            _END_MASTER( myThid )
373         ENDIF         ENDIF
374         IF ( krylovIterFail .GT. 0 ) THEN         IF ( krylovFails .GT. 0 ) THEN
375            _BEGIN_MASTER( myThid )
376          WRITE(msgBuf,'(A,I4,A,I10)')          WRITE(msgBuf,'(A,I4,A,I10)')
377       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
378       &       krylovIterFail, ' times in timestep ', myIter       &       krylovFails, ' times in timestep ', myIter+1
379          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
380       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
381            _END_MASTER( myThid )
382         ENDIF         ENDIF
383         WRITE(msgBuf,'(A,I6)')         _BEGIN_MASTER( myThid )
384           WRITE(msgBuf,'(A,I6,A,I10)')
385       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
386       &      totalKrylovIter       &      totalKrylovItersLoc, ' in timestep ', myIter+1
387          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
388       &       SQUEEZE_RIGHT, myThid )       &      SQUEEZE_RIGHT, myThid )
389                 _END_MASTER( myThid )
390        ENDIF        ENDIF
391    
392  #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */  #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22