C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_cg2d.F,v 1.42 2003/10/09 04:19:18 edhill Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: INI_CG2D C !INTERFACE: SUBROUTINE INI_CG2D( myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE INI_CG2D C | o Initialise 2d conjugate gradient solver operators. C *==========================================================* C | These arrays are purely a function of the basin geom. C | We set then here once and them use then repeatedly. C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" c#include "DYNVARS.h" #include "SURFACE.h" #include "CG2D.h" #ifdef ALLOW_OBCS #include "OBCS.h" #endif C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid - Thread no. that called this routine. INTEGER myThid C !LOCAL VARIABLES: C === Local variables === C xG, yG - Global coordinate location. C zG C iG, jG - Global coordinate index C bi,bj - Loop counters C faceArea - Temporary used to hold cell face areas. C I,J,K C myNorm - Work variable used in calculating normalisation factor C sumArea - Work variable used to compute the total Domain Area CHARACTER*(MAX_LEN_MBUF) msgBuf INTEGER bi, bj INTEGER I, J, K _RL faceArea _RS myNorm _RL sumArea _RL aC, aCw, aCs CEOP C-- Initialize arrays in common blocs (CG2D.h) ; not really necessary C but safer when EXCH do not fill all the overlap regions. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx aW2d(I,J,bi,bj) = 0. _d 0 aS2d(I,J,bi,bj) = 0. _d 0 pW(I,J,bi,bj) = 0. _d 0 pS(I,J,bi,bj) = 0. _d 0 pC(I,J,bi,bj) = 0. _d 0 cg2d_q(I,J,bi,bj) = 0. _d 0 ENDDO ENDDO DO J=1-1,sNy+1 DO I=1-1,sNx+1 cg2d_r(I,J,bi,bj) = 0. _d 0 cg2d_s(I,J,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO C-- Initialise laplace operator C aW2d: integral in Z Ax/dX C aS2d: integral in Z Ay/dY myNorm = 0. _d 0 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx aW2d(I,J,bi,bj) = 0. _d 0 aS2d(I,J,bi,bj) = 0. _d 0 ENDDO ENDDO DO K=1,Nr DO J=1,sNy DO I=1,sNx faceArea = _dyG(I,J,bi,bj)*drF(K) & *_hFacW(I,J,K,bi,bj) aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) + & implicSurfPress*implicDiv2DFlow & *faceArea*recip_dxC(I,J,bi,bj) faceArea = _dxG(I,J,bi,bj)*drF(K) & *_hFacS(I,J,K,bi,bj) aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) + & implicSurfPress*implicDiv2DFlow & *faceArea*recip_dyC(I,J,bi,bj) ENDDO ENDDO ENDDO #ifdef ALLOW_OBCS IF (useOBCS) THEN DO I=1,sNx IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0. IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0. IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0. IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0. ENDDO DO J=1,sNy IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0. IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0. IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0. IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0. ENDDO ENDIF #endif DO J=1,sNy DO I=1,sNx myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm) myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm) ENDDO ENDDO ENDDO ENDDO _GLOBAL_MAX_R4( myNorm, myThid ) IF ( myNorm .NE. 0. _d 0 ) THEN myNorm = 1. _d 0/myNorm ELSE myNorm = 1. _d 0 ENDIF cg2dNorm = myNorm _BEGIN_MASTER( myThid ) CcnhDebugStarts WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ', & cg2dNorm CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) WRITE(msgBuf,*) ' ' CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) CcnhDebugEnds _END_MASTER( myThid ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm ENDDO ENDDO ENDDO ENDDO C-- Update overlap regions CcnhDebugStarts C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid ) C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid ) CcnhDebugEnds c _EXCH_XY_R4(aW2d, myThid) c _EXCH_XY_R4(aS2d, myThid) CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid) CcnhDebugStarts C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid ) C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid ) CcnhDebugEnds C-- Define the solver tolerance in the appropriate Unit : cg2dNormaliseRHS = cg2dTargetResWunit.LE.0 IF (cg2dNormaliseRHS) THEN C- when using a normalisation of RHS, tolerance has no unit => no conversion cg2dTolerance = cg2dTargetResidual ELSE C- compute the total Area of the domain : sumArea = 0. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN sumArea = sumArea + rA(i,j,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO c WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea _GLOBAL_SUM_R8( sumArea, myThid ) C- convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2] cg2dTolerance = cg2dNorm * cg2dTargetResWunit & * sumArea / deltaTmom _BEGIN_MASTER( myThid ) WRITE(standardMessageUnit,'(2A,1P2E22.14)') ' ini_cg2d: ', & 'sumArea,cg2dTolerance =', sumArea,cg2dTolerance _END_MASTER( myThid ) ENDIF C-- Initialise preconditioner C Note. 20th May 1998 C I made a weird discovery! In the model paper we argue C for the form of the preconditioner used here ( see C A Finite-volume, Incompressible Navier-Stokes Model C ...., Marshall et. al ). The algebra gives a simple C 0.5 factor for the averaging of ac and aCw to get a C symmettric pre-conditioner. By using a factor of 0.51 C i.e. scaling the off-diagonal terms in the C preconditioner down slightly I managed to get the C number of iterations for convergence in a test case to C drop form 192 -> 134! Need to investigate this further! C For now I have introduced a parameter cg2dpcOffDFac which C defaults to 0.51 but can be set at runtime. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx pC(I,J,bi,bj) = 1. _d 0 aC = -( & aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj) & +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj) & +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)* & rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf & ) aCs = -( & aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj) & +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj) & +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)* & rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf & ) aCw = -( & aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj) & +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj) & +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)* & rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf & ) IF ( aC .EQ. 0. ) THEN pC(I,J,bi,bj) = 1. _d 0 ELSE pC(I,J,bi,bj) = 1. _d 0 / aC ENDIF IF ( aC + aCw .EQ. 0. ) THEN pW(I,J,bi,bj) = 0. ELSE pW(I,J,bi,bj) = & -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 ) ENDIF IF ( aC + aCs .EQ. 0. ) THEN pS(I,J,bi,bj) = 0. ELSE pS(I,J,bi,bj) = & -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 ) ENDIF C pC(I,J,bi,bj) = 1. C pW(I,J,bi,bj) = 0. C pS(I,J,bi,bj) = 0. ENDDO ENDDO ENDDO ENDDO C-- Update overlap regions _EXCH_XY_R4(pC, myThid) c _EXCH_XY_R4(pW, myThid) c _EXCH_XY_R4(pS, myThid) CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid) CcnhDebugStarts C CALL PLOT_FIELD_XYRS( pC, 'pC INI_CG2D.2' , 1, myThid ) C CALL PLOT_FIELD_XYRS( pW, 'pW INI_CG2D.2' , 1, myThid ) C CALL PLOT_FIELD_XYRS( pS, 'pS INI_CG2D.2' , 1, myThid ) CcnhDebugEnds RETURN END