pkg/diagnostics/diag_cg2d.F

d5ada4102b Jean*0001 #include "DIAG_OPTIONS.h"
e938313568 Jean*0002 #undef DEBUG_DIAG_CG2D
d5ada4102b Jean*0003
                0004 CBOP
                0005 C     !ROUTINE: DIAG_CG2D
                0006 C     !INTERFACE:
                0007       SUBROUTINE DIAG_CG2D(
                0008      I                aW2d, aS2d, b2d,
a63b8f5615 Jean*0009      I                offDiagFactor, residCriter,
d53cbeab8e Jean*0010      O                firstResidual, minResidual, lastResidual,
d5ada4102b Jean*0011      U                x2d, numIters,
d53cbeab8e Jean*0012      O                nIterMin,
                0013      I                printResidFrq, myThid )
d5ada4102b Jean*0014 C     !DESCRIPTION: \bv
                0015 C     *==========================================================*
                0016 C     | SUBROUTINE CG2D
                0017 C     | o Two-dimensional grid problem conjugate-gradient
                0018 C     |   inverter (with preconditioner).
                0019 C     *==========================================================*
                0020 C     | Con. grad is an iterative procedure for solving Ax = b.
                0021 C     | It requires the A be symmetric.
                0022 C     | This implementation assumes A is a five-diagonal
                0023 C     | matrix of the form that arises in the discrete
                0024 C     | representation of the del^2 operator in a
                0025 C     | two-dimensional space.
                0026 C     *==========================================================*
                0027 C     \ev
                0028
                0029 C     !USES:
                0030       IMPLICIT NONE
                0031 C     === Global data ===
                0032 #include "SIZE.h"
                0033 #include "EEPARAMS.h"
                0034 #include "PARAMS.h"
                0035
                0036 C     !INPUT/OUTPUT PARAMETERS:
                0037 C     === Routine arguments ===
d53cbeab8e Jean*0038 C     b2d           :: The source term or "right hand side"
                0039 C     x2d           :: The solution
a63b8f5615 Jean*0040 C     offDiagFactor :: preconditioner off-diagonal factor
                0041 C     residCriter   :: residual target for convergence
d5ada4102b Jean*0042 C     firstResidual :: the initial residual before any iterations
d53cbeab8e Jean*0043 C     minResidual   :: the lowest residual reached
d5ada4102b Jean*0044 C     lastResidual  :: the actual residual reached
                0045 C     numIters  :: Entry: the maximum number of iterations allowed
                0046 C                  Exit:  the actual number of iterations used
d53cbeab8e Jean*0047 C     nIterMin      :: iteration number corresponding to lowest residual
                0048 C     printResidFrq :: Frequency for printing residual in CG iterations
                0049 C     myThid        :: my Thread Id number
d5ada4102b Jean*0050       _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0051       _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052       _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053       _RL  residCriter
a63b8f5615 Jean*0054       _RL  offDiagFactor
d5ada4102b Jean*0055       _RL  firstResidual
d53cbeab8e Jean*0056       _RL  minResidual
d5ada4102b Jean*0057       _RL  lastResidual
                0058       _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059       INTEGER numIters
d53cbeab8e Jean*0060       INTEGER nIterMin
                0061       INTEGER printResidFrq
d5ada4102b Jean*0062       INTEGER myThid
                0063
                0064 C     !LOCAL VARIABLES:
                0065 C     === Local variables ====
d53cbeab8e Jean*0066 C     bi, bj     :: tile indices
d5ada4102b Jean*0067 C     eta_qrN    :: Used in computing search directions
                0068 C     eta_qrNM1     suffix N and NM1 denote current and
                0069 C     cgBeta        previous iterations respectively.
                0070 C     alpha
                0071 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a
                0072 C                   useful debuggin/trouble shooting diagnostic.
                0073 C                   For neumann problems sumRHS needs to be ~0.
                0074 C                   or they converge at a non-zero residual.
