/[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.15 by mlosch, Wed Jan 16 21:20:28 2013 UTC revision 1.19 by mlosch, Mon Feb 25 10:44:10 2013 UTC
# Line 60  C     !FUNCTIONS: Line 60  C     !FUNCTIONS:
60        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
61        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
62    
63    C     !LOCAL VARIABLES:
64    C     === Local variables ===
65  C     i,j,bi,bj :: loop indices  C     i,j,bi,bj :: loop indices
66        INTEGER i,j,bi,bj        INTEGER i,j,bi,bj
67  C     loop indices  C     loop indices
# Line 147  C     Start nonlinear Newton iteration: Line 149  C     Start nonlinear Newton iteration:
149         newtonIter = newtonIter + 1         newtonIter = newtonIter + 1
150  C     Compute initial residual F(u), (includes computation of global  C     Compute initial residual F(u), (includes computation of global
151  C     variables DWATN, zeta, and eta)  C     variables DWATN, zeta, and eta)
152  C     Update linear solution vector and return to Newton iteration         IF ( newtonIter .EQ. 1 ) CALL SEAICE_JFNK_UPDATE(
 C     Do the linesearch  
        CALL SEAICE_JFNK_UPDATE(  
153       I      duIce, dvIce,       I      duIce, dvIce,
154       U      uIce, vIce, JFNKresidual,       U      uIce, vIce, JFNKresidual,
155       O      uIceRes, vIceRes,       O      uIceRes, vIceRes,
156       I      newtonIter, myTime, myIter, myThid )       I      newtonIter, myTime, myIter, myThid )
 C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver  
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO J=1-Oly,sNy+Oly  
           DO I=1-Olx,sNx+Olx  
            duIce(I,J,bi,bj)= 0. _d 0  
            dvIce(I,J,bi,bj)= 0. _d 0  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
 CMLC     Do it again, Sam  
 CML       CALL SEAICE_CALC_RESIDUAL(  
 CML     I      uIce, vIce,  
 CML     O      uIceRes, vIceRes,  
 CML     I      newtonIter, 0, myTime, myIter, myThid )  
 CMLC     probably not necessary, will be removed later:  
 CML       CALL EXCH_UV_XY_RL( uIceRes, vIceRes,.TRUE.,myThid)  
 CMLC     Important: Compute the norm of the residual using the same scalar  
 CMLC     product that SEAICE_FGMRES does  
 CML       CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid)  
 CML       CALL SEAICE_SCALPROD(  
 CML     &      nVec,1,1,1,resTmp,resTmp,JFNKresidual,myThid)  
 CML       JFNKresidual = SQRT(JFNKresidual)  
