/[MITgcm]/MITgcm/model/src/cg2d_nsa.F
ViewVC logotype

Diff of /MITgcm/model/src/cg2d_nsa.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by heimbach, Wed Jun 7 01:45:42 2006 UTC revision 1.5 by jmc, Tue Jul 17 21:31:52 2012 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
 #ifdef ALLOW_USE_MPI  
 C HACK to avoid global_max  
 # define ALLOW_CONST_RHSMAX  
 #endif  
5    
6  CML THIS DOES NOT WORK +++++  CML THIS DOES NOT WORK +++++
7  #undef ALLOW_LOOP_DIRECTIVE  #undef ALLOW_LOOP_DIRECTIVE
8  CBOP  CBOP
9  C     !ROUTINE: CG2D_NSA  C     !ROUTINE: CG2D_NSA
10  C     !INTERFACE:  C     !INTERFACE:
11        SUBROUTINE CG2D_NSA(          SUBROUTINE CG2D_NSA(
12       I                cg2d_b,       U                cg2d_b, cg2d_x,
13       U                cg2d_x,       O                firstResidual, minResidualSq, lastResidual,
14       O                firstResidual,       U                numIters, nIterMin,
      O                lastResidual,  
      U                numIters,  
15       I                myThid )       I                myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | SUBROUTINE CG2D_NSA                                            C     | SUBROUTINE CG2D_NSA
19  C     | o Two-dimensional grid problem conjugate-gradient          C     | o Two-dimensional grid problem conjugate-gradient
20  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
21  C     | o This version is used only in the case when the matrix  C     | o This version is used only in the case when the matrix
22  C     |   operator is not "self-adjoint" (NSA). Any remaining  C     |   operator is not "self-adjoint" (NSA). Any remaining
23  C     |   residuals will immediately reported to the department  C     |   residuals will immediately reported to the department
24  C     |   of homeland security.  C     |   of homeland security.
25  C     *==========================================================*  C     *==========================================================*
26  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
27  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
28  C     | This implementation assumes A is a five-diagonal            C     | This implementation assumes A is a five-diagonal
29  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
30  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
31  C     | two-dimensional space.                                      C     | two-dimensional space.
32  C     | Notes:                                                      C     | Notes:
33  C     | ======                                                      C     | ======
34  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
35  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
36  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
37  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
38  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
39  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
40  C     | the threads to solve several different problems              C     | the threads to solve several different problems
41  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
42  C     *==========================================================*  C     *==========================================================*
43  C     \ev  C     \ev
44    
# Line 54  C     === Global data === Line 48  C     === Global data ===
48  #include "SIZE.h"  #include "SIZE.h"
49  #include "EEPARAMS.h"  #include "EEPARAMS.h"
50  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
51  #include "CG2D.h"  #include "CG2D.h"
 #include "SURFACE.h"  
52  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
53  # include "tamc.h"  # include "tamc.h"
54  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 64  C     === Global data === Line 56  C     === Global data ===
56    
57  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
58  C     === Routine arguments ===  C     === Routine arguments ===
59  C     myThid    - Thread on which I am working.  C     cg2d_b    :: The source term or "right hand side" (Output: normalised RHS)
60  C     cg2d_b    - The source term or "right hand side"  C     cg2d_x    :: The solution (Input: first guess)
61  C     cg2d_x    - The solution  C     firstResidual :: the initial residual before any iterations
62  C     firstResidual - the initial residual before any iterations  C     minResidualSq :: the lowest residual reached (squared)
63  C     lastResidual  - the actual residual reached  C     lastResidual  :: the actual residual reached
64  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  :: Inp: the maximum number of iterations allowed
65  C                 Exit:  the actual number of iterations used  C                  Out: the actual number of iterations used
66    C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
67    C                  Out: iteration number corresponding to lowest residual
68    C     myThid    :: Thread on which I am working.
69        _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)
70        _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)
71        _RL  firstResidual        _RL  firstResidual
72          _RL  minResidualSq
73        _RL  lastResidual        _RL  lastResidual
74        INTEGER numIters        INTEGER numIters
75          INTEGER nIterMin
76        INTEGER myThid        INTEGER myThid
77    
78  #ifdef ALLOW_CG2D_NSA  #ifdef ALLOW_CG2D_NSA
79  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
80  C     === Local variables ====  C     === Local variables ====
81  C     actualIts      - Number of iterations taken  C     bi, bj     :: tile index in X and Y.
82  C     actualResidual - residual  C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
83  C     bi          - Block index in X and Y.  C     actualIts  :: actual CG iteration number
84  C     bj  C     err_sq     :: Measure of the square of the residual of Ax - b.
85  C     eta_qrN     - Used in computing search directions  C     eta_qrN    :: Used in computing search directions; suffix N and NM1
86  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     denote current and previous iterations respectively.
87  C     cgBeta        previous iterations respectively.  C     recip_eta_qrNM1 :: reciprocal of eta_qrNM1
88  C     recip_eta_qrNM1 reciprocal of eta_qrNM1  C     cgBeta     :: coeff used to update conjugate direction vector "s".
89  C     alpha    C     alpha      :: coeff used to update solution & residual
90  C     alpha_aux   - to avoid the statement: alpha = 1./alpha (for TAMC/TAF)  C     alphaSum   :: to avoid the statement: alpha = 1./alpha (for TAMC/TAF)
91  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
92  C                   useful debuggin/trouble shooting diagnostic.  C                   debugging/trouble shooting diagnostic. For neumann problems
93  C                   For neumann problems sumRHS needs to be ~0.  C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
94  C                   or they converge at a non-zero residual.  C     cg2d_min   :: used to store solution corresponding to lowest residual.
95  C     err         - Measure of residual of Ax - b, usually the norm.  C     msgBuf     :: Informational/error message buffer
96  C     err_sq      - square of err (for TAMC/TAF)        INTEGER bi, bj
97  C     I, J, N     - Loop counters ( N counts CG iterations )        INTEGER i, j, it2d
98        INTEGER actualIts        INTEGER actualIts
       _RL    actualResidual  
       INTEGER bi, bj                
       INTEGER I, J, it2d  
   
       _RL    err  
       _RL    err_sq  
       _RL    eta_qrN  
       _RL    eta_qrNM1  
       _RL    recip_eta_qrNM1  
       _RL    cgBeta  
       _RL    alpha  
       _RL    alpha_aux  
       _RL    sumRHS  
       _RL    rhsMax, rhsMaxGlobal  
       _RL    rhsNorm  
