/[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.3 by mlosch, Tue Nov 6 12:53:14 2012 UTC revision 1.9 by mlosch, Mon Nov 12 09:46:38 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 69  C     parameters to compute convergence Line 74  C     parameters to compute convergence
74  C      C    
75        _RL     recip_deltaT        _RL     recip_deltaT
76        LOGICAL JFNKconverged, krylovConverged        LOGICAL JFNKconverged, krylovConverged
77          LOGICAL writeNow
78        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
79  C  C
80  C     u/vIceRes :: residual of sea-ice momentum equations  C     u/vIceRes :: residual of sea-ice momentum equations
# Line 81  C     precomputed (= constant per Newton Line 87  C     precomputed (= constant per Newton
87  C     zeta, eta, and DWATN, press  C     zeta, eta, and DWATN, press
88        _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
89        _RL etaPre  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL etaPre  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
90          _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
91        _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
92        _RL pressPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL pressPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
93  CEOP  CEOP
94    
95  C     Initialise  C     Initialise
96        newtonIter      = 0        newtonIter          = 0
97        newtonIterFail  = 0        krylovFails         = 0
98        krylovIterFail  = 0        totalKrylovItersLoc = 0
99        totalKrylovIter = 0        JFNKconverged       = .FALSE.
100        JFNKconverged   = .FALSE.        JFNKtol             = 0. _d 0
101        JFNKtol         = 0. _d 0        JFNKresidual        = 0. _d 0
102        JFNKresidual    = 0. _d 0        JFNKresidualKm1     = 0. _d 0
103        JFNKresidualKm1 = 0. _d 0        FGMRESeps           = 0. _d 0
104        FGMRESeps       = 0. _d 0        recip_deltaT        = 1. _d 0 / SEAICE_deltaTdyn
105        recip_deltaT    = 1. _d 0 / SEAICE_deltaTdyn  
106          iOutFGMRES=0
107    C     iOutFgmres=1 gives a little bit of output
108          IF ( debugLevel.GE.debLevA .AND.
109         &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
110         &     iOutFGMRES=1
111    
112  C      C    
113        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
114         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 141  C     constant for the preconditioner Line 154  C     constant for the preconditioner
154            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
155              zetaPre(I,J,bi,bj) =  zeta(I,J,bi,bj)              zetaPre(I,J,bi,bj) =  zeta(I,J,bi,bj)
156               etaPre(I,J,bi,bj) =   eta(I,J,bi,bj)               etaPre(I,J,bi,bj) =   eta(I,J,bi,bj)
157                etaZPre(I,J,bi,bj) =  etaZ(I,J,bi,bj)
158              dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)              dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)
159             pressPre(I,J,bi,bj) = press(I,J,bi,bj)             pressPre(I,J,bi,bj) = press(I,J,bi,bj)
160            ENDDO            ENDDO
# Line 194  C     in that routine Line 208  C     in that routine
208         krylovIter    = 0         krylovIter    = 0
209         iCode         = 0         iCode         = 0
210         IF ( debugLevel.GE.debLevA ) THEN           IF ( debugLevel.GE.debLevA ) THEN  
211            _BEGIN_MASTER( myThid )
212          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
213       &       ' S/R SEAICE_JFNK: newtonIter,',       &       ' S/R SEAICE_JFNK: newtonIter,',
214       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',
# Line 201  C     in that routine Line 216  C     in that routine
216       &       JFNKgamma_lin, JFNKresidual       &       JFNKgamma_lin, JFNKresidual
217          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
218       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
219            _END_MASTER( myThid )
220         ENDIF         ENDIF
221  C  C
222         JFNKconverged = JFNKresidual.LT.JFNKtol         JFNKconverged = JFNKresidual.LT.JFNKtol
# Line 221  C Line 237  C
237           CALL SEAICE_FGMRES_DRIVER(           CALL SEAICE_FGMRES_DRIVER(
238       I        uIceRes, vIceRes,       I        uIceRes, vIceRes,
239       U        duIce, dvIce, iCode,       U        duIce, dvIce, iCode,
240       I        FGMRESeps,         I        FGMRESeps, iOutFGMRES,
241       I        newtonIter, krylovIter, myTime, myIter, myThid )       I        newtonIter, krylovIter, myTime, myIter, myThid )
242  C     FGMRES returns iCode either asking for an new preconditioned vector  C     FGMRES returns iCode either asking for an new preconditioned vector
243  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate
244  C     iteration  C     iteration
245           IF (iCode.EQ.1) THEN           IF (iCode.EQ.1) THEN
246  C     Call preconditioner  C     Call preconditioner
247            CALL SEAICE_PRECONDITIONER(            IF ( SOLV_MAX_ITERS .GT. 0 )
248         &         CALL SEAICE_PRECONDITIONER(
249       U         duIce, dvIce,       U         duIce, dvIce,
250       I         zetaPre, etaPre, dwatPre, pressPre,       I         zetaPre, etaPre, etaZpre, dwatPre, pressPre,
251       I         newtonIter, krylovIter, myTime, myIter, myThid )       I         newtonIter, krylovIter, myTime, myIter, myThid )
252           ELSEIF (iCode.GE.2) THEN           ELSEIF (iCode.GE.2) THEN
253  C     Compute Jacobian times vector  C     Compute Jacobian times vector
# Line 242  C     Compute Jacobian times vector Line 259  C     Compute Jacobian times vector
259           krylovConverged = iCode.EQ.0           krylovConverged = iCode.EQ.0
260  C     End of Krylov iterate  C     End of Krylov iterate
261          ENDDO          ENDDO
262          totalKrylovIter = totalKrylovIter + krylovIter          totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
263  C     some output diagnostics  C     some output diagnostics
264          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
265             _BEGIN_MASTER( myThid )
266           WRITE(msgBuf,'(3(A,I6))')           WRITE(msgBuf,'(3(A,I6))')
267       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,
268       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
269       &        ', Nb. of FGMRES iterations = ', krylovIter       &        ', Nb. of FGMRES iterations = ', krylovIter
270           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
271       &        SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
272             _END_MASTER( myThid )
273          ENDIF          ENDIF
274          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
275           krylovIterFail = krylovIterFail + 1           krylovFails = krylovFails + 1
276          ENDIF          ENDIF
277  C     Update linear solution vector and return to Newton iteration  C     Update linear solution vector and return to Newton iteration
278          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 262  C     Update linear solution vector and Line 281  C     Update linear solution vector and
281             DO I=1-Olx,sNx+Olx             DO I=1-Olx,sNx+Olx
282              uIce(I,J,bi,bj) = uIce(I,J,bi,bj)+duIce(I,J,bi,bj)              uIce(I,J,bi,bj) = uIce(I,J,bi,bj)+duIce(I,J,bi,bj)
283              vIce(I,J,bi,bj) = vIce(I,J,bi,bj)+dvIce(I,J,bi,bj)              vIce(I,J,bi,bj) = vIce(I,J,bi,bj)+dvIce(I,J,bi,bj)
284    C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
285                duIce(I,J,bi,bj)= 0. _d 0
286                dvIce(I,J,bi,bj)= 0. _d 0
287             ENDDO             ENDDO
288            ENDDO            ENDDO
289           ENDDO           ENDDO
# Line 271  C     Set the stopping criterion for the Line 293  C     Set the stopping criterion for the
293         ENDIF         ENDIF
294  C     end of Newton iterate  C     end of Newton iterate
295        ENDDO        ENDDO
296  C     some output diagnostics  C
297        IF ( debugLevel.GE.debLevA ) THEN  C--   Output diagnostics
298    C
299          IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
300    C     Count iterations
301           totalJFNKtimeSteps = totalJFNKtimeSteps + 1
302           totalNewtonIters   = totalNewtonIters + newtonIter
303           totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
304  C     Record failure  C     Record failure
305           totalKrylovFails   = totalKrylovFails + krylovFails
306           IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
307            totalNewtonFails = totalNewtonFails + 1
308           ENDIF
309          ENDIF
310    C     Decide whether it is time to dump and reset the counter
311          writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
312         &     myTime+deltaTClock, deltaTClock)
313    #ifdef ALLOW_CAL
314          IF ( useCAL ) THEN
315           CALL CAL_TIME2DUMP(
316         I      zeroRL, SEAICE_monFreq,  deltaTClock,
317         U      writeNow,
318         I      myTime+deltaTclock, myIter+1, myThid )
319          ENDIF
320    #endif
321          IF ( writeNow ) THEN
322           _BEGIN_MASTER( myThid )
323           WRITE(msgBuf,'(A)')
324         &' // ======================================================='
325           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
326         &      SQUEEZE_RIGHT, myThid )
327           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
328           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
329         &      SQUEEZE_RIGHT, myThid )
330           WRITE(msgBuf,'(A)')
331         &' // ======================================================='
332           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
333         &      SQUEEZE_RIGHT, myThid )
334           WRITE(msgBuf,'(A,I10)')
335         &      ' %JFNK_MON: time step              = ', myIter+1
336           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
337         &      SQUEEZE_RIGHT, myThid )
338           WRITE(msgBuf,'(A,I10)')
339         &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
340           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
341         &      SQUEEZE_RIGHT, myThid )
342           WRITE(msgBuf,'(A,I10)')
343         &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
344           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
345         &      SQUEEZE_RIGHT, myThid )
346           WRITE(msgBuf,'(A,I10)')
347         &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
348           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
349         &      SQUEEZE_RIGHT, myThid )
350           WRITE(msgBuf,'(A,I10)')
351         &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
352           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
353         &      SQUEEZE_RIGHT, myThid )
354           WRITE(msgBuf,'(A,I10)')
355         &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
356           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
357         &      SQUEEZE_RIGHT, myThid )
358           WRITE(msgBuf,'(A)')
359         &' // ======================================================='
360           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
361         &      SQUEEZE_RIGHT, myThid )
362           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
363           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
364         &      SQUEEZE_RIGHT, myThid )
365           WRITE(msgBuf,'(A)')
366         &' // ======================================================='
367           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
368         &      SQUEEZE_RIGHT, myThid )
369           _END_MASTER( myThid )
370    C     reset and start again
371           totalJFNKtimeSteps = 0
372           totalNewtonIters   = 0
373           totalKrylovIters   = 0
374           totalKrylovFails   = 0
375           totalNewtonFails   = 0
376          ENDIF
377    
378    C     Print more debugging information
379          IF ( debugLevel.GE.debLevA ) THEN
380         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
381          newtonIterFail = newtonIterFail + 1          _BEGIN_MASTER( myThid )
382          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
383       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
384       &       myIter       &       myIter+1
385          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
386       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
387            _END_MASTER( myThid )
388         ENDIF         ENDIF
389         IF ( krylovIterFail .GT. 0 ) THEN         IF ( krylovFails .GT. 0 ) THEN
390            _BEGIN_MASTER( myThid )
391          WRITE(msgBuf,'(A,I4,A,I10)')          WRITE(msgBuf,'(A,I4,A,I10)')
392       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
393       &       krylovIterFail, ' times in timestep ', myIter       &       krylovFails, ' times in timestep ', myIter+1
394          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
395       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
396            _END_MASTER( myThid )
397         ENDIF         ENDIF
398         WRITE(msgBuf,'(A,I6)')         _BEGIN_MASTER( myThid )
399           WRITE(msgBuf,'(A,I6,A,I10)')
400       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
401       &      totalKrylovIter       &      totalKrylovItersLoc, ' in timestep ', myIter+1
402          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
403       &       SQUEEZE_RIGHT, myThid )       &      SQUEEZE_RIGHT, myThid )
404                 _END_MASTER( myThid )
405        ENDIF        ENDIF
406    
407  #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.3  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22