/[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.12 by mlosch, Mon Dec 17 10:08:16 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)
       _RL pressPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  
92  CEOP  CEOP
93    
94  C     Initialise  C     Initialise
95        newtonIter      = 0        newtonIter          = 0
96        newtonIterFail  = 0        krylovFails         = 0
97        krylovIterFail  = 0        totalKrylovItersLoc = 0
98        totalKrylovIter = 0        JFNKconverged       = .FALSE.
99        JFNKconverged   = .FALSE.        JFNKtol             = 0. _d 0
100        JFNKtol         = 0. _d 0        JFNKresidual        = 0. _d 0
101        JFNKresidual    = 0. _d 0        JFNKresidualKm1     = 0. _d 0
102        JFNKresidualKm1 = 0. _d 0        FGMRESeps           = 0. _d 0
103        FGMRESeps       = 0. _d 0        recip_deltaT        = 1. _d 0 / SEAICE_deltaTdyn
104        recip_deltaT    = 1. _d 0 / SEAICE_deltaTdyn  
105          iOutFGMRES=0
106    C     with iOutFgmres=1, seaice_fgmres prints the residual at each iteration
107          IF ( debugLevel.GE.debLevC .AND.
108         &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
109         &     iOutFGMRES=1
110    
111  C      C    
112        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
113         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 132  C     variables DWATN, zeta, and eta) Line 144  C     variables DWATN, zeta, and eta)
144       I      uIce, vIce,       I      uIce, vIce,
145       O      uIceRes, vIceRes,       O      uIceRes, vIceRes,
146       I      newtonIter, 0, myTime, myIter, myThid )       I      newtonIter, 0, myTime, myIter, myThid )
147    C     probably not necessary, will be removed later:
148         CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)         CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)
149  C     local copies of precomputed coefficients that are to stay  C     local copies of precomputed coefficients that are to stay
150  C     constant for the preconditioner  C     constant for the preconditioner
# Line 139  C     constant for the preconditioner Line 152  C     constant for the preconditioner
152          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
153           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
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              dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)             etaZPre(I,J,bi,bj) =  etaZ(I,J,bi,bj)
158             pressPre(I,J,bi,bj) = press(I,J,bi,bj)             dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)
159            ENDDO            ENDDO
160           ENDDO           ENDDO
161          ENDDO          ENDDO
# Line 194  C     in that routine Line 207  C     in that routine
207         krylovIter    = 0         krylovIter    = 0
208         iCode         = 0         iCode         = 0
209         IF ( debugLevel.GE.debLevA ) THEN           IF ( debugLevel.GE.debLevA ) THEN  
210            _BEGIN_MASTER( myThid )
211          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
212       &       ' S/R SEAICE_JFNK: newtonIter,',       &       ' S/R SEAICE_JFNK: newtonIter,',
213       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',
# Line 201  C     in that routine Line 215  C     in that routine
215       &       JFNKgamma_lin, JFNKresidual       &       JFNKgamma_lin, JFNKresidual
216          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
217       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
218            _END_MASTER( myThid )
219         ENDIF         ENDIF
220  C  C
221         JFNKconverged = JFNKresidual.LT.JFNKtol         JFNKconverged = JFNKresidual.LT.JFNKtol
# Line 221  C Line 236  C
236           CALL SEAICE_FGMRES_DRIVER(           CALL SEAICE_FGMRES_DRIVER(
237       I        uIceRes, vIceRes,       I        uIceRes, vIceRes,
238       U        duIce, dvIce, iCode,       U        duIce, dvIce, iCode,
239       I        FGMRESeps,         I        FGMRESeps, iOutFGMRES,
240       I        newtonIter, krylovIter, myTime, myIter, myThid )       I        newtonIter, krylovIter, myTime, myIter, myThid )
241  C     FGMRES returns iCode either asking for an new preconditioned vector  C     FGMRES returns iCode either asking for an new preconditioned vector
242  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate
243  C     iteration  C     iteration
244           IF (iCode.EQ.1) THEN           IF (iCode.EQ.1) THEN
245  C     Call preconditioner  C     Call preconditioner
246            CALL SEAICE_PRECONDITIONER(            IF ( SOLV_MAX_ITERS .GT. 0 )
247         &         CALL SEAICE_PRECONDITIONER(
248       U         duIce, dvIce,       U         duIce, dvIce,
249       I         zetaPre, etaPre, dwatPre, pressPre,       I         zetaPre, etaPre, etaZpre, dwatPre,
250       I         newtonIter, krylovIter, myTime, myIter, myThid )       I         newtonIter, krylovIter, myTime, myIter, myThid )
251           ELSEIF (iCode.GE.2) THEN           ELSEIF (iCode.GE.2) THEN
252  C     Compute Jacobian times vector  C     Compute Jacobian times vector
# Line 242  C     Compute Jacobian times vector Line 258  C     Compute Jacobian times vector
258           krylovConverged = iCode.EQ.0           krylovConverged = iCode.EQ.0
259  C     End of Krylov iterate  C     End of Krylov iterate
260          ENDDO          ENDDO
261          totalKrylovIter = totalKrylovIter + krylovIter          totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
262  C     some output diagnostics  C     some output diagnostics
263          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
264             _BEGIN_MASTER( myThid )
265           WRITE(msgBuf,'(3(A,I6))')           WRITE(msgBuf,'(3(A,I6))')
266       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,
267       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
268       &        ', Nb. of FGMRES iterations = ', krylovIter       &        ', Nb. of FGMRES iterations = ', krylovIter
269           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
270       &        SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
271             _END_MASTER( myThid )
272          ENDIF          ENDIF
273          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
274           krylovIterFail = krylovIterFail + 1           krylovFails = krylovFails + 1
275          ENDIF          ENDIF
276  C     Update linear solution vector and return to Newton iteration  C     Update linear solution vector and return to Newton iteration
277          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 274  C     Set the stopping criterion for the Line 292  C     Set the stopping criterion for the
292         ENDIF         ENDIF
293  C     end of Newton iterate  C     end of Newton iterate
294        ENDDO        ENDDO
295  C     some output diagnostics  C
296        IF ( debugLevel.GE.debLevA ) THEN  C--   Output diagnostics
297    C
298          IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
299    C     Count iterations
300           totalJFNKtimeSteps = totalJFNKtimeSteps + 1
301           totalNewtonIters   = totalNewtonIters + newtonIter
302           totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
303  C     Record failure  C     Record failure
304           totalKrylovFails   = totalKrylovFails + krylovFails
305           IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
306            totalNewtonFails = totalNewtonFails + 1
307           ENDIF
308          ENDIF
309    C     Decide whether it is time to dump and reset the counter
310          writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
311         &     myTime+deltaTClock, deltaTClock)
312    #ifdef ALLOW_CAL
313          IF ( useCAL ) THEN
314           CALL CAL_TIME2DUMP(
315         I      zeroRL, SEAICE_monFreq,  deltaTClock,
316         U      writeNow,
317         I      myTime+deltaTclock, myIter+1, myThid )
318          ENDIF
319    #endif
320          IF ( writeNow ) THEN
321           _BEGIN_MASTER( myThid )
322           WRITE(msgBuf,'(A)')
323         &' // ======================================================='
324           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
325         &      SQUEEZE_RIGHT, myThid )
326           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
327           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
328         &      SQUEEZE_RIGHT, myThid )
329           WRITE(msgBuf,'(A)')
330         &' // ======================================================='
331           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
332         &      SQUEEZE_RIGHT, myThid )
333           WRITE(msgBuf,'(A,I10)')
334         &      ' %JFNK_MON: time step              = ', myIter+1
335           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
336         &      SQUEEZE_RIGHT, myThid )
337           WRITE(msgBuf,'(A,I10)')
338         &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
339           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
340         &      SQUEEZE_RIGHT, myThid )
341           WRITE(msgBuf,'(A,I10)')
342         &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
343           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
344         &      SQUEEZE_RIGHT, myThid )
345           WRITE(msgBuf,'(A,I10)')
346         &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
347           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
348         &      SQUEEZE_RIGHT, myThid )
349           WRITE(msgBuf,'(A,I10)')
350         &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
351           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
352         &      SQUEEZE_RIGHT, myThid )
353           WRITE(msgBuf,'(A,I10)')
354         &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
355           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
356         &      SQUEEZE_RIGHT, myThid )
357           WRITE(msgBuf,'(A)')
358         &' // ======================================================='
359           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
360         &      SQUEEZE_RIGHT, myThid )
361           WRITE(msgBuf,'(A)') ' // End JFNK statistics'
362           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
363         &      SQUEEZE_RIGHT, myThid )
364           WRITE(msgBuf,'(A)')
365         &' // ======================================================='
366           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
367         &      SQUEEZE_RIGHT, myThid )
368           _END_MASTER( myThid )
369    C     reset and start again
370           totalJFNKtimeSteps = 0
371           totalNewtonIters   = 0
372           totalKrylovIters   = 0
373           totalKrylovFails   = 0
374           totalNewtonFails   = 0
375          ENDIF
376    
377    C     Print more debugging information
378          IF ( debugLevel.GE.debLevA ) THEN
379         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
380          newtonIterFail = newtonIterFail + 1          _BEGIN_MASTER( myThid )
381          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
382       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
383       &       myIter       &       myIter+1
384          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
385       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
386            _END_MASTER( myThid )
387         ENDIF         ENDIF
388         IF ( krylovIterFail .GT. 0 ) THEN         IF ( krylovFails .GT. 0 ) THEN
389            _BEGIN_MASTER( myThid )
390          WRITE(msgBuf,'(A,I4,A,I10)')          WRITE(msgBuf,'(A,I4,A,I10)')
391       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
392       &       krylovIterFail, ' times in timestep ', myIter       &       krylovFails, ' times in timestep ', myIter+1
393          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
394       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
395            _END_MASTER( myThid )
396         ENDIF         ENDIF
397         WRITE(msgBuf,'(A,I6)')         _BEGIN_MASTER( myThid )
398           WRITE(msgBuf,'(A,I6,A,I10)')
399       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
400       &      totalKrylovIter       &      totalKrylovItersLoc, ' in timestep ', myIter+1
401          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
402       &       SQUEEZE_RIGHT, myThid )       &      SQUEEZE_RIGHT, myThid )
403                 _END_MASTER( myThid )
404        ENDIF        ENDIF
405    
406  #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.12
  ViewVC Help
Powered by ViewVC 1.1.22