99        _RL    cg2dTolerance_sq        _RL    cg2dTolerance_sq
100          _RL    err_sq,  errTile(nSx,nSy)
101          _RL    eta_qrN, eta_qrNtile(nSx,nSy)
102          _RL    eta_qrNM1, recip_eta_qrNM1
103          _RL    cgBeta,  alpha
104          _RL    alphaSum,alphaTile(nSx,nSy)
105          _RL    sumRHS,  sumRHStile(nSx,nSy)
106          _RL    rhsMax,  rhsNorm
107          CHARACTER*(MAX_LEN_MBUF) msgBuf
108          LOGICAL printResidual
109  CEOP  CEOP
110    
111  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
112        IF ( numIters .GT. numItersMax ) THEN        IF ( numIters .GT. numItersMax ) THEN
113           WRITE(standardMessageUnit,'(A,I10)')          WRITE(msgBuf,'(A,I10)')
114       &        'CG2D_NSA: numIters > numItersMax = ', numItersMax       &    'CG2D_NSA: numIters > numItersMax =', numItersMax
115           STOP 'NON-NORMAL in CG2D_NSA'          CALL PRINT_ERROR( msgBuf, myThid )
116            STOP 'ABNORMAL END: S/R CG2D_NSA'
117          ENDIF
118          IF ( cg2dNormaliseRHS ) THEN
119            WRITE(msgBuf,'(A)') 'CG2D_NSA: cg2dNormaliseRHS is disabled'
120            CALL PRINT_ERROR( msgBuf, myThid )
121            WRITE(msgBuf,'(A)')
122         &    'set cg2dTargetResWunit (instead of cg2dTargetResidual)'
123            CALL PRINT_ERROR( msgBuf, myThid )
124            STOP 'ABNORMAL END: S/R CG2D_NSA'
125        ENDIF        ENDIF
126  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
127    
 CcnhDebugStarts  
 C     CHARACTER*(MAX_LEN_FNAM) suff  
 CcnhDebugEnds  
   
