summaryrefslogtreecommitdiff
path: root/Trivac/src/TRICH3.f
diff options
context:
space:
mode:
authorstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
committerstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
commit7dfcc480ba1e19bd3232349fc733caef94034292 (patch)
tree03ee104eb8846d5cc1a981d267687a729185d3f3 /Trivac/src/TRICH3.f
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Trivac/src/TRICH3.f')
-rwxr-xr-xTrivac/src/TRICH3.f257
1 files changed, 257 insertions, 0 deletions
diff --git a/Trivac/src/TRICH3.f b/Trivac/src/TRICH3.f
new file mode 100755
index 0000000..e96e2dd
--- /dev/null
+++ b/Trivac/src/TRICH3.f
@@ -0,0 +1,257 @@
+*DECK TRICH3
+ SUBROUTINE TRICH3(ISPLH,IPTRK,LX,LZ,L4,MAT,KN,MUW,MUX,MUY,MUZ,
+ 1 IPW,IPX,IPY,IPZ,IMPX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the compressed diagonal storage indices (MUW, MUX, MUY and
+* MUZ) an the permutation vectors (IPW, IPX, IPY and IPZ) for an ADI
+* splitting of a mesh corner finite difference discretization in
+* hexagonal geometry.
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Benaboud
+*
+*Parameters: input
+* ISPLH type of mesh-splitting: =1 for complete hexagons; =2 for
+* triangular mesh-splitting.
+* IPTRK L_TRACK pointer to the Trivac tracking information.
+* LX number of hexagons in a plane.
+* LZ number of axial planes.
+* L4 order of system matrices.
+* MAT mixture index assigned to each element.
+* KN element-ordered unknown list. Dimensionned to LL*LX*LZ
+* where LL=12 (hexagons) or 14 (triangles).
+* IMPX print parameter (equal to zero for no print).
+*
+*Parameters: output
+* MUW W-oriented compressed storage mode indices.
+* MUX X-oriented compressed storage mode indices.
+* MUY Y-oriented compressed storage mode indices.
+* MUZ Z-oriented compressed storage mode indices.
+* IPW W-oriented permutation matrices.
+* IPX X-oriented permutation matrices.
+* IPY Y-oriented permutation matrices.
+* IPZ Z-oriented permutation matrices.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK
+ INTEGER ISPLH,LX,LZ,L4,MAT(LX*LZ),KN(*),MUW(L4),MUX(L4),MUY(L4),
+ 1 MUZ(L4),IPW(L4),IPX(L4),IPY(L4),IPZ(L4),IMPX
+*----
+* LOCAL VARIABLES
+*----
+ REAL HW(14,14),HX(14,14),HY(14,14),HZ(14,14),HL(2,2),RFAC(28,7),
+ 1 RF6(24,6),RF7(28,7)
+ INTEGER NCODE(6),IJ1(14),IJ2(14),IJ27(14),IJ16(12),IJ26(12),
+ 1 IJ17(14)
+ INTEGER, DIMENSION(:), ALLOCATABLE :: IDX,IDY
+ COMMON /ELEMB/LC,T(5),TS(5),R(5,5),RS(5,5),Q(5,5),QS(5,5),V(5,4),
+ 1 E(5,5),RH(7,7),QH(7,7),RT(3,3),QT(3,3)
+ DATA HL / 1.0,2*0.0,1.0/
+ DATA IJ16,IJ26 /1,2,3,4,5,6,1,2,3,4,5,6,6*1,6*2/
+ DATA IJ17,IJ27 /1,2,3,4,5,6,7,1,2,3,4,5,6,7,7*1,7*2/
+ DATA RF6/
+ >1.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,1.0,0.5,
+ >1.0,0.0,1.0,1.0,0.0,0.5,1.0,0.0,0.0,0.0,0.0,0.0,
+ >0.0,1.0,1.0,1.0,0.5,0.0,1.0,1.0,0.0,0.0,0.5,1.0,
+ >0.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,
+ >0.0,1.0,1.0,0.5,1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,
+ >1.0,0.0,1.0,0.5,0.0,1.0,0.0,0.0,1.0,0.0,0.0,0.0,
+ >0.0,1.0,0.5,1.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,
+ >1.0,0.0,0.5,1.0,0.0,1.0,0.0,0.0,0.0,1.0,0.0,0.0,
+ >0.0,0.5,1.0,1.0,1.0,0.0,1.0,0.5,0.0,0.0,1.0,1.0,
+ >0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,
+ >0.0,0.0,0.0,0.0,0.0,1.0,0.5,1.0,0.0,0.0,1.0,1.0,
+ >0.5,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0/
+ DATA RF7/
+ >1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.5,0.0,1.0,0.5,
+ >1.0,0.0,1.0,0.5,1.0,0.0,0.5,1.0,0.0,0.0,0.0,0.0,0.0,0.0,
+ >0.0,1.0,1.0,0.5,1.0,0.5,0.0,1.0,1.0,0.0,0.5,0.0,0.5,1.0,
+ >0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,
+ >0.0,1.0,1.0,0.5,0.5,1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,
+ >1.0,0.0,1.0,0.5,0.5,0.0,1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,
+ >0.0,0.5,0.5,1.0,0.5,0.5,0.0,0.5,0.5,0.0,1.0,0.0,0.5,0.5,
+ >0.5,0.0,0.5,1.0,0.5,0.0,0.5,0.0,0.0,0.0,1.0,0.0,0.0,0.0,
+ >0.0,1.0,0.5,0.5,1.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,
+ >1.0,0.0,0.5,0.5,1.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,
+ >0.0,0.5,1.0,0.5,1.0,1.0,0.0,1.0,0.5,0.0,0.5,0.0,1.0,1.0,
+ >0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,
+ >0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.5,1.0,0.0,0.5,0.0,1.0,1.0,
+ >0.5,0.0,1.0,0.5,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0/
+*
+ IF(ISPLH.EQ.1) THEN
+ LC=6
+ DO 10 I=1,2*LC
+ IJ1(I)=IJ16(I)
+ IJ2(I)=IJ26(I)
+ 10 CONTINUE
+ DO 25 I=1,4*LC
+ DO 20 J=1,LC
+ RFAC(I,J)=RF6(I,J)
+ 20 CONTINUE
+ 25 CONTINUE
+ ELSE
+ LC=7
+ DO 30 I=1,2*LC
+ IJ1(I)=IJ17(I)
+ IJ2(I)=IJ27(I)
+ 30 CONTINUE
+ DO 45 I=1,4*LC
+ DO 40 J=1,LC
+ RFAC(I,J)=RF7(I,J)
+ 40 CONTINUE
+ 45 CONTINUE
+ ENDIF
+ LL=2*LC
+ DO 55 I=1,LL
+ I1=IJ1(I)
+ I2=IJ2(I)
+ DO 50 J=1,LL
+ J1=IJ1(J)
+ J2=IJ2(J)
+ HW(I,J) = RFAC(I1 ,J1) * HL(I2,J2)
+ HX(I,J) = RFAC(I1+LC ,J1) * HL(I2,J2)
+ HY(I,J) = RFAC(I1+2*LC,J1) * HL(I2,J2)
+ HZ(I,J) = RFAC(I1+3*LC,J1)
+ 50 CONTINUE
+ 55 CONTINUE
+*
+ DO 65 I=1,LL
+ I1 = IJ1(I)
+ I2 = IJ2(I)
+ DO 60 J=1,LL
+ J1 = IJ1(J)
+ J2 = IJ2(J)
+ HW(I,J) = RFAC(I1 ,J1) * HL(I2,J2)
+ HX(I,J) = RFAC(I1+LC ,J1) * HL(I2,J2)
+ HY(I,J) = RFAC(I1+2*LC,J1) * HL(I2,J2)
+ HZ(I,J) = RFAC(I1+3*LC,J1)
+ 60 CONTINUE
+ 65 CONTINUE
+*----
+* COMPUTE THE PERMUTATION VECTORS
+*----
+ DO 70 I=1,L4
+ IPW(I)=I
+ IPX(I)=0
+ IPY(I)=0
+ IPZ(I)=0
+ 70 CONTINUE
+ LT4 = L4
+ LPZ = LZ
+ IF(LZ.GT.1) THEN
+ LPZ = LZ+1
+ CALL LCMGET (IPTRK,'NCODE',NCODE)
+ IF((NCODE(5).EQ.7).OR.(NCODE(6).EQ.7)) LPZ = LZ
+ IF((NCODE(5).EQ.7).AND.(NCODE(6).EQ.7)) LPZ = LZ-1
+ LT4 = L4/LPZ
+ ENDIF
+ ALLOCATE(IDX(LT4),IDY(LT4))
+ CALL LCMGET (IPTRK,'ILX',IDX)
+ CALL LCMGET (IPTRK,'ILY',IDY)
+ DO 85 KZ = 1, LPZ
+ DO 80 KX = 1, LT4
+ IPX(KX+(KZ-1)*LT4) = IDX(KX) + (KZ-1)*LT4
+ IPY(KX+(KZ-1)*LT4) = IDY(KX) + (KZ-1)*LT4
+ 80 CONTINUE
+ 85 CONTINUE
+ DEALLOCATE(IDY,IDX)
+ KEL = 0
+ DO 95 KX = 1, LT4
+ DO 90 KZ = 1, LPZ
+ KEL = KEL + 1
+ IPZ(KX+(KZ-1)*LT4) = KEL
+ 90 CONTINUE
+ 95 CONTINUE
+*----
+* COMPUTE THE COMPRESSED DIAGONAL STORAGE INDICES
+*----
+ DO 100 I=1,L4
+ MUW(I)=1
+ MUX(I)=1
+ MUY(I)=1
+ MUZ(I)=1
+ 100 CONTINUE
+ NUM1=0
+ DO 130 K=1,LX*LZ
+ IF(MAT(K).LE.0) GO TO 130
+ DO 120 I=1,LL
+ INW1=KN(NUM1+I)
+ IF(INW1.EQ.0) GO TO 120
+ INX1=IPX(INW1)
+ INY1=IPY(INW1)
+ INZ1=IPZ(INW1)
+ DO 110 J=1,LL
+ INW2=KN(NUM1+J)
+ IF(INW2.EQ.0) GO TO 110
+ INX2=IPX(INW2)
+ INY2=IPY(INW2)
+ INZ2=IPZ(INW2)
+ IF((HW(I,J).NE.0.0).AND.(INW2.LT.INW1))
+ > MUW(INW1)=MAX0(MUW(INW1),INW1-INW2+1)
+ IF((HX(I,J).NE.0.0).AND.(INX2.LT.INX1))
+ > MUX(INX1)=MAX0(MUX(INX1),INX1-INX2+1)
+ IF((HY(I,J).NE.0.0).AND.(INY2.LT.INY1))
+ > MUY(INY1)=MAX0(MUY(INY1),INY1-INY2+1)
+ IF((HZ(I,J).NE.0.0).AND.(INZ2.LT.INZ1))
+ > MUZ(INZ1)=MAX0(MUZ(INZ1),INZ1-INZ2+1)
+ 110 CONTINUE
+ 120 CONTINUE
+ NUM1=NUM1+LL
+ 130 CONTINUE
+ IF(IMPX.GE.5) THEN
+ WRITE(6,510) 'IPW :',(IPW(I),I=1,L4)
+ WRITE(6,510) 'MUW :',(MUW(I),I=1,L4)
+ WRITE(6,510) 'IPX :',(IPX(I),I=1,L4)
+ WRITE(6,510) 'MUX :',(MUX(I),I=1,L4)
+ WRITE(6,510) 'IPY :',(IPY(I),I=1,L4)
+ WRITE(6,510) 'MUY :',(MUY(I),I=1,L4)
+ IF(LZ.GT.1) THEN
+ WRITE(6,510) 'IPZ :',(IPZ(I),I=1,L4)
+ WRITE(6,510) 'MUZ :',(MUZ(I),I=1,L4)
+ ENDIF
+ ENDIF
+*
+ MUWMAX=0
+ MUXMAX=0
+ MUYMAX=0
+ MUZMAX=0
+ IIMAWW=0
+ IIMAWX=0
+ IIMAWY=0
+ IIMAWZ=0
+ DO 140 I=1,L4
+ MUWMAX=MAX(MUWMAX,MUW(I))
+ MUXMAX=MAX(MUXMAX,MUX(I))
+ MUYMAX=MAX(MUYMAX,MUY(I))
+ MUZMAX=MAX(MUZMAX,MUZ(I))
+ IIMAWW=IIMAWW+MUW(I)
+ MUW(I)=IIMAWW
+ IIMAWX=IIMAWX+MUX(I)
+ MUX(I)=IIMAWX
+ IIMAWY=IIMAWY+MUY(I)
+ MUY(I)=IIMAWY
+ IIMAWZ=IIMAWZ+MUZ(I)
+ MUZ(I)=IIMAWZ
+ 140 CONTINUE
+ IF(IMPX.GT.0) WRITE (6,500) MUWMAX,MUXMAX,MUYMAX,MUZMAX
+ RETURN
+*
+ 500 FORMAT(/41H TRICH3: MAXIMUM BANDWIDTH ALONG W AXIS =,I5/
+ 1 27X,14HALONG X AXIS =,I5/27X,14HALONG Y AXIS =,I5/27X,
+ 2 14HALONG Z AXIS =,I5)
+ 510 FORMAT(/1X,A5/(1X,20I6))
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:19:05 +0000