/[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.18 by adcroft, Wed Dec 9 16:11:51 1998 UTC revision 1.30 by adcroft, Fri Mar 9 20:45:09 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6        SUBROUTINE CG2D(          SUBROUTINE CG2D(  
7       I                cg2d_b,       I                cg2d_b,
8       U                cg2d_x,       U                cg2d_x,
9         U                tolerance,
10         O                residual,
11         U                numIters,
12       I                myThid )       I                myThid )
13  C     /==========================================================\  C     /==========================================================\
14  C     | SUBROUTINE CG2D                                          |  C     | SUBROUTINE CG2D                                          |
# Line 36  C     === Global data === Line 40  C     === Global data ===
40  #include "PARAMS.h"  #include "PARAMS.h"
41  #include "GRID.h"  #include "GRID.h"
42  #include "CG2D_INTERNAL.h"  #include "CG2D_INTERNAL.h"
43    #include "SURFACE.h"
44    
45  C     === Routine arguments ===  C     === Routine arguments ===
46  C     myThid - Thread on which I am working.  C     myThid    - Thread on which I am working.
47        INTEGER myThid  C     cg2d_b    - The source term or "right hand side"
48        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     cg2d_x    - The solution
49    C     tolerance - Entry: the tolerance of accuracy to solve to
50    C                 Exit:  the initial residual before any iterations
51    C     residual  - Exit:  the actual residual reached
52    C     numIters  - Entry: the maximum number of iterations allowed
53    C                 Exit:  the actual number of iterations used
54        _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)
55          _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56          _RL  tolerance
57          _RL  residual
58          INTEGER numIters
59          INTEGER myThid
60    
61  C     === Local variables ====  C     === Local variables ====
62  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
63  C     actualResidual - residual  C     actualResidual - residual
64  C     bi          - Block index in X and Y.  C     bi          - Block index in X and Y.
65  C     bj  C     bj
66  C     etaN        - Used in computing search directions  C     eta_qrN     - Used in computing search directions
67  C     etaNM1        suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
68  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
69  C     alpha    C     alpha  
70  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS      - Sum of right-hand-side. Sometimes this is a
# Line 60  C                   or they converge at Line 74  C                   or they converge at
74  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
75  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
76        INTEGER actualIts        INTEGER actualIts
77        _RL    actualResidual        _RL    actualResidual,initialResidual
78        INTEGER bi, bj                      INTEGER bi, bj              
79        INTEGER I, J, it2d        INTEGER I, J, it2d
80        _RL    err        _RL    err
81        _RL    etaN        _RL    eta_qrN
82        _RL    etaNM1        _RL    eta_qrNM1
83        _RL    cgBeta        _RL    cgBeta
84        _RL    alpha        _RL    alpha
85        _RL    sumRHS        _RL    sumRHS
# Line 82  C     I, J, N     - Loop counters ( N co Line 96  C     I, J, N     - Loop counters ( N co
96    
97    
98  CcnhDebugStarts  CcnhDebugStarts
99        CHARACTER*(MAX_LEN_FNAM) suff  C     CHARACTER*(MAX_LEN_FNAM) suff
100  CcnhDebugEnds  CcnhDebugEnds
101    
102    
103  C--   Initialise inverter  C--   Initialise inverter
104        etaNBuf(1,myThid)   = 0. _d 0        eta_qrNM1 = 1. _d 0
       errBuf(1,myThid)    = 0. _d 0  
       sumRHSBuf(1,myThid) = 0. _d 0  
       etaNM1              = 1. _d 0  
