Initial commit from Polytechnique Montreal

1 files changed, 384 insertions, 0 deletions

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