--- MITgcm/pkg/seaice/seaice_jfnk.F 2013/05/30 14:07:19 1.23 +++ MITgcm/pkg/seaice/seaice_jfnk.F 2016/01/28 12:54:12 1.30 @@ -1,7 +1,10 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.23 2013/05/30 14:07:19 mlosch Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.30 2016/01/28 12:54:12 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 @@ -76,8 +85,8 @@ _RL JFNKgamma_lin _RL FGMRESeps _RL JFNKtol -C Adams-Bashforth extrapolation factors - _RL abFac, abAlpha +C backward differences extrapolation factors + _RL bdfFac, bdfAlpha C _RL recip_deltaT LOGICAL JFNKconverged, krylovConverged @@ -87,7 +96,7 @@ 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 Adams-Bashforth-2 time stepping +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 @@ -96,9 +105,14 @@ 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 @@ -118,16 +132,16 @@ & DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) ) & iOutFGMRES=1 -C Adams-Bashforth extrapolation factors - abFac = 0. _d 0 - IF ( SEAICEuseAB2 ) THEN - IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartAB.EQ.0 ) THEN - abFac = 0. _d 0 +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 - abFac = 0.5 _d 0 + SEAICE_abEps + bdfFac = 0.5 _d 0 ENDIF ENDIF - abAlpha = 1. _d 0 + abFac + bdfAlpha = 1. _d 0 + bdfFac DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) @@ -142,18 +156,19 @@ 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) * abAlpha - & + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * abFac - dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * abAlpha - & + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * abFac + 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 - IF ( .NOT.SEAICEuseIMEX ) THEN +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*(abA*u/vIceNm1+abB*(u/vIceNm1-u/vIceNm2))/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) @@ -162,11 +177,11 @@ & + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT ENDDO ENDDO - ENDIF +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 @@ -183,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) @@ -221,25 +237,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 @@ -256,7 +299,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 = ', @@ -272,13 +315,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 @@ -317,7 +360,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 @@ -391,7 +434,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 ',