diff options
Diffstat (limited to 'Trivac/src/TRIMTW.f')
| -rwxr-xr-x | Trivac/src/TRIMTW.f | 383 |
1 files changed, 383 insertions, 0 deletions
diff --git a/Trivac/src/TRIMTW.f b/Trivac/src/TRIMTW.f new file mode 100755 index 0000000..5fd6d9b --- /dev/null +++ b/Trivac/src/TRIMTW.f @@ -0,0 +1,383 @@ +*DECK TRIMTW + SUBROUTINE TRIMTW(ISPLH,IR,NEL,LL4,VOL,MAT,MATN,SGD,XSGD,SIDE, + 1 ZZ,KN,QFR,MUW,IPW,IPR,A11W) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of system matrices for a mesh centered finite difference +* discretization in hexagonal geometry (triangular sub meshs). +* 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. Benaboud +* +*Parameters: input +* ISPLH used to compute the number of triangles as 6*(ISPLH-1)**2. +* IR first dimension of matrix SGD. +* NEL total number of finite elements. +* LL4 order of system matrices. +* VOL volume of each element. +* MAT mixture index assigned to each hexagon. +* MATN mixture index assigned to each triangle. +* SGD nuclear properties per material mixtures: +* SGD(L,1): W-, X-, and Y-oriented diffusion coefficients; +* SGD(L,3): Z-oriented diffusion coefficients; +* SGD(L,4): removal macroscopic cross section. +* XSGD nuclear properties (IPR=0), derivatives (IPR=1) or first +* variations (IPR=2 or 3) of nuclear properties per material +* mixture. +* SIDE side of an hexagon. +* ZZ Z-directed mesh spacings. +* KN element-ordered unknown list. +* QFR element-ordered boundary conditions. +* 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. +* IPR type of calculation: +* =0: compute the system matrices; +* =1: compute the derivative of system matrices; +* =2 or =3: compute the variation of system matrices. +* +*Parameters: output +* A11W W-oriented matrices corresponding to the divergence (i.e +* leakage) and removal terms. Dimensionned to MUW(LL4). +* A11X X-oriented matrices corresponding to the divergence (i.e +* leakage) and removal terms. Dimensionned to MUX(LL4). +* A11Y Y-oriented matrices corresponding to the divergence (i.e +* leakage) and removal terms. Dimensionned to MUY(LL4). +* A11Z Z-oriented matrices corresponding to the divergence (i.e +* leakage) and removal terms. Dimensionned to MUZ(LL4). +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ISPLH,IR,NEL,LL4,MAT(NEL),MATN(LL4), + 1 KN((18*(ISPLH-1)**2+3)*NEL),MUW(LL4),IPW(LL4),IPR + REAL VOL(NEL),SGD(IR,4),XSGD(IR,4),SIDE,ZZ(NEL),QFR(8*NEL), + 1 A11W(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION A1(5),VAR1 + INTEGER, DIMENSION(:), ALLOCATABLE :: IWRK +*---- +* ASSEMBLY OF MATRIX A11W +*---- + NUM1 = 0 + NUM2 = 0 + NTPH = 6 * (ISPLH-1)**2 + NTPL = 1 + 2 * (ISPLH-1) + NVT1 = NTPL + 2 * (ISPLH-2) + NTPH / 2 + NVT2 = NTPH - NTPL - (ISPLH-4) * (NTPL+2) + NVT3 = NTPH - (ISPLH-4) * NTPL + IVAL = 3*NTPH+8 + IF(ISPLH.EQ.3) NVT2 = NTPH + IF(ISPLH.LE.3) ISAU = 2*(ISPLH-2) + IF(ISPLH.GE.4) ISAU = 6*(ISPLH-3) + ICR = ISAU*(1+2*(ISPLH-2)) + ALLOCATE(IWRK(NEL)) + MEL = 0 + DO 10 M=1,NEL + IF(MAT(M).LE.0) GO TO 10 + MEL = MEL + 1 + IWRK(MEL) = M + 10 CONTINUE + DO 40 K=1,NEL + L = MAT(K) + IF(L.EQ.0) GO TO 40 + VOL0 = VOL(K)/NTPH + IF(VOL0.EQ.0.0) GO TO 30 + KK4=KN(NUM1+3*NTPH+7) + KK5=KN(NUM1+3*NTPH+8) + IF(KK4.GT.0) KK4 = IWRK(KK4) + IF(KK5.GT.0) KK5 = IWRK(KK5) + DO 20 I = 1,NTPH +* + CALL TRINEI (3,1,1,ISPLH,ICR,I,KK1,KK2,KK3,KEL,IQF,NUM1, + > NTPH,NTPL,NVT1,NVT2,NVT3,IVAL,KN) +* + CALL TRITCO (NEL,LL4,ISPLH,IR,IQF,K,KK1,KK2,KK3,KK4,KK5, + > VOL0,MAT,MATN,SGD(1,1),XSGD(1,1),SIDE,ZZ,QFR(NUM2+1),IPR,A1) +* + INW1=IPW(KEL) + KEY0=MUW(INW1)-INW1 + IF(KK1.GT.0) THEN + INW2=IPW(KK1) + IF(INW2.LT.INW1) THEN + KEY=KEY0+INW2 + A11W(KEY)=A11W(KEY)-REAL(A1(1))/2. + ENDIF + ENDIF + IF(KK2.GT.0) THEN + INW2=IPW(KK2) + IF(INW2.LT.INW1) THEN + KEY=KEY0+INW2 + A11W(KEY)=A11W(KEY)-REAL(A1(2))/2. + ENDIF + ENDIF + KEY=KEY0+INW1 + VAR1 = A1(1)+A1(2)+A1(3)+A1(4)+A1(5) + A11W(KEY)=A11W(KEY)+REAL(VAR1)+XSGD(L,4)*VOL0 + 20 CONTINUE + 30 NUM1=NUM1+IVAL + NUM2=NUM2+8 + 40 CONTINUE + DEALLOCATE(IWRK) + RETURN + END +* + SUBROUTINE TRIMTX (ISPLH,IR,NEL,LL4,VOL,MAT,MATN,SGD,XSGD,SIDE, + 1 ZZ,KN,QFR,MUX,IPX,IPR,A11X) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ISPLH,IR,NEL,LL4,MAT(NEL),MATN(LL4), + 1 KN((18*(ISPLH-1)**2+3)*NEL),MUX(LL4),IPX(LL4),IPR + REAL VOL(NEL),SGD(IR,4),XSGD(IR,4),SIDE,ZZ(NEL),QFR(8*NEL), + 1 A11X(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION A1(5),VAR1 + INTEGER, DIMENSION(:), ALLOCATABLE :: IWRK +*---- +* ASSEMBLY OF MATRIX A11X +*---- + NUM1=0 + NUM2=0 + NTPH = 6*(ISPLH-1)**2 + NTPL = 1+2*(ISPLH-1) + NVT1 = NTPL + 2 * (ISPLH-2) + NTPH / 2 + NVT2 = NTPH - NTPL - (ISPLH-4) * (NTPL+2) + NVT3 = NTPH - (ISPLH-4) * NTPL + IVAL = 3*NTPH+8 + IF(ISPLH.EQ.3) NVT2 = NTPH + IF(ISPLH.LE.3) ISAU = 2*(ISPLH-2) + IF(ISPLH.GE.4) ISAU = 6*(ISPLH-3) + ICR = ISAU*(1+2*(ISPLH-2)) + ALLOCATE(IWRK(NEL)) + MEL = 0 + DO 105 M=1,NEL + IF(MAT(M).LE.0) GO TO 105 + MEL = MEL + 1 + IWRK(MEL) = M +105 CONTINUE + DO 130 K=1,NEL + L = MAT(K) + IF(L.EQ.0) GO TO 130 + VOL0 = VOL(K)/NTPH + IF(VOL0.EQ.0.0) GO TO 120 + KK4=KN(NUM1+3*NTPH+7) + KK5=KN(NUM1+3*NTPH+8) + IF(KK4.GT.0) KK4 = IWRK(KK4) + IF(KK5.GT.0) KK5 = IWRK(KK5) + DO 110 I = 1,NTPH +* + CALL TRINEI (3,2,1,ISPLH,ICR,I,KK1,KK2,KK3,KEL,IQF,NUM1, + > NTPH,NTPL,NVT1,NVT2,NVT3,IVAL,KN) +* + CALL TRITCO (NEL,LL4,ISPLH,IR,IQF,K,KK1,KK2,KK3,KK4,KK5, + > VOL0,MAT,MATN,SGD(1,1),XSGD(1,1),SIDE,ZZ,QFR(NUM2+1),IPR,A1) +* + INX1=IPX(KEL) + KEY0=MUX(INX1)-INX1 + IF(KK1.GT.0) THEN + INX2=IPX(KK1) + IF(INX2.LT.INX1) THEN + KEY=KEY0+INX2 + A11X(KEY)=A11X(KEY)-REAL(A1(1))/2. + ENDIF + ENDIF + IF(KK2.GT.0) THEN + INX2=IPX(KK2) + IF(INX2.LT.INX1) THEN + KEY=KEY0+INX2 + A11X(KEY)=A11X(KEY)-REAL(A1(2))/2. + ENDIF + ENDIF + KEY=KEY0+INX1 + VAR1 = A1(1)+A1(2)+A1(3)+A1(4)+A1(5) + A11X(KEY)=A11X(KEY)+REAL(VAR1)+XSGD(L,4)*VOL0 + 110 CONTINUE + 120 NUM1=NUM1+IVAL + NUM2=NUM2+8 + 130 CONTINUE + DEALLOCATE(IWRK) + RETURN + END +* + SUBROUTINE TRIMTY (ISPLH,IR,NEL,LL4,VOL,MAT,MATN,SGD,XSGD,SIDE, + 1 ZZ,KN,QFR,MUY,IPY,IPR,A11Y) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ISPLH,IR,NEL,LL4,MAT(NEL),MATN(LL4), + 1 KN((18*(ISPLH-1)**2+3)*NEL),MUY(LL4),IPY(LL4),IPR + REAL VOL(NEL),SGD(IR,4),XSGD(IR,4),SIDE,ZZ(NEL),QFR(8*NEL), + 1 A11Y(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION A1(5),VAR1 + INTEGER, DIMENSION(:), ALLOCATABLE :: IWRK +*---- +* ASSEMBLY OF MATRIX A11Y +*---- + NUM1=0 + NUM2=0 + NTPH = 6*(ISPLH-1)**2 + NTPL = 1+2*(ISPLH-1) + NVT1 = NTPL + 2 * (ISPLH-2) + NTPH / 2 + NVT2 = NTPH - NTPL - (ISPLH-4) * (NTPL+2) + NVT3 = NTPH - (ISPLH-4) * NTPL + IVAL = 3*NTPH+8 + IF(ISPLH.EQ.3) NVT2 = NTPH + IF(ISPLH.LE.3) ISAU = 2*(ISPLH-2) + IF(ISPLH.GE.4) ISAU = 6*(ISPLH-3) + ICR = ISAU*(1+2*(ISPLH-2)) + ALLOCATE(IWRK(NEL)) + MEL = 0 + DO 205 M=1,NEL + IF(MAT(M).LE.0) GO TO 205 + MEL = MEL + 1 + IWRK(MEL) = M +205 CONTINUE + DO 230 K=1,NEL + L = MAT(K) + IF(L.EQ.0) GO TO 230 + VOL0 = VOL(K)/NTPH + IF(VOL0.EQ.0.0) GO TO 220 + KK4=KN(NUM1+3*NTPH+7) + KK5=KN(NUM1+3*NTPH+8) + IF(KK4.GT.0) KK4 = IWRK(KK4) + IF(KK5.GT.0) KK5 = IWRK(KK5) + DO 210 I = 1,NTPH +* + CALL TRINEI (3,3,1,ISPLH,ICR,I,KK1,KK2,KK3,KEL,IQF,NUM1, + > NTPH,NTPL,NVT1,NVT2,NVT3,IVAL,KN) +* + CALL TRITCO (NEL,LL4,ISPLH,IR,IQF,K,KK1,KK2,KK3,KK4,KK5, + > VOL0,MAT,MATN,SGD(1,1),XSGD(1,1),SIDE,ZZ,QFR(NUM2+1),IPR,A1) +* + INY1=IPY(KEL) + KEY0=MUY(INY1)-INY1 + IF(KK1.GT.0) THEN + INY2=IPY(KK1) + IF(INY2.LT.INY1) THEN + KEY=KEY0+INY2 + A11Y(KEY)=A11Y(KEY)-REAL(A1(1))/2. + ENDIF + ENDIF + IF(KK2.GT.0) THEN + INY2=IPY(KK2) + IF(INY2.LT.INY1) THEN + KEY=KEY0+INY2 + A11Y(KEY)=A11Y(KEY)-REAL(A1(2))/2. + ENDIF + ENDIF + KEY=KEY0+INY1 + VAR1 = A1(1)+A1(2)+A1(3)+A1(4)+A1(5) + A11Y(KEY)=A11Y(KEY)+REAL(VAR1)+XSGD(L,4)*VOL0 + 210 CONTINUE + 220 NUM1=NUM1+IVAL + NUM2=NUM2+8 + 230 CONTINUE + DEALLOCATE(IWRK) + RETURN + END +* + SUBROUTINE TRIMTZ (ISPLH,IR,NEL,LL4,VOL,MAT,MATN,SGD,XSGD,SIDE, + 1 ZZ,KN,QFR,MUZ,IPZ,IPR,A11Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ISPLH,IR,NEL,LL4,MAT(NEL),MATN(LL4), + 1 KN((18*(ISPLH-1)**2+3)*NEL),MUZ(LL4),IPZ(LL4),IPR + REAL VOL(NEL),SGD(IR,4),XSGD(IR,4),SIDE,ZZ(NEL),QFR(8*NEL), + 1 A11Z(*) +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION A1(5),VAR1 + INTEGER, DIMENSION(:), ALLOCATABLE :: IWRK +*---- +* ASSEMBLY OF MATRIX A11Z +*---- + NUM1=0 + NUM2=0 + NTPH = 6*(ISPLH-1)**2 + NTPL = 1+2*(ISPLH-1) + NVT1 = NTPL + 2 * (ISPLH-2) + NTPH / 2 + NVT2 = NTPH - NTPL - (ISPLH-4) * (NTPL+2) + NVT3 = NTPH - (ISPLH-4) * NTPL + IVAL = 3*NTPH+8 + IF(ISPLH.EQ.3) NVT2 = NTPH + IF(ISPLH.LE.3) ISAU = 2*(ISPLH-2) + IF(ISPLH.GE.4) ISAU = 6*(ISPLH-3) + ICR = ISAU*(1+2*(ISPLH-2)) + ALLOCATE(IWRK(NEL)) + MEL = 0 + DO 305 M=1,NEL + IF(MAT(M).LE.0) GO TO 305 + MEL = MEL + 1 + IWRK(MEL) = M +305 CONTINUE + DO 330 K=1,NEL + L = MAT(K) + IF(L.EQ.0) GO TO 330 + VOL0 = VOL(K)/NTPH + IF(VOL0.EQ.0.0) GO TO 320 + DO 310 I = 1,NTPH +* + CALL TRINEI (3,1,1,ISPLH,ICR,I,KK1,KK2,KK3,KEL,IQF,NUM1, + > NTPH,NTPL,NVT1,NVT2,NVT3,IVAL,KN) + KK4 = KN(NUM1+NTPH+I) + KK5 = KN(NUM1+2*NTPH+I) + LK4 = KK4 + LK5 = KK5 + IF(LK4.GT.0) LK4 = IWRK(KN(NUM1+3*NTPH+7)) + IF(LK5.GT.0) LK5 = IWRK(KN(NUM1+3*NTPH+8)) +* + CALL TRITCO (NEL,LL4,ISPLH,IR,IQF,K,KK1,KK2,KK3,LK4,LK5, + > VOL0,MAT,MATN,SGD(1,1),XSGD(1,1),SIDE,ZZ,QFR(NUM2+1),IPR,A1) +* + INZ1=IPZ(KEL) + KEY0=MUZ(INZ1)-INZ1 + IF(KK4.GT.0) THEN + INZ2=IPZ(KK4) + IF(INZ2.LT.INZ1) THEN + KEY=KEY0+INZ2 + A11Z(KEY)=A11Z(KEY)-REAL(A1(4)) + ENDIF + ENDIF + IF(KK5.GT.0) THEN + INZ2=IPZ(KK5) + IF(INZ2.LT.INZ1) THEN + KEY=KEY0+INZ2 + A11Z(KEY)=A11Z(KEY)-REAL(A1(5)) + ENDIF + ENDIF + KEY=KEY0+INZ1 + VAR1 = A1(1)+A1(2)+A1(3)+A1(4)+A1(5) + A11Z(KEY)=A11Z(KEY)+REAL(VAR1)+XSGD(L,4)*VOL0 + 310 CONTINUE + 320 NUM1=NUM1+IVAL + NUM2=NUM2+8 + 330 CONTINUE + DEALLOCATE(IWRK) + RETURN + END |