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