--- MITgcm/pkg/seaice/seaice_jfnk.F 2013/04/04 07:02:51 1.21 +++ MITgcm/pkg/seaice/seaice_jfnk.F 2017/06/08 15:33:50 1.31 @@ -1,7 +1,10 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.21 2013/04/04 07:02:51 mlosch Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.31 2017/06/08 15:33:50 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif C-- File seaice_jfnk.F: seaice jfnk dynamical solver S/R: C-- Contents @@ -66,6 +69,12 @@ INTEGER newtonIter INTEGER krylovIter, krylovFails INTEGER totalKrylovItersLoc, totalNewtonItersLoc +C FGMRES parameters +C im :: size of Krylov space +C ifgmres :: interation counter + INTEGER im + PARAMETER ( im = 50 ) + INTEGER ifgmres C FGMRES flag that determines amount of output messages of fgmres INTEGER iOutFGMRES C FGMRES flag that indicates what fgmres wants us to do next @@ -73,10 +82,12 @@ _RL JFNKresidual _RL JFNKresidualKm1 C parameters to compute convergence criterion - _RL phi_e, alp_e, JFNKgamma_lin + _RL JFNKgamma_lin _RL FGMRESeps _RL JFNKtol - +C backward differences extrapolation factors + _RL bdfFac, bdfAlpha +C _RL recip_deltaT LOGICAL JFNKconverged, krylovConverged LOGICAL writeNow @@ -85,15 +96,23 @@ C u/vIceRes :: residual of sea-ice momentum equations _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 extra time level required for backward difference time stepping + _RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C du/vIce :: ice velocity increment to be added to u/vIce _RL duIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL dvIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C precomputed (= constant per Newton iteration) versions of C zeta, eta, and DWATN, press _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL zetaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +C work arrays + _RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy) + _RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy) + _RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy) CEOP C Initialise @@ -113,6 +132,17 @@ & DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) ) & iOutFGMRES=1 +C backward difference extrapolation factors + bdfFac = 0. _d 0 + IF ( SEAICEuseBDF2 ) THEN + IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN + bdfFac = 0. _d 0 + ELSE + bdfFac = 0.5 _d 0 + ENDIF + ENDIF + bdfAlpha = 1. _d 0 + bdfFac + DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy,sNy+OLy @@ -121,25 +151,37 @@ vIceRes(I,J,bi,bj) = 0. _d 0 duIce (I,J,bi,bj) = 0. _d 0 dvIce (I,J,bi,bj) = 0. _d 0 + ENDDO + ENDDO +C cycle ice velocities + DO J=1-OLy,sNy+OLy + DO I=1-OLx,sNx+OLx + duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha + & + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac + dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha + & + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj) vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj) ENDDO ENDDO +C As long as IMEX is not properly implemented leave this commented out +CML IF ( .NOT.SEAICEuseIMEX ) THEN C Compute things that do no change during the Newton iteration: C sea-surface tilt and wind stress: -C FORCEX/Y0 - mass*(u/vIceNm1)/deltaT +C FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj) - & + seaiceMassU(I,J,bi,bj)*uIceNm1(I,J,bi,bj)*recip_deltaT + & + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj) - & + seaiceMassV(I,J,bi,bj)*vIceNm1(I,J,bi,bj)*recip_deltaT + & + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT ENDDO ENDDO +CML ENDIF ENDDO ENDDO C Start nonlinear Newton iteration: outer loop iteration - DO WHILE ( newtonIter.LT.SEAICEnewtonIterMax .AND. + DO WHILE ( newtonIter.LT.SEAICEnonLinIterMax .AND. & .NOT.JFNKconverged ) newtonIter = newtonIter + 1 C Compute initial residual F(u), (includes computation of global @@ -156,6 +198,7 @@ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx zetaPre(I,J,bi,bj) = zeta(I,J,bi,bj) + zetaZPre(I,J,bi,bj)= zetaZ(I,J,bi,bj) etaPre(I,J,bi,bj) = eta(I,J,bi,bj) etaZPre(I,J,bi,bj) = etaZ(I,J,bi,bj) dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj) @@ -167,10 +210,9 @@ JFNKgamma_lin = JFNKgamma_lin_max 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 - alp_e = 1. _d 0 - JFNKgamma_lin = phi_e*( JFNKresidual/JFNKresidualKm1 )**alp_e +C Eisenstat and Walker (1996), eq.(2.6) + JFNKgamma_lin = SEAICE_JFNKphi + & *( JFNKresidual/JFNKresidualKm1 )**SEAICE_JFNKalpha JFNKgamma_lin = min(JFNKgamma_lin_max, JFNKgamma_lin) JFNKgamma_lin = max(JFNKgamma_lin_min, JFNKgamma_lin) ENDIF @@ -186,6 +228,7 @@ iCode = 0 JFNKconverged = JFNKresidual.LT.JFNKtol + & .OR.JFNKresidual.EQ.0. _d 0 C do Krylov loop only if convergence is not reached @@ -195,25 +238,52 @@ krylovConverged = .FALSE. FGMRESeps = JFNKgamma_lin * JFNKresidual +C map first guess sol; it is zero because the solution is a correction + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid) +C map rhs and change its sign because we are solving J*u = -F + CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,rhs,.TRUE.,myThid) + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1,nVec + rhs(j,bi,bj) = - rhs(j,bi,bj) + ENDDO + ENDDO + ENDDO DO WHILE ( .NOT.krylovConverged ) C solution vector sol = du/vIce C residual vector (rhs) Fu = u/vIceRes C output work vectors wk1, -> input work vector wk2 - CALL SEAICE_FGMRES_DRIVER( - I uIceRes, vIceRes, - U duIce, dvIce, iCode, - I FGMRESeps, iOutFGMRES, - I newtonIter, krylovIter, myTime, myIter, myThid ) +C map preconditioner results or Jacobian times vector, +C stored in du/vIce to wk2, for iCode=0, wk2 is set to zero, +C because du/vIce = 0 + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk2,.TRUE.,myThid) +C + CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter, + U vv,w,wk1,wk2, + I FGMRESeps,SEAICElinearIterMax,iOutFGMRES, + U iCode, + I myThid) +C + IF ( iCode .EQ. 0 ) THEN +C map sol(ution) vector to du/vIce + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.FALSE.,myThid) + ELSE +C map work vector to du/vIce to either compute a preconditioner +C solution (wk1=rhs) or a Jacobian times wk1 + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk1,.FALSE.,myThid) + ENDIF +C Fill overlaps in updated fields + CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid) C FGMRES returns iCode either asking for an new preconditioned vector C or product of matrix (Jacobian) times vector. For iCode = 0, terminate C iteration IF (iCode.EQ.1) THEN C Call preconditioner - IF ( SOLV_MAX_ITERS .GT. 0 ) + IF ( SEAICEpreconLinIter .GT. 0 ) & CALL SEAICE_PRECONDITIONER( U duIce, dvIce, - I zetaPre, etaPre, etaZpre, dwatPre, + I zetaPre, etaPre, etaZpre, zetaZpre, dwatPre, I newtonIter, krylovIter, myTime, myIter, myThid ) ELSEIF (iCode.GE.2) THEN C Compute Jacobian times vector @@ -230,7 +300,7 @@ IF ( debugLevel.GE.debLevA ) THEN _BEGIN_MASTER( myThid ) totalNewtonItersLoc = - & SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter + & SEAICEnonLinIterMax*(myIter-nIter0)+newtonIter WRITE(msgBuf,'(2A,2(1XI6),2E12.5)') & ' S/R SEAICE_JFNK: Newton iterate / total, ', & 'JFNKgamma_lin, initial norm = ', @@ -246,13 +316,13 @@ & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) ENDIF - IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN + IF ( krylovIter.EQ.SEAICElinearIterMax ) THEN krylovFails = krylovFails + 1 ENDIF 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 + JFNKtol=SEAICEnonLinTol*JFNKresidual IF ( JFNKres_tFac .NE. UNSET_RL ) & JFNKres_t = JFNKresidual * JFNKres_tFac ENDIF @@ -291,7 +361,7 @@ totalKrylovIters = totalKrylovIters + totalKrylovItersLoc C Record failure totalKrylovFails = totalKrylovFails + krylovFails - IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN + IF ( newtonIter .EQ. SEAICEnonLinIterMax ) THEN totalNewtonFails = totalNewtonFails + 1 ENDIF ENDIF @@ -365,7 +435,7 @@ C Print more debugging information IF ( debugLevel.GE.debLevA ) THEN - IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN + IF ( newtonIter .EQ. SEAICEnonLinIterMax ) THEN _BEGIN_MASTER( myThid ) WRITE(msgBuf,'(A,I10)') & ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',