/[MITgcm]/MITgcm/model/src/cg2d.F
ViewVC logotype

Diff of /MITgcm/model/src/cg2d.F

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

revision 1.44 by jmc, Fri Nov 4 01:19:25 2005 UTC revision 1.45 by jmc, Tue Dec 5 05:25:08 2006 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: CG2D  C     !ROUTINE: CG2D
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE CG2D(          SUBROUTINE CG2D(
10       I                cg2d_b,       I                cg2d_b,
11       U                cg2d_x,       U                cg2d_x,
12       O                firstResidual,       O                firstResidual,
# Line 15  C     !INTERFACE: Line 15  C     !INTERFACE:
15       I                myThid )       I                myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | SUBROUTINE CG2D                                            C     | SUBROUTINE CG2D
19  C     | o Two-dimensional grid problem conjugate-gradient          C     | o Two-dimensional grid problem conjugate-gradient
20  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
21  C     *==========================================================*  C     *==========================================================*
22  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
23  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
24  C     | This implementation assumes A is a five-diagonal            C     | This implementation assumes A is a five-diagonal
25  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
26  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
27  C     | two-dimensional space.                                      C     | two-dimensional space.
28  C     | Notes:                                                      C     | Notes:
29  C     | ======                                                      C     | ======
30  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
31  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
32  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
33  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
34  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
35  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
36  C     | the threads to solve several different problems              C     | the threads to solve several different problems
37  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
38  C     *==========================================================*  C     *==========================================================*
39  C     \ev  C     \ev
40    
# Line 44  C     === Global data === Line 44  C     === Global data ===
44  #include "SIZE.h"  #include "SIZE.h"
45  #include "EEPARAMS.h"  #include "EEPARAMS.h"
46  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
47  #include "CG2D.h"  #include "CG2D.h"
48    #include "GRID.h"
49  #include "SURFACE.h"  #include "SURFACE.h"
50    
51  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
52  C     === Routine arguments ===  C     === Routine arguments ===
53  C     myThid    - Thread on which I am working.  C     myThid    :: Thread on which I am working.
54  C     cg2d_b    - The source term or "right hand side"  C     cg2d_b    :: The source term or "right hand side"
55  C     cg2d_x    - The solution  C     cg2d_x    :: The solution
56  C     firstResidual - the initial residual before any iterations  C     firstResidual :: the initial residual before any iterations
57  C     lastResidual  - the actual residual reached  C     lastResidual  :: the actual residual reached
58  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  :: Entry: the maximum number of iterations allowed
59  C                 Exit:  the actual number of iterations used  C                  Exit:  the actual number of iterations used
60        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
61        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62        _RL  firstResidual        _RL  firstResidual
# Line 66  C                 Exit:  the actual numb Line 66  C                 Exit:  the actual numb
66    
67  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
68  C     === Local variables ====  C     === Local variables ====
69  C     actualIts      - Number of iterations taken  C     actualIts      :: Number of iterations taken
70  C     actualResidual - residual  C     actualResidual :: residual
71  C     bi          - Block index in X and Y.  C     bi, bj     :: Block index in X and Y.
72  C     bj  C     eta_qrN    :: Used in computing search directions
 C     eta_qrN     - Used in computing search directions  
73  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
74  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
75  C     alpha    C     alpha
76  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a
77  C                   useful debuggin/trouble shooting diagnostic.  C                   useful debuggin/trouble shooting diagnostic.
78  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
79  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
80  C     err         - Measure of residual of Ax - b, usually the norm.  C     err        :: Measure of residual of Ax - b, usually the norm.
81  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, it2d :: Loop counters ( it2d counts CG iterations )
82        INTEGER actualIts        INTEGER actualIts
83        _RL    actualResidual        _RL    actualResidual
84        INTEGER bi, bj                      INTEGER bi, bj
85        INTEGER I, J, it2d        INTEGER I, J, it2d
86          INTEGER ks
87        _RL    err,errTile        _RL    err,errTile
88        _RL    eta_qrN,eta_qrNtile        _RL    eta_qrN,eta_qrNtile
89        _RL    eta_qrNM1        _RL    eta_qrNM1
# Line 131  C     _GLOBAL_MAX_R8( rhsMax, myThid ) Line 131  C     _GLOBAL_MAX_R8( rhsMax, myThid )
131        rhsMax=1.        rhsMax=1.
132  #else  #else
133        _GLOBAL_MAX_R8( rhsMax, myThid )        _GLOBAL_MAX_R8( rhsMax, myThid )
 Catm  rhsMax=1.  
134  #endif  #endif
135        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
136        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
# Line 166  C--   Initial residual calculation Line 165  C--   Initial residual calculation
165          errTile    = 0. _d 0          errTile    = 0. _d 0
166          DO J=1,sNy          DO J=1,sNy
167           DO I=1,sNx           DO I=1,sNx
168              ks = ksurfC(I,J,bi,bj)
169            cg2d_s(I,J,bi,bj) = 0.            cg2d_s(I,J,bi,bj) = 0.
170            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -
171       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)
# Line 176  C--   Initial residual calculation Line 176  C--   Initial residual calculation
176       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
177       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
178       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
179       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*       &    -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)*recip_Bo(i,j,bi,bj)
180       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm       &                        *cg2d_x(I  ,J  ,bi,bj)
181         &                        /deltaTMom/deltaTfreesurf*cg2dNorm
182    c    &    +aC2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
183       &    )       &    )
184            errTile        = errTile        +            errTile    = errTile + cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
185       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)            sumRHStile = sumRHStile + cg2d_b(I,J,bi,bj)
           sumRHStile        = sumRHStile        +  
      &     cg2d_b(I,J,bi,bj)  
