--- MITgcm/pkg/seaice/seaice_jfnk.F 2013/01/16 21:20:28 1.15 +++ MITgcm/pkg/seaice/seaice_jfnk.F 2013/02/25 10:44:10 1.19 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.15 2013/01/16 21:20:28 mlosch Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.19 2013/02/25 10:44:10 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" @@ -60,6 +60,8 @@ LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE +C !LOCAL VARIABLES: +C === Local variables === C i,j,bi,bj :: loop indices INTEGER i,j,bi,bj C loop indices @@ -147,37 +149,11 @@ newtonIter = newtonIter + 1 C Compute initial residual F(u), (includes computation of global C variables DWATN, zeta, and eta) -C Update linear solution vector and return to Newton iteration -C Do the linesearch - CALL SEAICE_JFNK_UPDATE( + IF ( newtonIter .EQ. 1 ) CALL SEAICE_JFNK_UPDATE( I duIce, dvIce, U uIce, vIce, JFNKresidual, O uIceRes, vIceRes, 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) C local copies of precomputed coefficients that are to stay C constant for the preconditioner DO bj=myByLo(myThid),myByHi(myThid) @@ -194,7 +170,7 @@ ENDDO C compute convergence criterion for linear preconditioned FGMRES JFNKgamma_lin = JFNKgamma_lin_max - IF ( newtonIter.GT.1.AND.newtonIter.LE.100 + IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter & .AND.JFNKresidual.LT.JFNKres_t ) THEN C Eisenstat, 1996, equ.(2.6) phi_e = 1. _d 0 @@ -278,8 +254,13 @@ IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN krylovFails = krylovFails + 1 ENDIF -C Set the stopping criterion for the Newton iteration - IF ( newtonIter .EQ. 1 ) JFNKtol=JFNKgamma_nonlin*JFNKresidual +C Set the stopping criterion for the Newton iteration and the +C criterion for the transition from accurate to approximate FGMRES + IF ( newtonIter .EQ. 1 ) THEN + JFNKtol=JFNKgamma_nonlin*JFNKresidual + IF ( JFNKres_tFac .NE. UNSET_RL ) + & JFNKres_t = JFNKresidual * JFNKres_tFac + ENDIF C Update linear solution vector and return to Newton iteration C Do a linesearch if necessary, and compute a new residual. C Note that it should be possible to do the following operations @@ -419,6 +400,7 @@ RETURN END +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_JFNK_UPDATE C !INTERFACE: @@ -475,7 +457,8 @@ _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) -C Local variables: +C !LOCAL VARIABLES: +C === Local variables === C i,j,bi,bj :: loop indices INTEGER i,j,bi,bj INTEGER l @@ -496,17 +479,6 @@ facLS = 1. _d 0 doLineSearch = .TRUE. 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 C Create update DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) @@ -529,6 +501,18 @@ CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid) CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,resLoc,myThid) resLoc = SQRT(resLoc) +C Determine, if we need more iterations + doLineSearch = resLoc .GE. JFNKresidual +C Limit the maximum number of iterations arbitrarily to four + doLineSearch = doLineSearch .AND. l .LT. 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 C some output diagnostics IF ( debugLevel.GE.debLevA .AND. doLineSearch ) THEN _BEGIN_MASTER( myThid )