diff options
Diffstat (limited to 'Trivac/src/TRIHWW.f')
| -rwxr-xr-x | Trivac/src/TRIHWW.f | 418 |
1 files changed, 418 insertions, 0 deletions
diff --git a/Trivac/src/TRIHWW.f b/Trivac/src/TRIHWW.f new file mode 100755 index 0000000..9f049d1 --- /dev/null +++ b/Trivac/src/TRIHWW.f @@ -0,0 +1,418 @@ +*DECK TRIHWW + SUBROUTINE TRIHWW(NBMIX,NBLOS,IELEM,LL4F,LL4W,MAT,SIDE,ZZ,FRZ, + 1 QFR,IPERT,KN,SGD,XSGD,MUW,IPBBW,LC,R,V,BBW,TTF,AW,C11W) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of system matrices for a Thomas-Raviart-Schneider (dual) +* finite element method in hexagonal 3-D diffusion approximation. +* Note: system matrices should be initialized by the calling program. +* +*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 +* NBMIX number of mixtures. +* NBLOS number of lozenges per direction, taking into account +* mesh-splitting. +* IELEM degree of the Lagrangian finite elements: =1 (linear); +* =2 (parabolic); =3 (cubic). +* ICOL type of quadrature: =1 (analytical integration); +* =2 (Gauss-Lobatto); =3 (Gauss-Legendre). +* ISPLH mesh-splitting index. Each hexagon is splitted into 3*ISPLH**2 +* lozenges. +* LL4F number of flux components. +* LL4W number of W-directed currents. +* LL4X number of X-directed currents. +* LL4Y number of Y-directed currents. +* LL4Z number of Z-directed currents. +* MAT mixture index assigned to each element. +* SIDE side of an hexagon. +* ZZ Z-directed mesh spacings. +* FRZ volume fractions for the axial SYME boundary condition. +* QFR element-ordered boundary conditions. +* IPERT mixture permutation index. +* KN ADI permutation indices for the volumes and currents. +* SGD nuclear properties by material mixture: +* SGD(L,1)= X-oriented diffusion coefficients; +* SGD(L,2)= Y-oriented diffusion coefficients; +* SGD(L,3)= Z-oriented diffusion coefficients; +* SGD(L,4)= removal macroscopic cross section. +* XSGD one over nuclear properties. +* MUW W-directed compressed storage mode indices. +* MUX X-directed compressed storage mode indices. +* MUY Y-directed compressed storage mode indices. +* MUZ Z-directed compressed storage 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. +* LC order of the unit matrices. +* R unit matrix. +* V unit matrix. +* BBW W-directed flux-current matrices. +* BBX X-directed flux-current matrices. +* BBY Y-directed flux-current matrices. +* BBZ Z-directed flux-current matrices. +* +*Parameters: output +* TTF flux-flux matrices. +* AW W-directed main current-current matrices. Dimensionned to +* MUW(LL4W). +* AX X-directed main current-current matrices. Dimensionned to +* MUX(LL4X). +* AY Y-directed main current-current matrices. Dimensionned to +* MUY(LL4Y). +* AZ Z-directed main current-current matrices. Dimensionned to +* MUZ(LL4Z). +* C11W W-directed main current-current matrices to be factorized. +* Dimensionned to MUW(LL4W). +* C11X X-directed main current-current matrices to be factorized. +* Dimensionned to MUX(LL4X). +* C11Y Y-directed main current-current matrices to be factorized. +* Dimensionned to MUY(LL4Y). +* C11Z Z-directed main current-current matrices to be factorized. +* Dimensionned to MUZ(LL4Z). +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,LL4F,LL4W,MAT(3,NBLOS),IPERT(NBLOS), + 1 KN(NBLOS,3+6*(IELEM+2)*IELEM**2),MUW(LL4W),IPBBW(2*IELEM,LL4W),LC + REAL SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8),SGD(NBMIX,4), + 1 XSGD(NBMIX,4),R(LC,LC),V(LC,LC-1),TTF(LL4F),BBW(2*IELEM,LL4W), + 2 AW(*),C11W(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION FFF,TTTT + REAL QQ(5,5) +*---- +* W-ORIENTED COUPLINGS +*---- + DO 25 I0=1,IELEM + DO 20 J0=1,IELEM + FFF=0.0D0 + DO 10 K0=2,IELEM + FFF=FFF+V(K0,I0)*V(K0,J0)/R(K0,K0) + 10 CONTINUE + IF(ABS(FFF).LE.1.0E-6) FFF=0.0D0 + QQ(I0,J0)=REAL(FFF) + 20 CONTINUE + 25 CONTINUE +* + NELEH=(IELEM+1)*IELEM**2 + IIMAW=MUW(LL4W) + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 50 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 50 + NUM=NUM+1 + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 50 + DZ=ZZ(1,IPERT(KEL)) + VOL0=REAL(TTTT*DZ*FRZ(KEL)) + DINV=XSGD(IBM,1) + SIG3=SGD(IBM,3)/(DZ*DZ) + SIG4=SGD(IBM,4) + DO 34 K5=0,1 + DO 33 K4=0,IELEM-1 + DO 32 K3=0,IELEM-1 + DO 31 K2=1,IELEM+1 + KNW1=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INW1=ABS(KNW1) + DO 30 K1=1,IELEM+1 + KNW2=KN(NUM,3+K5*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INW2=ABS(KNW2) + IF((KNW2.NE.0).AND.(KNW1.NE.0)) THEN + L=MUW(INW1)-INW1+INW2 + SG=REAL(SIGN(1,KNW1)*SIGN(1,KNW2)) + IF(K1.LE.K2) AW(L)=AW(L)-(4./3.)*SG*VOL0*DINV*R(K2,K1) + IF(K1.EQ.K2) THEN + IF((K1.EQ.1).AND.(K5.EQ.0)) AW(L)=AW(L)-QFR(NUM,1) + IF((K1.EQ.IELEM+1).AND.(K5.EQ.1)) AW(L)=AW(L)-QFR(NUM,2) + ENDIF + ENDIF + 30 CONTINUE + 31 CONTINUE + 32 CONTINUE + 33 CONTINUE + 34 CONTINUE + DO 42 K3=0,IELEM-1 + DO 41 K2=0,IELEM-1 + DO 40 K1=0,IELEM-1 + JND1=(NUM-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + JND2=(KN(NUM,1)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + JND3=(KN(NUM,2)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + TTF(JND1)=TTF(JND1)+VOL0*SIG4+VOL0*QQ(K3+1,K3+1)*SIG3 + TTF(JND2)=TTF(JND2)+VOL0*SIG4+VOL0*QQ(K3+1,K3+1)*SIG3 + TTF(JND3)=TTF(JND3)+VOL0*SIG4+VOL0*QQ(K3+1,K3+1)*SIG3 + 40 CONTINUE + 41 CONTINUE + 42 CONTINUE + 50 CONTINUE +*---- +* COMPUTE THE W-ORIENTED SYSTEM MATRIX AFTER FLUX ELIMINATION +*---- + DO 60 I0=1,IIMAW + C11W(I0)=-AW(I0) + 60 CONTINUE + MUIM1=0 + DO 90 I=1,LL4W + MUI=MUW(I) + DO 80 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 75 I0=1,2*IELEM + II=IPBBW(I0,I) + IF(II.EQ.0) GO TO 80 + DO 70 J0=1,2*IELEM + JJ=IPBBW(J0,J) + IF(II.EQ.JJ) C11W(KEY)=C11W(KEY)+BBW(I0,I)*BBW(J0,J)/TTF(II) + 70 CONTINUE + 75 CONTINUE + 80 CONTINUE + MUIM1=MUI + 90 CONTINUE + RETURN + END +* + SUBROUTINE TRIHWX(NBMIX,NBLOS,IELEM,LL4F,LL4W,LL4X,MAT,SIDE,ZZ, + 1 FRZ,QFR,IPERT,KN,XSGD,MUX,IPBBX,LC,R,BBX,TTF,AX,C11X) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,LL4F,LL4W,LL4X,MAT(3,NBLOS), + 1 MUX(LL4X),IPBBX(2*IELEM,LL4X),LC,IPERT(NBLOS), + 2 KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8),XSGD(NBMIX,4), + 1 R(LC,LC),TTF(LL4F),BBX(2*IELEM,LL4X),AX(*),C11X(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT +*---- +* X-ORIENTED COUPLINGS +*---- + NELEH=(IELEM+1)*IELEM**2 + IIMAX=MUX(LL4X) + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 120 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 120 + NUM=NUM+1 + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 120 + VOL0=REAL(TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL)) + DINV=XSGD(IBM,1) + DO 114 K5=0,1 + DO 113 K4=0,IELEM-1 + DO 112 K3=0,IELEM-1 + DO 111 K2=1,IELEM+1 + KNX1=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INX1=ABS(KNX1)-LL4W + DO 110 K1=1,IELEM+1 + KNX2=KN(NUM,3+(K5+2)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INX2=ABS(KNX2)-LL4W + IF((KNX2.NE.0).AND.(KNX1.NE.0)) THEN + L=MUX(INX1)-INX1+INX2 + SG=REAL(SIGN(1,KNX1)*SIGN(1,KNX2)) + IF(K1.LE.K2) AX(L)=AX(L)-(4./3.)*SG*VOL0*DINV*R(K2,K1) + IF(K1.EQ.K2) THEN + IF((K1.EQ.1).AND.(K5.EQ.0)) AX(L)=AX(L)-QFR(NUM,3) + IF((K1.EQ.IELEM+1).AND.(K5.EQ.1)) AX(L)=AX(L)-QFR(NUM,4) + ENDIF + ENDIF + 110 CONTINUE + 111 CONTINUE + 112 CONTINUE + 113 CONTINUE + 114 CONTINUE + 120 CONTINUE +*---- +* COMPUTE THE X-ORIENTED SYSTEM MATRIX AFTER FLUX ELIMINATION +*---- + DO 130 I0=1,IIMAX + C11X(I0)=-AX(I0) + 130 CONTINUE + MUIM1=0 + DO 160 I=1,LL4X + MUI=MUX(I) + DO 150 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 145 I0=1,2*IELEM + II=IPBBX(I0,I) + IF(II.EQ.0) GO TO 150 + DO 140 J0=1,2*IELEM + JJ=IPBBX(J0,J) + IF(II.EQ.JJ) C11X(KEY)=C11X(KEY)+BBX(I0,I)*BBX(J0,J)/TTF(II) + 140 CONTINUE + 145 CONTINUE + 150 CONTINUE + MUIM1=MUI + 160 CONTINUE + RETURN + END +* + SUBROUTINE TRIHWY(NBMIX,NBLOS,IELEM,LL4F,LL4W,LL4X,LL4Y,MAT, + 1 SIDE,ZZ,FRZ,QFR,IPERT,KN,XSGD,MUY,IPBBY,LC,R,BBY,TTF,AY,C11Y) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,LL4F,LL4W,LL4X,LL4Y,MAT(3,NBLOS), + 1 MUY(LL4Y),IPBBY(2*IELEM,LL4Y),LC,IPERT(NBLOS), + 2 KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8),XSGD(NBMIX,4), + 1 R(LC,LC),TTF(LL4F),BBY(2*IELEM,LL4Y),AY(*),C11Y(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT +*---- +* Y-ORIENTED COUPLINGS +*---- + NELEH=(IELEM+1)*IELEM**2 + IIMAY=MUY(LL4Y) + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 220 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 220 + NUM=NUM+1 + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 220 + VOL0=REAL(TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL)) + DINV=XSGD(IBM,1) + DO 214 K5=0,1 + DO 213 K4=0,IELEM-1 + DO 212 K3=0,IELEM-1 + DO 211 K2=1,IELEM+1 + KNY1=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K2) + INY1=ABS(KNY1)-LL4W-LL4X + DO 210 K1=1,IELEM+1 + KNY2=KN(NUM,3+(K5+4)*NELEH+(K4*IELEM+K3)*(IELEM+1)+K1) + INY2=ABS(KNY2)-LL4W-LL4X + IF((KNY2.NE.0).AND.(KNY1.NE.0)) THEN + L=MUY(INY1)-INY1+INY2 + SG=REAL(SIGN(1,KNY1)*SIGN(1,KNY2)) + IF(K1.LE.K2) AY(L)=AY(L)-(4./3.)*SG*VOL0*DINV*R(K2,K1) + IF(K1.EQ.K2) THEN + IF((K1.EQ.1).AND.(K5.EQ.0)) AY(L)=AY(L)-QFR(NUM,5) + IF((K1.EQ.IELEM+1).AND.(K5.EQ.1)) AY(L)=AY(L)-QFR(NUM,6) + ENDIF + ENDIF + 210 CONTINUE + 211 CONTINUE + 212 CONTINUE + 213 CONTINUE + 214 CONTINUE + 220 CONTINUE +*---- +* COMPUTE THE Y-ORIENTED SYSTEM MATRIX AFTER FLUX ELIMINATION +*---- + DO 230 I0=1,IIMAY + C11Y(I0)=-AY(I0) + 230 CONTINUE + MUIM1=0 + DO 260 I=1,LL4Y + MUI=MUY(I) + DO 250 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 245 I0=1,2*IELEM + II=IPBBY(I0,I) + IF(II.EQ.0) GO TO 250 + DO 240 J0=1,2*IELEM + JJ=IPBBY(J0,J) + IF(II.EQ.JJ) C11Y(KEY)=C11Y(KEY)+BBY(I0,I)*BBY(J0,J)/TTF(II) + 240 CONTINUE + 245 CONTINUE + 250 CONTINUE + MUIM1=MUI + 260 CONTINUE + RETURN + END +* + SUBROUTINE TRIHWZ(NBMIX,NBLOS,IELEM,ICOL,LL4F,LL4W,LL4X,LL4Y, + 1 LL4Z,MAT,SIDE,ZZ,FRZ,QFR,IPERT,KN,XSGD,MUZ,IPBBZ,LC,R,BBZ,TTF, + 2 AZ,C11Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,NBLOS,IELEM,ICOL,LL4F,LL4W,LL4X,LL4Y,LL4Z, + 1 MAT(3,NBLOS),MUZ(LL4Z),IPBBZ(2*IELEM,LL4Z),LC,IPERT(NBLOS), + 2 KN(NBLOS,3+6*(IELEM+2)*IELEM**2) + REAL SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8),XSGD(NBMIX,4), + 1 R(LC,LC),TTF(LL4F),BBZ(2*IELEM,LL4Z),AZ(*),C11Z(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION TTTT +*---- +* Z-ORIENTED COUPLINGS +*---- + NELEH=(IELEM+1)*IELEM**2 + IIMAZ=MUZ(LL4Z) + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 340 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 340 + NUM=NUM+1 + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 340 + VOL0=REAL(TTTT*ZZ(1,IPERT(KEL))*FRZ(KEL)) + DINV=XSGD(IBM,3) + DO 292 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 291 K2=0,IELEM-1 + DO 290 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) THEN + KEY=MUZ(INZ1) + AZ(KEY)=AZ(KEY)-VOL0*R(1,1)*DINV-QFR(NUM,7) + ENDIF + IF(KNZ2.NE.0) THEN + KEY=MUZ(INZ2) + AZ(KEY)=AZ(KEY)-VOL0*R(IELEM+1,IELEM+1)*DINV-QFR(NUM,8) + ENDIF + IF((ICOL.NE.2).AND.(KNZ1.NE.0).AND.(KNZ2.NE.0)) THEN + IF(INZ2.GT.INZ1) KEY=MUZ(INZ2)-INZ2+INZ1 + IF(INZ2.LE.INZ1) KEY=MUZ(INZ1)-INZ1+INZ2 + SG=REAL(SIGN(1,KNZ1)*SIGN(1,KNZ2)) + IF(INZ1.EQ.INZ2) SG=2.0*SG + AZ(KEY)=AZ(KEY)-SG*VOL0*R(IELEM+1,1)*DINV + ENDIF + 290 CONTINUE + 291 CONTINUE + 292 CONTINUE + 340 CONTINUE +* + DO 350 I0=1,IIMAZ + C11Z(I0)=-AZ(I0) + 350 CONTINUE + MUIM1=0 + DO 380 I=1,LL4Z + MUI=MUZ(I) + DO 370 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 365 I0=1,2*IELEM + II=IPBBZ(I0,I) + IF(II.EQ.0) GO TO 370 + DO 360 J0=1,2*IELEM + JJ=IPBBZ(J0,J) + IF(II.EQ.JJ) C11Z(KEY)=C11Z(KEY)+BBZ(I0,I)*BBZ(J0,J)/TTF(II) + 360 CONTINUE + 365 CONTINUE + 370 CONTINUE + MUIM1=MUI + 380 CONTINUE + RETURN + END |