| 63 |
|
|
| 64 |
C !LOCAL VARIABLES: |
C !LOCAL VARIABLES: |
| 65 |
C === Local variables === |
C === Local variables === |
| 66 |
C i,j,bi,bj :: loop indices |
C i,j,k,bi,bj :: loop indices |
| 67 |
INTEGER i,j,bi,bj |
INTEGER i,j,k,bi,bj |
| 68 |
C loop indices |
C loop indices |
| 69 |
INTEGER newtonIter |
INTEGER newtonIter |
| 70 |
INTEGER krylovIter, krylovFails |
INTEGER krylovIter, krylovFails |
| 71 |
INTEGER totalKrylovItersLoc, totalNewtonItersLoc |
INTEGER totalKrylovItersLoc, totalNewtonItersLoc |
| 72 |
|
C FGMRES parameters |
| 73 |
|
C im :: size of Krylov space |
| 74 |
|
C ifgmres :: interation counter |
| 75 |
|
INTEGER im |
| 76 |
|
PARAMETER ( im = 50 ) |
| 77 |
|
INTEGER ifgmres |
| 78 |
C FGMRES flag that determines amount of output messages of fgmres |
C FGMRES flag that determines amount of output messages of fgmres |
| 79 |
INTEGER iOutFGMRES |
INTEGER iOutFGMRES |
| 80 |
C FGMRES flag that indicates what fgmres wants us to do next |
C FGMRES flag that indicates what fgmres wants us to do next |
| 109 |
_RL etaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
_RL etaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
| 110 |
_RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
_RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
| 111 |
_RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
_RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
| 112 |
|
C work arrays |
| 113 |
|
_RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy) |
| 114 |
|
_RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy) |
| 115 |
|
_RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy) |
| 116 |
CEOP |
CEOP |
| 117 |
|
|
| 118 |
C Initialise |
C Initialise |
| 237 |
|
|
| 238 |
krylovConverged = .FALSE. |
krylovConverged = .FALSE. |
| 239 |
FGMRESeps = JFNKgamma_lin * JFNKresidual |
FGMRESeps = JFNKgamma_lin * JFNKresidual |
| 240 |
|
C map first guess sol; it is zero because the solution is a correction |
| 241 |
|
CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid) |
| 242 |
|
C map rhs and change its sign because we are solving J*u = -F |
| 243 |
|
CALL SEAICE_MAP2VEC(nVec,-uIceRes,-vIceRes,rhs,.TRUE.,myThid) |
| 244 |
DO WHILE ( .NOT.krylovConverged ) |
DO WHILE ( .NOT.krylovConverged ) |
| 245 |
C solution vector sol = du/vIce |
C solution vector sol = du/vIce |
| 246 |
C residual vector (rhs) Fu = u/vIceRes |
C residual vector (rhs) Fu = u/vIceRes |
| 247 |
C output work vectors wk1, -> input work vector wk2 |
C output work vectors wk1, -> input work vector wk2 |
| 248 |
|
|
| 249 |
CALL SEAICE_FGMRES_DRIVER( |
C map preconditioner results or Jacobian times vector, |
| 250 |
I uIceRes, vIceRes, |
C stored in du/vIce to wk2, for iCode=0, wk2 is set to zero, |
| 251 |
U duIce, dvIce, iCode, |
C because du/vIce = 0 |
| 252 |
I FGMRESeps, iOutFGMRES, |
CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk2,.TRUE.,myThid) |
| 253 |
I newtonIter, krylovIter, myTime, myIter, myThid ) |
C |
| 254 |
|
CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter, |
| 255 |
|
U vv,w,wk1,wk2, |
| 256 |
|
I FGMRESeps,SEAICEkrylovIterMax,iOutFGMRES, |
| 257 |
|
U iCode, |
| 258 |
|
I myThid) |
| 259 |
|
C |
| 260 |
|
IF ( iCode .EQ. 0 ) THEN |
| 261 |
|
C map sol(ution) vector to du/vIce |
| 262 |
|
CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.FALSE.,myThid) |
| 263 |
|
ELSE |
| 264 |
|
C map work vector to du/vIce to either compute a preconditioner |
| 265 |
|
C solution (wk1=rhs) or a Jacobian times wk1 |
| 266 |
|
CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk1,.FALSE.,myThid) |
| 267 |
|
ENDIF |
| 268 |
|
C Fill overlaps in updated fields |
| 269 |
|
CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid) |
| 270 |
C FGMRES returns iCode either asking for an new preconditioned vector |
C FGMRES returns iCode either asking for an new preconditioned vector |
| 271 |
C or product of matrix (Jacobian) times vector. For iCode = 0, terminate |
C or product of matrix (Jacobian) times vector. For iCode = 0, terminate |
| 272 |
C iteration |
C iteration |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch