/[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.2 by jmc, Tue Apr 28 18:01:14 2009 UTC revision 1.3 by jmc, Sat Jul 11 21:45:44 2009 UTC
# Line 5  C $Name$ Line 5  C $Name$
5  #ifdef ALLOW_USE_MPI  #ifdef ALLOW_USE_MPI
6  C HACK to avoid global_max  C HACK to avoid global_max
7  # define ALLOW_CONST_RHSMAX  # define ALLOW_CONST_RHSMAX
8  #endif  #endif
9    
10  CML THIS DOES NOT WORK +++++  CML THIS DOES NOT WORK +++++
11  #undef ALLOW_LOOP_DIRECTIVE  #undef ALLOW_LOOP_DIRECTIVE
12  CBOP  CBOP
13  C     !ROUTINE: CG2D_NSA  C     !ROUTINE: CG2D_NSA
14  C     !INTERFACE:  C     !INTERFACE:
15        SUBROUTINE CG2D_NSA(          SUBROUTINE CG2D_NSA(
16       I                cg2d_b,       I                cg2d_b,
17       U                cg2d_x,       U                cg2d_x,
18       O                firstResidual,       O                firstResidual,
# Line 21  C     !INTERFACE: Line 21  C     !INTERFACE:
21       I                myThid )       I                myThid )
22  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
23  C     *==========================================================*  C     *==========================================================*
24  C     | SUBROUTINE CG2D_NSA                                            C     | SUBROUTINE CG2D_NSA
25  C     | o Two-dimensional grid problem conjugate-gradient          C     | o Two-dimensional grid problem conjugate-gradient
26  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
27  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
28  C     |   operator is not "self-adjoint" (NSA). Any remaining  C     |   operator is not "self-adjoint" (NSA). Any remaining
29  C     |   residuals will immediately reported to the department  C     |   residuals will immediately reported to the department
30  C     |   of homeland security.  C     |   of homeland security.
31  C     *==========================================================*  C     *==========================================================*
32  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
33  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
34  C     | This implementation assumes A is a five-diagonal            C     | This implementation assumes A is a five-diagonal
35  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
36  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
37  C     | two-dimensional space.                                      C     | two-dimensional space.
38  C     | Notes:                                                      C     | Notes:
39  C     | ======                                                      C     | ======
40  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
41  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
42  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
43  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
44  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
45  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
46  C     | the threads to solve several different problems              C     | the threads to solve several different problems
47  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
48  C     *==========================================================*  C     *==========================================================*
49  C     \ev  C     \ev
50    
# Line 54  C     === Global data === Line 54  C     === Global data ===
54  #include "SIZE.h"  #include "SIZE.h"
55  #include "EEPARAMS.h"  #include "EEPARAMS.h"
56  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
57  #include "CG2D.h"  #include "CG2D.h"
58  #include "SURFACE.h"  c#include "GRID.h"
59    c#include "SURFACE.h"
60  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
61  # include "tamc.h"  # include "tamc.h"
62  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 64  C     === Global data === Line 64  C     === Global data ===
64    
65  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
66  C     === Routine arguments ===  C     === Routine arguments ===
67  C     myThid    - Thread on which I am working.  C     cg2d_b    :: The source term or "right hand side"
68  C     cg2d_b    - The source term or "right hand side"  C     cg2d_x    :: The solution
69  C     cg2d_x    - The solution  C     firstResidual :: the initial residual before any iterations
70  C     firstResidual - the initial residual before any iterations  C     lastResidual  :: the actual residual reached
71  C     lastResidual  - the actual residual reached  C     numIters  :: Entry: the maximum number of iterations allowed
72  C     numIters  - Entry: the maximum number of iterations allowed  C                  Exit:  the actual number of iterations used
73  C                 Exit:  the actual number of iterations used  C     myThid    :: Thread on which I am working.
74        _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)
75        _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)
76        _RL  firstResidual        _RL  firstResidual
# Line 81  C                 Exit:  the actual numb Line 81  C                 Exit:  the actual numb
81  #ifdef ALLOW_CG2D_NSA  #ifdef ALLOW_CG2D_NSA
82  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
83  C     === Local variables ====  C     === Local variables ====
84  C     actualIts      - Number of iterations taken  C     actualIts      :: Number of iterations taken
85  C     actualResidual - residual  C     actualResidual :: residual
86  C     bi          - Block index in X and Y.  C     bi, bj     :: Block index in X and Y.
87  C     bj  C     eta_qrN    :: Used in computing search directions
 C     eta_qrN     - Used in computing search directions  