128  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
129            act1 = myThid - 1            act1 = myThid - 1
130            max1 = nTx*nTy            max1 = nTx*nTy
# Line 136  CcnhDebugEnds Line 132  CcnhDebugEnds
132            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
133  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
134    
135  C--   Initialise inverter  C--   Initialise auxiliary constant, some output variable and inverter
136        eta_qrNM1 = 1. _d 0        cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
137        recip_eta_qrNM1 = 1./eta_qrNM1        minResidualSq = -1. _d 0
138          eta_qrNM1     =  1. _d 0
139  CcnhDebugStarts        recip_eta_qrNM1= 1. _d 0
 C     _EXCH_XY_R8( cg2d_b, myThid )  
 C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )  
 C     suff = 'unnormalised'  
 C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)  
 C     STOP  
 CcnhDebugEnds  
140    
141  C--   Normalise RHS  C--   Normalise RHS
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
142        rhsMax = 0. _d 0        rhsMax = 0. _d 0
143        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
144         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
145          DO J=1,sNy          DO j=1,sNy
146           DO I=1,sNx           DO i=1,sNx
147            cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*cg2dNorm            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
148            rhsMax = MAX(ABS(cg2d_b(I,J,bi,bj)),rhsMax)            rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
149           ENDDO           ENDDO
150          ENDDO          ENDDO
151         ENDDO         ENDDO
152        ENDDO        ENDDO
153    
154    #ifndef ALLOW_AUTODIFF_TAMC
155        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
156  C     -  Normalise RHS :  C-    Normalise RHS :
157  #ifdef LETS_MAKE_JAM         _GLOBAL_MAX_RL( rhsMax, myThid )
 C     _GLOBAL_MAX_R8( rhsMax, myThid )  
       rhsMaxGlobal=1.  
 #else  
 #ifdef ALLOW_CONST_RHSMAX  
       rhsMaxGlobal=1.  
 #else  
       rhsMaxGlobal=rhsMax  
       _GLOBAL_MAX_R8( rhsMaxGlobal, myThid )  
 #endif /* ALLOW_CONST_RHSMAX */  
 #endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
       IF ( rhsMaxGlobal .NE. 0. ) THEN  
        rhsNorm = 1. _d 0 / rhsMaxGlobal  
       ELSE  
158         rhsNorm = 1. _d 0         rhsNorm = 1. _d 0
159        ENDIF         IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
160           DO bj=myByLo(myThid),myByHi(myThid)
161  #ifdef ALLOW_AUTODIFF_TAMC          DO bi=myBxLo(myThid),myBxHi(myThid)
162  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte           DO j=1,sNy
163  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte            DO i=1,sNx
164  #endif /* ALLOW_AUTODIFF_TAMC */             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
165        DO bj=myByLo(myThid),myByHi(myThid)             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
166         DO bi=myBxLo(myThid),myBxHi(myThid)            ENDDO
         DO J=1,sNy  
          DO I=1,sNx  
           cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*rhsNorm  
           cg2d_x(I,J,bi,bj) = cg2d_x(I,J,bi,bj)*rhsNorm  
167           ENDDO           ENDDO
168          ENDDO          ENDDO
169         ENDDO         ENDDO
170        ENDDO  C-    end Normalise RHS
 C- end Normalise RHS  
171        ENDIF        ENDIF
172    #endif /* ndef ALLOW_AUTODIFF_TAMC */
173    
174  C--   Update overlaps  C--   Update overlaps
175        _EXCH_XY_R8( cg2d_b, myThid )        CALL EXCH_XY_RL( cg2d_x, myThid )
       _EXCH_XY_R8( cg2d_x, myThid )  
 CcnhDebugStarts  
 C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )  
 C     suff = 'normalised'  
 C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)  
 CcnhDebugEnds  
