pkg/exch2/w2_set_f2f_index.F

017b6b2289 Jean*0001 #include "CPP_EEOPTIONS.h"
                0002 #include "W2_OPTIONS.h"
                0003
                0004 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0005 CBOP 0
                0006 C !ROUTINE: W2_SET_F2F_INDEX
                0007
                0008 C !INTERFACE:
                0009       SUBROUTINE W2_SET_F2F_INDEX( myThid )
                0010
                0011 C     !DESCRIPTION:
                0012 C     Set-up index correspondence matrix for connected Facet-Edges
                0013
                0014 C     !USES:
                0015       IMPLICIT NONE
                0016
d6ea3164dc Jean*0017 C      Tile topology settings data structures
017b6b2289 Jean*0018 #include "SIZE.h"
                0019 #include "EEPARAMS.h"
                0020 #include "W2_EXCH2_SIZE.h"
                0021 #include "W2_EXCH2_PARAMS.h"
                0022 #include "W2_EXCH2_TOPOLOGY.h"
                0023
                0024 C     !INPUT PARAMETERS:
                0025 C     myThid  :: my Thread Id number
                0026 C               (Note: not relevant since threading has not yet started)
                0027       INTEGER myThid
                0028
                0029 C     !LOCAL VARIABLES:
                0030 C     === Local variables ===