105    
106  CcnhDebugStarts  CcnhDebugStarts
107  C     _EXCH_XY_R8( cg2d_b, myThid )  C     _EXCH_XY_R8( cg2d_b, myThid )
# Line 112  C--   Normalise RHS Line 123  C--   Normalise RHS
123          ENDDO          ENDDO
124         ENDDO         ENDDO
125        ENDDO        ENDDO
126        rhsMaxBuf(1,myThid) = rhsMax  #ifdef LETS_MAKE_JAM
127        _GLOBAL_MAX_R8( rhsMaxbuf, rhsMax, myThid )  C     _GLOBAL_MAX_R8( rhsMax, myThid )
128        rhsMax =  rhsMaxBuf(1,1)        rhsMax=1.
129    #else
130          _GLOBAL_MAX_R8( rhsMax, myThid )
131    Catm  rhsMax=1.
132    #endif
133        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
134        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
135        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 154  C--   Initial residual calculation Line 169  C--   Initial residual calculation
169       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
170       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
171       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
172       &    -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
173       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
174       &    )       &    )
175            err            = err            +            err            = err            +
# Line 166  C--   Initial residual calculation Line 181  C--   Initial residual calculation
181         ENDDO         ENDDO
182        ENDDO        ENDDO
183  C     _EXCH_XY_R8( cg2d_r, myThid )  C     _EXCH_XY_R8( cg2d_r, myThid )
184    #ifdef LETS_MAKE_JAM
185          CALL EXCH_XY_O1_R8_JAM( cg2d_r )
186    #else
187         OLw        = 1         OLw        = 1
188         OLe        = 1         OLe        = 1
189         OLn        = 1         OLn        = 1
# Line 177  C     _EXCH_XY_R8( cg2d_r, myThid ) Line 195  C     _EXCH_XY_R8( cg2d_r, myThid )
195       I            OLw, OLe, OLs, OLn, myNz,       I            OLw, OLe, OLs, OLn, myNz,
196       I            exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
197       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
198    #endif
199  C     _EXCH_XY_R8( cg2d_s, myThid )  C     _EXCH_XY_R8( cg2d_s, myThid )
200    #ifdef LETS_MAKE_JAM
201          CALL EXCH_XY_O1_R8_JAM( cg2d_s )
202    #else
203         OLw        = 1         OLw        = 1
204         OLe        = 1         OLe        = 1
205         OLn        = 1         OLn        = 1
# Line 189  C     _EXCH_XY_R8( cg2d_s, myThid ) Line 211  C     _EXCH_XY_R8( cg2d_s, myThid )
211       I            OLw, OLe, OLs, OLn, myNz,       I            OLw, OLe, OLs, OLn, myNz,
212       I            exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
213       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
214        sumRHSBuf(1,myThid) = sumRHS  #endif
215        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHS, myThid )
       sumRHS = sumRHSBuf(1,1)  
       errBuf(1,myThid)    = err  
216  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)
217        _GLOBAL_SUM_R8( errBuf    , err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
       err    = errBuf(1,1)  
218                
219        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
220        write(0,'(A,1PE30.14)') 'cg2d: Sum(rhs) = ',sumRHS        write(0,'(A,1PE30.14)') ' cg2d: Sum(rhs) = ',sumRHS
221        _END_MASTER( )        _END_MASTER( )
222    
223        actualIts      = 0        actualIts      = 0
224        actualResidual = SQRT(err)        actualResidual = SQRT(err)
225  C     _BARRIER  C     _BARRIER
226        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
227         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
228       & actualIts, actualResidual  c    & actualIts, actualResidual
229        _END_MASTER( )  c     _END_MASTER( )
230          initialResidual=actualResidual
231    
232  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
233        DO 10 it2d=1, cg2dMaxIters        DO 10 it2d=1, numIters
234    
235  CcnhDebugStarts  CcnhDebugStarts
236  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',
237  C    &  actualResidual  C    &  actualResidual
238  CcnhDebugEnds  CcnhDebugEnds
239         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. tolerance ) GOTO 11
240  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
241  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
242         etaN = 0. _d 0         eta_qrN = 0. _d 0
243         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
244          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
245           DO J=1,sNy           DO J=1,sNy
# Line 233  C--    conjugate direction vector "s". Line 253  C--    conjugate direction vector "s".
253  CcnhDebugStarts  CcnhDebugStarts
254  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)
255  CcnhDebugEnds  CcnhDebugEnds
256             etaN = etaN             eta_qrN = eta_qrN
257       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
258            ENDDO            ENDDO
259           ENDDO           ENDDO
260          ENDDO          ENDDO
261         ENDDO         ENDDO
262    
263         etanBuf(1,myThid) = etaN         _GLOBAL_SUM_R8(eta_qrN, myThid)
        _GLOBAL_SUM_R8(etaNbuf,etaN, myThid)  
        etaN = etaNBuf(1,1)  