176    
177  C--   Initial residual calculation  C--   Initial residual calculation
       err = 0. _d 0  
       err_sq = 0. _d 0  
       sumRHS = 0. _d 0  
178  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
179  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte
180  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte
181  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
182        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
183         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
184          DO J=1,sNy          errTile(bi,bj)    = 0. _d 0
185           DO I=1,sNx          sumRHStile(bi,bj) = 0. _d 0
186            cg2d_s(I,J,bi,bj) = 0.          DO j=0,sNy+1
187            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -           DO i=0,sNx+1
188       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)            cg2d_s(i,j,bi,bj) = 0.
189       &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)           ENDDO
190       &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)          ENDDO
191       &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)          DO j=1,sNy
192       &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)           DO i=1,sNx
193       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)            cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
194       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
195       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
196       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*       &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
197       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm       &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
198  CML     &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
199       &    )       &    )
200            err_sq         = err_sq         +            errTile(bi,bj)    = errTile(bi,bj)
201       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
202            sumRHS            = sumRHS            +            sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
      &     cg2d_b(I,J,bi,bj)  
203           ENDDO           ENDDO
204          ENDDO          ENDDO
205         ENDDO         ENDDO
206        ENDDO        ENDDO
207    
208  #ifdef LETS_MAKE_JAM  c     CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
209        CALL EXCH_XY_O1_R8_JAM( cg2d_r )        CALL EXCH_XY_RL ( cg2d_r, myThid )
210        CALL EXCH_XY_O1_R8_JAM( cg2d_s )        CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
211  #else        CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
212        _EXCH_XY_R8( cg2d_r, myThid )        actualIts = 0
213        _EXCH_XY_R8( cg2d_s, myThid )        IF ( err_sq .NE. 0. ) THEN
214  #endif          firstResidual = SQRT(err_sq)
215         _GLOBAL_SUM_R8( sumRHS, myThid )        ELSE
216         _GLOBAL_SUM_R8( err_sq, myThid )          firstResidual = 0.
217         if ( err_sq .ne. 0. ) then        ENDIF
218            err = SQRT(err_sq)  
219         else        printResidual = .FALSE.
220            err = 0.        IF ( debugLevel .GE. debLevZero ) THEN
221         end if          _BEGIN_MASTER( myThid )
222         actualIts      = 0          printResidual = printResidualFreq.GE.1
223         actualResidual = err          WRITE(standardmessageunit,'(A,1P2E22.14)')
224               &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
225        _BEGIN_MASTER( myThid )          _END_MASTER( myThid )
226        write(standardMessageUnit,'(A,1P2E22.14)')        ENDIF
      &     ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMaxGlobal  
       _END_MASTER( )  
 C     _BARRIER  
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',  
 c    & actualIts, actualResidual  
 c     _END_MASTER( )  
       firstResidual=actualResidual  
       cg2dTolerance_sq = cg2dTolerance**2  
227    
228  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
229  Cml   begin main solver loop  Cml   begin main solver loop
# Line 284  CADJ LOOP = iteration, cg2d_x = comlev_c Line 233  CADJ LOOP = iteration, cg2d_x = comlev_c
233        DO it2d=1, numIters        DO it2d=1, numIters
234  #ifdef ALLOW_LOOP_DIRECTIVE  #ifdef ALLOW_LOOP_DIRECTIVE
235  CML      it2d = 0  CML      it2d = 0
236  CML      DO WHILE ( err_sq .GT. cg2dTolerance_sq .and. it2d .LT. numIters )  CML      DO WHILE ( err_sq .GT. cg2dTolerance_sq .and. it2d .LT. numIters )
237  CML      it2d = it2d+1  CML      it2d = it2d+1
238  #endif /* ALLOW_LOOP_DIRECTIVE */  #endif /* ALLOW_LOOP_DIRECTIVE */
239    
240  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
241         icg2dkey = (ikey-1)*numItersMax + it2d         icg2dkey = (ikey-1)*numItersMax + it2d
242  CMLCADJ STORE err = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE err_sq = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
 CADJ STORE err_sq = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  
