/[MITgcm]/MITgcm/pkg/diagnostics/diag_cg2d.F
ViewVC logotype

Diff of /MITgcm/pkg/diagnostics/diag_cg2d.F

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

revision 1.1 by jmc, Tue Jun 14 00:18:36 2011 UTC revision 1.2 by jmc, Fri Jul 22 19:53:39 2011 UTC
# Line 9  C     !INTERFACE: Line 9  C     !INTERFACE:
9        SUBROUTINE DIAG_CG2D(        SUBROUTINE DIAG_CG2D(
10       I                aW2d, aS2d, b2d,       I                aW2d, aS2d, b2d,
11       I                residCriter,       I                residCriter,
12       O                firstResidual, lastResidual,       O                firstResidual, minResidual, lastResidual,
13       U                x2d, numIters,       U                x2d, numIters,
14       I                myThid )       O                nIterMin,
15         I                printResidFrq, myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | SUBROUTINE CG2D  C     | SUBROUTINE CG2D
# Line 39  c#include "SURFACE.h" Line 40  c#include "SURFACE.h"
40    
41  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
42  C     === Routine arguments ===  C     === Routine arguments ===
43  C     b2d    :: The source term or "right hand side"  C     b2d           :: The source term or "right hand side"
44  C     x2d    :: The solution  C     x2d           :: The solution
45  C     firstResidual :: the initial residual before any iterations  C     firstResidual :: the initial residual before any iterations
46    C     minResidual   :: the lowest residual reached
47  C     lastResidual  :: the actual residual reached  C     lastResidual  :: the actual residual reached
48  C     numIters  :: Entry: the maximum number of iterations allowed  C     numIters  :: Entry: the maximum number of iterations allowed
49  C                  Exit:  the actual number of iterations used  C                  Exit:  the actual number of iterations used
50  C     myThid    :: Thread on which I am working.  C     nIterMin      :: iteration number corresponding to lowest residual
51    C     printResidFrq :: Frequency for printing residual in CG iterations
52    C     myThid        :: my Thread Id number
53        _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54        _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55        _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56        _RL  residCriter        _RL  residCriter
57        _RL  firstResidual        _RL  firstResidual
58          _RL  minResidual
59        _RL  lastResidual        _RL  lastResidual
60        _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
61        INTEGER numIters        INTEGER numIters
62          INTEGER nIterMin
63          INTEGER printResidFrq
64        INTEGER myThid        INTEGER myThid
65    
66  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
67  C     === Local variables ====  C     === Local variables ====
68  C     actualIts      :: Number of iterations taken  C     bi, bj     :: tile indices
 C     actualResidual :: residual  
 C     bi, bj     :: Block index in X and Y.  
69  C     eta_qrN    :: Used in computing search directions  C     eta_qrN    :: Used in computing search directions
70  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
71  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
# Line 69  C     sumRHS     :: Sum of right-hand-si Line 74  C     sumRHS     :: Sum of right-hand-si
74  C                   useful debuggin/trouble shooting diagnostic.  C                   useful debuggin/trouble shooting diagnostic.
75  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
76  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
77  C     err        :: Measure of residual of Ax - b, usually the norm.  C     err        :: Measure of current residual of Ax - b, usually the norm.
78  C     i, j, it2d :: Loop counters ( it2d counts CG iterations )  C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
       INTEGER actualIts  
       _RL    actualResidual  
79        INTEGER bi, bj        INTEGER bi, bj
80        INTEGER i, j, it2d        INTEGER i, j, it2d
81        _RL    err,     errTile(nSx,nSy)        _RL    err,     errTile(nSx,nSy)
# Line 92  CEOP Line 95  CEOP
95        _RL  q2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)        _RL  q2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
96        _RL  r2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)        _RL  r2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
97        _RL  s2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)        _RL  s2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
98          _RL  x2dm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
99    
100  C--   Set matrice main diagnonal:  C--   Set matrice main diagnonal:
101        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 147  C--   Initial residual calculation Line 151  C--   Initial residual calculation
151          DO j=1-1,sNy+1          DO j=1-1,sNy+1
152           DO i=1-1,sNx+1           DO i=1-1,sNx+1
153            s2d(i,j,bi,bj) = 0.            s2d(i,j,bi,bj) = 0.
154              x2dm(i,j,bi,bj) = x2d(i,j,bi,bj)
155           ENDDO           ENDDO
156          ENDDO          ENDDO
157          sumRHStile(bi,bj) = 0. _d 0          sumRHStile(bi,bj) = 0. _d 0
# Line 171  C--   Initial residual calculation Line 176  C--   Initial residual calculation
176        CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )        CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
177        CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )        CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
178        err = SQRT(err)        err = SQRT(err)
179        actualIts      = 0        it2d = 0
180        actualResidual = err        firstResidual = err
181        firstResidual  = err        minResidual   = err
182          nIterMin = it2d
183    
184        printResidual = .FALSE.        printResidual = .FALSE.
185        IF ( debugLevel .GE. debLevZero ) THEN        IF ( debugLevel .GE. debLevZero ) THEN
186          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
187          printResidual = printResidualFreq.GE.1          printResidual = printResidFrq.GE.1
188  c       WRITE(standardmessageunit,'(A,1P2E22.14)')  c       WRITE(standardmessageunit,'(A,1P2E22.14)')
189  c    &  ' diag_cg2d: Sum(rhs),rhsMax = ', sumRHS, rhsMax  c    &  ' diag_cg2d: Sum(rhs),rhsMax = ', sumRHS, rhsMax
190          _END_MASTER( myThid )          _END_MASTER( myThid )
# Line 267  C      Now compute "interior" points. Line 273  C      Now compute "interior" points.
273    
274         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
275         err = SQRT(err)         err = SQRT(err)
        actualIts      = it2d  
        actualResidual = err  
