C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diag_cg2d.F,v 1.2 2011/07/22 19:53:39 jmc Exp $ C $Name: $ #include "DIAG_OPTIONS.h" CBOP C !ROUTINE: DIAG_CG2D C !INTERFACE: SUBROUTINE DIAG_CG2D( I aW2d, aS2d, b2d, I residCriter, O firstResidual, minResidual, lastResidual, U x2d, numIters, O nIterMin, I printResidFrq, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE CG2D C | o Two-dimensional grid problem conjugate-gradient C | inverter (with preconditioner). C *==========================================================* C | Con. grad is an iterative procedure for solving Ax = b. C | It requires the A be symmetric. C | This implementation assumes A is a five-diagonal C | matrix of the form that arises in the discrete C | representation of the del^2 operator in a C | two-dimensional space. C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global data === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" c#include "CG2D.h" c#include "GRID.h" c#include "SURFACE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C b2d :: The source term or "right hand side" C x2d :: The solution C firstResidual :: the initial residual before any iterations C minResidual :: the lowest residual reached C lastResidual :: the actual residual reached C numIters :: Entry: the maximum number of iterations allowed C Exit: the actual number of iterations used C nIterMin :: iteration number corresponding to lowest residual C printResidFrq :: Frequency for printing residual in CG iterations C myThid :: my Thread Id number _RS aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL residCriter _RL firstResidual _RL minResidual _RL lastResidual _RL x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) INTEGER numIters INTEGER nIterMin INTEGER printResidFrq INTEGER myThid C !LOCAL VARIABLES: C === Local variables ==== C bi, bj :: tile indices C eta_qrN :: Used in computing search directions C eta_qrNM1 suffix N and NM1 denote current and C cgBeta previous iterations respectively. C alpha C sumRHS :: Sum of right-hand-side. Sometimes this is a C useful debuggin/trouble shooting diagnostic. C For neumann problems sumRHS needs to be ~0. C or they converge at a non-zero residual. C err :: Measure of current residual of Ax - b, usually the norm. C i, j, it2d :: Loop counters ( it2d counts CG iterations ) INTEGER bi, bj INTEGER i, j, it2d _RL err, errTile(nSx,nSy) _RL eta_qrN, eta_qrNtile(nSx,nSy) _RL eta_qrNM1 _RL cgBeta _RL alpha, alphaTile(nSx,nSy) _RL sumRHS, sumRHStile(nSx,nSy) _RL pW_tmp, pS_tmp CHARACTER*(MAX_LEN_MBUF) msgBuf LOGICAL printResidual CEOP _RS aC2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS pW (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS pS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS pC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL q2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy) _RL r2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy) _RL s2d(1-1:sNx+1,1-1:sNy+1,nSx,nSy) _RL x2dm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C-- Set matrice main diagnonal: DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx aC2d(i,j,bi,bj) = -( ( aW2d(i,j,bi,bj)+aW2d(i+1,j,bi,bj) ) & +( aS2d(i,j,bi,bj)+aS2d(i,j+1,bi,bj) ) & ) ENDDO ENDDO ENDDO ENDDO CALL EXCH_XY_RS(aC2d, myThid) C-- Initialise preconditioner DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy+1 DO i=1,sNx+1 IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN pC(i,j,bi,bj) = 1. _d 0 ELSE pC(i,j,bi,bj) = 1. _d 0 / aC2d(i,j,bi,bj) ENDIF pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj) IF ( pW_tmp .EQ. 0. ) THEN pW(i,j,bi,bj) = 0. ELSE pW(i,j,bi,bj) = -aW2d(i,j,bi,bj) & /( (cg2dpcOffDFac*pW_tmp)**2 ) ENDIF pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj) IF ( pS_tmp .EQ. 0. ) THEN pS(i,j,bi,bj) = 0. ELSE pS(i,j,bi,bj) = -aS2d(i,j,bi,bj) & /( (cg2dpcOffDFac*pS_tmp)**2 ) ENDIF ENDDO ENDDO ENDDO ENDDO C-- Initialise inverter eta_qrNM1 = 1. _d 0 CALL EXCH_XY_RL( x2d, myThid ) C-- Initial residual calculation DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-1,sNy+1 DO i=1-1,sNx+1 s2d(i,j,bi,bj) = 0. x2dm(i,j,bi,bj) = x2d(i,j,bi,bj) ENDDO ENDDO sumRHStile(bi,bj) = 0. _d 0 errTile(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx r2d(i,j,bi,bj) = b2d(i,j,bi,bj) - & (aW2d(i ,j ,bi,bj)*x2d(i-1,j ,bi,bj) & +aW2d(i+1,j ,bi,bj)*x2d(i+1,j ,bi,bj) & +aS2d(i ,j ,bi,bj)*x2d(i ,j-1,bi,bj) & +aS2d(i ,j+1,bi,bj)*x2d(i ,j+1,bi,bj) & +aC2d(i ,j ,bi,bj)*x2d(i ,j ,bi,bj) & ) errTile(bi,bj) = errTile(bi,bj) & + r2d(i,j,bi,bj)*r2d(i,j,bi,bj) sumRHStile(bi,bj) = sumRHStile(bi,bj) + b2d(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO CALL EXCH_S3D_RL( r2d, 1, myThid ) CALL GLOBAL_SUM_TILE_RL( errTile, err, myThid ) CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid ) err = SQRT(err) it2d = 0 firstResidual = err minResidual = err nIterMin = it2d printResidual = .FALSE. IF ( debugLevel .GE. debLevZero ) THEN _BEGIN_MASTER( myThid ) printResidual = printResidFrq.GE.1 c WRITE(standardmessageunit,'(A,1P2E22.14)') c & ' diag_cg2d: Sum(rhs),rhsMax = ', sumRHS, rhsMax _END_MASTER( myThid ) ENDIF IF ( err .LT. residCriter ) GOTO 11 C >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DO 10 it2d=1, numIters C-- Solve preconditioning equation and update C-- conjugate direction vector "s". DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) eta_qrNtile(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx q2d(i,j,bi,bj) = & pC(i ,j ,bi,bj)*r2d(i ,j ,bi,bj) & +pW(i ,j ,bi,bj)*r2d(i-1,j ,bi,bj) & +pW(i+1,j ,bi,bj)*r2d(i+1,j ,bi,bj) & +pS(i ,j ,bi,bj)*r2d(i ,j-1,bi,bj) & +pS(i ,j+1,bi,bj)*r2d(i ,j+1,bi,bj) eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj) & +q2d(i,j,bi,bj)*r2d(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid ) cgBeta = eta_qrN/eta_qrNM1 eta_qrNM1 = eta_qrN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx s2d(i,j,bi,bj) = q2d(i,j,bi,bj) & + cgBeta*s2d(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO C-- Do exchanges that require messages i.e. between processes. CALL EXCH_S3D_RL( s2d, 1, myThid ) C== Evaluate laplace operator on conjugate gradient vector C== q = A.s DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) alphaTile(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx q2d(i,j,bi,bj) = & aW2d(i ,j ,bi,bj)*s2d(i-1,j ,bi,bj) & +aW2d(i+1,j ,bi,bj)*s2d(i+1,j ,bi,bj) & +aS2d(i ,j ,bi,bj)*s2d(i ,j-1,bi,bj) & +aS2d(i ,j+1,bi,bj)*s2d(i ,j+1,bi,bj) & +aC2d(i ,j ,bi,bj)*s2d(i ,j ,bi,bj) alphaTile(bi,bj) = alphaTile(bi,bj) & + s2d(i,j,bi,bj)*q2d(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( alphaTile, alpha, myThid ) alpha = eta_qrN/alpha C== Update solution and residual vectors C Now compute "interior" points. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) errTile(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx x2d(i,j,bi,bj)=x2d(i,j,bi,bj)+alpha*s2d(i,j,bi,bj) r2d(i,j,bi,bj)=r2d(i,j,bi,bj)-alpha*q2d(i,j,bi,bj) errTile(bi,bj) = errTile(bi,bj) & + r2d(i,j,bi,bj)*r2d(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( errTile, err, myThid ) err = SQRT(err) IF ( printResidual ) THEN IF ( MOD( it2d-1, printResidFrq ).EQ.0 ) THEN WRITE(msgBuf,'(A,I6,A,1PE21.14)') & ' diag_cg2d: iter=', it2d, ' ; resid.= ', err CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) ENDIF ENDIF IF ( err .LT. residCriter ) GOTO 11 IF ( err .LT. minResidual ) THEN C- Store lowest residual solution minResidual = err nIterMin = it2d DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx x2dm(i,j,bi,bj) = x2d(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDIF CALL EXCH_S3D_RL( r2d, 1, myThid ) 10 CONTINUE it2d = numIters 11 CONTINUE C-- Return parameters to caller lastResidual = err numIters = it2d IF ( err .GT. minResidual ) THEN C- use the lowest residual solution (instead of current one <-> last residual) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx x2d(i,j,bi,bj) = x2dm(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDIF c CALL EXCH_XY_RL( x2d, myThid ) RETURN END