157  C     local copies of precomputed coefficients that are to stay  C     local copies of precomputed coefficients that are to stay
158  C     constant for the preconditioner  C     constant for the preconditioner
159         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
# Line 194  C     constant for the preconditioner Line 170  C     constant for the preconditioner
170         ENDDO         ENDDO
171  C     compute convergence criterion for linear preconditioned FGMRES  C     compute convergence criterion for linear preconditioned FGMRES
172         JFNKgamma_lin = JFNKgamma_lin_max         JFNKgamma_lin = JFNKgamma_lin_max
173         IF ( newtonIter.GT.1.AND.newtonIter.LE.100         IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter
174       &      .AND.JFNKresidual.LT.JFNKres_t ) THEN       &      .AND.JFNKresidual.LT.JFNKres_t ) THEN
175  C     Eisenstat, 1996, equ.(2.6)        C     Eisenstat, 1996, equ.(2.6)      
176          phi_e = 1. _d 0          phi_e = 1. _d 0
# Line 278  C     some output diagnostics Line 254  C     some output diagnostics
254          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN          IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
255           krylovFails = krylovFails + 1           krylovFails = krylovFails + 1
256          ENDIF          ENDIF
257  C     Set the stopping criterion for the Newton iteration  C     Set the stopping criterion for the Newton iteration and the
258          IF ( newtonIter .EQ. 1 ) JFNKtol=JFNKgamma_nonlin*JFNKresidual  C     criterion for the transition from accurate to approximate FGMRES
259            IF ( newtonIter .EQ. 1 ) THEN
260             JFNKtol=JFNKgamma_nonlin*JFNKresidual
261             IF ( JFNKres_tFac .NE. UNSET_RL )
262         &        JFNKres_t = JFNKresidual * JFNKres_tFac
263            ENDIF
264  C     Update linear solution vector and return to Newton iteration  C     Update linear solution vector and return to Newton iteration
265  C     Do a linesearch if necessary, and compute a new residual.  C     Do a linesearch if necessary, and compute a new residual.
266  C     Note that it should be possible to do the following operations  C     Note that it should be possible to do the following operations
# Line 419  C     Print more debugging information Line 400  C     Print more debugging information
400        RETURN        RETURN
401        END        END
402    
403    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
404  CBOP  CBOP
405  C     !ROUTINE: SEAICE_JFNK_UPDATE  C     !ROUTINE: SEAICE_JFNK_UPDATE
406  C     !INTERFACE:  C     !INTERFACE:
# Line 475  C     u/vIceRes :: residual of sea-ice m Line 457  C     u/vIceRes :: residual of sea-ice m
457        _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
458        _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
459    
460  C     Local variables:  C     !LOCAL VARIABLES:
461    C     === Local variables ===
462  C     i,j,bi,bj :: loop indices  C     i,j,bi,bj :: loop indices
463        INTEGER i,j,bi,bj        INTEGER i,j,bi,bj
464        INTEGER l        INTEGER l
# Line 496  C     Initialise some local variables Line 479  C     Initialise some local variables
479        facLS = 1. _d 0        facLS = 1. _d 0
480        doLineSearch = .TRUE.        doLineSearch = .TRUE.
481        DO WHILE ( doLineSearch )        DO WHILE ( doLineSearch )
 C     Determine, if we need more iterations  
        doLineSearch = resLoc .GE. JFNKresidual  
        doLineSearch = doLineSearch .AND. l .LE. 4  
 C     For the first iteration du/vIce = 0 and there will be no  
 C     improvement of the residual possible, so we do only the first  
 C     iteration  
        IF ( newtonIter .EQ. 1 ) doLineSearch = .FALSE.  
 C     Only start a linesearch after some Newton iterations  
        IF ( newtonIter .LE. SEAICE_JFNK_lsIter ) doLineSearch = .FALSE.  
 C     Increment counter  
        l = l + 1  
482  C     Create update  C     Create update
483         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
484          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 529  C     product that SEAICE_FGMRES does Line 501  C     product that SEAICE_FGMRES does
501         CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid)         CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid)
502         CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,resLoc,myThid)         CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,resLoc,myThid)
503         resLoc = SQRT(resLoc)         resLoc = SQRT(resLoc)
504    C     Determine, if we need more iterations
505           doLineSearch = resLoc .GE. JFNKresidual
506    C     Limit the maximum number of iterations arbitrarily to four
507           doLineSearch = doLineSearch .AND. l .LT. 4
508    C     For the first iteration du/vIce = 0 and there will be no
509    C     improvement of the residual possible, so we do only the first
510    C     iteration
511           IF ( newtonIter .EQ. 1 ) doLineSearch = .FALSE.
512    C     Only start a linesearch after some Newton iterations
513           IF ( newtonIter .LE. SEAICE_JFNK_lsIter ) doLineSearch = .FALSE.
514    C     Increment counter
515           l = l + 1
516  C     some output diagnostics  C     some output diagnostics
517         IF ( debugLevel.GE.debLevA .AND. doLineSearch ) THEN         IF ( debugLevel.GE.debLevA .AND. doLineSearch ) THEN
518          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22