d53cbeab8e Jean*0075 C     err        :: Measure of current residual of Ax - b, usually the norm.
d5ada4102b Jean*0076 C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
                0077       INTEGER bi, bj
                0078       INTEGER i, j, it2d
                0079       _RL    err,     errTile(nSx,nSy)
                0080       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
                0081       _RL    eta_qrNM1
                0082       _RL    cgBeta
                0083       _RL    alpha,  alphaTile(nSx,nSy)
                0084       _RL    sumRHS, sumRHStile(nSx,nSy)
                0085       _RL    pW_tmp, pS_tmp
                0086       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0087       LOGICAL printResidual
                0088 CEOP
                0089       _RS  aC2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0090       _RS  pW  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0091       _RS  pS  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0092       _RS  pC  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
11c6dae13e Jean*0093       _RL  q2d(0:sNx+1,0:sNy+1,nSx,nSy)
e938313568 Jean*0094 #ifdef DEBUG_DIAG_CG2D
                0095       CHARACTER*(10) sufx
                0096       _RL  r2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0097 #else
11c6dae13e Jean*0098       _RL  r2d(0:sNx+1,0:sNy+1,nSx,nSy)
e938313568 Jean*0099 #endif
11c6dae13e Jean*0100       _RL  s2d(0:sNx+1,0:sNy+1,nSx,nSy)
d53cbeab8e Jean*0101       _RL  x2dm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
d5ada4102b Jean*0102
ce1941842b Jean*0103 #ifdef ALLOW_DEBUG
                0104       IF (debugMode) CALL DEBUG_ENTER('DIAG_CG2D',myThid)
                0105 #endif
                0106
d5ada4102b Jean*0107 C--   Set matrice main diagnonal:
                0108       DO bj=myByLo(myThid),myByHi(myThid)
                0109        DO bi=myBxLo(myThid),myBxHi(myThid)
ce1941842b Jean*0110 C-    Initialise overlap regions (in case EXCH call do not fill them all)
11c6dae13e Jean*0111         DO j = 1-OLy,sNy+OLy
                0112          DO i = 1-OLx,sNx+OLx
ce1941842b Jean*0113            aC2d(i,j,bi,bj) = 0.
                0114          ENDDO
                0115         ENDDO
d5ada4102b Jean*0116         DO j=1,sNy
                0117          DO i=1,sNx
                0118            aC2d(i,j,bi,bj) = -( ( aW2d(i,j,bi,bj)+aW2d(i+1,j,bi,bj) )
                0119      &                         +( aS2d(i,j,bi,bj)+aS2d(i,j+1,bi,bj) )
                0120      &                        )
                0121          ENDDO
                0122         ENDDO
                0123        ENDDO
                0124       ENDDO
                0125       CALL EXCH_XY_RS(aC2d, myThid)
                0126
                0127 C--   Initialise preconditioner
                0128       DO bj=myByLo(myThid),myByHi(myThid)
                0129        DO bi=myBxLo(myThid),myBxHi(myThid)
                0130         DO j=1,sNy+1
                0131          DO i=1,sNx+1
                0132           IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
                0133            pC(i,j,bi,bj) = 1. _d 0
                0134           ELSE
                0135            pC(i,j,bi,bj) = 1. _d 0 / aC2d(i,j,bi,bj)
                0136           ENDIF
                0137           pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj)
                0138           IF ( pW_tmp .EQ. 0. ) THEN
                0139            pW(i,j,bi,bj) = 0.
                0140           ELSE
a63b8f5615 Jean*0141            pW(i,j,bi,bj) = -aW2d(i,j,bi,bj)*offDiagFactor
                0142      &                     *4. _d 0/( pW_tmp*pW_tmp )
d5ada4102b Jean*0143           ENDIF
                0144           pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
                0145           IF ( pS_tmp .EQ. 0. ) THEN
                0146            pS(i,j,bi,bj) = 0.
                0147           ELSE
a63b8f5615 Jean*0148            pS(i,j,bi,bj) = -aS2d(i,j,bi,bj)*offDiagFactor
                0149      &                     *4. _d 0/( pS_tmp*pS_tmp )