243  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
244  CML      IF ( err .LT. cg2dTolerance ) THEN        IF ( err_sq .GE. cg2dTolerance_sq ) THEN
       IF ( err_sq .LT. cg2dTolerance_sq ) THEN  
 Cml      DO NOTHING  
       ELSE  
245    
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' residual = ',  
 C    &  actualResidual  
 CcnhDebugEnds  
246  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
247  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
        eta_qrN = 0. _d 0  
248  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
249  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
250  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
251  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
252         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
253          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
254           DO J=1,sNy           eta_qrNtile(bi,bj) = 0. _d 0
255            DO I=1,sNx           DO j=1,sNy
256             cg2d_z(I,J,bi,bj) =            DO i=1,sNx
257       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)             cg2d_z(i,j,bi,bj) =
258       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)       &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
259       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)       &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
260       &     +pS(I  ,J  ,bi,bj)*cg2d_r(I  ,J-1,bi,bj)       &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
261       &     +pS(I  ,J+1,bi,bj)*cg2d_r(I  ,J+1,bi,bj)       &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
262  CcnhDebugStarts       &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
263  C          cg2d_z(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)             eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
264  CcnhDebugEnds       &     +cg2d_z(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
            eta_qrN = eta_qrN  
      &     +cg2d_z(I,J,bi,bj)*cg2d_r(I,J,bi,bj)  
265            ENDDO            ENDDO
266           ENDDO           ENDDO
267          ENDDO          ENDDO
268         ENDDO         ENDDO
269    
270         _GLOBAL_SUM_R8(eta_qrN, myThid)         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' eta_qrN = ',eta_qrN  
 CcnhDebugEnds  
271  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
272  CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
273  CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
274  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
275  CML       cgBeta   = eta_qrN/eta_qrNM1  CML       cgBeta   = eta_qrN/eta_qrNM1
276         cgBeta   = eta_qrN*recip_eta_qrNM1         cgBeta   = eta_qrN*recip_eta_qrNM1
277  CcnhDebugStarts  Cml    store normalisation factor for the next interation (in case there is one).
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' beta = ',cgBeta  
 CcnhDebugEnds  
 Cml    store normalisation factor for the next interation  
 Cml    (in case there is one).  
278  CML    store the inverse of the normalization factor for higher precision  CML    store the inverse of the normalization factor for higher precision
279  CML       eta_qrNM1 = eta_qrN  CML       eta_qrNM1 = eta_qrN
280         recip_eta_qrNM1 = 1./eta_qrN         recip_eta_qrNM1 = 1. _d 0/eta_qrN
281    
282         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
283          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
284           DO J=1,sNy           DO j=1,sNy
285            DO I=1,sNx            DO i=1,sNx
286             cg2d_s(I,J,bi,bj) = cg2d_z(I,J,bi,bj)             cg2d_s(i,j,bi,bj) = cg2d_z(i,j,bi,bj)
287       &                       + cgBeta*cg2d_s(I,J,bi,bj)       &                       + cgBeta*cg2d_s(i,j,bi,bj)
288            ENDDO            ENDDO
289           ENDDO           ENDDO
290          ENDDO          ENDDO
# Line 361  CML       eta_qrNM1 = eta_qrN Line 292  CML       eta_qrNM1 = eta_qrN
292    
293  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
294  C--    processes.  C--    processes.
295  #ifdef LETS_MAKE_JAM  c      CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
296        CALL EXCH_XY_O1_R8_JAM( cg2d_s )         CALL EXCH_XY_RL ( cg2d_s, myThid )
 #else  
        _EXCH_XY_R8( cg2d_s, myThid )  
 #endif  
297    
298  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
299  C==    q = A.s  C==    q = A.s
        alpha     = 0. _d 0  
        alpha_aux = 0. _d 0  
300  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
301  #ifndef ALLOW_LOOP_DIRECTIVE  #ifndef ALLOW_LOOP_DIRECTIVE
302  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
303  #endif /* not ALLOW_LOOP_DIRECTIVE */  #endif /* not ALLOW_LOOP_DIRECTIVE */
304  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
305         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
306          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
307           DO J=1,sNy           alphaTile(bi,bj) = 0. _d 0
308            DO I=1,sNx           DO j=1,sNy
309             cg2d_q(I,J,bi,bj) =            DO i=1,sNx
310       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)             cg2d_q(i,j,bi,bj) =
311       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)       &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
312       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)       &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
313       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)       &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
314       &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
315       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
316       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)             alphaTile(bi,bj) = alphaTile(bi,bj)
317       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &                      + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
      &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  
      &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  
 CML     &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm  
           alpha_aux = alpha_aux+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)  
