diff options
Diffstat (limited to 'Trivac/src/TRIDXX.f')
| -rwxr-xr-x | Trivac/src/TRIDXX.f | 322 |
1 files changed, 322 insertions, 0 deletions
diff --git a/Trivac/src/TRIDXX.f b/Trivac/src/TRIDXX.f new file mode 100755 index 0000000..86c6280 --- /dev/null +++ b/Trivac/src/TRIDXX.f @@ -0,0 +1,322 @@ +*DECK TRIDXX + SUBROUTINE TRIDXX(NBMIX,CYLIND,IELEM,ICOL,NEL,LL4F,LL4X,MAT,VOL, + 1 XX,YY,ZZ,DD,KN,QFR,SGD,XSGD,MUX,IPBBX,LC,R,V,BBX,TTF,AX,C11X) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of system matrices for a Thomas-Raviart (dual) finite element +* method in Cartesian 3-D diffusion approximation. +* Note: system matrices should be initialized by the calling program. +* +*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. Hebert +* +*Parameters: input +* NBMIX number of mixtures. +* CYLIND cylindrical geometry flag (set with CYLIND=.true.). +* 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). +* NEL total number of finite elements. +* LL4F number of flux components. +* 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. +* VOL volume of each element. +* XX X-directed mesh spacings. +* YY Y-directed mesh spacings. +* ZZ Z-directed mesh spacings. +* DD used with cylindrical geometry. +* KN element-ordered unknown list. +* QFR element-ordered boundary conditions. +* 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. +* MUX X-directed compressed storage mode indices. +* MUY Y-directed compressed storage mode indices. +* MUZ Z-directed compressed storage mode 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. +* BBX X-directed flux-current matrices. +* BBY Y-directed flux-current matrices. +* BBZ Z-directed flux-current matrices. +* +*Parameters: output +* TTF flux-flux matrices. +* 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). +* 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,IELEM,ICOL,NEL,LL4F,LL4X,MAT(NEL), + 1 KN(NEL*(1+6*IELEM**2)),MUX(LL4X),IPBBX(2*IELEM,LL4X),LC + REAL VOL(NEL),XX(NEL),YY(NEL),ZZ(NEL),DD(NEL),QFR(6*NEL), + 1 SGD(NBMIX,4),XSGD(NBMIX,4),R(LC,LC),V(LC,LC-1),TTF(LL4F), + 2 BBX(2*IELEM,LL4X),AX(*),C11X(*) + LOGICAL CYLIND +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION FFF + REAL QQ(5,5) +*---- +* X-ORIENTED COUPLINGS +*---- + IF((CYLIND).AND.((IELEM.GT.1).OR.(ICOL.NE.2))) + 1 CALL XABORT('TRIDXX: TYPE OF DISCRETIZATION NOT IMPLEMENTED.') + 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 +* + NUM1=0 + NUM2=0 + DO 60 IE=1,NEL + L=MAT(IE) + IF(L.EQ.0) GO TO 60 + VOL0=VOL(IE) + IF(VOL0.EQ.0.0) GO TO 50 + DX=XX(IE) + DY=YY(IE) + DZ=ZZ(IE) + IF(CYLIND) THEN + DIN=1.0-0.5*DX/DD(IE) + DOT=1.0+0.5*DX/DD(IE) + ELSE + DIN=1.0 + DOT=1.0 + ENDIF +* + DO 45 K3=0,IELEM-1 + DO 40 K2=0,IELEM-1 + KN1=KN(NUM1+2+K3*IELEM+K2) + KN2=KN(NUM1+2+IELEM**2+K3*IELEM+K2) + INX1=ABS(KN1)-LL4F + INX2=ABS(KN2)-LL4F + DO 30 K1=0,IELEM-1 + JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1 + TTF(JND1)=TTF(JND1)+VOL0*SGD(L,4)+VOL0*QQ(K1+1,K1+1)*SGD(L,1)/ + 1 (DX*DX) + TTF(JND1)=TTF(JND1)+VOL0*QQ(K2+1,K2+1)*SGD(L,2)/(DY*DY) + TTF(JND1)=TTF(JND1)+VOL0*QQ(K3+1,K3+1)*SGD(L,3)/(DZ*DZ) + 30 CONTINUE + IF(KN1.NE.0) THEN + KEY=MUX(INX1) + AX(KEY)=AX(KEY)-DIN*(VOL0*R(1,1)*XSGD(L,1)+QFR(NUM2+1)) + ENDIF + IF(KN2.NE.0) THEN + KEY=MUX(INX2) + AX(KEY)=AX(KEY)-DOT*(VOL0*R(IELEM+1,IELEM+1)*XSGD(L,1) + 1 +QFR(NUM2+2)) + ENDIF + IF((ICOL.NE.2).AND.(KN1.NE.0).AND.(KN2.NE.0)) THEN + IF(INX2.GT.INX1) KEY=MUX(INX2)-INX2+INX1 + IF(INX2.LE.INX1) KEY=MUX(INX1)-INX1+INX2 + SG=REAL(SIGN(1,KN1)*SIGN(1,KN2)) + IF(INX1.EQ.INX2) SG=2.0*SG + AX(KEY)=AX(KEY)-SG*VOL0*R(IELEM+1,1)*XSGD(L,1) + ENDIF + 40 CONTINUE + 45 CONTINUE + 50 NUM1=NUM1+1+6*IELEM**2 + NUM2=NUM2+6 + 60 CONTINUE +* + DO 121 I0=1,MUX(LL4X) + C11X(I0)=-AX(I0) + 121 CONTINUE + MUIM1=0 + DO 716 I=1,LL4X + MUI=MUX(I) + DO 715 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 714 I0=1,2*IELEM + II=IPBBX(I0,I) + IF(II.EQ.0) GO TO 715 + DO 713 J0=1,2*IELEM + JJ=IPBBX(J0,J) + IF(II.EQ.JJ) C11X(KEY)=C11X(KEY)+BBX(I0,I)*BBX(J0,J)/TTF(II) + 713 CONTINUE + 714 CONTINUE + 715 CONTINUE + MUIM1=MUI + 716 CONTINUE + RETURN + END +* + SUBROUTINE TRIDXY(NBMIX,IELEM,ICOL,NEL,LL4F,LL4X,LL4Y,MAT,VOL,YY, + 1 KN,QFR,XSGD,MUY,IPBBY,LC,R,BBY,TTF,AY,C11Y) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,IELEM,ICOL,NEL,LL4F,LL4X,LL4Y,MAT(NEL), + 1 KN(NEL*(1+6*IELEM**2)),MUY(LL4Y),IPBBY(2*IELEM,LL4Y),LC + REAL VOL(NEL),YY(NEL),QFR(6*NEL),XSGD(NBMIX,4),R(LC,LC),TTF(LL4F), + 1 BBY(2*IELEM,LL4Y),AY(*),C11Y(*) +*---- +* Y-ORIENTED COUPLINGS +*---- + NUM1=0 + NUM2=0 + DO 240 IE=1,NEL + L=MAT(IE) + IF(L.EQ.0) GO TO 240 + VOL0=VOL(IE) + IF(VOL0.EQ.0.0) GO TO 230 + DY=YY(IE) +* + DO 195 K3=0,IELEM-1 + DO 190 K1=0,IELEM-1 + KN1=KN(NUM1+2+2*IELEM**2+K3*IELEM+K1) + KN2=KN(NUM1+2+3*IELEM**2+K3*IELEM+K1) + INY1=ABS(KN1)-LL4F-LL4X + INY2=ABS(KN2)-LL4F-LL4X + IF(KN1.NE.0) THEN + KEY=MUY(INY1) + AY(KEY)=AY(KEY)-VOL0*R(1,1)*XSGD(L,2)-QFR(NUM2+3) + ENDIF + IF(KN2.NE.0) THEN + KEY=MUY(INY2) + AY(KEY)=AY(KEY)-VOL0*R(IELEM+1,IELEM+1)*XSGD(L,2) + 1 -QFR(NUM2+4) + ENDIF + IF((ICOL.NE.2).AND.(KN1.NE.0).AND.(KN2.NE.0)) THEN + IF(INY2.GT.INY1) KEY=MUY(INY2)-INY2+INY1 + IF(INY2.LE.INY1) KEY=MUY(INY1)-INY1+INY2 + SG=REAL(SIGN(1,KN1)*SIGN(1,KN2)) + IF(INY1.EQ.INY2) SG=2.0*SG + AY(KEY)=AY(KEY)-SG*VOL0*R(IELEM+1,1)*XSGD(L,2) + ENDIF + 190 CONTINUE + 195 CONTINUE + 230 NUM1=NUM1+1+6*IELEM**2 + NUM2=NUM2+6 + 240 CONTINUE +* + DO 212 I0=1,MUY(LL4Y) + C11Y(I0)=-AY(I0) + 212 CONTINUE + MUIM1=0 + DO 216 I=1,LL4Y + MUI=MUY(I) + DO 215 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 214 I0=1,2*IELEM + II=IPBBY(I0,I) + IF(II.EQ.0) GO TO 215 + DO 213 J0=1,2*IELEM + JJ=IPBBY(J0,J) + IF(II.EQ.JJ) C11Y(KEY)=C11Y(KEY)+BBY(I0,I)*BBY(J0,J)/TTF(II) + 213 CONTINUE + 214 CONTINUE + 215 CONTINUE + MUIM1=MUI + 216 CONTINUE + RETURN + END +* + SUBROUTINE TRIDXZ(NBMIX,IELEM,ICOL,NEL,LL4F,LL4X,LL4Y,LL4Z,MAT, + 1 VOL,ZZ,KN,QFR,XSGD,MUZ,IPBBZ,LC,R,BBZ,TTF,AZ,C11Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,IELEM,ICOL,NEL,LL4F,LL4X,LL4Y,LL4Z,MAT(NEL), + 1 KN(NEL*(1+6*IELEM**2)),MUZ(LL4Z),IPBBZ(2*IELEM,LL4Z),LC + REAL VOL(NEL),ZZ(NEL),QFR(6*NEL),XSGD(NBMIX,4),R(LC,LC),TTF(LL4F), + 1 BBZ(2*IELEM,LL4Z),AZ(*),C11Z(*) +*---- +* Z-ORIENTED COUPLINGS +*---- + NUM1=0 + NUM2=0 + DO 340 IE=1,NEL + L=MAT(IE) + IF(L.EQ.0) GO TO 340 + VOL0=VOL(IE) + IF(VOL0.EQ.0.0) GO TO 330 + DZ=ZZ(IE) +* + DO 295 K2=0,IELEM-1 + DO 290 K1=0,IELEM-1 + KN1=KN(NUM1+2+4*IELEM**2+K2*IELEM+K1) + KN2=KN(NUM1+2+5*IELEM**2+K2*IELEM+K1) + INZ1=ABS(KN1)-LL4F-LL4X-LL4Y + INZ2=ABS(KN2)-LL4F-LL4X-LL4Y + IF(KN1.NE.0) THEN + KEY=MUZ(INZ1) + AZ(KEY)=AZ(KEY)-VOL0*R(1,1)*XSGD(L,3)-QFR(NUM2+5) + ENDIF + IF(KN2.NE.0) THEN + KEY=MUZ(INZ2) + AZ(KEY)=AZ(KEY)-VOL0*R(IELEM+1,IELEM+1)*XSGD(L,3) + 1 -QFR(NUM2+6) + ENDIF + IF((ICOL.NE.2).AND.(KN1.NE.0).AND.(KN2.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,KN1)*SIGN(1,KN2)) + IF(INZ1.EQ.INZ2) SG=2.0*SG + AZ(KEY)=AZ(KEY)-SG*VOL0*R(IELEM+1,1)*XSGD(L,3) + ENDIF + 290 CONTINUE + 295 CONTINUE + 330 NUM1=NUM1+1+6*IELEM**2 + NUM2=NUM2+6 + 340 CONTINUE +* + DO 312 I0=1,MUZ(LL4Z) + C11Z(I0)=-AZ(I0) + 312 CONTINUE + MUIM1=0 + DO 316 I=1,LL4Z + MUI=MUZ(I) + DO 315 J=I-(MUI-MUIM1)+1,I + KEY=MUI-I+J + DO 314 I0=1,2*IELEM + II=IPBBZ(I0,I) + IF(II.EQ.0) GO TO 315 + DO 313 J0=1,2*IELEM + JJ=IPBBZ(J0,J) + IF(II.EQ.JJ) C11Z(KEY)=C11Z(KEY)+BBZ(I0,I)*BBZ(J0,J)/TTF(II) + 313 CONTINUE + 314 CONTINUE + 315 CONTINUE + MUIM1=MUI + 316 CONTINUE + RETURN + END |