d5ada4102b Jean*0150           ENDIF
e938313568 Jean*0151 c         pC(i,j,bi,bj) = 1.
                0152 c         pW(i,j,bi,bj) = 0.
                0153 c         pS(i,j,bi,bj) = 0.
d5ada4102b Jean*0154          ENDDO
                0155         ENDDO
                0156        ENDDO
                0157       ENDDO
                0158
                0159 C--   Initialise inverter
                0160       eta_qrNM1 = 1. _d 0
                0161
                0162       CALL EXCH_XY_RL( x2d, myThid )
                0163
                0164 C--   Initial residual calculation
                0165       DO bj=myByLo(myThid),myByHi(myThid)
                0166        DO bi=myBxLo(myThid),myBxHi(myThid)
11c6dae13e Jean*0167         DO j=0,sNy+1
                0168          DO i=0,sNx+1
                0169           r2d(i,j,bi,bj) = 0.
d5ada4102b Jean*0170           s2d(i,j,bi,bj) = 0.
d53cbeab8e Jean*0171           x2dm(i,j,bi,bj) = x2d(i,j,bi,bj)
d5ada4102b Jean*0172          ENDDO
                0173         ENDDO
                0174         sumRHStile(bi,bj) = 0. _d 0
                0175         errTile(bi,bj)    = 0. _d 0
                0176         DO j=1,sNy
                0177          DO i=1,sNx
                0178           r2d(i,j,bi,bj) = b2d(i,j,bi,bj) -
                0179      &      (aW2d(i  ,j  ,bi,bj)*x2d(i-1,j  ,bi,bj)
                0180      &      +aW2d(i+1,j  ,bi,bj)*x2d(i+1,j  ,bi,bj)
                0181      &      +aS2d(i  ,j  ,bi,bj)*x2d(i  ,j-1,bi,bj)
                0182      &      +aS2d(i  ,j+1,bi,bj)*x2d(i  ,j+1,bi,bj)
                0183      &      +aC2d(i  ,j  ,bi,bj)*x2d(i  ,j  ,bi,bj)
                0184      &      )
                0185           errTile(bi,bj) = errTile(bi,bj)
                0186      &                  + r2d(i,j,bi,bj)*r2d(i,j,bi,bj)
                0187           sumRHStile(bi,bj) = sumRHStile(bi,bj) + b2d(i,j,bi,bj)
                0188          ENDDO
                0189         ENDDO
                0190        ENDDO
                0191       ENDDO
e938313568 Jean*0192 #ifdef DEBUG_DIAG_CG2D
                0193       CALL EXCH_XY_RL ( r2d, myThid )
                0194 #else
d5ada4102b Jean*0195       CALL EXCH_S3D_RL( r2d, 1, myThid )
e938313568 Jean*0196 #endif
d5ada4102b Jean*0197       CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
e938313568 Jean*0198       IF ( printResidFrq.GE.1 )
                0199      &  CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
d5ada4102b Jean*0200       err = SQRT(err)
d53cbeab8e Jean*0201       it2d = 0
                0202       firstResidual = err
                0203       minResidual   = err
                0204       nIterMin = it2d
d5ada4102b Jean*0205
                0206       printResidual = .FALSE.
                0207       IF ( debugLevel .GE. debLevZero ) THEN
                0208         _BEGIN_MASTER( myThid )
d53cbeab8e Jean*0209         printResidual = printResidFrq.GE.1
e938313568 Jean*0210         IF ( printResidual ) THEN
                0211          WRITE(msgBuf,'(2A,I6,A,1PE17.9,A,1PE14.6)')' diag_cg2d:',
                0212      &    ' iter=', it2d, ' ; resid.=', err, ' ; sumRHS=', sumRHS
                0213          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0214      &                       SQUEEZE_RIGHT, myThid )
                0215         ENDIF
d5ada4102b Jean*0216         _END_MASTER( myThid )
                0217       ENDIF
