/[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.19 by adcroft, Mon Mar 22 15:54:03 1999 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          CHARACTER*(MAX_LEN_FNAM) suff
86    CcnhDebugEnds
87    
88    
89  C--   Initialise inverter  C--   Initialise inverter
90        etaNBuf(1,myThid)   = 0. D0        etaNBuf(1,myThid)   = 0. _d 0
91        errBuf(1,myThid)    = 0. D0        errBuf(1,myThid)    = 0. _d 0
92        sumRHSBuf(1,myThid) = 0. D0        sumRHSBuf(1,myThid) = 0. _d 0
93        etaNM1              = 1. D0        etaNM1              = 1. _d 0
94    
95  CcnhDebugStarts  CcnhDebugStarts
96        _EXCH_XY_R8( cg2d_b, myThid )  C     _EXCH_XY_R8( cg2d_b, myThid )
97        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 )
98    C     suff = 'unnormalised'
99    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
100    C     STOP
101  CcnhDebugEnds  CcnhDebugEnds
102    
103  C--   Normalise RHS  C--   Normalise RHS
# 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 )         OLw        = 1
170           OLe        = 1
171           OLn        = 1
172           OLs        = 1
173           exchWidthX = 1
174           exchWidthY = 1
175           myNz       = 1
176           CALL EXCH_RL( cg2d_r,
177         I            OLw, OLe, OLs, OLn, myNz,
178         I            exchWidthX, exchWidthY,
179         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
180    C     _EXCH_XY_R8( cg2d_s, myThid )
181           OLw        = 1
182           OLe        = 1
183           OLn        = 1
184           OLs        = 1
185           exchWidthX = 1
186           exchWidthY = 1
187           myNz       = 1
188           CALL EXCH_RL( cg2d_s,
189         I            OLw, OLe, OLs, OLn, myNz,
190         I            exchWidthX, exchWidthY,
191         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
192        sumRHSBuf(1,myThid) = sumRHS        sumRHSBuf(1,myThid) = sumRHS
193        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )
194        sumRHS = sumRHSBuf(1,1)        sumRHS = sumRHSBuf(1,1)
# Line 149  C--   Initial residual calculation Line 196  C--   Initial residual calculation
196  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)
197        _GLOBAL_SUM_R8( errBuf    , err   , myThid )        _GLOBAL_SUM_R8( errBuf    , err   , myThid )
198        err    = errBuf(1,1)        err    = errBuf(1,1)
199        write(0,*) 'cg2d: Sum(rhs) = ',sumRHS        
200          _BEGIN_MASTER( myThid )
201          write(0,'(A,1PE30.14)') ' cg2d: Sum(rhs) = ',sumRHS
202          _END_MASTER( )
203    
204        actualIts      = 0        actualIts      = 0
205        actualResidual = SQRT(err)        actualResidual = SQRT(err)
206  C     _BARRIER  C     _BARRIER
207        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
208         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
209         & actualIts, actualResidual
210        _END_MASTER( )        _END_MASTER( )
211    
212  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
213        DO 10 it2d=1, cg2dMaxIters        DO 10 it2d=1, cg2dMaxIters
214    
215  CcnhDebugStarts  CcnhDebugStarts
216  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',actualResidual  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',
217    C    &  actualResidual
218  CcnhDebugEnds  CcnhDebugEnds
219         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
220  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
# Line 204  CcnhDebugEnds Line 256  CcnhDebugEnds
256          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
257           DO J=1,sNy           DO J=1,sNy
258            DO I=1,sNx            DO I=1,sNx
259             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)
260         &                       + cgBeta*cg2d_s(I,J,bi,bj)
261            ENDDO            ENDDO
262           ENDDO           ENDDO
263          ENDDO          ENDDO
# Line 212  CcnhDebugEnds Line 265  CcnhDebugEnds
265    
266  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
267  C--    processes.  C--    processes.
268         _EXCH_XY_R8( cg2d_s, myThid )  C      _EXCH_XY_R8( cg2d_s, myThid )
269           OLw        = 1
270           OLe        = 1
271           OLn        = 1
272           OLs        = 1
273           exchWidthX = 1
274           exchWidthY = 1
275           myNz       = 1
276           CALL EXCH_RL( cg2d_s,
277         I            OLw, OLe, OLs, OLn, myNz,
278         I            exchWidthX, exchWidthY,
279         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
280    
281    
282  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
283  C==    q = A.s  C==    q = A.s
# Line 270  C      Now compute "interior" points. Line 335  C      Now compute "interior" points.
335         actualIts      = it2d         actualIts      = it2d
336         actualResidual = err         actualResidual = err
337         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
338         _EXCH_XY_R8(cg2d_r, myThid )  C      _EXCH_XY_R8(cg2d_r, myThid )
339           OLw        = 1
340           OLe        = 1
341           OLn        = 1
342           OLs        = 1
343           exchWidthX = 1
344           exchWidthY = 1
345           myNz       = 1
346           CALL EXCH_RL( cg2d_r,
347         I             OLw, OLe, OLs, OLn, myNz,
348         I             exchWidthX, exchWidthY,
349         I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
350    
351     10 CONTINUE     10 CONTINUE
352     11 CONTINUE     11 CONTINUE
353    
# Line 287  C--   Un-normalise the answer Line 364  C--   Un-normalise the answer
364    
365        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
366        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
367         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
368         & actualIts, actualResidual
369        _END_MASTER( )        _END_MASTER( )
370    
371  CcnhDebugStarts  CcnhDebugStarts
# Line 318  C     err    = errBuf(1,1) Line 396  C     err    = errBuf(1,1)
396  C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)  C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)
397  CcnhDebugEnds  CcnhDebugEnds
398    
399          RETURN
400        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22