/[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.1 by mlosch, Tue Oct 16 07:00:21 2012 UTC revision 1.11 by mlosch, Mon Dec 3 15:49:17 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 78  C     du/vIce   :: ice velocity incremen Line 84  C     du/vIce   :: ice velocity incremen
84        _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85        _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
86  C     precomputed (= constant per Newton iteration) versions of  C     precomputed (= constant per Newton iteration) versions of
87  C     zeta, eta, and DWATN  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 dwatPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _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)
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     iOutFgmres=1 gives a little bit of output
107          IF ( debugLevel.GE.debLevA .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 131  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           CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)
148  C     local copies of precomputed coefficients that are to stay  C     local copies of precomputed coefficients that are to stay
149  C     constant for the preconditioner  C     constant for the preconditioner
150         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
# Line 139  C     constant for the preconditioner Line 153  C     constant for the preconditioner
153            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
154             zetaPre(I,J,bi,bj) =  zeta(I,J,bi,bj)             zetaPre(I,J,bi,bj) =  zeta(I,J,bi,bj)
155              etaPre(I,J,bi,bj) =   eta(I,J,bi,bj)              etaPre(I,J,bi,bj) =   eta(I,J,bi,bj)
156               etaZPre(I,J,bi,bj) =  etaZ(I,J,bi,bj)
157             dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)             dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)
158            ENDDO            ENDDO
159           ENDDO           ENDDO
# Line 191  C     in that routine Line 206  C     in that routine
206         krylovIter    = 0         krylovIter    = 0
207         iCode         = 0         iCode         = 0
208         IF ( debugLevel.GE.debLevA ) THEN           IF ( debugLevel.GE.debLevA ) THEN  
209            _BEGIN_MASTER( myThid )
210          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
211       &       ' S/R SEAICE_JFNK: newtonIter,',       &       ' S/R SEAICE_JFNK: newtonIter,',
212       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',       &       ' total newtonIter, JFNKgamma_lin, initial norm = ',
213       &       newtonIter, SEAICEnewtonIterMax*myIter+newtonIter,       &       newtonIter,SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
214       &       JFNKgamma_lin, JFNKresidual       &       JFNKgamma_lin, JFNKresidual
215          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
216       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
217            _END_MASTER( myThid )
218         ENDIF         ENDIF
219  C  C
220         JFNKconverged = JFNKresidual.LT.JFNKtol         JFNKconverged = JFNKresidual.LT.JFNKtol
# Line 215  C     solution vector sol = du/vIce Line 232  C     solution vector sol = du/vIce
232  C     residual vector (rhs) Fu = u/vIceRes  C     residual vector (rhs) Fu = u/vIceRes
233  C     output work vectors wk1, -> input work vector wk2  C     output work vectors wk1, -> input work vector wk2
234  C      C    
          CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)  
          CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)  