264  CcnhDebugStarts  CcnhDebugStarts
265  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' etaN = ',etaN  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' eta_qrN = ',eta_qrN
266  CcnhDebugEnds  CcnhDebugEnds
267         cgBeta   = etaN/etaNM1         cgBeta   = eta_qrN/eta_qrNM1
268  CcnhDebugStarts  CcnhDebugStarts
269  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta
270  CcnhDebugEnds  CcnhDebugEnds
271         etaNM1 = etaN         eta_qrNM1 = eta_qrN
272    
273         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
274          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 266  CcnhDebugEnds Line 284  CcnhDebugEnds
284  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
285  C--    processes.  C--    processes.
286  C      _EXCH_XY_R8( cg2d_s, myThid )  C      _EXCH_XY_R8( cg2d_s, myThid )
287    #ifdef LETS_MAKE_JAM
288          CALL EXCH_XY_O1_R8_JAM( cg2d_s )
289    #else
290         OLw        = 1         OLw        = 1
291         OLe        = 1         OLe        = 1
292         OLn        = 1         OLn        = 1
# Line 277  C      _EXCH_XY_R8( cg2d_s, myThid ) Line 298  C      _EXCH_XY_R8( cg2d_s, myThid )
298       I            OLw, OLe, OLs, OLn, myNz,       I            OLw, OLe, OLs, OLn, myNz,
299       I            exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
300       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
301    #endif
302    
303  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
304  C==    q = A.s  C==    q = A.s
# Line 295  C==    q = A.s Line 316  C==    q = A.s
316       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
317       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
318       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
319       &    -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
320       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
321            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
322            ENDDO            ENDDO
323           ENDDO           ENDDO
324          ENDDO          ENDDO
325         ENDDO         ENDDO
326         alphaBuf(1,myThid) = alpha         _GLOBAL_SUM_R8(alpha,myThid)
        _GLOBAL_SUM_R8(alphaBuf,alpha,myThid)  
        alpha = alphaBuf(1,1)  
327  CcnhDebugStarts  CcnhDebugStarts
328  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha
329  CcnhDebugEnds  CcnhDebugEnds
330         alpha = etaN/alpha         alpha = eta_qrN/alpha
331  CcnhDebugStarts  CcnhDebugStarts
332  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
333  CcnhDebugEnds  CcnhDebugEnds
# Line 328  C      Now compute "interior" points. Line 347  C      Now compute "interior" points.
347          ENDDO          ENDDO
348         ENDDO         ENDDO
349    
350         errBuf(1,myThid) = err         _GLOBAL_SUM_R8( err   , myThid )
        _GLOBAL_SUM_R8( errBuf    , err   , myThid )  
        err = errBuf(1,1)  
351         err = SQRT(err)         err = SQRT(err)
352         actualIts      = it2d         actualIts      = it2d
353         actualResidual = err         actualResidual = err
354         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
355  C      _EXCH_XY_R8(cg2d_r, myThid )  C      _EXCH_XY_R8(cg2d_r, myThid )
356    #ifdef LETS_MAKE_JAM
357          CALL EXCH_XY_O1_R8_JAM( cg2d_r )
358    #else
359         OLw        = 1         OLw        = 1
360         OLe        = 1         OLe        = 1
361         OLn        = 1         OLn        = 1
# Line 347  C      _EXCH_XY_R8(cg2d_r, myThid ) Line 367  C      _EXCH_XY_R8(cg2d_r, myThid )
367       I             OLw, OLe, OLs, OLn, myNz,       I             OLw, OLe, OLs, OLn, myNz,
368       I             exchWidthX, exchWidthY,       I             exchWidthX, exchWidthY,
369       I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
370    #endif
371    
372     10 CONTINUE     10 CONTINUE
373     11 CONTINUE     11 CONTINUE
# Line 362  C--   Un-normalise the answer Line 383  C--   Un-normalise the answer
383         ENDDO         ENDDO
384        ENDDO        ENDDO
385    
386        _EXCH_XY_R8(cg2d_x, myThid )  C     The following exchange was moved up to solve_for_pressure
387        _BEGIN_MASTER( myThid )  C     for compatibility with TAMC.
388         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  C     _EXCH_XY_R8(cg2d_x, myThid )
389       & actualIts, actualResidual  c     _BEGIN_MASTER( myThid )
390        _END_MASTER( )  c      WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
391    c    & actualIts, actualResidual
392    c     _END_MASTER( )
393    
394    C--   Return parameters to caller
395          tolerance=initialResidual
396          residual=actualResidual
397          numIters=actualIts
398    
399  CcnhDebugStarts  CcnhDebugStarts
400  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )
# Line 390  C        ENDDO Line 418  C        ENDDO
418  C       ENDDO  C       ENDDO
419  C      ENDDO  C      ENDDO
420  C     ENDDO  C     ENDDO
421  C     errBuf(1,myThid)    = err  C     _GLOBAL_SUM_R8( err   , myThid )
 C     _GLOBAL_SUM_R8( errBuf    , err   , myThid )  
 C     err    = errBuf(1,1)  
422  C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)  C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)
423  CcnhDebugEnds  CcnhDebugEnds
424    
425          RETURN
426        END        END
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22