88  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
89  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
90  C     recip_eta_qrNM1 reciprocal of eta_qrNM1  C     recip_eta_qrNM1 :: reciprocal of eta_qrNM1
91  C     alpha    C     alpha
92  C     alpha_aux   - to avoid the statement: alpha = 1./alpha (for TAMC/TAF)  C     alpha_aux  :: to avoid the statement: alpha = 1./alpha (for TAMC/TAF)
93  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a
94  C                   useful debuggin/trouble shooting diagnostic.  C                   useful debuggin/trouble shooting diagnostic.
95  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
96  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
97  C     err         - Measure of residual of Ax - b, usually the norm.  C     err        :: Measure of residual of Ax - b, usually the norm.
98  C     err_sq      - square of err (for TAMC/TAF)  C     err_sq     :: square of err (for TAMC/TAF)
99  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, it2d :: Loop counters ( it2d counts CG iterations )
100        INTEGER actualIts        INTEGER actualIts
101        _RL    actualResidual        _RL    actualResidual
102        INTEGER bi, bj                      INTEGER bi, bj
103        INTEGER I, J, it2d        INTEGER I, J, it2d
104    
105        _RL    err        _RL    err
# Line 119  CEOP Line 118  CEOP
118    
119  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
120        IF ( numIters .GT. numItersMax ) THEN        IF ( numIters .GT. numItersMax ) THEN
121           WRITE(standardMessageUnit,'(A,I10)')           WRITE(standardMessageUnit,'(A,I10)')
122       &        'CG2D_NSA: numIters > numItersMax = ', numItersMax       &        'CG2D_NSA: numIters > numItersMax = ', numItersMax
123           STOP 'NON-NORMAL in CG2D_NSA'           STOP 'NON-NORMAL in CG2D_NSA'
124        ENDIF        ENDIF
# Line 150  CcnhDebugEnds Line 149  CcnhDebugEnds
149    
150  C--   Normalise RHS  C--   Normalise RHS
151  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
152  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte
153  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
154        rhsMax = 0. _d 0        rhsMax = 0. _d 0
155        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 178  C     _GLOBAL_MAX_RL( rhsMax, myThid ) Line 177  C     _GLOBAL_MAX_RL( rhsMax, myThid )
177  #endif /* ALLOW_CONST_RHSMAX */  #endif /* ALLOW_CONST_RHSMAX */
178  #endif  #endif
179  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
180  CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte
181  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
182        IF ( rhsMaxGlobal .NE. 0. ) THEN        IF ( rhsMaxGlobal .NE. 0. ) THEN
183         rhsNorm = 1. _d 0 / rhsMaxGlobal         rhsNorm = 1. _d 0 / rhsMaxGlobal
# Line 187  CADJ STORE rhsNorm = comlev1_cg2d, key = Line 186  CADJ STORE rhsNorm = comlev1_cg2d, key =
186        ENDIF        ENDIF
187    
188  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
189  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte
190  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte
191  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
192        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
193         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 204  C- end Normalise RHS Line 203  C- end Normalise RHS
203        ENDIF        ENDIF
204    
205  C--   Update overlaps  C--   Update overlaps
206        _EXCH_XY_RL( cg2d_b, myThid )  c     CALL EXCH_XY_RL( cg2d_b, myThid )
207        _EXCH_XY_RL( cg2d_x, myThid )        CALL EXCH_XY_RL( cg2d_x, myThid )
208  CcnhDebugStarts  CcnhDebugStarts
209  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )
210  C     suff = 'normalised'  C     suff = 'normalised'
# Line 217  C--   Initial residual calculation Line 216  C--   Initial residual calculation
216        err_sq = 0. _d 0        err_sq = 0. _d 0
217        sumRHS = 0. _d 0        sumRHS = 0. _d 0
218  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
219  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_b = comlev1_cg2d, key = ikey, byte = isbyte
220  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte
221  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
222        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
223         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
224            DO J=1-1,sNy+1
225             DO I=1-1,sNx+1
226              cg2d_s(I,J,bi,bj) = 0.
227             ENDDO
228            ENDDO
229          DO J=1,sNy          DO J=1,sNy
230           DO I=1,sNx           DO I=1,sNx
           cg2d_s(I,J,bi,bj) = 0.  