e938313568 Jean*0218 #ifdef DEBUG_DIAG_CG2D
                0219       IF ( printResidFrq.GE.1 ) THEN
                0220         WRITE(sufx,'(I10.10)') 0
                0221         CALL WRITE_FLD_XY_RL( 'r2d.',sufx, r2d, 1, myThid )
                0222       ENDIF
                0223 #endif
d5ada4102b Jean*0224
                0225       IF ( err .LT. residCriter ) GOTO 11
                0226
                0227 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                0228       DO 10 it2d=1, numIters
                0229
                0230 C--    Solve preconditioning equation and update
                0231 C--    conjugate direction vector "s".
                0232        DO bj=myByLo(myThid),myByHi(myThid)
                0233         DO bi=myBxLo(myThid),myBxHi(myThid)
                0234          eta_qrNtile(bi,bj) = 0. _d 0
                0235          DO j=1,sNy
                0236           DO i=1,sNx
                0237            q2d(i,j,bi,bj) =
                0238      &        pC(i  ,j  ,bi,bj)*r2d(i  ,j  ,bi,bj)
                0239      &       +pW(i  ,j  ,bi,bj)*r2d(i-1,j  ,bi,bj)
                0240      &       +pW(i+1,j  ,bi,bj)*r2d(i+1,j  ,bi,bj)
                0241      &       +pS(i  ,j  ,bi,bj)*r2d(i  ,j-1,bi,bj)
                0242      &       +pS(i  ,j+1,bi,bj)*r2d(i  ,j+1,bi,bj)
                0243            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
                0244      &       +q2d(i,j,bi,bj)*r2d(i,j,bi,bj)
                0245           ENDDO
                0246          ENDDO
                0247         ENDDO
                0248        ENDDO
                0249
                0250        CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
                0251        cgBeta   = eta_qrN/eta_qrNM1
                0252        eta_qrNM1 = eta_qrN
                0253
                0254        DO bj=myByLo(myThid),myByHi(myThid)
                0255         DO bi=myBxLo(myThid),myBxHi(myThid)
                0256          DO j=1,sNy
                0257           DO i=1,sNx
                0258            s2d(i,j,bi,bj) = q2d(i,j,bi,bj)
                0259      &                    + cgBeta*s2d(i,j,bi,bj)
                0260           ENDDO
                0261          ENDDO
                0262         ENDDO
                0263        ENDDO
                0264
                0265 C--    Do exchanges that require messages i.e. between processes.
                0266        CALL EXCH_S3D_RL( s2d, 1, myThid )
                0267
                0268 C==    Evaluate laplace operator on conjugate gradient vector
                0269 C==    q = A.s
                0270        DO bj=myByLo(myThid),myByHi(myThid)
                0271         DO bi=myBxLo(myThid),myBxHi(myThid)
                0272          alphaTile(bi,bj) = 0. _d 0
                0273          DO j=1,sNy
                0274           DO i=1,sNx
                0275            q2d(i,j,bi,bj) =
                0276      &       aW2d(i  ,j  ,bi,bj)*s2d(i-1,j  ,bi,bj)
                0277      &      +aW2d(i+1,j  ,bi,bj)*s2d(i+1,j  ,bi,bj)
                0278      &      +aS2d(i  ,j  ,bi,bj)*s2d(i  ,j-1,bi,bj)
                0279      &      +aS2d(i  ,j+1,bi,bj)*s2d(i  ,j+1,bi,bj)
                0280      &      +aC2d(i  ,j  ,bi,bj)*s2d(i  ,j  ,bi,bj)
                0281            alphaTile(bi,bj) = alphaTile(bi,bj)
                0282      &                      + s2d(i,j,bi,bj)*q2d(i,j,bi,bj)
                0283           ENDDO
                0284          ENDDO
                0285         ENDDO
                0286        ENDDO
                0287        CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
                0288        alpha = eta_qrN/alpha
                0289
                0290 C==    Update solution and residual vectors
                0291 C      Now compute "interior" points.
                0292        DO bj=myByLo(myThid),myByHi(myThid)
                0293         DO bi=myBxLo(myThid),myBxHi(myThid)
                0294          errTile(bi,bj) = 0. _d 0
                0295          DO j=1,sNy
                0296           DO i=1,sNx
                0297            x2d(i,j,bi,bj)=x2d(i,j,bi,bj)+alpha*s2d(i,j,bi,bj)
                0298            r2d(i,j,bi,bj)=r2d(i,j,bi,bj)-alpha*q2d(i,j,bi,bj)
                0299            errTile(bi,bj) = errTile(bi,bj)
                0300      &                    + r2d(i,j,bi,bj)*r2d(i,j,bi,bj)
                0301           ENDDO
                0302          ENDDO
                0303         ENDDO
                0304        ENDDO
                0305
                0306        CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
                0307        err = SQRT(err)
                0308        IF ( printResidual ) THEN
