/[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.10 by cnh, Sun Sep 6 17:35:19 1998 UTC revision 1.24 by adcroft, Fri Mar 24 18:31:31 2000 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5        SUBROUTINE CG2D(          SUBROUTINE CG2D(  
6         I                cg2d_b,
7         U                cg2d_x,
8       I                myThid )       I                myThid )
9  C     /==========================================================\  C     /==========================================================\
10  C     | SUBROUTINE CG2D                                          |  C     | SUBROUTINE CG2D                                          |
# Line 26  C     | solving the same problem. On the Line 28  C     | solving the same problem. On the
28  C     | the threads to solve several different problems          |  C     | the threads to solve several different problems          |
29  C     | concurrently this implementation will not work.          |  C     | concurrently this implementation will not work.          |
30  C     \==========================================================/  C     \==========================================================/
31          IMPLICIT NONE
32    
33  C     === Global data ===  C     === Global data ===
34  #include "SIZE.h"  #include "SIZE.h"
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37  #include "GRID.h"  #include "GRID.h"
38  #include "CG2D.h"  #include "CG2D_INTERNAL.h"
39    
40  C     === Routine arguments ===  C     === Routine arguments ===
41  C     myThid - Thread on which I am working.  C     myThid - Thread on which I am working.
42        INTEGER myThid        INTEGER myThid
43          _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44          _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45    
46    
47  C     === Local variables ====  C     === Local variables ====
48  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
# Line 54  C                   or they converge at Line 60  C                   or they converge at
60  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
61  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
62        INTEGER actualIts        INTEGER actualIts
63        REAL    actualResidual        _RL    actualResidual
64        INTEGER bi, bj                      INTEGER bi, bj              
65        INTEGER I, J, it2d        INTEGER I, J, it2d
66        REAL    err        _RL    err
67        REAL    etaN        _RL    etaN
68        REAL    etaNM1        _RL    etaNM1
69        REAL    cgBeta        _RL    cgBeta
70        REAL    alpha        _RL    alpha
71        REAL    sumRHS        _RL    sumRHS
72        REAL    rhsMax        _RL    rhsMax
73        REAL    rhsNorm        _RL    rhsNorm
74    
75          INTEGER OLw
76          INTEGER OLe
77          INTEGER OLn
78          INTEGER OLs
79          INTEGER exchWidthX
80          INTEGER exchWidthY
81          INTEGER myNz
82    
83    
84    CcnhDebugStarts
85    C     CHARACTER*(MAX_LEN_FNAM) suff
86    CcnhDebugEnds
87    
88    
89  C--   Initialise inverter  C--   Initialise inverter
90        etaNBuf(1,myThid)   = 0. D0        etaNM1              = 1. _d 0
       errBuf(1,myThid)    = 0. D0  
       sumRHSBuf(1,myThid) = 0. D0  
       etaNM1              = 1. D0  
91    
92  CcnhDebugStarts  CcnhDebugStarts
93        _EXCH_XY_R8( cg2d_b, myThid )  C     _EXCH_XY_R8( cg2d_b, myThid )
94        CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )
95    C     suff = 'unnormalised'
96    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
97    C     STOP
98  CcnhDebugEnds  CcnhDebugEnds
99    
100  C--   Normalise RHS  C--   Normalise RHS
# Line 89  C--   Normalise RHS Line 109  C--   Normalise RHS
109          ENDDO          ENDDO
110         ENDDO         ENDDO
111        ENDDO        ENDDO
112        rhsMaxBuf(1,myThid) = rhsMax  #ifdef LETS_MAKE_JAM
113        _GLOBAL_MAX_R8( rhsMaxbuf, rhsMax, myThid )        _GLOBAL_MAX_R8( rhsMax, myThid )
114        rhsMax =  rhsMaxBuf(1,1)  C     rhsMax=1.
115    #else
116          _GLOBAL_MAX_R8( rhsMax, myThid )
117    #endif
118        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
119        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
120        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 109  C--   Update overlaps Line 132  C--   Update overlaps
132        _EXCH_XY_R8( cg2d_b, myThid )        _EXCH_XY_R8( cg2d_b, myThid )
133        _EXCH_XY_R8( cg2d_x, myThid )        _EXCH_XY_R8( cg2d_x, myThid )
134  CcnhDebugStarts  CcnhDebugStarts
135        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 )
136    C     suff = 'normalised'
137    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
138  CcnhDebugEnds  CcnhDebugEnds
139    
140  C--   Initial residual calculation  C--   Initial residual calculation
# Line 140  C--   Initial residual calculation Line 165  C--   Initial residual calculation
165          ENDDO          ENDDO
166         ENDDO         ENDDO
167        ENDDO        ENDDO
168        _EXCH_XY_R8( cg2d_r, myThid )  C     _EXCH_XY_R8( cg2d_r, myThid )
169        _EXCH_XY_R8( cg2d_s, myThid )  #ifdef LETS_MAKE_JAM
170        sumRHSBuf(1,myThid) = sumRHS        CALL EXCH_XY_O1_R8_JAM( cg2d_r )
171        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )  #else
172        sumRHS = sumRHSBuf(1,1)         OLw        = 1
173        errBuf(1,myThid)    = err         OLe        = 1
174           OLn        = 1
175           OLs        = 1
176           exchWidthX = 1
177           exchWidthY = 1
178           myNz       = 1
179           CALL EXCH_RL( cg2d_r,
180         I            OLw, OLe, OLs, OLn, myNz,
181         I            exchWidthX, exchWidthY,
182         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
183    #endif
184    C     _EXCH_XY_R8( cg2d_s, myThid )
185    #ifdef LETS_MAKE_JAM
186          CALL EXCH_XY_O1_R8_JAM( cg2d_s )
187    #else
188           OLw        = 1
189           OLe        = 1
190           OLn        = 1
191           OLs        = 1
192           exchWidthX = 1
193           exchWidthY = 1
194           myNz       = 1
195           CALL EXCH_RL( cg2d_s,
196         I            OLw, OLe, OLs, OLn, myNz,
197         I            exchWidthX, exchWidthY,
198         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
199    #endif
200          _GLOBAL_SUM_R8( sumRHS, myThid )
201  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)
202        _GLOBAL_SUM_R8( errBuf    , err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
203        err    = errBuf(1,1)        
204        write(0,*) 'cg2d: Sum(rhs) = ',sumRHS        _BEGIN_MASTER( myThid )
205          write(0,'(A,1PE30.14)') ' cg2d: Sum(rhs) = ',sumRHS
206          _END_MASTER( )
207    
208        actualIts      = 0        actualIts      = 0
209        actualResidual = SQRT(err)        actualResidual = SQRT(err)
210  C     _BARRIER  C     _BARRIER
211        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
212         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
213         & actualIts, actualResidual
214        _END_MASTER( )        _END_MASTER( )
215    
216  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
217        DO 10 it2d=1, cg2dMaxIters        DO 10 it2d=1, cg2dMaxIters
218    
219  CcnhDebugStarts  CcnhDebugStarts
220  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',actualResidual  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',
221    C    &  actualResidual
222  CcnhDebugEnds  CcnhDebugEnds
223         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
224  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
# Line 188  CcnhDebugEnds Line 244  CcnhDebugEnds
244          ENDDO          ENDDO
245         ENDDO         ENDDO
246    
247         etanBuf(1,myThid) = etaN         _GLOBAL_SUM_R8(etaN, myThid)
        _GLOBAL_SUM_R8(etaNbuf,etaN, myThid)  
        etaN = etaNBuf(1,1)  
248  CcnhDebugStarts  CcnhDebugStarts
249  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' etaN = ',etaN  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' etaN = ',etaN
250  CcnhDebugEnds  CcnhDebugEnds
# Line 204  CcnhDebugEnds Line 258  CcnhDebugEnds
258          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
259           DO J=1,sNy           DO J=1,sNy
260            DO I=1,sNx            DO I=1,sNx
261             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)
262         &                       + cgBeta*cg2d_s(I,J,bi,bj)
263            ENDDO            ENDDO
264           ENDDO           ENDDO
265          ENDDO          ENDDO
# Line 212  CcnhDebugEnds Line 267  CcnhDebugEnds
267    
268  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
269  C--    processes.  C--    processes.
270         _EXCH_XY_R8( cg2d_s, myThid )  C      _EXCH_XY_R8( cg2d_s, myThid )
271    #ifdef LETS_MAKE_JAM
272          CALL EXCH_XY_O1_R8_JAM( cg2d_s )
273    #else
274           OLw        = 1
275           OLe        = 1
276           OLn        = 1
277           OLs        = 1
278           exchWidthX = 1
279           exchWidthY = 1
280           myNz       = 1
281           CALL EXCH_RL( cg2d_s,
282         I            OLw, OLe, OLs, OLn, myNz,
283         I            exchWidthX, exchWidthY,
284         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
285    #endif
286    
287  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
288  C==    q = A.s  C==    q = A.s
# Line 237  C==    q = A.s Line 307  C==    q = A.s
307           ENDDO           ENDDO
308          ENDDO          ENDDO
309         ENDDO         ENDDO
310         alphaBuf(1,myThid) = alpha         _GLOBAL_SUM_R8(alpha,myThid)
        _GLOBAL_SUM_R8(alphaBuf,alpha,myThid)  
        alpha = alphaBuf(1,1)  
311  CcnhDebugStarts  CcnhDebugStarts
312  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha
313  CcnhDebugEnds  CcnhDebugEnds
# Line 263  C      Now compute "interior" points. Line 331  C      Now compute "interior" points.
331          ENDDO          ENDDO
332         ENDDO         ENDDO
333    
334         errBuf(1,myThid) = err         _GLOBAL_SUM_R8( err   , myThid )
        _GLOBAL_SUM_R8( errBuf    , err   , myThid )  
        err = errBuf(1,1)  
335         err = SQRT(err)         err = SQRT(err)
336         actualIts      = it2d         actualIts      = it2d
337         actualResidual = err         actualResidual = err
338         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
339         _EXCH_XY_R8(cg2d_r, myThid )  C      _EXCH_XY_R8(cg2d_r, myThid )
340    #ifdef LETS_MAKE_JAM
341          CALL EXCH_XY_O1_R8_JAM( cg2d_r )
342    #else
343           OLw        = 1
344           OLe        = 1
345           OLn        = 1
346           OLs        = 1
347           exchWidthX = 1
348           exchWidthY = 1
349           myNz       = 1
350           CALL EXCH_RL( cg2d_r,
351         I             OLw, OLe, OLs, OLn, myNz,
352         I             exchWidthX, exchWidthY,
353         I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
354    #endif
355    
356     10 CONTINUE     10 CONTINUE
357     11 CONTINUE     11 CONTINUE
358    
# Line 285  C--   Un-normalise the answer Line 367  C--   Un-normalise the answer
367         ENDDO         ENDDO
368        ENDDO        ENDDO
369    
370        _EXCH_XY_R8(cg2d_x, myThid )  C     The following exchange was moved up to solve_for_pressure
371    C     for compatibility with TAMC.
372    C     _EXCH_XY_R8(cg2d_x, myThid )
373        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
374         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
375         & actualIts, actualResidual
376        _END_MASTER( )        _END_MASTER( )
377    
378  CcnhDebugStarts  CcnhDebugStarts
# Line 312  C        ENDDO Line 397  C        ENDDO
397  C       ENDDO  C       ENDDO
398  C      ENDDO  C      ENDDO
399  C     ENDDO  C     ENDDO
400  C     errBuf(1,myThid)    = err  C     _GLOBAL_SUM_R8( err   , myThid )
 C     _GLOBAL_SUM_R8( errBuf    , err   , myThid )  
 C     err    = errBuf(1,1)  
401  C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)  C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)
402  CcnhDebugEnds  CcnhDebugEnds
403    
404          RETURN
405        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22