/[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.29 by mlosch, Wed Jan 27 14:03:34 2016 UTC revision 1.30 by mlosch, Thu Jan 28 12:54:12 2016 UTC
# Line 63  C     !FUNCTIONS: Line 63  C     !FUNCTIONS:
63    
64  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
65  C     === Local variables ===  C     === Local variables ===
66  C     i,j,k,bi,bj :: loop indices  C     i,j,bi,bj :: loop indices
67        INTEGER i,j,k,bi,bj        INTEGER i,j,bi,bj
68  C     loop indices  C     loop indices
69        INTEGER newtonIter        INTEGER newtonIter
70        INTEGER krylovIter, krylovFails        INTEGER krylovIter, krylovFails
# Line 181  CML        ENDIF Line 181  CML        ENDIF
181         ENDDO         ENDDO
182        ENDDO        ENDDO
183  C     Start nonlinear Newton iteration: outer loop iteration  C     Start nonlinear Newton iteration: outer loop iteration
184        DO WHILE ( newtonIter.LT.SEAICEnewtonIterMax .AND.        DO WHILE ( newtonIter.LT.SEAICEnonLinIterMax .AND.
185       &     .NOT.JFNKconverged )       &     .NOT.JFNKconverged )
186         newtonIter = newtonIter + 1         newtonIter = newtonIter + 1
187  C     Compute initial residual F(u), (includes computation of global  C     Compute initial residual F(u), (includes computation of global
# Line 240  C     start Krylov iteration (FGMRES) Line 240  C     start Krylov iteration (FGMRES)
240  C     map first guess sol; it is zero because the solution is a correction  C     map first guess sol; it is zero because the solution is a correction
241         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid)         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid)
242  C     map rhs and change its sign because we are solving J*u = -F  C     map rhs and change its sign because we are solving J*u = -F
243          CALL SEAICE_MAP2VEC(nVec,-uIceRes,-vIceRes,rhs,.TRUE.,myThid)          CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,rhs,.TRUE.,myThid)
244            DO bj=myByLo(myThid),myByHi(myThid)
245             DO bi=myBxLo(myThid),myBxHi(myThid)
246              DO j=1,nVec
247               rhs(j,bi,bj) = - rhs(j,bi,bj)
248              ENDDO
249             ENDDO
250            ENDDO
251          DO WHILE ( .NOT.krylovConverged )          DO WHILE ( .NOT.krylovConverged )
252  C     solution vector sol = du/vIce  C     solution vector sol = du/vIce
253  C     residual vector (rhs) Fu = u/vIceRes  C     residual vector (rhs) Fu = u/vIceRes
# Line 253  C     because du/vIce = 0 Line 260  C     because du/vIce = 0
260  C  C
261           CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter,           CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter,
262       U        vv,w,wk1,wk2,       U        vv,w,wk1,wk2,
263       I        FGMRESeps,SEAICEkrylovIterMax,iOutFGMRES,       I        FGMRESeps,SEAICElinearIterMax,iOutFGMRES,
264       U        iCode,       U        iCode,
265       I        myThid)       I        myThid)
266  C      C    
# Line 272  C     or product of matrix (Jacobian) ti Line 279  C     or product of matrix (Jacobian) ti
279  C     iteration  C     iteration
280           IF (iCode.EQ.1) THEN           IF (iCode.EQ.1) THEN
281  C     Call preconditioner  C     Call preconditioner
282            IF ( SOLV_MAX_ITERS .GT. 0 )            IF ( SEAICEpreconLinIter .GT. 0 )
283       &         CALL SEAICE_PRECONDITIONER(       &         CALL SEAICE_PRECONDITIONER(
284       U         duIce, dvIce,       U         duIce, dvIce,
285       I         zetaPre, etaPre, etaZpre, zetaZpre, dwatPre,       I         zetaPre, etaPre, etaZpre, zetaZpre, dwatPre,
# Line 292  C     some output diagnostics Line 299  C     some output diagnostics
299          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
300           _BEGIN_MASTER( myThid )           _BEGIN_MASTER( myThid )
301           totalNewtonItersLoc =           totalNewtonItersLoc =
302       &        SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter       &        SEAICEnonLinIterMax*(myIter-nIter0)+newtonIter
303           WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')           WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
304       &        ' S/R SEAICE_JFNK: Newton iterate / total, ',       &        ' S/R SEAICE_JFNK: Newton iterate / total, ',
305       &        'JFNKgamma_lin, initial norm = ',       &        'JFNKgamma_lin, initial norm = ',
# Line 308  C     some output diagnostics Line 315  C     some output diagnostics
315       &        SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
316           _END_MASTER( myThid )           _END_MASTER( myThid )
317          ENDIF          ENDIF
318          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICElinearIterMax ) THEN
319           krylovFails = krylovFails + 1           krylovFails = krylovFails + 1
320          ENDIF          ENDIF
321  C     Set the stopping criterion for the Newton iteration and the  C     Set the stopping criterion for the Newton iteration and the
322  C     criterion for the transition from accurate to approximate FGMRES  C     criterion for the transition from accurate to approximate FGMRES
323          IF ( newtonIter .EQ. 1 ) THEN          IF ( newtonIter .EQ. 1 ) THEN
324           JFNKtol=JFNKgamma_nonlin*JFNKresidual           JFNKtol=SEAICEnonLinTol*JFNKresidual
325           IF ( JFNKres_tFac .NE. UNSET_RL )           IF ( JFNKres_tFac .NE. UNSET_RL )
326       &        JFNKres_t = JFNKresidual * JFNKres_tFac       &        JFNKres_t = JFNKresidual * JFNKres_tFac
327          ENDIF          ENDIF
# Line 353  C     Count iterations Line 360  C     Count iterations
360         totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc         totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
361  C     Record failure  C     Record failure
362         totalKrylovFails   = totalKrylovFails + krylovFails         totalKrylovFails   = totalKrylovFails + krylovFails
363         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnonLinIterMax ) THEN
364          totalNewtonFails = totalNewtonFails + 1          totalNewtonFails = totalNewtonFails + 1
365         ENDIF         ENDIF
366        ENDIF        ENDIF
# Line 427  C     reset and start again Line 434  C     reset and start again
434    
435  C     Print more debugging information  C     Print more debugging information
436        IF ( debugLevel.GE.debLevA ) THEN        IF ( debugLevel.GE.debLevA ) THEN
437         IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN         IF ( newtonIter .EQ. SEAICEnonLinIterMax ) THEN
438          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
439          WRITE(msgBuf,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
440       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',       &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22