/[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.2 by mlosch, Wed Oct 17 14:53:51 2012 UTC revision 1.7 by mlosch, Fri Nov 9 12:56:00 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 132  C     variables DWATN, zeta, and eta) Line 143  C     variables DWATN, zeta, and eta)
143       I      uIce, vIce,       I      uIce, vIce,
144       O      uIceRes, vIceRes,       O      uIceRes, vIceRes,
145       I      newtonIter, 0, myTime, myIter, myThid )       I      newtonIter, 0, myTime, myIter, myThid )
146           CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)
147  C     local copies of precomputed coefficients that are to stay  C     local copies of precomputed coefficients that are to stay
148  C     constant for the preconditioner  C     constant for the preconditioner
149         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
# Line 193  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 = ',
212       &       newtonIter, SEAICEnewtonIterMax*myIter+newtonIter,       &       newtonIter,SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
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 217  C     solution vector sol = du/vIce Line 231  C     solution vector sol = du/vIce
231  C     residual vector (rhs) Fu = u/vIceRes  C     residual vector (rhs) Fu = u/vIceRes
232  C     output work vectors wk1, -> input work vector wk2  C     output work vectors wk1, -> input work vector wk2
233  C      C    
          CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)  
          CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)  
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
241  C     iteration  C     iteration
242           IF (iCode.EQ.1) THEN           IF (iCode.EQ.1) THEN
243  C     Call preconditioner  C     Call preconditioner
244            CALL SEAICE_PRECONDITIONER(            IF ( SOLV_MAX_ITERS .GT. 0 )
245         &         CALL SEAICE_PRECONDITIONER(
246       U         duIce, dvIce,       U         duIce, dvIce,
247       I         zetaPre, etaPre, dwatPre, pressPre,       I         zetaPre, etaPre, dwatPre, pressPre,
248       I         newtonIter, krylovIter, myTime, myIter, myThid )       I         newtonIter, krylovIter, myTime, myIter, myThid )
# Line 243  C     Compute Jacobian times vector Line 256  C     Compute Jacobian times vector
256           krylovConverged = iCode.EQ.0           krylovConverged = iCode.EQ.0
257  C     End of Krylov iterate  C     End of Krylov iterate
258          ENDDO          ENDDO
259          totalKrylovIter = totalKrylovIter + krylovIter          totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
260  C     some output diagnostics  C     some output diagnostics
261          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
262             _BEGIN_MASTER( myThid )
263           WRITE(msgBuf,'(3(A,I6))')           WRITE(msgBuf,'(3(A,I6))')
264       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,
265       &        ' / ', SEAICEnewtonIterMax*myIter+newtonIter,       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
266       &        ', Nb. of FGMRES iterations = ', krylovIter       &        ', Nb. of FGMRES iterations = ', krylovIter
267           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
268       &        SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
269             _END_MASTER( myThid )
270          ENDIF          ENDIF
271          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
272           krylovIterFail = krylovIterFail + 1           krylovFails = krylovFails + 1
273          ENDIF          ENDIF
274  C     Update linear solution vector and return to Newton iteration  C     Update linear solution vector and return to Newton iteration
275          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 263  C     Update linear solution vector and Line 278  C     Update linear solution vector and
278             DO I=1-Olx,sNx+Olx             DO I=1-Olx,sNx+Olx
279              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)
280              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)
281    C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
282                duIce(I,J,bi,bj)= 0. _d 0
283                dvIce(I,J,bi,bj)= 0. _d 0
284             ENDDO             ENDDO
285            ENDDO            ENDDO
286           ENDDO           ENDDO
# Line 272  C     Set the stopping criterion for the Line 290  C     Set the stopping criterion for the
290         ENDIF         ENDIF
291  C     end of Newton iterate  C     end of Newton iterate
292        ENDDO        ENDDO
293  C     some output diagnostics  C
294        IF ( debugLevel.GE.debLevA ) THEN  C--   Output diagnostics
295    C
296          IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
297    C     Count iterations
298           totalJFNKtimeSteps = totalJFNKtimeSteps + 1
299           totalNewtonIters   = totalNewtonIters + newtonIter
300           totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
301  C     Record failure  C     Record failure
302           totalKrylovFails   = totalKrylovFails + krylovFails
303           IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
304            totalNewtonFails = totalNewtonFails + 1
305           ENDIF
306          ENDIF
307    C     Decide whether it is time to dump and reset the counter
308          IF ( DIFFERENT_MULTIPLE(SEAICE_monFreq,myTime+deltaTClock,
309         &     deltaTClock) ) THEN
310           _BEGIN_MASTER( myThid )
311           WRITE(msgBuf,'(A)')
312         &' // ======================================================='
313           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
314         &      SQUEEZE_RIGHT, myThid )
315           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
316           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
317         &      SQUEEZE_RIGHT, myThid )
318           WRITE(msgBuf,'(A)')
319         &' // ======================================================='
320           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
321         &      SQUEEZE_RIGHT, myThid )
322           WRITE(msgBuf,'(A,I10)')
323         &      ' %JFNK_MON: time step              = ', myIter+1
324           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
325         &      SQUEEZE_RIGHT, myThid )
326           WRITE(msgBuf,'(A,I10)')
327         &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
328           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
329         &      SQUEEZE_RIGHT, myThid )
330           WRITE(msgBuf,'(A,I10)')
331         &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
332           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
333         &      SQUEEZE_RIGHT, myThid )
334           WRITE(msgBuf,'(A,I10)')
335         &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
336           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
337         &      SQUEEZE_RIGHT, myThid )
338           WRITE(msgBuf,'(A,I10)')
339         &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
340           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
341         &      SQUEEZE_RIGHT, myThid )
342           WRITE(msgBuf,'(A,I10)')
343         &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
344           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
345         &      SQUEEZE_RIGHT, myThid )
346           WRITE(msgBuf,'(A)')
347         &' // ======================================================='
348           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
349         &      SQUEEZE_RIGHT, myThid )
350           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
351           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
352         &      SQUEEZE_RIGHT, myThid )
353           WRITE(msgBuf,'(A)')
354         &' // ======================================================='
355           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
356         &      SQUEEZE_RIGHT, myThid )
357           _END_MASTER( myThid )
358    C     reset and start again
359           totalJFNKtimeSteps = 0
360           totalNewtonIters   = 0
361           totalKrylovIters   = 0
362           totalKrylovFails   = 0
363           totalNewtonFails   = 0
364          ENDIF
365    
366    C     Print more debugging information
367          IF ( debugLevel.GE.debLevA ) THEN
368         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
369          newtonIterFail = newtonIterFail + 1          _BEGIN_MASTER( myThid )
370          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
371       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
372       &       myIter       &       myIter+1
373          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
374       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
375            _END_MASTER( myThid )
376         ENDIF         ENDIF
377         IF ( krylovIterFail .GT. 0 ) THEN         IF ( krylovFails .GT. 0 ) THEN
378            _BEGIN_MASTER( myThid )
379          WRITE(msgBuf,'(A,I4,A,I10)')          WRITE(msgBuf,'(A,I4,A,I10)')
380       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
381       &       krylovIterFail, ' times in timestep ', myIter       &       krylovFails, ' times in timestep ', myIter+1
382          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
383       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
384            _END_MASTER( myThid )
385         ENDIF         ENDIF
386         WRITE(msgBuf,'(A,I6)')         _BEGIN_MASTER( myThid )
387           WRITE(msgBuf,'(A,I6,A,I10)')
388       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
389       &      totalKrylovIter       &      totalKrylovItersLoc, ' in timestep ', myIter+1
390          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
391       &       SQUEEZE_RIGHT, myThid )       &      SQUEEZE_RIGHT, myThid )
392                 _END_MASTER( myThid )
393        ENDIF        ENDIF
394    
395  #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.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22