235           CALL SEAICE_FGMRES_DRIVER(           CALL SEAICE_FGMRES_DRIVER(
236       I        uIceRes, vIceRes,       I        uIceRes, vIceRes,
237       U        duIce, dvIce, iCode,       U        duIce, dvIce, iCode,
238       I        FGMRESeps,         I        FGMRESeps, iOutFGMRES,
239       I        newtonIter, krylovIter, myTime, myIter, myThid )       I        newtonIter, krylovIter, myTime, myIter, myThid )
240  C     FGMRES returns iCode either asking for an new preconditioned vector  C     FGMRES returns iCode either asking for an new preconditioned vector
241  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate  C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate
242  C     iteration  C     iteration
243           IF (iCode.EQ.1) THEN           IF (iCode.EQ.1) THEN
244  C     Call preconditioner  C     Call preconditioner
245            CALL SEAICE_PRECONDITIONER(            IF ( SOLV_MAX_ITERS .GT. 0 )
246         &         CALL SEAICE_PRECONDITIONER(
247       U         duIce, dvIce,       U         duIce, dvIce,
248       I         zetaPre, etaPre, dwatPre,       I         zetaPre, etaPre, etaZpre, dwatPre,
249       I         newtonIter, krylovIter, myTime, myIter, myThid )       I         newtonIter, krylovIter, myTime, myIter, myThid )
250           ELSEIF (iCode.GE.2) THEN           ELSEIF (iCode.GE.2) THEN
251  C     Compute Jacobian times vector  C     Compute Jacobian times vector
# Line 241  C     Compute Jacobian times vector Line 257  C     Compute Jacobian times vector
257           krylovConverged = iCode.EQ.0           krylovConverged = iCode.EQ.0
258  C     End of Krylov iterate  C     End of Krylov iterate
259          ENDDO          ENDDO
260          totalKrylovIter = totalKrylovIter + krylovIter          totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
261  C     some output diagnostics  C     some output diagnostics
262          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
263             _BEGIN_MASTER( myThid )
264           WRITE(msgBuf,'(3(A,I6))')           WRITE(msgBuf,'(3(A,I6))')
265       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,       &        ' S/R SEAICE_JFNK: Newton iterate / total = ', newtonIter,
266       &        ' / ', SEAICEnewtonIterMax*myIter+newtonIter,       &        ' / ', SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
267       &        ', Nb. of FGMRES iterations = ', krylovIter       &        ', Nb. of FGMRES iterations = ', krylovIter
268           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
269       &        SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
270             _END_MASTER( myThid )
271          ENDIF          ENDIF
272          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
273           krylovIterFail = krylovIterFail + 1           krylovFails = krylovFails + 1
274          ENDIF          ENDIF
275  C     Update linear solution vector and return to Newton iteration  C     Update linear solution vector and return to Newton iteration
276          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 261  C     Update linear solution vector and Line 279  C     Update linear solution vector and
279             DO I=1-Olx,sNx+Olx             DO I=1-Olx,sNx+Olx
280              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)
281              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)
282    C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
283                duIce(I,J,bi,bj)= 0. _d 0
284                dvIce(I,J,bi,bj)= 0. _d 0
285             ENDDO             ENDDO
286            ENDDO            ENDDO
287           ENDDO           ENDDO
# Line 270  C     Set the stopping criterion for the Line 291  C     Set the stopping criterion for the
291         ENDIF         ENDIF
292  C     end of Newton iterate  C     end of Newton iterate
293        ENDDO        ENDDO
294  C     some output diagnostics  C
295        IF ( debugLevel.GE.debLevA ) THEN  C--   Output diagnostics
296    C
297          IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
298    C     Count iterations
299           totalJFNKtimeSteps = totalJFNKtimeSteps + 1
300           totalNewtonIters   = totalNewtonIters + newtonIter
301           totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
302  C     Record failure  C     Record failure
303           totalKrylovFails   = totalKrylovFails + krylovFails
304           IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
305            totalNewtonFails = totalNewtonFails + 1
306           ENDIF
307          ENDIF
308    C     Decide whether it is time to dump and reset the counter
309          writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
310         &     myTime+deltaTClock, deltaTClock)
311    #ifdef ALLOW_CAL
312          IF ( useCAL ) THEN
313           CALL CAL_TIME2DUMP(
314         I      zeroRL, SEAICE_monFreq,  deltaTClock,
315         U      writeNow,
316         I      myTime+deltaTclock, myIter+1, myThid )
317          ENDIF
318    #endif
319          IF ( writeNow ) THEN
320           _BEGIN_MASTER( myThid )
321           WRITE(msgBuf,'(A)')
322         &' // ======================================================='
323           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
324         &      SQUEEZE_RIGHT, myThid )
325           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
326           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
327         &      SQUEEZE_RIGHT, myThid )
328           WRITE(msgBuf,'(A)')
329         &' // ======================================================='
330           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
331         &      SQUEEZE_RIGHT, myThid )
332           WRITE(msgBuf,'(A,I10)')
333         &      ' %JFNK_MON: time step              = ', myIter+1
334           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
335         &      SQUEEZE_RIGHT, myThid )
336           WRITE(msgBuf,'(A,I10)')
337         &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
338           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
339         &      SQUEEZE_RIGHT, myThid )
340           WRITE(msgBuf,'(A,I10)')
341         &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
342           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
343         &      SQUEEZE_RIGHT, myThid )
344           WRITE(msgBuf,'(A,I10)')
345         &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
346           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
347         &      SQUEEZE_RIGHT, myThid )
348           WRITE(msgBuf,'(A,I10)')
349         &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
350           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
351         &      SQUEEZE_RIGHT, myThid )
352           WRITE(msgBuf,'(A,I10)')
353         &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
354           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
355         &      SQUEEZE_RIGHT, myThid )
356           WRITE(msgBuf,'(A)')
357         &' // ======================================================='
358           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
359         &      SQUEEZE_RIGHT, myThid )
360           WRITE(msgBuf,'(A)') ' // End JFNK statistics'
361           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
362         &      SQUEEZE_RIGHT, myThid )
363           WRITE(msgBuf,'(A)')
364         &' // ======================================================='
365           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
366         &      SQUEEZE_RIGHT, myThid )
367           _END_MASTER( myThid )
368    C     reset and start again
369           totalJFNKtimeSteps = 0
370           totalNewtonIters   = 0
371           totalKrylovIters   = 0
372           totalKrylovFails   = 0
373           totalNewtonFails   = 0
374          ENDIF
375    
376    C     Print more debugging information
377          IF ( debugLevel.GE.debLevA ) THEN
378         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
379          newtonIterFail = newtonIterFail + 1          _BEGIN_MASTER( myThid )
380          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
381       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
382       &       myIter       &       myIter+1
383          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
384       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
385            _END_MASTER( myThid )
386         ENDIF         ENDIF
387         IF ( krylovIterFail .GT. 0 ) THEN         IF ( krylovFails .GT. 0 ) THEN
388            _BEGIN_MASTER( myThid )
389          WRITE(msgBuf,'(A,I4,A,I10)')          WRITE(msgBuf,'(A,I4,A,I10)')
390       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
391       &       krylovIterFail, ' times in timestep ', myIter       &       krylovFails, ' times in timestep ', myIter+1
392          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
393       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
394            _END_MASTER( myThid )
395         ENDIF         ENDIF
396         WRITE(msgBuf,'(A,I6)')         _BEGIN_MASTER( myThid )
397           WRITE(msgBuf,'(A,I6,A,I10)')
398       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
399       &      totalKrylovIter       &      totalKrylovItersLoc, ' in timestep ', myIter+1
400          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
401       &       SQUEEZE_RIGHT, myThid )       &      SQUEEZE_RIGHT, myThid )
402                 _END_MASTER( myThid )
403        ENDIF        ENDIF
404    
405  #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.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22