231            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -
232       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)
233       &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)       &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)
234       &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)       &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)
235       &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)       &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)
236       &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    +aC2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
      &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
      &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
      &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
      &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  
      &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  
 CML     &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm  
237       &    )       &    )
238            err_sq         = err_sq         +  c    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
239    c    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
240    c    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
241    c    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
242    c    &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
243    c    &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm
244    cML     &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
245    c    &    )
246              err_sq         = err_sq         +
247       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
248            sumRHS            = sumRHS            +            sumRHS            = sumRHS            +
249       &     cg2d_b(I,J,bi,bj)       &     cg2d_b(I,J,bi,bj)
# Line 247  CML     &     cg2d_x(I  ,J  ,bi,bj)/delt Line 252  CML     &     cg2d_x(I  ,J  ,bi,bj)/delt
252         ENDDO         ENDDO
253        ENDDO        ENDDO
254    
255  #ifdef LETS_MAKE_JAM  c     CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
256        CALL EXCH_XY_O1_R8_JAM( cg2d_r )        CALL EXCH_XY_RL ( cg2d_r, myThid )
257        CALL EXCH_XY_O1_R8_JAM( cg2d_s )        _GLOBAL_SUM_RL( sumRHS, myThid )
258  #else        _GLOBAL_SUM_RL( err_sq, myThid )
259        _EXCH_XY_RL( cg2d_r, myThid )         IF ( err_sq .NE. 0. ) THEN
       _EXCH_XY_RL( cg2d_s, myThid )  
 #endif  
        _GLOBAL_SUM_RL( sumRHS, myThid )  
        _GLOBAL_SUM_RL( err_sq, myThid )  
        if ( err_sq .ne. 0. ) then  
260            err = SQRT(err_sq)            err = SQRT(err_sq)
261         else         ELSE
262            err = 0.            err = 0.
263         end if         ENDIF
264         actualIts      = 0         actualIts      = 0
265         actualResidual = err         actualResidual = err
266          
267        _BEGIN_MASTER( myThid )        IF ( debugLevel .GE. debLevZero ) THEN
268        write(standardMessageUnit,'(A,1P2E22.14)')          _BEGIN_MASTER( myThid )
269       &     ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMaxGlobal          WRITE(standardmessageunit,'(A,1P2E22.14)')
270        _END_MASTER( )       &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMaxGlobal
271            _END_MASTER( myThid )
272          ENDIF
273  C     _BARRIER  C     _BARRIER
274  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
275  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',
276  c    & actualIts, actualResidual  c    & actualIts, actualResidual
277  c     _END_MASTER( )  c     _END_MASTER( myThid )
278        firstResidual=actualResidual        firstResidual=actualResidual
279        cg2dTolerance_sq = cg2dTolerance**2        cg2dTolerance_sq = cg2dTolerance**2
280    
# Line 284  CADJ LOOP = iteration, cg2d_x = comlev_c Line 286  CADJ LOOP = iteration, cg2d_x = comlev_c
286        DO it2d=1, numIters        DO it2d=1, numIters
287  #ifdef ALLOW_LOOP_DIRECTIVE  #ifdef ALLOW_LOOP_DIRECTIVE
288  CML      it2d = 0  CML      it2d = 0
289  CML      DO WHILE ( err_sq .GT. cg2dTolerance_sq .and. it2d .LT. numIters )  CML      DO WHILE ( err_sq .GT. cg2dTolerance_sq .and. it2d .LT. numIters )
290  CML      it2d = it2d+1  CML      it2d = it2d+1
291  #endif /* ALLOW_LOOP_DIRECTIVE */  #endif /* ALLOW_LOOP_DIRECTIVE */
292    
293  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
294         icg2dkey = (ikey-1)*numItersMax + it2d         icg2dkey = (ikey-1)*numItersMax + it2d
295  CMLCADJ STORE err = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CMLCADJ STORE err = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
296  CADJ STORE err_sq = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE err_sq = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
297  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
298  CML      IF ( err .LT. cg2dTolerance ) THEN  CML      IF ( err .LT. cg2dTolerance ) THEN
299        IF ( err_sq .LT. cg2dTolerance_sq ) THEN        IF ( err_sq .LT. cg2dTolerance_sq ) THEN
# Line 306  C--    Solve preconditioning equation an Line 308  C--    Solve preconditioning equation an
308  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
309         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
310  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
311  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
312  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
313  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
314         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
315          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
316           DO J=1,sNy           DO J=1,sNy
317            DO I=1,sNx            DO I=1,sNx
318             cg2d_z(I,J,bi,bj) =             cg2d_z(I,J,bi,bj) =
319       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)
320       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)
321       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)
# Line 334  CcnhDebugStarts Line 336  CcnhDebugStarts
336  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' eta_qrN = ',eta_qrN  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' eta_qrN = ',eta_qrN
337  CcnhDebugEnds  CcnhDebugEnds
338  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
339  CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CMLCADJ STORE eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
340  CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE recip_eta_qrNM1 = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
341  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
342  CML       cgBeta   = eta_qrN/eta_qrNM1  CML       cgBeta   = eta_qrN/eta_qrNM1
343         cgBeta   = eta_qrN*recip_eta_qrNM1         cgBeta   = eta_qrN*recip_eta_qrNM1
344  CcnhDebugStarts  CcnhDebugStarts
345  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' beta = ',cgBeta  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' beta = ',cgBeta
346  CcnhDebugEnds  CcnhDebugEnds
347  Cml    store normalisation factor for the next interation  Cml    store normalisation factor for the next interation
348  Cml    (in case there is one).  Cml    (in case there is one).
349  CML    store the inverse of the normalization factor for higher precision  CML    store the inverse of the normalization factor for higher precision
350  CML       eta_qrNM1 = eta_qrN  CML       eta_qrNM1 = eta_qrN
# Line 361  CML       eta_qrNM1 = eta_qrN Line 363  CML       eta_qrNM1 = eta_qrN
363    
364  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
365  C--    processes.  C--    processes.
366  #ifdef LETS_MAKE_JAM  c      CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
367        CALL EXCH_XY_O1_R8_JAM( cg2d_s )         CALL EXCH_XY_RL ( cg2d_s, myThid )
 #else  
        _EXCH_XY_RL( cg2d_s, myThid )  
 #endif  
