/[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.8 by mlosch, Fri Nov 9 14:11:43 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 81  C     precomputed (= constant per Newton Line 86  C     precomputed (= constant per Newton
86  C     zeta, eta, and DWATN, press  C     zeta, eta, and DWATN, press
87        _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
88        _RL etaPre  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL etaPre  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
89          _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
90        _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
91        _RL pressPre(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     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 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           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 140  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             pressPre(I,J,bi,bj) = press(I,J,bi,bj)             pressPre(I,J,bi,bj) = press(I,J,bi,bj)
159            ENDDO            ENDDO
# Line 193  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 = ',
214       &       newtonIter, SEAICEnewtonIterMax*myIter+newtonIter,       &       newtonIter,SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter,
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 217  C     solution vector sol = du/vIce Line 233  C     solution vector sol = du/vIce
233  C     residual vector (rhs) Fu = u/vIceRes  C     residual vector (rhs) Fu = u/vIceRes
234  C     output work vectors wk1, -> input work vector wk2  C     output work vectors wk1, -> input work vector wk2
235  C      C    
          CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)  
          CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)  
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, pressPre,
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 243  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+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 263  C     Update linear solution vector and Line 280  C     Update linear solution vector and
280             DO I=1-Olx,sNx+Olx             DO I=1-Olx,sNx+Olx
281              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)
282              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)
283    C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
284                duIce(I,J,bi,bj)= 0. _d 0
285                dvIce(I,J,bi,bj)= 0. _d 0
286             ENDDO             ENDDO
287            ENDDO            ENDDO
288           ENDDO           ENDDO
# Line 272  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          IF ( DIFFERENT_MULTIPLE(SEAICE_monFreq,myTime+deltaTClock,
311         &     deltaTClock) ) THEN
312           _BEGIN_MASTER( myThid )
313           WRITE(msgBuf,'(A)')
314         &' // ======================================================='
315           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
316         &      SQUEEZE_RIGHT, myThid )
317           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
318           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
319         &      SQUEEZE_RIGHT, myThid )
320           WRITE(msgBuf,'(A)')
321         &' // ======================================================='
322           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
323         &      SQUEEZE_RIGHT, myThid )
324           WRITE(msgBuf,'(A,I10)')
325         &      ' %JFNK_MON: time step              = ', myIter+1
326           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
327         &      SQUEEZE_RIGHT, myThid )
328           WRITE(msgBuf,'(A,I10)')
329         &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
330           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
331         &      SQUEEZE_RIGHT, myThid )
332           WRITE(msgBuf,'(A,I10)')
333         &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
334           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
335         &      SQUEEZE_RIGHT, myThid )
336           WRITE(msgBuf,'(A,I10)')
337         &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
338           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
339         &      SQUEEZE_RIGHT, myThid )
340           WRITE(msgBuf,'(A,I10)')
341         &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
342           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
343         &      SQUEEZE_RIGHT, myThid )
344           WRITE(msgBuf,'(A,I10)')
345         &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
346           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
347         &      SQUEEZE_RIGHT, myThid )
348           WRITE(msgBuf,'(A)')
349         &' // ======================================================='
350           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
351         &      SQUEEZE_RIGHT, myThid )
352           WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
353           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
354         &      SQUEEZE_RIGHT, myThid )
355           WRITE(msgBuf,'(A)')
356         &' // ======================================================='
357           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
358         &      SQUEEZE_RIGHT, myThid )
359           _END_MASTER( myThid )
360    C     reset and start again
361           totalJFNKtimeSteps = 0
362           totalNewtonIters   = 0
363           totalKrylovIters   = 0
364           totalKrylovFails   = 0
365           totalNewtonFails   = 0
366          ENDIF
367    
368    C     Print more debugging information
369          IF ( debugLevel.GE.debLevA ) THEN
370         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
371          newtonIterFail = newtonIterFail + 1          _BEGIN_MASTER( myThid )
372          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
373       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
374       &       myIter       &       myIter+1
375          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
376       &       SQUEEZE_RIGHT, myThid )       &       SQUEEZE_RIGHT, myThid )
377            _END_MASTER( myThid )
378         ENDIF         ENDIF
379         IF ( krylovIterFail .GT. 0 ) THEN         IF ( krylovFails .GT. 0 ) THEN
380            _BEGIN_MASTER( myThid )
381          WRITE(msgBuf,'(A,I4,A,I10)')          WRITE(msgBuf,'(A,I4,A,I10)')
382       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',       &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
383       &       krylovIterFail, ' times in timestep ', myIter       &       krylovFails, ' times in timestep ', 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         WRITE(msgBuf,'(A,I6)')         _BEGIN_MASTER( myThid )
389           WRITE(msgBuf,'(A,I6,A,I10)')
390       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',       &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
391       &      totalKrylovIter       &      totalKrylovItersLoc, ' 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    
397  #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.8
  ViewVC Help
Powered by ViewVC 1.1.22