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