/[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.32 by heimbach, Mon May 14 21:33:29 2001 UTC revision 1.33 by adcroft, Tue May 29 14:01:36 2001 UTC
# Line 6  C $Name$ Line 6  C $Name$
6        SUBROUTINE CG2D(          SUBROUTINE CG2D(  
7       I                cg2d_b,       I                cg2d_b,
8       U                cg2d_x,       U                cg2d_x,
9       U                tolerance,       O                firstResidual,
10       O                residual,       O                lastResidual,
11       U                numIters,       U                numIters,
12       I                myThid )       I                myThid )
13  C     /==========================================================\  C     /==========================================================\
# Line 39  C     === Global data === Line 39  C     === Global data ===
39  #include "EEPARAMS.h"  #include "EEPARAMS.h"
40  #include "PARAMS.h"  #include "PARAMS.h"
41  #include "GRID.h"  #include "GRID.h"
42  #include "CG2D_INTERNAL.h"  #include "CG2D.h"
43  #include "SURFACE.h"  #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"  C     cg2d_b    - The source term or "right hand side"
48  C     cg2d_x    - The solution  C     cg2d_x    - The solution
49  C     tolerance - Entry: the tolerance of accuracy to solve to  C     firstResidual - the initial residual before any iterations
50  C                 Exit:  the initial residual before any iterations  C     lastResidual  - the actual residual reached
 C     residual  - Exit:  the actual residual reached  
51  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  - Entry: the maximum number of iterations allowed
52  C                 Exit:  the actual number of iterations used  C                 Exit:  the actual number of iterations used
53        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55        _RL  tolerance        _RL  firstResidual
56        _RL  residual        _RL  lastResidual
57        INTEGER numIters        INTEGER numIters
58        INTEGER myThid        INTEGER myThid
59    
# Line 74  C                   or they converge at Line 73  C                   or they converge at
73  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
74  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
75        INTEGER actualIts        INTEGER actualIts
76        _RL    actualResidual,initialResidual        _RL    actualResidual
77        INTEGER bi, bj                      INTEGER bi, bj              
78        INTEGER I, J, it2d        INTEGER I, J, it2d
79        _RL    err        _RL    err
# Line 123  C--   Normalise RHS Line 122  C--   Normalise RHS
122          ENDDO          ENDDO
123         ENDDO         ENDDO
124        ENDDO        ENDDO
125    
126          IF (cg2dNormaliseRHS) THEN
127    C-  Normalise RHS :
128  #ifdef LETS_MAKE_JAM  #ifdef LETS_MAKE_JAM
129  C     _GLOBAL_MAX_R8( rhsMax, myThid )  C     _GLOBAL_MAX_R8( rhsMax, myThid )
130        rhsMax=1.        rhsMax=1.
# Line 142  Catm  rhsMax=1. Line 144  Catm  rhsMax=1.
144          ENDDO          ENDDO
145         ENDDO         ENDDO
146        ENDDO        ENDDO
147    C- end Normalise RHS
148          ENDIF
149    
150  C--   Update overlaps  C--   Update overlaps
151        _EXCH_XY_R8( cg2d_b, myThid )        _EXCH_XY_R8( cg2d_b, myThid )
# Line 191  C     _EXCH_XY_R8( cg2d_r, myThid ) Line 195  C     _EXCH_XY_R8( cg2d_r, myThid )
195         exchWidthX = 1         exchWidthX = 1
196         exchWidthY = 1         exchWidthY = 1
197         myNz       = 1         myNz       = 1
198          IF (useCubedSphereExchange) THEN
199           CALL EXCH_RL_CUBE( cg2d_r,
200         I            OLw, OLe, OLs, OLn, myNz,
201         I            exchWidthX, exchWidthY,
202         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
203          ELSE
204         CALL EXCH_RL( cg2d_r,         CALL EXCH_RL( cg2d_r,
205       I            OLw, OLe, OLs, OLn, sNx, sNy, myNz,       I            OLw, OLe, OLs, OLn, myNz,
206       I            exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
207       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
208          ENDIF
209  #endif  #endif
210  C     _EXCH_XY_R8( cg2d_s, myThid )  C     _EXCH_XY_R8( cg2d_s, myThid )
211  #ifdef LETS_MAKE_JAM  #ifdef LETS_MAKE_JAM
# Line 207  C     _EXCH_XY_R8( cg2d_s, myThid ) Line 218  C     _EXCH_XY_R8( cg2d_s, myThid )
218         exchWidthX = 1         exchWidthX = 1
219         exchWidthY = 1         exchWidthY = 1
220         myNz       = 1         myNz       = 1
221          IF (useCubedSphereExchange) THEN
222           CALL EXCH_RL_CUBE( cg2d_s,
223         I            OLw, OLe, OLs, OLn, myNz,
224         I            exchWidthX, exchWidthY,
225         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
226          ELSE
227         CALL EXCH_RL( cg2d_s,         CALL EXCH_RL( cg2d_s,
228       I            OLw, OLe, OLs, OLn, sNx, sNy, myNz,       I            OLw, OLe, OLs, OLn, myNz,
229       I            exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
230       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
231          ENDIF
232  #endif  #endif
233        _GLOBAL_SUM_R8( sumRHS, myThid )         _GLOBAL_SUM_R8( sumRHS, myThid )
234  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)         _GLOBAL_SUM_R8( err   , myThid )
235        _GLOBAL_SUM_R8( err   , myThid )         err = SQRT(err)
236           actualIts      = 0
237           actualResidual = err
238                
239        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
240        write(*,'(A,1PE30.14)') ' cg2d: Sum(rhs) = ',sumRHS        write(*,'(A,1P2E22.14)')' cg2d: Sum(rhs),rhsMax = ',
241         &                                  sumRHS,rhsMax
242        _END_MASTER( )        _END_MASTER( )
   
       actualIts      = 0  
       actualResidual = SQRT(err)  
243  C     _BARRIER  C     _BARRIER
244  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
245  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
246  c    & actualIts, actualResidual  c    & actualIts, actualResidual
247  c     _END_MASTER( )  c     _END_MASTER( )
248        initialResidual=actualResidual        firstResidual=actualResidual
249    
250          IF ( err .LT. cg2dTolerance ) GOTO 11
251    
252  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
253        DO 10 it2d=1, numIters        DO 10 it2d=1, numIters
# Line 236  CcnhDebugStarts Line 256  CcnhDebugStarts
256  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',
257  C    &  actualResidual  C    &  actualResidual
258  CcnhDebugEnds  CcnhDebugEnds
        IF ( err .LT. tolerance ) GOTO 11  
259  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
260  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
261         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
# Line 294  C      _EXCH_XY_R8( cg2d_s, myThid ) Line 313  C      _EXCH_XY_R8( cg2d_s, myThid )
313         exchWidthX = 1         exchWidthX = 1
314         exchWidthY = 1         exchWidthY = 1
315         myNz       = 1         myNz       = 1
316          IF (useCubedSphereExchange) THEN
317           CALL EXCH_RL_CUBE( cg2d_s,
318         I            OLw, OLe, OLs, OLn, myNz,
319         I            exchWidthX, exchWidthY,
320         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
321          ELSE
322         CALL EXCH_RL( cg2d_s,         CALL EXCH_RL( cg2d_s,
323       I            OLw, OLe, OLs, OLn, sNx, sNy, myNz,       I            OLw, OLe, OLs, OLn, myNz,
324       I            exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
325       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
326          ENDIF
327  #endif  #endif
328    
329  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
# Line 351  C      Now compute "interior" points. Line 377  C      Now compute "interior" points.
377         err = SQRT(err)         err = SQRT(err)
378         actualIts      = it2d         actualIts      = it2d
379         actualResidual = err         actualResidual = err
380         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTolerance ) GOTO 11
381  C      _EXCH_XY_R8(cg2d_r, myThid )  C      _EXCH_XY_R8(cg2d_r, myThid )
382  #ifdef LETS_MAKE_JAM  #ifdef LETS_MAKE_JAM
383        CALL EXCH_XY_O1_R8_JAM( cg2d_r )        CALL EXCH_XY_O1_R8_JAM( cg2d_r )
# Line 363  C      _EXCH_XY_R8(cg2d_r, myThid ) Line 389  C      _EXCH_XY_R8(cg2d_r, myThid )
389         exchWidthX = 1         exchWidthX = 1
390         exchWidthY = 1         exchWidthY = 1
391         myNz       = 1         myNz       = 1
392          IF (useCubedSphereExchange) THEN
393           CALL EXCH_RL_CUBE( cg2d_r,
394         I            OLw, OLe, OLs, OLn, myNz,
395         I            exchWidthX, exchWidthY,
396         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
397          ELSE
398         CALL EXCH_RL( cg2d_r,         CALL EXCH_RL( cg2d_r,
399       I             OLw, OLe, OLs, OLn, sNx, sNy, myNz,       I            OLw, OLe, OLs, OLn, myNz,
400       I             exchWidthX, exchWidthY,       I            exchWidthX, exchWidthY,
401       I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )       I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
402          ENDIF
403  #endif  #endif
404    
405     10 CONTINUE     10 CONTINUE
406     11 CONTINUE     11 CONTINUE
407    
408          IF (cg2dNormaliseRHS) THEN
409  C--   Un-normalise the answer  C--   Un-normalise the answer
410        DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
411         DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
412          DO J=1,sNy            DO J=1,sNy
413           DO I=1,sNx             DO I=1,sNx
414            cg2d_x(I  ,J  ,bi,bj) = cg2d_x(I  ,J  ,bi,bj)/rhsNorm              cg2d_x(I  ,J  ,bi,bj) = cg2d_x(I  ,J  ,bi,bj)/rhsNorm
415               ENDDO
416              ENDDO
417           ENDDO           ENDDO
418          ENDDO          ENDDO
419         ENDDO        ENDIF
       ENDDO  
420    
421  C     The following exchange was moved up to solve_for_pressure  C     The following exchange was moved up to solve_for_pressure
422  C     for compatibility with TAMC.  C     for compatibility with TAMC.
# Line 392  c    & actualIts, actualResidual Line 427  c    & actualIts, actualResidual
427  c     _END_MASTER( )  c     _END_MASTER( )
428    
429  C--   Return parameters to caller  C--   Return parameters to caller
430        tolerance=initialResidual        lastResidual=actualResidual
       residual=actualResidual  
431        numIters=actualIts        numIters=actualIts
432    
433  CcnhDebugStarts  CcnhDebugStarts
Legend:
Removed from v.1.32  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22