/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_fgmres.F
ViewVC logotype

Diff of /MITgcm_contrib/torge/itd/code/seaice_fgmres.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by torge, Fri Nov 2 16:52:08 2012 UTC revision 1.3 by torge, Mon Dec 10 22:19:49 2012 UTC
# Line 18  C     !INTERFACE: Line 18  C     !INTERFACE:
18       I     uIceRes, vIceRes,       I     uIceRes, vIceRes,
19       U     duIce, dvIce,       U     duIce, dvIce,
20       U     iCode,       U     iCode,
21       I     FGMRESeps,       I     FGMRESeps, iOutFGMRES,
22       I     newtonIter, krylovIter, myTime, myIter, myThid )       I     newtonIter, krylovIter, myTime, myIter, myThid )
23    
24  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
# Line 48  C     myIter :: Simulation timestep numb Line 48  C     myIter :: Simulation timestep numb
48  C     myThid :: my Thread Id. number  C     myThid :: my Thread Id. number
49  C     newtonIter :: current iterate of Newton iteration  C     newtonIter :: current iterate of Newton iteration
50  C     krylovIter :: current iterate of Newton iteration  C     krylovIter :: current iterate of Newton iteration
51  C     iCode  :: FGMRES parameter to determine next step  C     iCode      :: FGMRES parameter to determine next step
52    C     iOutFGMRES :: control output of fgmres
53        _RL     myTime        _RL     myTime
54        INTEGER myIter        INTEGER myIter
55        INTEGER myThid        INTEGER myThid
56        INTEGER newtonIter        INTEGER newtonIter
57        INTEGER krylovIter        INTEGER krylovIter
58          INTEGER iOutFGMRES
59        INTEGER iCode        INTEGER iCode
60  C     FGMRESeps :: tolerance for FGMRES  C     FGMRESeps :: tolerance for FGMRES
61        _RL FGMRESeps        _RL FGMRESeps
# Line 74  C     FGMRES parameters Line 76  C     FGMRES parameters
76  C     n  :: size of the input vector(s)  C     n  :: size of the input vector(s)
77  C     im :: size of Krylov space  C     im :: size of Krylov space
78  C     ifgmres :: interation counter  C     ifgmres :: interation counter
 C     iout :: control output of fgmres  
79        INTEGER n        INTEGER n
80        PARAMETER ( n  = 2*sNx*sNy*nSx*nSy )        PARAMETER ( n  = 2*sNx*sNy*nSx*nSy )
81        INTEGER im        INTEGER im
82        PARAMETER ( im = 50 )        PARAMETER ( im = 50 )
83        INTEGER ifgmres, iout        INTEGER ifgmres
84  C     work arrays  C     work arrays
85        _RL rhs(n), sol(n)        _RL rhs(n), sol(n)
86        _RL vv(n,im+1), w(n,im)        _RL vv(n,im+1), w(n,im)
# Line 90  C     they are not forgotten between Kry Line 91  C     they are not forgotten between Kry
91        COMMON /FGMRES_RL/ sol, rhs, vv, w        COMMON /FGMRES_RL/ sol, rhs, vv, w
92  CEOP  CEOP
93    
 C     iout=1 give a little bit of output  
       iout=1  
   
94  C     For now, let only the master thread do all the work  C     For now, let only the master thread do all the work
95  C     - copy from 2D arrays to 1D-vector  C     - copy from 2D arrays to 1D-vector
96  C     - perform fgmres step (including global sums)  C     - perform fgmres step (including global sums)
# Line 101  C     not sure if this works properly Line 99  C     not sure if this works properly
99    
100        _BEGIN_MASTER ( myThid )        _BEGIN_MASTER ( myThid )
101        IF ( iCode .EQ. 0 ) THEN        IF ( iCode .EQ. 0 ) THEN
102  C     first guess is zero because it is a correction  C     The first guess is zero because it is a correction, but this
103    C     is implemented by setting du/vIce=0 outside of this routine;
104    C     this make it possible to restart FGMRES with a nonzero sol
105           CALL SEAICE_MAP2VEC(n,duIce,dvIce,sol,.TRUE.,myThid)
106  C     wk2 needs to be reset for iCode = 0, because it may contain  C     wk2 needs to be reset for iCode = 0, because it may contain
107  C     remains of the previous Krylov iteration  C     remains of the previous Krylov iteration
108         DO k=1,n         DO k=1,n
         sol(k) = 0. _d 0  
109          wk2(k) = 0. _d 0          wk2(k) = 0. _d 0
110         ENDDO         ENDDO
111        ELSEIF ( iCode .EQ. 3 ) THEN        ELSEIF ( iCode .EQ. 3 ) THEN
# Line 125  C     stored in du/vIce to wk2 Line 125  C     stored in du/vIce to wk2
125  C  C
126        CALL SEAICE_FGMRES (n,im,rhs,sol,ifgmres,vv,w,wk1,wk2,        CALL SEAICE_FGMRES (n,im,rhs,sol,ifgmres,vv,w,wk1,wk2,
127       &     FGMRESeps,SEAICEkrylovIterMax,       &     FGMRESeps,SEAICEkrylovIterMax,
128       &     iout,icode,krylovIter,myThid)       &     iOutFGMRES,iCode,krylovIter,myThid)
129  C  C
130        IF ( iCode .EQ. 0 ) THEN        IF ( iCode .EQ. 0 ) THEN
131  C     map sol(ution) vector to du/vIce  C     map sol(ution) vector to du/vIce
# Line 511  C         call daxpy(n, t, w(1,j), 1, so Line 511  C         call daxpy(n, t, w(1,j), 1, so
511  C  C
512  C     test for return  C     test for return
513  C  C
       print *, 'ml-fgmres: its, maxits: ', its, maxits, ro, '<', eps1  
514        if (ro .le. eps1 .or. its .ge. maxits) goto 999        if (ro .le. eps1 .or. its .ge. maxits) goto 999
515  C  C
516  C     else compute residual vector and continue..  C     else compute residual vector and continue..
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22