276         IF ( printResidual ) THEN         IF ( printResidual ) THEN
277          IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN          IF ( MOD( it2d-1, printResidFrq ).EQ.0 ) THEN
278           WRITE(msgBuf,'(A,I6,A,1PE21.14)')           WRITE(msgBuf,'(A,I6,A,1PE21.14)')
279       &    ' diag_cg2d: iter=', actualIts, ' ; resid.= ', actualResidual       &    ' diag_cg2d: iter=', it2d, ' ; resid.= ', err
280           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
281       &                       SQUEEZE_RIGHT, myThid )       &                       SQUEEZE_RIGHT, myThid )
282          ENDIF          ENDIF
283         ENDIF         ENDIF
284         IF ( err .LT. residCriter ) GOTO 11         IF ( err .LT. residCriter ) GOTO 11
285           IF ( err .LT. minResidual ) THEN
286    C-     Store lowest residual solution
287             minResidual = err
288             nIterMin = it2d
289             DO bj=myByLo(myThid),myByHi(myThid)
290              DO bi=myBxLo(myThid),myBxHi(myThid)
291               DO j=1,sNy
292                DO i=1,sNx
293                 x2dm(i,j,bi,bj) = x2d(i,j,bi,bj)
294                ENDDO
295               ENDDO
296              ENDDO
297             ENDDO
298           ENDIF
299    
300         CALL EXCH_S3D_RL( r2d, 1, myThid )         CALL EXCH_S3D_RL( r2d, 1, myThid )
301    
302     10 CONTINUE     10 CONTINUE
303          it2d = numIters
304     11 CONTINUE     11 CONTINUE
305    
 c     CALL EXCH_XY_RL( x2d, myThid )  
   
306  C--   Return parameters to caller  C--   Return parameters to caller
307        lastResidual = actualResidual        lastResidual = err
308        numIters = actualIts        numIters = it2d
309    
310          IF ( err .GT. minResidual ) THEN
311    C-    use the lowest residual solution (instead of current one <-> last residual)
312            DO bj=myByLo(myThid),myByHi(myThid)
313             DO bi=myBxLo(myThid),myBxHi(myThid)
314              DO j=1,sNy
315               DO i=1,sNx
316                 x2d(i,j,bi,bj) = x2dm(i,j,bi,bj)
317               ENDDO
318              ENDDO
319             ENDDO
320            ENDDO
321          ENDIF
322    c     CALL EXCH_XY_RL( x2d, myThid )
323    
324        RETURN        RETURN
325        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22