d6ea3164dc Jean*0031 C     msgBuf     :: Informational/error message buffer
017b6b2289 Jean*0032       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0033       CHARACTER*1 edge(4)
                0034       INTEGER i, j, ii, jj, i1, j1, k, lo, ll
                0035       INTEGER errCnt
                0036       INTEGER chk1, chk2, chk3, chk4, chk5, chk6
                0037       LOGICAL prtFlag
                0038 CEOP
                0039       DATA edge / 'N' , 'S' , 'E' , 'W' /
                0040
                0041       WRITE(msgBuf,'(2A)') 'W2_SET_F2F_INDEX:',
                0042      &       ' index matrix for connected Facet-Edges:'
                0043       CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
                0044       prtFlag = ABS(W2_printMsg).GE.2
                0045      &       .OR. ( W2_printMsg .NE.0 .AND. myProcId.EQ.0 )
                0046
                0047 C--   Check Topology
                0048       errCnt = 0
                0049       DO j=1,nFacets
                0050        DO i=1,4
                0051 C-    connected to:
                0052          jj = INT(facet_link(i,j))
                0053          ii = MOD( NINT(facet_link(i,j)*10.), 10 )
                0054          IF ( facet_link(i,j).EQ. 0. ) THEN
                0055           WRITE(msgBuf,'(3A,I3,A,F6.2,A)')
                0056      &           '** WARNING ** ',  edge(i), '.Edge of facet #',
                0057      &        j, ' disconnected (facet_link=',facet_link(i,j),')'
                0058           CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
                0059           CALL PRINT_MESSAGE( msgBuf,errorMessageUnit,SQUEEZE_RIGHT,1 )
                0060          ELSEIF ( jj.LT.1 .OR. jj.GT.nFacets
                0061      &            .OR. ii.LT.1 .OR. ii.GT.4 ) THEN
                0062           WRITE(msgBuf,'(2A,I3,A,F6.2,A)') edge(i), '.Edge of facet #',
                0063      &      j, ' : bad connection (facet_link=',facet_link(i,j),')'
                0064           CALL PRINT_ERROR( msgBuf, myThid )
                0065           errCnt = errCnt + 1
                0066          ELSE
                0067 C-    check connection is set-up both ways:
                0068           j1 = INT(facet_link(ii,jj))
                0069           i1 = MOD( NINT(facet_link(ii,jj)*10.), 10 )
                0070           IF ( j1.NE.j .OR. i1.NE.i ) THEN
                0071            WRITE(msgBuf,'(2(2A,I3,A),F5.2,A)')
                0072      &      edge(i), '.Edge facet #', j,' connect to: ',
                0073      &      edge(ii),'.Edge facet #',jj,' (',facet_link(i,j),' )'
                0074            CALL PRINT_ERROR( msgBuf, myThid )
                0075           IF ( i1.GE.1 .AND. i1.LE.4 ) THEN
                0076            WRITE(msgBuf,'(A,2(2A,I3,A),F5.2,A)') ' but ',
                0077      &      edge(ii),'.Edge facet #',jj,' connect to: ',
                0078      &      edge(i1),'.Edge facet #',j1,' (',facet_link(ii,jj),' )'
                0079           ELSE
                0080            WRITE(msgBuf,'(A,2(2A,I3,A),F5.2,A)') ' but ',
                0081      &      edge(ii),'.Edge facet #',jj,' connect to: ',
                0082      &           '?' ,'.Edge facet #',j1,' (',facet_link(ii,jj),' )'
                0083           ENDIF
                0084            CALL PRINT_ERROR( msgBuf, myThid )
                0085            errCnt = errCnt + 1
                0086           ENDIF
                0087          ENDIF
                0088        ENDDO
                0089       ENDDO
                0090       IF ( errCnt.GT.0 ) THEN
                0091         WRITE(msgBuf,'(A,I3,A)')
                0092      &   ' W2_SET_F2F_INDEX: found', errCnt, ' Topology errors'
                0093         CALL PRINT_ERROR( msgBuf, myThid )
                0094         STOP 'ABNORMAL END: S/R W2_SET_F2F_INDEX'
                0095       ENDIF
                0096
                0097 C--   Check length of connection (facet size) between connected facet-edges
                0098       errCnt = 0
                0099       DO j=1,nFacets
                0100        DO i=1,4
                0101 C-    Length of N or S Edge = x-size, E or W Edge = y-size
                0102         lo = 2*(j-1) + (i+1)/2
                0103         lo = facet_dims( lo )
                0104 C-    connected to:
                0105         jj = INT(facet_link(i,j))
                0106         ii = MOD( NINT(facet_link(i,j)*10.), 10 )
                0107         IF ( jj.GE.1 ) THEN
                0108          ll = 2*(jj-1)+(ii+1)/2
                0109          ll = facet_dims( ll )
                0110          IF ( lo.NE.ll ) THEN
                0111            WRITE(msgBuf,'(3A,I3,A,I8,A)') 'Length of connection: ',
                0112      &      edge(i), '.Edge facet #', j,  ' (=',lo,')'
                0113            CALL PRINT_ERROR( msgBuf, myThid )
                0114            WRITE(msgBuf,'(3A,I3,A,I8,A)') '  to: ',
                0115      &      edge(ii),'.Edge facet #', jj, ' (=',ll,') are different'
                0116            CALL PRINT_ERROR( msgBuf, myThid )
                0117            errCnt = errCnt + 1
                0118          ENDIF
                0119 C--   Set indices correspondence matrix (facet_pij):
                0120 C     suffix "so" for indices of source facet j ;
                0121 C     suffix "tg" for indices of target facet jj= INT(facet_link(i,j))
                0122 C     pij(:,i,j) : matrix which gives so indices when applied to tg indices
                0123 C        iso = pij(1)*itg + pij(2)*jtg + oi
                0124 C        jso = pij(3)*itg + pij(4)*jtg + oj
                0125
                0126 C-    Default: Same orientation (works for N <-> S or E <-> W)
                0127          facet_pij(1,i,j) = 1
                0128          facet_pij(2,i,j) = 0
                0129          facet_pij(3,i,j) = 0
                0130          facet_pij(4,i,j) = 1
                0131 C--   1rst: cases with same orientation
                0132          IF ( i.EQ.1 .AND. ii.EQ.2 ) THEN
                0133 C-    N <-- S ; [:,jHi]_so <-- [:,0]_tg
                0134            facet_oi(i,j) = 0
                0135            facet_oj(i,j) = +facet_dims(2*j)
                0136          ELSEIF ( i.EQ.2 .AND. ii.EQ.1 ) THEN
                0137 C-    S <-- N ; [:, 1 ]_so <-- [:,1+jHi]_tg
                0138            facet_oi(i,j) = 0
                0139            facet_oj(i,j) = -facet_dims(2*jj)
                0140          ELSEIF ( i.EQ.3 .AND. ii.EQ.4 ) THEN
                0141 C-    E <-- W ; [iHi,:]_so <-- [0,:]_tg
                0142            facet_oi(i,j) = +facet_dims(2*j-1)
                0143            facet_oj(i,j) = 0
                0144          ELSEIF ( i.EQ.4 .AND. ii.EQ.3 ) THEN
                0145 C-    W <-- E  (i=4 & ii=3); [ 1 ,:]_so <-- [iHi,:]_tg
                0146            facet_oi(i,j) = -facet_dims(2*jj-1)
                0147            facet_oj(i,j) = 0
                0148 C--   2nd : cases where orientation changes
                0149          ELSEIF ( i.EQ.1 .AND. ii.EQ.4 ) THEN
                0150 C-    N <-- W ; mapping W.face target indices to N.face source indices
                0151 C-         [:,jHi]_so <-- [0,:]_tg
                0152            facet_pij(1,i,j) = 0
                0153            facet_pij(2,i,j) =-1
                0154            facet_pij(3,i,j) = 1
                0155            facet_pij(4,i,j) = 0
                0156            facet_oi(i,j) = lo+1
                0157            facet_oj(i,j) = +facet_dims(2*j)
                0158          ELSEIF ( i.EQ.2 .AND. ii.EQ.3 ) THEN
                0159 C-    S <-- E ; mapping E.face target indices to S.face source indices
                0160 C-         [:,1]_so <-- [1+iHi,:]_tg
                0161            facet_pij(1,i,j) = 0
                0162            facet_pij(2,i,j) =-1
                0163            facet_pij(3,i,j) = 1
                0164            facet_pij(4,i,j) = 0
                0165            facet_oi(i,j) = lo+1
                0166            facet_oj(i,j) = -facet_dims(2*jj-1)
                0167          ELSEIF ( i.EQ.3 .AND. ii.EQ.2 ) THEN
                0168 C-    E <-- S ; mapping S.face target indices to E.face source indices
                0169 C-         [iHi,:]_so <-- [:,0]_tg
                0170            facet_pij(1,i,j) = 0
                0171            facet_pij(2,i,j) = 1
                0172            facet_pij(3,i,j) =-1
                0173            facet_pij(4,i,j) = 0
                0174            facet_oi(i,j) = +facet_dims(2*j-1)
                0175            facet_oj(i,j) = lo+1
                0176          ELSEIF ( i.EQ.4 .AND. ii.EQ.1 ) THEN
                0177 C-    W <-- N ; mapping N.face target indices to W.face source indices
                0178 C-         [ 1 ,:]_so <-- [:,1+jHi]_tg
                0179            facet_pij(1,i,j) = 0
                0180            facet_pij(2,i,j) = 1
                0181            facet_pij(3,i,j) =-1
                0182            facet_pij(4,i,j) = 0
                0183            facet_oi(i,j) = -facet_dims(2*jj)
                0184            facet_oj(i,j) = lo+1
                0185          ELSE
                0186            WRITE(msgBuf,'(2(3A,I3),A)') ' connect ',
                0187      &      edge(i), '.Edge (facet#', j, ' ) to: ',
                0188      &      edge(ii),'.Edge (facet#', jj,' )'
                0189            CALL PRINT_ERROR( msgBuf, myThid )
                0190            WRITE(msgBuf,'(A)')
                0191      &       ' W2_SET_F2F_INDEX: Connection not supported'
                0192            CALL PRINT_ERROR( msgBuf, myThid )
                0193            errCnt = errCnt + 1
                0194          ENDIF
                0195 C--   Print resulting index matrix:
                0196          IF ( prtFlag ) WRITE(W2_oUnit,'(2(3A,I3),A,4I3,A,2I6)')
                0197      &    '  ', edge(i), '.Edge Facet', j, ' <-- ',
                0198      &          edge(ii),'.Edge Facet', jj,
                0199      &    ' : pij=', (facet_pij(k,i,j),k=1,4),
                0200      &    ' ; oi,oj=', facet_oi(i,j), facet_oj(i,j)
                0201         ENDIF
                0202        ENDDO
                0203       ENDDO
                0204       IF ( errCnt.GT.0 ) THEN
                0205         WRITE(msgBuf,'(A,I3,A)')
                0206      &   ' W2_SET_F2F_INDEX: found', errCnt, ' Connection errors'
                0207         CALL PRINT_ERROR( msgBuf, myThid )
                0208         STOP 'ABNORMAL END: S/R W2_SET_F2F_INDEX'
                0209       ENDIF
                0210
                0211 C--   Check indices correspondence matrix reciprocity :
                0212 C     This is not necessary (if code above is right) ; However:
                0213 C     a) reciprocity is used later-on => good to check;
                0214 C     b) easy to add an error in setting offset => worth to check.
                0215       errCnt = 0
                0216       DO j=1,nFacets
                0217        DO i=1,4
                0218 C-    connected to:
                0219         jj = INT(facet_link(i,j))
                0220         ii = MOD( NINT(facet_link(i,j)*10.), 10 )
                0221         IF ( jj.GE.1 ) THEN
                0222 C        iso = pij(1)*( Pij(1)*iso + Pij(2)*jso + Oi )
                0223 C            + pij(2)*( Pij(3)*iso + Pij(4)*jso + Oj )
                0224 C            + oi
                0225 C        jso = pij(3)*( Pij(1)*iso + Pij(2)*jso + Oi )
                0226 C            + pij(4)*( Pij(3)*iso + Pij(4)*jso + Oj )
                0227 C            + oj
                0228 C-      Matrix product:
                0229          chk1 = facet_pij(1,i,j)*facet_pij(1,ii,jj)
                0230      &        + facet_pij(2,i,j)*facet_pij(3,ii,jj)
                0231          chk2 = facet_pij(1,i,j)*facet_pij(2,ii,jj)
                0232      &        + facet_pij(2,i,j)*facet_pij(4,ii,jj)
                0233          chk3 = facet_pij(3,i,j)*facet_pij(1,ii,jj)
                0234      &        + facet_pij(4,i,j)*facet_pij(3,ii,jj)
                0235          chk4 = facet_pij(3,i,j)*facet_pij(2,ii,jj)
                0236      &        + facet_pij(4,i,j)*facet_pij(4,ii,jj)
                0237 C-      Offsets:
                0238          chk5 = facet_pij(1,i,j)*facet_oi(ii,jj)
                0239      &        + facet_pij(2,i,j)*facet_oj(ii,jj)
                0240      &        + facet_oi(i,j)
                0241          chk6 = facet_pij(3,i,j)*facet_oi(ii,jj)
                0242      &        + facet_pij(4,i,j)*facet_oj(ii,jj)
                0243      &        + facet_oj(i,j)
                0244          IF ( chk1.NE.1 .OR. chk2.NE.0 .OR. chk5.NE.0 .OR.
                0245      &        chk3.NE.0 .OR. chk4.NE.1 .OR. chk6.NE.0 ) THEN
                0246            WRITE(msgBuf,'(2(3A,I3),A)') ' connect ',
                0247      &      edge(i), '.Edge (facet#', j, ' ) to: ',
                0248      &      edge(ii),'.Edge (facet#', jj,' )'
                0249            CALL PRINT_ERROR( msgBuf, myThid )
                0250            WRITE(msgBuf,'(A,4I4,2I8)')
                0251      &      ' Bug in Matrix Product:',chk1,chk2,chk3,chk4,chk5,chk6
                0252            CALL PRINT_ERROR( msgBuf, myThid )
                0253            errCnt = errCnt + 1
                0254          ENDIF
                0255         ENDIF
                0256        ENDDO
                0257       ENDDO
                0258       IF ( errCnt.GT.0 ) THEN
                0259         WRITE(msgBuf,'(A,I3,A)')
                0260      &   ' W2_SET_F2F_INDEX: found', errCnt, ' bugs in Matrix product'
                0261         CALL PRINT_ERROR( msgBuf, myThid )
                0262         STOP 'ABNORMAL END: S/R W2_SET_F2F_INDEX'
                0263       ENDIF
                0264
                0265       RETURN
                0266       END