/[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.26 by adcroft, Fri Feb 2 21:04:47 2001 UTC revision 1.30 by adcroft, Fri Mar 9 20:45:09 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6        SUBROUTINE CG2D(          SUBROUTINE CG2D(  
7       I                cg2d_b,       I                cg2d_b,
8       U                cg2d_x,       U                cg2d_x,
9         U                tolerance,
10         O                residual,
11         U                numIters,
12       I                myThid )       I                myThid )
13  C     /==========================================================\  C     /==========================================================\
14  C     | SUBROUTINE CG2D                                          |  C     | SUBROUTINE CG2D                                          |
# Line 36  C     === Global data === Line 40  C     === Global data ===
40  #include "PARAMS.h"  #include "PARAMS.h"
41  #include "GRID.h"  #include "GRID.h"
42  #include "CG2D_INTERNAL.h"  #include "CG2D_INTERNAL.h"
43    #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        INTEGER myThid  C     cg2d_b    - The source term or "right hand side"
48        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     cg2d_x    - The solution
49    C     tolerance - Entry: the tolerance of accuracy to solve to
50    C                 Exit:  the initial residual before any iterations
51    C     residual  - Exit:  the actual residual reached
52    C     numIters  - Entry: the maximum number of iterations allowed
53    C                 Exit:  the actual number of iterations used
54        _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)
55          _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56          _RL  tolerance
57          _RL  residual
58          INTEGER numIters
59          INTEGER myThid
60    
61  C     === Local variables ====  C     === Local variables ====
62  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
63  C     actualResidual - residual  C     actualResidual - residual
64  C     bi          - Block index in X and Y.  C     bi          - Block index in X and Y.
65  C     bj  C     bj
66  C     etaN        - Used in computing search directions  C     eta_qrN     - Used in computing search directions
67  C     etaNM1        suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
68  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
69  C     alpha    C     alpha  
70  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS      - Sum of right-hand-side. Sometimes this is a
# Line 60  C                   or they converge at Line 74  C                   or they converge at
74  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
75  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
76        INTEGER actualIts        INTEGER actualIts
77        _RL    actualResidual        _RL    actualResidual,initialResidual
78        INTEGER bi, bj                      INTEGER bi, bj              
79        INTEGER I, J, it2d        INTEGER I, J, it2d
80        _RL    err        _RL    err
81        _RL    etaN        _RL    eta_qrN
82        _RL    etaNM1        _RL    eta_qrNM1
83        _RL    cgBeta        _RL    cgBeta
84        _RL    alpha        _RL    alpha
85        _RL    sumRHS        _RL    sumRHS
# Line 87  CcnhDebugEnds Line 101  CcnhDebugEnds
101    
102    
103  C--   Initialise inverter  C--   Initialise inverter
104        etaNM1              = 1. _d 0        eta_qrNM1 = 1. _d 0
105    
106  CcnhDebugStarts  CcnhDebugStarts
107  C     _EXCH_XY_R8( cg2d_b, myThid )  C     _EXCH_XY_R8( cg2d_b, myThid )
# Line 155  C--   Initial residual calculation Line 169  C--   Initial residual calculation
169       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
170       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
171       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
172       &    -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
173       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
174       &    )       &    )
175            err            = err            +            err            = err            +
# Line 209  C     WRITE(6,*) ' mythid, err = ', myth Line 223  C     WRITE(6,*) ' mythid, err = ', myth
223        actualIts      = 0        actualIts      = 0
224        actualResidual = SQRT(err)        actualResidual = SQRT(err)
225  C     _BARRIER  C     _BARRIER
226        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
227         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
228       & actualIts, actualResidual  c    & actualIts, actualResidual
229        _END_MASTER( )  c     _END_MASTER( )
230          initialResidual=actualResidual
231    
232  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
233        DO 10 it2d=1, cg2dMaxIters        DO 10 it2d=1, numIters
234    
235  CcnhDebugStarts  CcnhDebugStarts
236  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',
237  C    &  actualResidual  C    &  actualResidual
238  CcnhDebugEnds  CcnhDebugEnds
239         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. tolerance ) GOTO 11
240  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
241  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
242         etaN = 0. _d 0         eta_qrN = 0. _d 0
243         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
244          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
245           DO J=1,sNy           DO J=1,sNy
# Line 238  C--    conjugate direction vector "s". Line 253  C--    conjugate direction vector "s".
253  CcnhDebugStarts  CcnhDebugStarts
254  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)
255  CcnhDebugEnds  CcnhDebugEnds
256             etaN = etaN             eta_qrN = eta_qrN
257       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
258            ENDDO            ENDDO
259           ENDDO           ENDDO
260          ENDDO          ENDDO
261         ENDDO         ENDDO
262    
263         _GLOBAL_SUM_R8(etaN, myThid)         _GLOBAL_SUM_R8(eta_qrN, myThid)
264  CcnhDebugStarts  CcnhDebugStarts
265  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' etaN = ',etaN  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' eta_qrN = ',eta_qrN
266  CcnhDebugEnds  CcnhDebugEnds
267         cgBeta   = etaN/etaNM1         cgBeta   = eta_qrN/eta_qrNM1
268  CcnhDebugStarts  CcnhDebugStarts
269  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta
270  CcnhDebugEnds  CcnhDebugEnds
271         etaNM1 = etaN         eta_qrNM1 = eta_qrN
272    
273         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
274          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 301  C==    q = A.s Line 316  C==    q = A.s
316       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
317       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
318       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
319       &    -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
320       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
321            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
322            ENDDO            ENDDO
# Line 312  C==    q = A.s Line 327  C==    q = A.s
327  CcnhDebugStarts  CcnhDebugStarts
328  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha
329  CcnhDebugEnds  CcnhDebugEnds
330         alpha = etaN/alpha         alpha = eta_qrN/alpha
331  CcnhDebugStarts  CcnhDebugStarts
332  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
333  CcnhDebugEnds  CcnhDebugEnds
# Line 371  C--   Un-normalise the answer Line 386  C--   Un-normalise the answer
386  C     The following exchange was moved up to solve_for_pressure  C     The following exchange was moved up to solve_for_pressure
387  C     for compatibility with TAMC.  C     for compatibility with TAMC.
388  C     _EXCH_XY_R8(cg2d_x, myThid )  C     _EXCH_XY_R8(cg2d_x, myThid )
389        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
390         WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
391       & actualIts, actualResidual  c    & actualIts, actualResidual
392        _END_MASTER( )  c     _END_MASTER( )
393    
394    C--   Return parameters to caller
395          tolerance=initialResidual
396          residual=actualResidual
397          numIters=actualIts
398    
399  CcnhDebugStarts  CcnhDebugStarts
400  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22