d53cbeab8e Jean*0309         IF ( MOD( it2d-1, printResidFrq ).EQ.0 ) THEN
e938313568 Jean*0310          WRITE(msgBuf,'(A,I6,A,1PE17.9)')
                0311      &    ' diag_cg2d: iter=', it2d, ' ; resid.=', err
d5ada4102b Jean*0312          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0313      &                       SQUEEZE_RIGHT, myThid )
                0314         ENDIF
                0315        ENDIF
                0316        IF ( err .LT. residCriter ) GOTO 11
d53cbeab8e Jean*0317        IF ( err .LT. minResidual ) THEN
                0318 C-     Store lowest residual solution
                0319          minResidual = err
                0320          nIterMin = it2d
                0321          DO bj=myByLo(myThid),myByHi(myThid)
                0322           DO bi=myBxLo(myThid),myBxHi(myThid)
                0323            DO j=1,sNy
                0324             DO i=1,sNx
                0325              x2dm(i,j,bi,bj) = x2d(i,j,bi,bj)
                0326             ENDDO
                0327            ENDDO
                0328           ENDDO
                0329          ENDDO
                0330        ENDIF
d5ada4102b Jean*0331
e938313568 Jean*0332 #ifdef DEBUG_DIAG_CG2D
                0333        CALL EXCH_XY_RL( r2d, myThid )
                0334        IF ( printResidFrq.GE.1 ) THEN
                0335         WRITE(sufx,'(I10.10)') it2d
                0336         CALL WRITE_FLD_XY_RL( 'r2d.',sufx, r2d, 1, myThid )
                0337         CALL WRITE_FLD_XY_RL( 'x2d.',sufx, x2d, 1, myThid )
                0338        ENDIF
                0339 #else
d5ada4102b Jean*0340        CALL EXCH_S3D_RL( r2d, 1, myThid )
e938313568 Jean*0341 #endif
d5ada4102b Jean*0342
                0343    10 CONTINUE
d53cbeab8e Jean*0344       it2d = numIters
d5ada4102b Jean*0345    11 CONTINUE
                0346
                0347 C--   Return parameters to caller
d53cbeab8e Jean*0348       lastResidual = err
                0349       numIters = it2d
                0350
                0351       IF ( err .GT. minResidual ) THEN
                0352 C-    use the lowest residual solution (instead of current one <-> last residual)
                0353         DO bj=myByLo(myThid),myByHi(myThid)
                0354          DO bi=myBxLo(myThid),myBxHi(myThid)
                0355           DO j=1,sNy
                0356            DO i=1,sNx
                0357              x2d(i,j,bi,bj) = x2dm(i,j,bi,bj)
                0358            ENDDO
                0359           ENDDO
                0360          ENDDO
                0361         ENDDO
                0362       ENDIF
                0363 c     CALL EXCH_XY_RL( x2d, myThid )
d5ada4102b Jean*0364
ce1941842b Jean*0365 #ifdef ALLOW_DEBUG
                0366       IF (debugMode) CALL DEBUG_LEAVE('DIAG_CG2D',myThid)
                0367 #endif
                0368
d5ada4102b Jean*0369       RETURN
                0370       END