/[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.7 by cnh, Mon Jun 15 05:13:55 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_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6        SUBROUTINE CG2D(          SUBROUTINE CG2D(  
7         I                cg2d_b,
8         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 26  C     | solving the same problem. On the Line 32  C     | solving the same problem. On the
32  C     | the threads to solve several different problems          |  C     | the threads to solve several different problems          |
33  C     | concurrently this implementation will not work.          |  C     | concurrently this implementation will not work.          |
34  C     \==========================================================/  C     \==========================================================/
35          IMPLICIT NONE
36    
37  C     === Global data ===  C     === Global data ===
38  #include "SIZE.h"  #include "SIZE.h"
39  #include "EEPARAMS.h"  #include "EEPARAMS.h"
40  #include "PARAMS.h"  #include "PARAMS.h"
41  #include "GRID.h"  #include "GRID.h"
42  #include "CG2D.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    C     cg2d_b    - The source term or "right hand side"
48    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)
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        INTEGER myThid
60    
61  C     === Local variables ====  C     === Local variables ====
# Line 43  C     actualIts      - Number of iterati Line 63  C     actualIts      - Number of iterati
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 54  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        REAL    actualResidual        _RL    actualResidual,initialResidual
78        INTEGER bi, bj                      INTEGER bi, bj              
79        INTEGER I, J, it2d        INTEGER I, J, it2d
80        REAL    err        _RL    err
81        REAL    etaN        _RL    eta_qrN
82        REAL    etaNM1        _RL    eta_qrNM1
83        REAL    cgBeta        _RL    cgBeta
84        REAL    alpha        _RL    alpha
85        REAL    sumRHS        _RL    sumRHS
86        REAL    rhsMax        _RL    rhsMax
87        REAL    rhsNorm        _RL    rhsNorm
88    
89          INTEGER OLw
90          INTEGER OLe
91          INTEGER OLn
92          INTEGER OLs
93          INTEGER exchWidthX
94          INTEGER exchWidthY
95          INTEGER myNz
96    
97    
98    CcnhDebugStarts
99    C     CHARACTER*(MAX_LEN_FNAM) suff
100    CcnhDebugEnds
101    
102    
103  C--   Initialise inverter  C--   Initialise inverter
104        etaNBuf(1,myThid)   = 0. D0        eta_qrNM1 = 1. _d 0
105        errBuf(1,myThid)    = 0. D0  
106        sumRHSBuf(1,myThid) = 0. D0  CcnhDebugStarts
107        etaNM1              = 1. D0  C     _EXCH_XY_R8( cg2d_b, myThid )
108    C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )
109    C     suff = 'unnormalised'
110    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
111    C     STOP
112    CcnhDebugEnds
113    
114  C--   Normalise RHS  C--   Normalise RHS
115        rhsMax = 0. _d 0        rhsMax = 0. _d 0
# Line 84  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 105  C--   Update overlaps Line 148  C--   Update overlaps
148        _EXCH_XY_R8( cg2d_x, myThid )        _EXCH_XY_R8( cg2d_x, myThid )
149  CcnhDebugStarts  CcnhDebugStarts
150  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )
151    C     suff = 'normalised'
152    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
153  CcnhDebugEnds  CcnhDebugEnds
154    
155  C--   Initial residual calculation  C--   Initial residual calculation
# Line 124  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*_zA(i,j,bi,bj)*       &    -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 135  C--   Initial residual calculation Line 180  C--   Initial residual calculation
180          ENDDO          ENDDO
181         ENDDO         ENDDO
182        ENDDO        ENDDO
183        _EXCH_XY_R8( cg2d_r, myThid )  C     _EXCH_XY_R8( cg2d_r, myThid )
184        _EXCH_XY_R8( cg2d_s, myThid )  #ifdef LETS_MAKE_JAM
185        sumRHSBuf(1,myThid) = sumRHS        CALL EXCH_XY_O1_R8_JAM( cg2d_r )
186        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )  #else
187        sumRHS = sumRHSBuf(1,1)         OLw        = 1
188        errBuf(1,myThid)    = err         OLe        = 1
189           OLn        = 1
190           OLs        = 1
191           exchWidthX = 1
192           exchWidthY = 1
193           myNz       = 1
194           CALL EXCH_RL( cg2d_r,
195         I            OLw, OLe, OLs, OLn, myNz,
196         I            exchWidthX, exchWidthY,
197         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
198    #endif
199    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
204           OLe        = 1
205           OLn        = 1
206           OLs        = 1
207           exchWidthX = 1
208           exchWidthY = 1
209           myNz       = 1
210           CALL EXCH_RL( cg2d_s,
211         I            OLw, OLe, OLs, OLn, myNz,
212         I            exchWidthX, exchWidthY,
213         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
214    #endif
215          _GLOBAL_SUM_R8( sumRHS, myThid )
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 )
218        err    = errBuf(1,1)        
219  C     write(0,*) 'cg2d: Sum(rhs) = ',sumRHS        _BEGIN_MASTER( myThid )
220          write(0,'(A,1PE30.14)') ' cg2d: Sum(rhs) = ',sumRHS
221          _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,*) ' CG2D iters, err = ', actualIts, actualResidual  c      WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
228        _END_MASTER( )  c    & actualIts, actualResidual
229    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 = ',actualResidual  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',
237    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 176  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)
275           DO J=1,sNy           DO J=1,sNy
276            DO I=1,sNx            DO I=1,sNx
277             cg2d_s(I,J,bi,bj) = cg2d_q(I,J,bi,bj) + cgBeta*cg2d_s(I,J,bi,bj)             cg2d_s(I,J,bi,bj) = cg2d_q(I,J,bi,bj)
278         &                       + cgBeta*cg2d_s(I,J,bi,bj)
279            ENDDO            ENDDO
280           ENDDO           ENDDO
281          ENDDO          ENDDO
# Line 207  CcnhDebugEnds Line 283  CcnhDebugEnds
283    
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         _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
291           OLe        = 1
292           OLn        = 1
293           OLs        = 1
294           exchWidthX = 1
295           exchWidthY = 1
296           myNz       = 1
297           CALL EXCH_RL( cg2d_s,
298         I            OLw, OLe, OLs, OLn, myNz,
299         I            exchWidthX, exchWidthY,
300         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 225  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*_zA(i,j,bi,bj)*       &    -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 258  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         _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
360           OLe        = 1
361           OLn        = 1
362           OLs        = 1
363           exchWidthX = 1
364           exchWidthY = 1
365           myNz       = 1
366           CALL EXCH_RL( cg2d_r,
367         I             OLw, OLe, OLs, OLn, myNz,
368         I             exchWidthX, exchWidthY,
369         I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
370    #endif
371    
372     10 CONTINUE     10 CONTINUE
373     11 CONTINUE     11 CONTINUE
374    
# Line 280  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,*) ' CG2D iters, err = ', actualIts, actualResidual  C     _EXCH_XY_R8(cg2d_x, myThid )
389        _END_MASTER( )  c     _BEGIN_MASTER( myThid )
390    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 307  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.7  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22