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/TRICHH.f | 364 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 364 insertions(+) create mode 100755 Trivac/src/TRICHH.f (limited to 'Trivac/src/TRICHH.f') diff --git a/Trivac/src/TRICHH.f b/Trivac/src/TRICHH.f new file mode 100755 index 0000000..5f46445 --- /dev/null +++ b/Trivac/src/TRICHH.f @@ -0,0 +1,364 @@ +*DECK TRICHH + SUBROUTINE TRICHH(IMPX,MAXKN,NBLOS,LXH,LZ,IELEM,ISPLH,L4,LL4F, + 1 LL4W,LL4X,LL4Y,LL4Z,SIDE,ZZ,FRZ,IPERT,KN,V,H,MUW,MUX,MUY,MUZ, + 2 IPBBW,IPBBX,IPBBY,IPBBZ,BBW,BBX,BBY,BBZ,CTRAN) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Thomas-Raviart-Schneider (dual) finite element unknown numbering for +* ADI solution in a 3D hexagonal domain. Compute the storage info for +* ADI matrices in compressed diagonal storage mode. Compute the ADI +* permutation vectors. Compute the group-independent WB, XB, YB and ZB +* matrices. +* +*Copyright: +* Copyright (C) 2006 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. Hebert +* +*Parameters: input +* IMPX print parameter. +* MAXKN number of components in KN. +* NBLOS number of lozenges per direction in 3D with mesh-splitting. +* LXH number of hexagons in a plane. +* LZ number of elements along the Z axis. +* IELEM degree of the Lagrangian finite elements: =1 (linear); +* =2 (parabolic); =3 (cubic). +* ISPLH mesh-splitting in 3*ISPLH**2 lozenges per hexagon. +* L4 total number of unknown (variational coefficients) per +* energy group (order of system matrices). +* LL4F exact number of flux unknowns. +* LL4W exact number of W-directed current unknowns. +* LL4X exact number of X-directed current unknowns. +* LL4Y exact number of Y-directed current unknowns. +* LL4Z exact number of Z-directed current unknowns. +* SIDE side of an hexagon. +* ZZ Z-directed mesh spacings. +* FRZ volume fractions for the axial SYME boundary condition. +* IPERT mixture permutation index. +* KN ADI permutation indices for the volumes and currents. +* V nodal coupling matrix matrix. +* H Piolat (hexagonal) coupling matrix. +* +*Parameters: output +* MUW W-directed compressed diagonal mode indices. +* MUX X-directed compressed diagonal mode indices. +* MUY Y-directed compressed diagonal mode indices. +* MUZ Z-directed compressed diagonal mode indices. +* IPBBW W-directed perdue storage indices. +* IPBBX X-directed perdue storage indices. +* IPBBY Y-directed perdue storage indices. +* IPBBZ Z-directed perdue storage indices. +* BBW W-directed flux-current matrices. +* BBX X-directed flux-current matrices. +* BBY Y-directed flux-current matrices. +* BBZ Z-directed flux-current matrices. +* CTRAN tranverse coupling Piolat unit matrix. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IMPX,MAXKN,NBLOS,LXH,LZ,IELEM,ISPLH,L4,IPERT(NBLOS), + 1 KN(NBLOS,MAXKN/NBLOS),LL4F,LL4W,LL4X,LL4Y,LL4Z,MUW(L4), + 2 MUX(L4),MUY(L4),MUZ(L4),IPBBW(2*IELEM,LL4W),IPBBX(2*IELEM,LL4X), + 3 IPBBY(2*IELEM,LL4Y),IPBBZ(2*IELEM,LL4Z) + REAL SIDE,ZZ(3,NBLOS),FRZ(NBLOS),V(IELEM+1,IELEM), + 1 H(IELEM+1,IELEM),BBW(2*IELEM,LL4W),BBX(2*IELEM,LL4X), + 2 BBY(2*IELEM,LL4Y),BBZ(2*IELEM,LL4Z) + DOUBLE PRECISION CTRAN((IELEM+1)*IELEM,(IELEM+1)*IELEM) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT,DENOM,VOL0 +* + NELEH=(IELEM+1)*IELEM**2 + NELEZ=6*IELEM**2 + NBC=INT((SQRT(REAL((4*LXH-1)/3))+1.)/2.) + IF(LL4F.GT.3*NBLOS*IELEM**3) CALL XABORT('TRICHH: BUG1.') + IF(LL4W.GT.(2*NBLOS*IELEM+(2*NBC-1)*ISPLH*LZ)*IELEM**2) + 1 CALL XABORT('TRICHH: BUG2.') +*---- +* COMPUTE THE TRANVERSE COUPLING PIOLAT UNIT MATRIX +*---- + CTRAN(:(IELEM+1)*IELEM,:(IELEM+1)*IELEM)=0.0D0 + CNORM=SIDE*SIDE/SQRT(3.) + I=0 + DO 22 JS=1,IELEM + DO 21 JT=1,IELEM+1 + J=0 + I=I+1 + SSS=1.0 + DO 20 IT=1,IELEM + DO 10 IS=1,IELEM+1 + J=J+1 + CTRAN(I,J)=SSS*CNORM*H(IS,JS)*H(JT,IT) + 10 CONTINUE + SSS=-SSS + 20 CONTINUE + 21 CONTINUE + 22 CONTINUE + IF(IMPX.GT.1) THEN + WRITE(6,*) 'TRICHH: MATRIX CTRAN' + DO 30 I=1,(IELEM+1)*IELEM + WRITE(6,'(10(1X,1P,E12.4))') (CTRAN(I,J),J=1,(IELEM+1)*IELEM) + 30 CONTINUE + WRITE(6,*) ' ' + ENDIF +*---- +* COMPUTE THE W-, X- ,Y- AND Z-ORIENTED SYSTEM BANDWIDTH VECTORS +*---- + MUW(:L4)=1 + MUX(:L4)=1 + MUY(:L4)=1 + MUZ(:L4)=1 + IPBBW(:2*IELEM,:LL4W)=0 + IPBBX(:2*IELEM,:LL4X)=0 + IPBBY(:2*IELEM,:LL4Y)=0 + IPBBZ(:2*IELEM,:LL4Z)=0 + NUM=0 + DO 80 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 80 + NUM=NUM+1 + DO 64 K5=0,1 ! TWO LOZENGES PER HEXAGON + DO 63 K4=0,IELEM-1 + DO 62 K3=0,IELEM-1 + DO 61 K2=1,IELEM+1 + KNW1=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + KNX1=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + KNY1=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INW1=ABS(KNW1) + INX1=ABS(KNX1)-LL4W + INY1=ABS(KNY1)-LL4W-LL4X + DO 40 K1=1,IELEM+1 + KNW2=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + KNX2=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + KNY2=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INW2=ABS(KNW2) + INX2=ABS(KNX2)-LL4W + INY2=ABS(KNY2)-LL4W-LL4X + IF((KNW2.NE.0).AND.(KNW1.NE.0)) THEN + MUW(INW1)=MAX(MUW(INW1),INW1-INW2+1) + MUW(INW2)=MAX(MUW(INW2),INW2-INW1+1) + ENDIF + IF((KNX2.NE.0).AND.(KNX1.NE.0)) THEN + MUX(INX1)=MAX(MUX(INX1),INX1-INX2+1) + MUX(INX2)=MAX(MUX(INX2),INX2-INX1+1) + ENDIF + IF((KNY2.NE.0).AND.(KNY1.NE.0)) THEN + MUY(INY1)=MAX(MUY(INY1),INY1-INY2+1) + MUY(INY2)=MAX(MUY(INY2),INY2-INY1+1) + ENDIF + 40 CONTINUE + DO 60 K1=0,IELEM-1 + IF(V(K2,K1+1).EQ.0.0) GO TO 60 + IF(K5.EQ.0) THEN + JND1=(NUM-1)*IELEM**3+K4*IELEM**2+K3*IELEM+K1+1 + JND2=(KN(NUM,1)-1)*IELEM**3+K4*IELEM**2+K3*IELEM+K1+1 + JND3=(KN(NUM,2)-1)*IELEM**3+K4*IELEM**2+K3*IELEM+K1+1 + ELSE + JND1=(KN(NUM,1)-1)*IELEM**3+K4*IELEM**2+K1*IELEM+K3+1 + JND2=(KN(NUM,2)-1)*IELEM**3+K4*IELEM**2+K1*IELEM+K3+1 + JND3=(KN(NUM,3)-1)*IELEM**3+K4*IELEM**2+K1*IELEM+K3+1 + ENDIF + IF(KNW1.NE.0) CALL TRINDX(JND1,IPBBW(1,INW1),2*IELEM) + IF(KNX1.NE.0) CALL TRINDX(JND2,IPBBX(1,INX1),2*IELEM) + IF(KNY1.NE.0) CALL TRINDX(JND3,IPBBY(1,INY1),2*IELEM) + 60 CONTINUE + 61 CONTINUE + 62 CONTINUE + 63 CONTINUE + 64 CONTINUE + DO 73 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 72 K2=0,IELEM-1 + DO 71 K1=0,IELEM-1 + KNZ1=KN(NUM,3+6*NELEH+2*K5*IELEM**2+K2*IELEM+K1+1) + KNZ2=KN(NUM,3+6*NELEH+(2*K5+1)*IELEM**2+K2*IELEM+K1+1) + INZ1=ABS(KNZ1)-LL4W-LL4X-LL4Y + INZ2=ABS(KNZ2)-LL4W-LL4X-LL4Y + IF((KNZ1.NE.0).AND.(KNZ2.NE.0)) THEN + MUZ(INZ2)=MAX(MUZ(INZ2),INZ2-INZ1+1) + MUZ(INZ1)=MAX(MUZ(INZ1),INZ1-INZ2+1) + ENDIF + DO 70 K3=0,IELEM-1 + IF(K5.EQ.0) THEN + JND1=(NUM-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + ELSE + JND1=(KN(NUM,K5)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + ENDIF + IF(KNZ1.NE.0) CALL TRINDX(JND1,IPBBZ(1,INZ1),2*IELEM) + IF(KNZ2.NE.0) CALL TRINDX(JND1,IPBBZ(1,INZ2),2*IELEM) + 70 CONTINUE + 71 CONTINUE + 72 CONTINUE + 73 CONTINUE + 80 CONTINUE +* + MUWMAX=0 + IIMAXW=0 + DO 90 I=1,LL4W + MUWMAX=MAX(MUWMAX,MUW(I)) + IIMAXW=IIMAXW+MUW(I) + MUW(I)=IIMAXW + 90 CONTINUE + MUXMAX=0 + IIMAXX=0 + DO 100 I=1,LL4X + MUXMAX=MAX(MUXMAX,MUX(I)) + IIMAXX=IIMAXX+MUX(I) + MUX(I)=IIMAXX + 100 CONTINUE + MUYMAX=0 + IIMAXY=0 + DO 110 I=1,LL4Y + MUYMAX=MAX(MUYMAX,MUY(I)) + IIMAXY=IIMAXY+MUY(I) + MUY(I)=IIMAXY + 110 CONTINUE + MUZMAX=0 + IIMAXZ=0 + DO 120 I=1,LL4Z + MUZMAX=MAX(MUZMAX,MUZ(I)) + IIMAXZ=IIMAXZ+MUZ(I) + MUZ(I)=IIMAXZ + 120 CONTINUE + IF(IMPX.GT.0) THEN + WRITE (6,600) MUWMAX,MUXMAX,MUYMAX,MUZMAX + WRITE (6,610) IIMAXW,IIMAXX,IIMAXY,IIMAXZ + ENDIF +*---- +* COMPUTE THE FLUX-CURRENT COUPLING MATRICES WB, XB, YB AND ZB. +*---- + BBW(:2*IELEM,:LL4W)=0.0 + BBX(:2*IELEM,:LL4X)=0.0 + BBY(:2*IELEM,:LL4Y)=0.0 + BBZ(:2*IELEM,:LL4Z)=0.0 + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + DENOM=0.5D0*SQRT(3.D00)*SIDE + NUM=0 + DO 260 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 260 + NUM=NUM+1 + DZ=ZZ(1,IPERT(KEL)) + VOL0=TTTT*DZ*FRZ(KEL) + DO 194 K5=0,1 + DO 193 K4=0,IELEM-1 + DO 192 K3=0,IELEM-1 + DO 191 K2=1,IELEM+1 + KNW1=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + KNX1=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + KNY1=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INW1=ABS(KNW1) + INX1=ABS(KNX1)-LL4W + INY1=ABS(KNY1)-LL4W-LL4X + DO 190 K1=0,IELEM-1 + IF(V(K2,K1+1).EQ.0.0) GO TO 190 + IF(K5.EQ.0) THEN + SSS=(-1.0)**K1 + JND1=(NUM-1)*IELEM**3+K4*IELEM**2+K3*IELEM+K1+1 + JND2=(KN(NUM,1)-1)*IELEM**3+K4*IELEM**2+K3*IELEM+K1+1 + JND3=(KN(NUM,2)-1)*IELEM**3+K4*IELEM**2+K3*IELEM+K1+1 + ELSE + SSS=1.0 + JND1=(KN(NUM,1)-1)*IELEM**3+K4*IELEM**2+K1*IELEM+K3+1 + JND2=(KN(NUM,2)-1)*IELEM**3+K4*IELEM**2+K1*IELEM+K3+1 + JND3=(KN(NUM,3)-1)*IELEM**3+K4*IELEM**2+K1*IELEM+K3+1 + ENDIF + IF(KNW1.NE.0.0) THEN + KK=0 + DO 130 I=1,2*IELEM + IF(IPBBW(I,INW1).EQ.JND1) THEN + KK=I + GO TO 140 + ENDIF + 130 CONTINUE + CALL XABORT('TRICHH: BUG3.') + 140 SG=REAL(SIGN(1,KNW1)) + BBW(KK,INW1)=BBW(KK,INW1)+SG*SSS*REAL(VOL0/DENOM)*V(K2,K1+1) + ENDIF + IF(KNX1.NE.0.0) THEN + KK=0 + DO 150 I=1,2*IELEM + IF(IPBBX(I,INX1).EQ.JND2) THEN + KK=I + GO TO 160 + ENDIF + 150 CONTINUE + CALL XABORT('TRICHH: BUG4.') + 160 SG=REAL(SIGN(1,KNX1)) + BBX(KK,INX1)=BBX(KK,INX1)+SG*SSS*REAL(VOL0/DENOM)*V(K2,K1+1) + ENDIF + IF(KNY1.NE.0.0) THEN + KK=0 + DO 170 I=1,2*IELEM + IF(IPBBY(I,INY1).EQ.JND3) THEN + KK=I + GO TO 180 + ENDIF + 170 CONTINUE + CALL XABORT('TRICHH: BUG5.') + 180 SG=REAL(SIGN(1,KNY1)) + BBY(KK,INY1)=BBY(KK,INY1)+SG*SSS*REAL(VOL0/DENOM)*V(K2,K1+1) + ENDIF + 190 CONTINUE + 191 CONTINUE + 192 CONTINUE + 193 CONTINUE + 194 CONTINUE + DO 253 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 252 K2=0,IELEM-1 + DO 251 K1=0,IELEM-1 + KNZ1=KN(NUM,3+6*NELEH+2*K5*IELEM**2+K2*IELEM+K1+1) + KNZ2=KN(NUM,3+6*NELEH+(2*K5+1)*IELEM**2+K2*IELEM+K1+1) + INZ1=ABS(KNZ1)-LL4W-LL4X-LL4Y + INZ2=ABS(KNZ2)-LL4W-LL4X-LL4Y + DO 250 K3=0,IELEM-1 + IF(K5.EQ.0) THEN + JND1=(NUM-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + ELSE + JND1=(KN(NUM,K5)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + ENDIF + IF(KNZ1.NE.0) THEN + KK=0 + DO 210 I=1,2*IELEM + IF(IPBBZ(I,INZ1).EQ.JND1) THEN + KK=I + GO TO 220 + ENDIF + 210 CONTINUE + CALL XABORT('TRICHH: BUG6.') + 220 SG=REAL(SIGN(1,KNZ1)) + BBZ(KK,INZ1)=BBZ(KK,INZ1)+SG*REAL(VOL0/DZ)*V(1,K3+1) + ENDIF + IF(KNZ2.NE.0) THEN + KK=0 + DO 230 I=1,2*IELEM + IF(IPBBZ(I,INZ2).EQ.JND1) THEN + KK=I + GO TO 240 + ENDIF + 230 CONTINUE + CALL XABORT('TRICHH: BUG7.') + 240 SG=REAL(SIGN(1,KNZ2)) + BBZ(KK,INZ2)=BBZ(KK,INZ2)+SG*REAL(VOL0/DZ)*V(IELEM+1,K3+1) + ENDIF + 250 CONTINUE + 251 CONTINUE + 252 CONTINUE + 253 CONTINUE + 260 CONTINUE + RETURN +* + 600 FORMAT(/52H TRICHH: MAXIMUM BANDWIDTH FOR W-ORIENTED MATRICES =, + 1 I4/27X,25HFOR X-ORIENTED MATRICES =,I4/27X,16HFOR Y-ORIENTED M, + 2 9HATRICES =,I4/27X,25HFOR Z-ORIENTED MATRICES =,I4) + 610 FORMAT(/40H TRICHH: LENGTH OF W-ORIENTED MATRICES =,I10/16X, + 1 24HOF X-ORIENTED MATRICES =,I10/16X,24HOF Y-ORIENTED MATRICES =, + 2 I10/16X,24HOF Z-ORIENTED MATRICES =,I10) + END -- cgit v1.2.3