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/PNFH3E.f | 384 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 384 insertions(+) create mode 100755 Trivac/src/PNFH3E.f (limited to 'Trivac/src/PNFH3E.f') diff --git a/Trivac/src/PNFH3E.f b/Trivac/src/PNFH3E.f new file mode 100755 index 0000000..622d895 --- /dev/null +++ b/Trivac/src/PNFH3E.f @@ -0,0 +1,384 @@ +*DECK PNFH3E + SUBROUTINE PNFH3E (IL,NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,L4,LL4F, + 1 MAT,SIGTI,SIDE,ZZ,FRZ,QFR,IPERT,KN,LC,R,V,SUNKNO,FUNKNO) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Perform a one-group SPN flux iteration in hexagonal 3D geometry. +* Raviart-Thomas-Schneider method in hexagonal geometry. +* +*Copyright: +* Copyright (C) 2009 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 +* IL current Legendre order. +* 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); =4 (quartic). +* ICOL type of quadrature: =1 (analytical integration); +* =2 (Gauss-Lobatto); =3 (Gauss-Legendre). +* NLF number of Legendre orders for the flux (even number). +* NVD type of void boundary condition if NLF>0 and ICOL=3. +* NAN number of Legendre orders for the cross sections. +* L4 number of unknowns per energy group and per set of two +* Legendre orders. +* LL4F number of flux components. +* MAT index-number of the mixture type assigned to each volume. +* SIGTI inverse macroscopic cross sections ordered by mixture. +* SIGTI(:,NAN) generally contains the inverse total cross +* section only. +* 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. +* LC order of the unit matrices. +* R unit Cartesian mass matrix. +* V unit nodal coupling matrix. +* SUNKNO sources. +* FUNKNO initial fluxes. +* +*Parameters: output +* FUNKNO right-hand-side of the linear system. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IL,NBMIX,NBLOS,IELEM,ICOL,NLF,NVD,NAN,L4,LL4F, + 1 MAT(3,NBLOS),IPERT(NBLOS),KN(NBLOS,3+6*(IELEM+2)*IELEM**2),LC + REAL SIGTI(NBMIX,NAN),SIDE,ZZ(3,NBLOS),FRZ(NBLOS),QFR(NBLOS,8), + 1 R(LC,LC),V(LC,LC-1),SUNKNO(L4*NLF/2),FUNKNO(L4*NLF/2) +*---- +* LOCAL VARIABLES +*---- + REAL QQ(5,5) + DOUBLE PRECISION FFF,TTTT,UUUU,VOL0,GARSI,FACT,VAR1 +* + IF(ICOL.EQ.3) THEN + IF(NVD.EQ.0) THEN + NZMAR=NLF+1 + ELSE IF(NVD.EQ.1) THEN + NZMAR=NLF + ELSE IF(NVD.EQ.2) THEN + NZMAR=65 + ENDIF + ELSE + NZMAR=65 + ENDIF + DO 16 I0=1,IELEM + DO 15 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) + 15 CONTINUE + 16 CONTINUE + JOFF=(IL/2)*L4 + FACT=REAL(2*IL+1) + IF(MOD(IL,2).EQ.0) THEN + DO 20 I=1,L4 + FUNKNO(JOFF+I)=SUNKNO(JOFF+I) + 20 CONTINUE + ENDIF +*---- +* COMPUTE THE SOLUTION AT ORDER IL. +*---- + NELEH=(IELEM+1)*IELEM**2 + TTTT=0.5D0*SQRT(3.D00)*SIDE*SIDE + NUM=0 + DO 150 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 150 + NUM=NUM+1 + DZ=ZZ(1,IPERT(KEL)) + VOL0=TTTT*DZ*FRZ(KEL) + UUUU=SIDE*DZ*FRZ(KEL) + IF(MOD(IL,2).EQ.0) THEN +* EVEN PARITY EQUATION + IF(IL.GE.2) THEN + DO 34 K5=0,1 ! TWO LOZENGES PER HEXAGON + 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) + 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=JOFF+LL4F+ABS(KNW1) + INX1=JOFF+LL4F+ABS(KNX1) + INY1=JOFF+LL4F+ABS(KNY1) + DO 30 K1=0,IELEM-1 + IF(V(K2,K1+1).EQ.0.0) GO TO 30 + IF(K5.EQ.0) THEN + SSS=(-1.0)**K1 + JND1=JOFF+(((NUM-1)*IELEM+K4)*IELEM+K3)*IELEM+K1+1 + JND2=JOFF+(((KN(NUM,1)-1)*IELEM+K4)*IELEM+K3)*IELEM+K1+1 + JND3=JOFF+(((KN(NUM,2)-1)*IELEM+K4)*IELEM+K3)*IELEM+K1+1 + ELSE + SSS=1.0 + JND1=JOFF+(((KN(NUM,1)-1)*IELEM+K4)*IELEM+K1)*IELEM+K3+1 + JND2=JOFF+(((KN(NUM,2)-1)*IELEM+K4)*IELEM+K1)*IELEM+K3+1 + JND3=JOFF+(((KN(NUM,3)-1)*IELEM+K4)*IELEM+K1)*IELEM+K3+1 + ENDIF + VAR1=SSS*REAL(IL)*UUUU*V(K2,K1+1) + IF(KNW1.NE.0) THEN + SG=REAL(SIGN(1,KNW1)) + FUNKNO(JND1)=FUNKNO(JND1)-SG*REAL(VAR1)*FUNKNO(INW1-L4) + ENDIF + IF(KNX1.NE.0) THEN + SG=REAL(SIGN(1,KNX1)) + FUNKNO(JND2)=FUNKNO(JND2)-SG*REAL(VAR1)*FUNKNO(INX1-L4) + ENDIF + IF(KNY1.NE.0) THEN + SG=REAL(SIGN(1,KNY1)) + FUNKNO(JND3)=FUNKNO(JND3)-SG*REAL(VAR1)*FUNKNO(INY1-L4) + ENDIF + 30 CONTINUE + 31 CONTINUE + 32 CONTINUE + 33 CONTINUE + 34 CONTINUE + DO 43 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 42 K2=0,IELEM-1 + DO 41 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=JOFF+LL4F+ABS(KNZ1) + INZ2=JOFF+LL4F+ABS(KNZ2) + DO 40 K3=0,IELEM-1 + IF(K5.EQ.0) THEN + JND1=JOFF+((((NUM-1)*IELEM)+K3)*IELEM+K2)*IELEM+K1+1 + ELSE + JND1=JOFF+(((KN(NUM,K5)-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + ENDIF + IF(KNZ1.NE.0) THEN + SG=REAL(SIGN(1,KNZ1)) + VAR1=SG*(VOL0/DZ)*REAL(IL)*V(1,K3+1)*FUNKNO(INZ1-L4) + FUNKNO(JND1)=FUNKNO(JND1)-REAL(VAR1) + ENDIF + IF(KNZ2.NE.0) THEN + SG=REAL(SIGN(1,KNZ2)) + VAR1=SG*(VOL0/DZ)*REAL(IL)*V(IELEM+1,K3+1)* + 1 FUNKNO(INZ2-L4) + FUNKNO(JND1)=FUNKNO(JND1)-REAL(VAR1) + ENDIF + 40 CONTINUE + 41 CONTINUE + 42 CONTINUE + 43 CONTINUE + ENDIF + ELSE +* PARTIAL INVERSION OF THE ODD PARITY EQUATION. MODIFICATION +* OF THE EVEN PARITY EQUATION. + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 150 + IF(IELEM.GT.1) THEN + DO 52 K3=0,IELEM-1 + DO 51 K2=0,IELEM-1 + DO 50 K1=0,IELEM-1 + IF(QQ(K3+1,K3+1).EQ.0.0) GO TO 50 + JND1=JOFF+(NUM-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + JND2=JOFF+(KN(NUM,1)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + JND3=JOFF+(KN(NUM,2)-1)*IELEM**3+K3*IELEM**2+K2*IELEM+K1+1 + IF(IL.GE.3) THEN + GARSI=SIGTI(IBM,MIN(IL-1,NAN)) + KND1=JND1-L4 + KND2=JND2-L4 + KND3=JND3-L4 + VAR1=(REAL(IL-1)*REAL(IL-2))*VOL0*QQ(K3+1,K3+1)*GARSI + 1 /(REAL(2*IL-3)*DZ*DZ) + FUNKNO(JND1)=FUNKNO(JND1)-REAL(VAR1)*FUNKNO(KND1) + FUNKNO(JND2)=FUNKNO(JND2)-REAL(VAR1)*FUNKNO(KND2) + FUNKNO(JND3)=FUNKNO(JND3)-REAL(VAR1)*FUNKNO(KND3) + ENDIF + IF(IL.LE.NLF-3) THEN + GARSI=SIGTI(IBM,MIN(IL+1,NAN)) + KND1=JND1+L4 + KND2=JND2+L4 + KND3=JND3+L4 + VAR1=(REAL(IL)*REAL(IL+1))*VOL0*QQ(K3+1,K3+1)*GARSI + 1 /(FACT*DZ*DZ) + FUNKNO(JND1)=FUNKNO(JND1)-REAL(VAR1)*FUNKNO(KND1) + FUNKNO(JND2)=FUNKNO(JND2)-REAL(VAR1)*FUNKNO(KND2) + FUNKNO(JND3)=FUNKNO(JND3)-REAL(VAR1)*FUNKNO(KND3) + ENDIF + 50 CONTINUE + 51 CONTINUE + 52 CONTINUE + ENDIF +* +* ODD PARITY EQUATION + DO 93 K5=0,1 ! TWO LOZENGES PER HEXAGON + DO 92 K4=0,IELEM-1 + DO 91 K3=0,IELEM-1 + DO 90 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=JOFF+LL4F+ABS(KNW1) + INX1=JOFF+LL4F+ABS(KNX1) + INY1=JOFF+LL4F+ABS(KNY1) + IF(KNW1.NE.0) THEN + DO 60 IL2=1,NLF-1,2 + IF(IL2.EQ.IL) GO TO 60 + ZMARS=PNMAR2(NZMAR,IL2,IL) + INW2=(IL2/2)*L4+LL4F+ABS(KNW1) + IF((K2.EQ.1).AND.(K5.EQ.0)) THEN + VAR1=0.5*FACT*QFR(NUM,1)*ZMARS*FUNKNO(INW2) + FUNKNO(INW1)=FUNKNO(INW1)+REAL(VAR1) + ELSE IF((K2.EQ.IELEM+1).AND.(K5.EQ.1)) THEN + VAR1=0.5*FACT*QFR(NUM,2)*ZMARS*FUNKNO(INW2) + FUNKNO(INW1)=FUNKNO(INW1)+REAL(VAR1) + ENDIF + 60 CONTINUE + ENDIF + IF(KNX1.NE.0) THEN + DO 70 IL2=1,NLF-1,2 + IF(IL2.EQ.IL) GO TO 70 + ZMARS=PNMAR2(NZMAR,IL2,IL) + INX2=(IL2/2)*L4+LL4F+ABS(KNX1) + IF((K2.EQ.1).AND.(K5.EQ.0)) THEN + VAR1=0.5*FACT*QFR(NUM,3)*ZMARS*FUNKNO(INX2) + FUNKNO(INX1)=FUNKNO(INX1)+REAL(VAR1) + ELSE IF((K2.EQ.IELEM+1).AND.(K5.EQ.1)) THEN + VAR1=0.5*FACT*QFR(NUM,4)*ZMARS*FUNKNO(INX2) + FUNKNO(INX1)=FUNKNO(INX1)+REAL(VAR1) + ENDIF + 70 CONTINUE + ENDIF + IF(KNY1.NE.0) THEN + DO 80 IL2=1,NLF-1,2 + IF(IL2.EQ.IL) GO TO 80 + ZMARS=PNMAR2(NZMAR,IL2,IL) + INY2=(IL2/2)*L4+LL4F+ABS(KNY1) + IF((K2.EQ.1).AND.(K5.EQ.0)) THEN + VAR1=0.5*FACT*QFR(NUM,5)*ZMARS*FUNKNO(INY2) + FUNKNO(INY1)=FUNKNO(INY1)+REAL(VAR1) + ELSE IF((K2.EQ.IELEM+1).AND.(K5.EQ.1)) THEN + VAR1=0.5*FACT*QFR(NUM,6)*ZMARS*FUNKNO(INY2) + FUNKNO(INY1)=FUNKNO(INY1)+REAL(VAR1) + ENDIF + 80 CONTINUE + ENDIF + 90 CONTINUE + 91 CONTINUE + 92 CONTINUE + 93 CONTINUE + DO 122 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 121 K2=0,IELEM-1 + DO 120 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=JOFF+LL4F+ABS(KNZ1) + INZ2=JOFF+LL4F+ABS(KNZ2) + IF((QFR(NUM,7).NE.0.0).AND.(KNZ1.NE.0)) THEN +* ZINF SIDE. + DO 100 IL2=1,NLF-1,2 + IF(IL2.EQ.IL) GO TO 100 + ZMARS=PNMAR2(NZMAR,IL2,IL) + INDL=(IL2/2)*L4+LL4F+ABS(KNZ1) + VAR1=0.5*FACT*QFR(NUM,7)*ZMARS*FUNKNO(INDL) + FUNKNO(INZ1)=FUNKNO(INZ1)+REAL(VAR1) + 100 CONTINUE + ENDIF + IF((QFR(NUM,8).NE.0.0).AND.(KNZ2.NE.0)) THEN +* ZSUP SIDE. + DO 110 IL2=1,NLF-1,2 + IF(IL2.EQ.IL) GO TO 110 + ZMARS=PNMAR2(NZMAR,IL2,IL) + INDL=(IL2/2)*L4+LL4F+ABS(KNZ2) + VAR1=0.5*FACT*QFR(NUM,8)*ZMARS*FUNKNO(INDL) + FUNKNO(INZ2)=FUNKNO(INZ2)+REAL(VAR1) + 110 CONTINUE + ENDIF + 120 CONTINUE + 121 CONTINUE + 122 CONTINUE +* + IF(IL.LE.NLF-3) THEN + DO 134 K5=0,1 ! TWO LOZENGES PER HEXAGON + DO 133 K4=0,IELEM-1 + DO 132 K3=0,IELEM-1 + DO 131 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=JOFF+LL4F+ABS(KNW1) + INX1=JOFF+LL4F+ABS(KNX1) + INY1=JOFF+LL4F+ABS(KNY1) + DO 130 K1=0,IELEM-1 + IF(V(K2,K1+1).EQ.0.0) GO TO 130 + IF(K5.EQ.0) THEN + SSS=(-1.0)**K1 + JND1=JOFF+(((NUM-1)*IELEM+K4)*IELEM+K3)*IELEM+K1+1 + JND2=JOFF+(((KN(NUM,1)-1)*IELEM+K4)*IELEM+K3)*IELEM+K1+1 + JND3=JOFF+(((KN(NUM,2)-1)*IELEM+K4)*IELEM+K3)*IELEM+K1+1 + ELSE + SSS=1.0 + JND1=JOFF+(((KN(NUM,1)-1)*IELEM+K4)*IELEM+K1)*IELEM+K3+1 + JND2=JOFF+(((KN(NUM,2)-1)*IELEM+K4)*IELEM+K1)*IELEM+K3+1 + JND3=JOFF+(((KN(NUM,3)-1)*IELEM+K4)*IELEM+K1)*IELEM+K3+1 + ENDIF + VAR1=SSS*REAL(IL+1)*UUUU*V(K2,K1+1) + IF(KNW1.NE.0) THEN + SG=REAL(SIGN(1,KNW1)) + FUNKNO(INW1)=FUNKNO(INW1)-SG*REAL(VAR1)*FUNKNO(JND1+L4) + ENDIF + IF(KNX1.NE.0) THEN + SG=REAL(SIGN(1,KNX1)) + FUNKNO(INX1)=FUNKNO(INX1)-SG*REAL(VAR1)*FUNKNO(JND2+L4) + ENDIF + IF(KNY1.NE.0) THEN + SG=REAL(SIGN(1,KNY1)) + FUNKNO(INY1)=FUNKNO(INY1)-SG*REAL(VAR1)*FUNKNO(JND3+L4) + ENDIF + 130 CONTINUE + 131 CONTINUE + 132 CONTINUE + 133 CONTINUE + 134 CONTINUE + DO 143 K5=0,2 ! THREE LOZENGES PER HEXAGON + DO 142 K2=0,IELEM-1 + DO 141 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=JOFF+LL4F+ABS(KNZ1) + INZ2=JOFF+LL4F+ABS(KNZ2) + DO 140 K3=0,IELEM-1 + IF(K5.EQ.0) THEN + JND1=JOFF+((((NUM-1)*IELEM)+K3)*IELEM+K2)*IELEM+K1+1 + ELSE + JND1=JOFF+(((KN(NUM,K5)-1)*IELEM+K3)*IELEM+K2)*IELEM+K1+1 + ENDIF + IF(KNZ1.NE.0) THEN + SG=REAL(SIGN(1,KNZ1)) + VAR1=SG*(VOL0/DZ)*REAL(IL+1)*V(1,K3+1)*FUNKNO(JND1+L4) + FUNKNO(INZ1)=FUNKNO(INZ1)-REAL(VAR1) + ENDIF + IF(KNZ2.NE.0) THEN + SG=REAL(SIGN(1,KNZ2)) + VAR1=SG*(VOL0/DZ)*REAL(IL+1)*V(IELEM+1,K3+1)* + 1 FUNKNO(JND1+L4) + FUNKNO(INZ2)=FUNKNO(INZ2)-REAL(VAR1) + ENDIF + 140 CONTINUE + 141 CONTINUE + 142 CONTINUE + 143 CONTINUE + ENDIF + ENDIF + 150 CONTINUE + RETURN + END -- cgit v1.2.3