368    
369  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
370  C==    q = A.s  C==    q = A.s
# Line 373  C==    q = A.s Line 372  C==    q = A.s
372         alpha_aux = 0. _d 0         alpha_aux = 0. _d 0
373  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
374  #ifndef ALLOW_LOOP_DIRECTIVE  #ifndef ALLOW_LOOP_DIRECTIVE
375  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_s = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
376  #endif /* not ALLOW_LOOP_DIRECTIVE */  #endif /* not ALLOW_LOOP_DIRECTIVE */
377  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
378         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
379          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
380           DO J=1,sNy           DO J=1,sNy
381            DO I=1,sNx            DO I=1,sNx
382             cg2d_q(I,J,bi,bj) =             cg2d_q(I,J,bi,bj) =
383       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)
384       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)
385       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)
386       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)
387       &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aC2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
388       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  c    &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
389       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  c    &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
390       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  c    &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
391       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  c    &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
392       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  c    &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
393  CML     &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm  c    &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm
394            alpha_aux = alpha_aux+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)  cML     &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
395               alpha_aux = alpha_aux+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
396            ENDDO            ENDDO
397           ENDDO           ENDDO
398          ENDDO          ENDDO
# Line 405  CcnhDebugEnds Line 405  CcnhDebugEnds
405  CcnhDebugStarts  CcnhDebugStarts
406  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG2D_NSA: Iteration ',it2d-1,' alpha= ',alpha
407  CcnhDebugEnds  CcnhDebugEnds
408        
409  C==    Update solution and residual vectors  C==    Update solution and residual vectors
410  C      Now compute "interior" points.  C      Now compute "interior" points.
411         err    = 0. _d 0         err    = 0. _d 0
412         err_sq = 0. _d 0         err_sq = 0. _d 0
413  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
414  #ifndef ALLOW_LOOP_DIRECTIVE  #ifndef ALLOW_LOOP_DIRECTIVE
415  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte  CADJ STORE cg2d_r = comlev1_cg2d_iter, key = icg2dkey, byte = isbyte
416  #endif /* ALLOW_LOOP_DIRECTIVE */  #endif /* ALLOW_LOOP_DIRECTIVE */
417  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
418         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
# Line 436  CADJ STORE cg2d_r = comlev1_cg2d_iter, k Line 436  CADJ STORE cg2d_r = comlev1_cg2d_iter, k
436         actualIts      = it2d         actualIts      = it2d
437         actualResidual = err         actualResidual = err
438    
439  #ifdef LETS_MAKE_JAM  c      CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
440        CALL EXCH_XY_O1_R8_JAM( cg2d_r )         CALL EXCH_XY_RL ( cg2d_r, myThid )
 #else  
       _EXCH_XY_RL( cg2d_r, myThid )  
       _EXCH_XY_RL( cg2d_x, myThid )  
 #endif  