318            ENDDO            ENDDO
319           ENDDO           ENDDO
320          ENDDO          ENDDO
321         ENDDO         ENDDO
322         _GLOBAL_SUM_R8(alpha_aux,myThid)         CALL GLOBAL_SUM_TILE_RL( alphaTile, alphaSum, myThid )
323  CcnhDebugStarts         alpha = eta_qrN/alphaSum
324  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' SUM(s*q)= ',alpha_aux  
325  CcnhDebugEnds  C==    Update simultaneously solution and residual vectors (and Iter number)
        alpha = eta_qrN/alpha_aux  
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' alpha= ',alpha  
 CcnhDebugEnds  
       
 C==    Update solution and residual vectors  
326  C      Now compute "interior" points.  C      Now compute "interior" points.
        err    = 0. _d 0  
        err_sq = 0. _d 0  
327  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
328  #ifndef ALLOW_LOOP_DIRECTIVE  #ifndef ALLOW_LOOP_DIRECTIVE
329  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
330  #endif /* ALLOW_LOOP_DIRECTIVE */  #endif /* ALLOW_LOOP_DIRECTIVE */
331  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
332         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
333          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
334           DO J=1,sNy           errTile(bi,bj) = 0. _d 0
335            DO I=1,sNx           DO j=1,sNy
336             cg2d_x(I,J,bi,bj)=cg2d_x(I,J,bi,bj)+alpha*cg2d_s(I,J,bi,bj)            DO i=1,sNx
337             cg2d_r(I,J,bi,bj)=cg2d_r(I,J,bi,bj)-alpha*cg2d_q(I,J,bi,bj)             cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
338             err_sq = err_sq+cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)             cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
339               errTile(bi,bj) = errTile(bi,bj)
340         &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
341            ENDDO            ENDDO
342           ENDDO           ENDDO
343          ENDDO          ENDDO
344         ENDDO         ENDDO
345           actualIts = it2d
346    
347         _GLOBAL_SUM_R8( err_sq   , myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq,    myThid )
348         if ( err_sq .ne. 0. ) then         IF ( printResidual ) THEN
349            err = SQRT(err_sq)          IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
350         else           WRITE(msgBuf,'(A,I6,A,1PE21.14)')
351            err = 0.       &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
352         end if           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
353         actualIts      = it2d       &                       SQUEEZE_RIGHT, myThid )
354         actualResidual = err          ENDIF
355           ENDIF
356    
357  #ifdef LETS_MAKE_JAM  c      CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
358        CALL EXCH_XY_O1_R8_JAM( cg2d_r )         CALL EXCH_XY_RL ( cg2d_r, myThid )
 #else  
       _EXCH_XY_R8( cg2d_r, myThid )  
       _EXCH_XY_R8( cg2d_x, myThid )  
 #endif  
359    
360  Cml   end of IF ( err .LT. cg2dTolerance ) THEN; ELSE  Cml   end of if "err >= cg2dTolerance" block ; end main solver loop
361        ENDIF        ENDIF
 Cml     end main solver loop  
362        ENDDO        ENDDO
363    
364    #ifndef ALLOW_AUTODIFF_TAMC
365        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte  
 CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