186           ENDDO           ENDDO
187          ENDDO          ENDDO
188          sumRHS = sumRHS + sumRHStile          sumRHS = sumRHS + sumRHStile
# Line 209  C     _EXCH_XY_R8( cg2d_s, myThid ) Line 209  C     _EXCH_XY_R8( cg2d_s, myThid )
209    
210         IF ( debugLevel .GE. debLevZero ) THEN         IF ( debugLevel .GE. debLevZero ) THEN
211          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
212          write(standardmessageunit,'(A,1P2E22.14)')          WRITE(standardmessageunit,'(A,1P2E22.14)')
213       &  ' cg2d: Sum(rhs),rhsMax = ',       &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
      &                                  sumRHS,rhsMax  
214          _END_MASTER( myThid )          _END_MASTER( myThid )
215         ENDIF         ENDIF
216  C     _BARRIER  C     _BARRIER
217  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
218  c      WRITE(standardmessageunit,'(A,I6,1PE30.14)')  c      WRITE(standardmessageunit,'(A,I6,1PE30.14)')
219  c    & ' CG2D iters, err = ',  c    & ' CG2D iters, err = ',
220  c    & actualIts, actualResidual  c    & actualIts, actualResidual
221  c     _END_MASTER( myThid )  c     _END_MASTER( myThid )
222        firstResidual=actualResidual        firstResidual=actualResidual
# Line 239  C--    conjugate direction vector "s". Line 238  C--    conjugate direction vector "s".
238           eta_qrNtile = 0. _d 0           eta_qrNtile = 0. _d 0
239           DO J=1,sNy           DO J=1,sNy
240            DO I=1,sNx            DO I=1,sNx
241             cg2d_q(I,J,bi,bj) =             cg2d_q(I,J,bi,bj) =
242       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)
243       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)
244       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)
# Line 294  C==    q = A.s Line 293  C==    q = A.s
293           alphaTile = 0. _d 0           alphaTile = 0. _d 0
294           DO J=1,sNy           DO J=1,sNy
295            DO I=1,sNx            DO I=1,sNx
296             cg2d_q(I,J,bi,bj) =             ks = ksurfC(I,J,bi,bj)
297               cg2d_q(I,J,bi,bj) =
298       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)
299       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)
300       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)
# Line 303  C==    q = A.s Line 303  C==    q = A.s
303       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
304       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
305       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
306       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*       &    -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)*recip_Bo(i,j,bi,bj)
307       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm       &                        *cg2d_s(I  ,J  ,bi,bj)
308         &                        /deltaTMom/deltaTfreesurf*cg2dNorm
309    c    &    +aC2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
310            alphaTile = alphaTile+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)            alphaTile = alphaTile+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
311            ENDDO            ENDDO
312           ENDDO           ENDDO
# Line 319  CcnhDebugEnds Line 321  CcnhDebugEnds
321  CcnhDebugStarts  CcnhDebugStarts
322  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
323  CcnhDebugEnds  CcnhDebugEnds
324        
325  C==    Update solution and residual vectors  C==    Update solution and residual vectors
326  C      Now compute "interior" points.  C      Now compute "interior" points.
327         err = 0. _d 0         err = 0. _d 0
# Line 369  C     The following exchange was moved u Line 371  C     The following exchange was moved u
371  C     for compatibility with TAMC.  C     for compatibility with TAMC.
372  C     _EXCH_XY_R8(cg2d_x, myThid )  C     _EXCH_XY_R8(cg2d_x, myThid )
373  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
374  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
375  c    & actualIts, actualResidual  c    & actualIts, actualResidual
376  c     _END_MASTER( myThid )  c     _END_MASTER( myThid )
377    
# Line 393  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I Line 395  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I
395  C    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  C    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
396  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
397  C    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj))  C    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj))
398  C         err            = err            +  C         err            = err            +
399  C    &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)  C    &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
400  C        ENDDO  C        ENDDO
401  C       ENDDO  C       ENDDO
Legend:
Removed from v.1.44  
changed lines
  Added in v.1.45
  ViewVC Help
Powered by ViewVC 1.1.22