441    
442  Cml   end of IF ( err .LT. cg2dTolerance ) THEN; ELSE  Cml   end of IF ( err .LT. cg2dTolerance ) THEN; ELSE
443        ENDIF        ENDIF
# Line 450  Cml     end main solver loop Line 446  Cml     end main solver loop
446    
447        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
448  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
449  CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE rhsNorm = comlev1_cg2d, key = ikey, byte = isbyte
450  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte  CADJ STORE cg2d_x = comlev1_cg2d, key = ikey, byte = isbyte
451  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
452  C--   Un-normalise the answer  C--   Un-normalise the answer
453          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 469  C     The following exchange was moved u Line 465  C     The following exchange was moved u
465  C     for compatibility with TAMC.  C     for compatibility with TAMC.
466  C     _EXCH_XY_RL(cg2d_x, myThid )  C     _EXCH_XY_RL(cg2d_x, myThid )
467  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
468  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D_NSA iters, err = ',
469  c    & actualIts, actualResidual  c    & actualIts, actualResidual
470  c     _END_MASTER( )  c     _END_MASTER( )
471    
# Line 496  C Line 492  C
492        character*(*) chardum        character*(*) chardum
493        integer int1, int2, int3, idow, icount        integer int1, int2, int3, idow, icount
494    
495  C     the length of this vector must be greater or equal  C     the length of this vector must be greater or equal
496  C     twice the number of timesteps  C     twice the number of timesteps
497        integer nidow        integer nidow
498  #ifdef ALLOW_TAMC_CHECKPOINTING  #ifdef ALLOW_TAMC_CHECKPOINTING
# Line 507  C     twice the number of timesteps Line 503  C     twice the number of timesteps
503        integer istoreidow(nidow)        integer istoreidow(nidow)
504        common /istorecommon/ istoreidow        common /istorecommon/ istoreidow
505    
506        print *, 'adstore: ', chardum, int1, idow, int2, int3, icount        print *, 'adstore: ', chardum, int1, idow, int2, int3, icount
507    
508        if ( icount .gt. nidow ) then        if ( icount .gt. nidow ) then
509           print *, 'adstore: error: icount > nidow = ', nidow           print *, 'adstore: error: icount > nidow = ', nidow
# Line 530  C     twice the number of timesteps Line 526  C     twice the number of timesteps
526        integer int1, int2, int3, idow, icount        integer int1, int2, int3, idow, icount
527    
528    
529  C     the length of this vector must be greater or equal  C     the length of this vector must be greater or equal
530  C     twice the number of timesteps  C     twice the number of timesteps
531        integer nidow        integer nidow
532  #ifdef ALLOW_TAMC_CHECKPOINTING  #ifdef ALLOW_TAMC_CHECKPOINTING
# Line 541  C     twice the number of timesteps Line 537  C     twice the number of timesteps
537        integer istoreidow(nidow)        integer istoreidow(nidow)
538        common /istorecommon/ istoreidow        common /istorecommon/ istoreidow
539    
540        print *, 'adresto: ', chardum, int1, idow, int2, int3, icount        print *, 'adresto: ', chardum, int1, idow, int2, int3, icount
541    
542        if ( icount .gt. nidow ) then        if ( icount .gt. nidow ) then
543           print *, 'adstore: error: icount > nidow = ', nidow           print *, 'adstore: error: icount > nidow = ', nidow
# Line 554  C     twice the number of timesteps Line 550  C     twice the number of timesteps
550        end        end
551  # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */  # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */
552    
   
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22