366  C--   Un-normalise the answer  C--   Un-normalise the answer
367          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
368           DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
369            DO J=1,sNy            DO j=1,sNy
370             DO I=1,sNx             DO i=1,sNx
371              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
372             ENDDO             ENDDO
373            ENDDO            ENDDO
374           ENDDO           ENDDO
375          ENDDO          ENDDO
376        ENDIF        ENDIF
377    #endif /* ndef ALLOW_AUTODIFF_TAMC */
 C     The following exchange was moved up to solve_for_pressure  
 C     for compatibility with TAMC.  
 C     _EXCH_XY_R8(cg2d_x, myThid )  
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',  
 c    & actualIts, actualResidual  
 c     _END_MASTER( )  
378    
379  C--   Return parameters to caller  C--   Return parameters to caller
380        lastResidual=actualResidual        IF ( err_sq .NE. 0. ) THEN
381        numIters=actualIts          lastResidual = SQRT(err_sq)
382          ELSE
383            lastResidual = 0.
384          ENDIF
385          numIters = actualIts
386    
387  #endif /* ALLOW_CG2D_NSA */  #endif /* ALLOW_CG2D_NSA */
388        RETURN        RETURN
389        END        END
390    
391  # if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))  #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
392  C  
393  C     These routines are routinely part of the TAMC/TAF library that is  C     These routines are routinely part of the TAMC/TAF library that is
394  C     not included in the MITcgm, therefore they are mimicked here.  C     not included in the MITcgm, therefore they are mimicked here.
395  C  
396        subroutine adstore(chardum,int1,idow,int2,int3,icount)        subroutine adstore(chardum,int1,idow,int2,int3,icount)
397    
398        implicit none        implicit none
# Line 496  C Line 403  C
403        character*(*) chardum        character*(*) chardum
404        integer int1, int2, int3, idow, icount        integer int1, int2, int3, idow, icount
405    
406  C     the length of this vector must be greater or equal  C     the length of this vector must be greater or equal
407  C     twice the number of timesteps  C     twice the number of timesteps
408        integer nidow        integer nidow
409  #ifdef ALLOW_TAMC_CHECKPOINTING  #ifdef ALLOW_TAMC_CHECKPOINTING
# Line 507  C     twice the number of timesteps Line 414  C     twice the number of timesteps
414        integer istoreidow(nidow)        integer istoreidow(nidow)
415        common /istorecommon/ istoreidow        common /istorecommon/ istoreidow
416    
417        print *, 'adstore: ', chardum, int1, idow, int2, int3, icount        print *, 'adstore: ', chardum, int1, idow, int2, int3, icount
418    
419        if ( icount .gt. nidow ) then        if ( icount .gt. nidow ) then
420           print *, 'adstore: error: icount > nidow = ', nidow           print *, 'adstore: error: icount > nidow = ', nidow
# Line 529  C     twice the number of timesteps Line 436  C     twice the number of timesteps
436        character*(*) chardum        character*(*) chardum
437        integer int1, int2, int3, idow, icount        integer int1, int2, int3, idow, icount
438    
439    C     the length of this vector must be greater or equal
 C     the length of this vector must be greater or equal  
440  C     twice the number of timesteps  C     twice the number of timesteps
441        integer nidow        integer nidow
442  #ifdef ALLOW_TAMC_CHECKPOINTING  #ifdef ALLOW_TAMC_CHECKPOINTING
# Line 541  C     twice the number of timesteps Line 447  C     twice the number of timesteps
447        integer istoreidow(nidow)        integer istoreidow(nidow)
448        common /istorecommon/ istoreidow        common /istorecommon/ istoreidow
449    
450        print *, 'adresto: ', chardum, int1, idow, int2, int3, icount        print *, 'adresto: ', chardum, int1, idow, int2, int3, icount
451    
452        if ( icount .gt. nidow ) then        if ( icount .gt. nidow ) then
453           print *, 'adstore: error: icount > nidow = ', nidow           print *, 'adstore: error: icount > nidow = ', nidow
# Line 552  C     twice the number of timesteps Line 458  C     twice the number of timesteps
458    
459        return        return
460        end        end
461  # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */  #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */
   
   
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22