From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Trivac/src/TRICH3.f | 257 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 257 insertions(+) create mode 100755 Trivac/src/TRICH3.f (limited to 'Trivac/src/TRICH3.f') 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 -- cgit v1.2.3