/[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.14 by cnh, Wed Oct 28 03:11:36 1998 UTC
# Line 3  C $Header$ Line 3  C $Header$
3  #include "CPP_EEOPTIONS.h"  #include "CPP_EEOPTIONS.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 32  C     === Global data === Line 34  C     === Global data ===
34  #include "EEPARAMS.h"  #include "EEPARAMS.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36  #include "GRID.h"  #include "GRID.h"
37  #include "CG2D.h"  #include "CG2D_INTERNAL.h"
38    
39  C     === Routine arguments ===  C     === Routine arguments ===
40  C     myThid - Thread on which I am working.  C     myThid - Thread on which I am working.
41        INTEGER myThid        INTEGER myThid
42          _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
43          _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44    
45    
46  C     === Local variables ====  C     === Local variables ====
47  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
# Line 54  C                   or they converge at Line 59  C                   or they converge at
59  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
60  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
61        INTEGER actualIts        INTEGER actualIts
62        REAL    actualResidual        _RL    actualResidual
63        INTEGER bi, bj                      INTEGER bi, bj              
64        INTEGER I, J, it2d        INTEGER I, J, it2d
65        REAL    err        _RL    err
66        REAL    etaN        _RL    etaN
67        REAL    etaNM1        _RL    etaNM1
68        REAL    cgBeta        _RL    cgBeta
69        REAL    alpha        _RL    alpha
70        REAL    sumRHS        _RL    sumRHS
71        REAL    rhsMax        _RL    rhsMax
72        REAL    rhsNorm        _RL    rhsNorm
73    
74          INTEGER OLw
75          INTEGER OLe
76          INTEGER OLn
77          INTEGER OLs
78          INTEGER exchWidthX
79          INTEGER exchWidthY
80          INTEGER myNz
81    
82    
83    CcnhDebugStarts
84          CHARACTER*(MAX_LEN_FNAM) suff
85    CcnhDebugEnds
86    
87    
88  C--   Initialise inverter  C--   Initialise inverter
89        etaNBuf(1,myThid)   = 0. D0        etaNBuf(1,myThid)   = 0. D0
# Line 73  C--   Initialise inverter Line 92  C--   Initialise inverter
92        etaNM1              = 1. D0        etaNM1              = 1. D0
93    
94  CcnhDebugStarts  CcnhDebugStarts
95        _EXCH_XY_R8( cg2d_b, myThid )  C     _EXCH_XY_R8( cg2d_b, myThid )
96        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 )
97    C     suff = 'unnormalised'
98    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
99    C     STOP
100  CcnhDebugEnds  CcnhDebugEnds
101    
102  C--   Normalise RHS  C--   Normalise RHS
# Line 109  C--   Update overlaps Line 131  C--   Update overlaps
131        _EXCH_XY_R8( cg2d_b, myThid )        _EXCH_XY_R8( cg2d_b, myThid )
132        _EXCH_XY_R8( cg2d_x, myThid )        _EXCH_XY_R8( cg2d_x, myThid )
133  CcnhDebugStarts  CcnhDebugStarts
134        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 )
135    C     suff = 'normalised'
136    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
137  CcnhDebugEnds  CcnhDebugEnds
138    
139  C--   Initial residual calculation  C--   Initial residual calculation
# Line 140  C--   Initial residual calculation Line 164  C--   Initial residual calculation
164          ENDDO          ENDDO
165         ENDDO         ENDDO
166        ENDDO        ENDDO
167        _EXCH_XY_R8( cg2d_r, myThid )  C     _EXCH_XY_R8( cg2d_r, myThid )
168        _EXCH_XY_R8( cg2d_s, myThid )         OLw        = 1
169           OLe        = 1
170           OLn        = 1
171           OLs        = 1
172           exchWidthX = 1
173           exchWidthY = 1
174           myNz       = 1
175           CALL EXCH_RL( cg2d_r,
176         I            OLw, OLe, OLs, OLn, myNz,
177         I            exchWidthX, exchWidthY,
178         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
179    C     _EXCH_XY_R8( cg2d_s, myThid )
180           OLw        = 1
181           OLe        = 1
182           OLn        = 1
183           OLs        = 1
184           exchWidthX = 1
185           exchWidthY = 1
186           myNz       = 1
187           CALL EXCH_RL( cg2d_s,
188         I            OLw, OLe, OLs, OLn, myNz,
189         I            exchWidthX, exchWidthY,
190         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
191        sumRHSBuf(1,myThid) = sumRHS        sumRHSBuf(1,myThid) = sumRHS
192        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )
193        sumRHS = sumRHSBuf(1,1)        sumRHS = sumRHSBuf(1,1)
# Line 149  C--   Initial residual calculation Line 195  C--   Initial residual calculation
195  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)
196        _GLOBAL_SUM_R8( errBuf    , err   , myThid )        _GLOBAL_SUM_R8( errBuf    , err   , myThid )
197        err    = errBuf(1,1)        err    = errBuf(1,1)
198        write(0,*) 'cg2d: Sum(rhs) = ',sumRHS        
199          _BEGIN_MASTER( myThid )
200          write(0,'(A,1PE30.20)') 'cg2d: Sum(rhs) = ',sumRHS
201          _END_MASTER( )
202    
203        actualIts      = 0        actualIts      = 0
204        actualResidual = SQRT(err)        actualResidual = SQRT(err)
205  C     _BARRIER  C     _BARRIER
206        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
207         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(0,'(A,I6,1PE30.20)') ' CG2D iters, err = ',
208         & actualIts, actualResidual
209        _END_MASTER( )        _END_MASTER( )
210    
211  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
212        DO 10 it2d=1, cg2dMaxIters        DO 10 it2d=1, cg2dMaxIters
213    
214  CcnhDebugStarts  CcnhDebugStarts
215  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',actualResidual  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',
216    C    &  actualResidual
217  CcnhDebugEnds  CcnhDebugEnds
218         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
219  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
# Line 204  CcnhDebugEnds Line 255  CcnhDebugEnds
255          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
256           DO J=1,sNy           DO J=1,sNy
257            DO I=1,sNx            DO I=1,sNx
258             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)
259         &                       + cgBeta*cg2d_s(I,J,bi,bj)
260            ENDDO            ENDDO
261           ENDDO           ENDDO
262          ENDDO          ENDDO
# Line 212  CcnhDebugEnds Line 264  CcnhDebugEnds
264    
265  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
266  C--    processes.  C--    processes.
267         _EXCH_XY_R8( cg2d_s, myThid )  C      _EXCH_XY_R8( cg2d_s, myThid )
268           OLw        = 1
269           OLe        = 1
270           OLn        = 1
271           OLs        = 1
272           exchWidthX = 1
273           exchWidthY = 1
274           myNz       = 1
275           CALL EXCH_RL( cg2d_s,
276         I            OLw, OLe, OLs, OLn, myNz,
277         I            exchWidthX, exchWidthY,
278         I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
279    
280    
281  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
282  C==    q = A.s  C==    q = A.s
# Line 270  C      Now compute "interior" points. Line 334  C      Now compute "interior" points.
334         actualIts      = it2d         actualIts      = it2d
335         actualResidual = err         actualResidual = err
336         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
337         _EXCH_XY_R8(cg2d_r, myThid )  C      _EXCH_XY_R8(cg2d_r, myThid )
338           OLw        = 1
339           OLe        = 1
340           OLn        = 1
341           OLs        = 1
342           exchWidthX = 1
343           exchWidthY = 1
344           myNz       = 1
345           CALL EXCH_RL( cg2d_r,
346         I             OLw, OLe, OLs, OLn, myNz,
347         I             exchWidthX, exchWidthY,
348         I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
349    
350     10 CONTINUE     10 CONTINUE
351     11 CONTINUE     11 CONTINUE
352    
# Line 287  C--   Un-normalise the answer Line 363  C--   Un-normalise the answer
363    
364        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
365        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
366         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(0,'(A,I6,1PE30.20)') ' CG2D iters, err = ',
367         & actualIts, actualResidual
368        _END_MASTER( )        _END_MASTER( )
369    
370  CcnhDebugStarts  CcnhDebugStarts
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22