2 |
C $Name$ |
C $Name$ |
3 |
|
|
4 |
#include "SEAICE_OPTIONS.h" |
#include "SEAICE_OPTIONS.h" |
5 |
|
#ifdef ALLOW_AUTODIFF |
6 |
|
# include "AUTODIFF_OPTIONS.h" |
7 |
|
#endif |
8 |
|
|
9 |
C-- File seaice_jfnk.F: seaice jfnk dynamical solver S/R: |
C-- File seaice_jfnk.F: seaice jfnk dynamical solver S/R: |
10 |
C-- Contents |
C-- Contents |
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 |
105 |
C precomputed (= constant per Newton iteration) versions of |
C precomputed (= constant per Newton iteration) versions of |
106 |
C zeta, eta, and DWATN, press |
C zeta, eta, and DWATN, press |
107 |
_RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
_RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
108 |
|
_RL zetaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
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 |
198 |
DO j=1-OLy,sNy+OLy |
DO j=1-OLy,sNy+OLy |
199 |
DO i=1-OLx,sNx+OLx |
DO i=1-OLx,sNx+OLx |
200 |
zetaPre(I,J,bi,bj) = zeta(I,J,bi,bj) |
zetaPre(I,J,bi,bj) = zeta(I,J,bi,bj) |
201 |
|
zetaZPre(I,J,bi,bj)= zetaZ(I,J,bi,bj) |
202 |
etaPre(I,J,bi,bj) = eta(I,J,bi,bj) |
etaPre(I,J,bi,bj) = eta(I,J,bi,bj) |
203 |
etaZPre(I,J,bi,bj) = etaZ(I,J,bi,bj) |
etaZPre(I,J,bi,bj) = etaZ(I,J,bi,bj) |
204 |
dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj) |
dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj) |
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 |
275 |
IF ( SOLV_MAX_ITERS .GT. 0 ) |
IF ( SOLV_MAX_ITERS .GT. 0 ) |
276 |
& CALL SEAICE_PRECONDITIONER( |
& CALL SEAICE_PRECONDITIONER( |
277 |
U duIce, dvIce, |
U duIce, dvIce, |
278 |
I zetaPre, etaPre, etaZpre, dwatPre, |
I zetaPre, etaPre, etaZpre, zetaZpre, dwatPre, |
279 |
I newtonIter, krylovIter, myTime, myIter, myThid ) |
I newtonIter, krylovIter, myTime, myIter, myThid ) |
280 |
ELSEIF (iCode.GE.2) THEN |
ELSEIF (iCode.GE.2) THEN |
281 |
C Compute Jacobian times vector |
C Compute Jacobian times vector |