Initial commit from Polytechnique Montreal

1 files changed, 225 insertions, 0 deletions

diff --git a/Trivac/src/PNFH2E.f b/Trivac/src/PNFH2E.f
new file mode 100755
index 0000000..ee0998e
--- /dev/null
+++ b/Trivac/src/PNFH2E.f
@@ -0,0 +1,225 @@
+*DECK PNFH2E
+      SUBROUTINE PNFH2E (IELEM,ICOL,NBLOS,SIDE,NLF,NVD,L4,IPERT,KN,
+     1 QFR,MU,IIMAX,LC,V,SYS,SUNKNO,FUNKNO,NADI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform a one-group SPN flux iteration in hexagonal 2D 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
+* 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).
+* NBLOS   number of lozenges per direction, taking into account
+*         mesh-splitting.
+* SIDE    side of the hexagons.
+* NLF     number of Legendre orders for the flux (even number).
+* NVD     type of void boundary condition if NLF>0 and ICOL=3.
+* L4      number of unknowns per energy group and per set of two
+*         Legendre orders.
+* IPERT   mixture permutation index.
+* KN      element-ordered unknown list.
+* QFR     element-ordered boundary conditions.
+* MU      profiled storage indices for matrix SYS.
+* IIMAX   dimension of SYS.
+* LC      order of the unit matrices.
+* V       unit nodal coupling matrix.
+* SYS     LU factors of the system matrix.
+* SUNKNO  sources.
+* FUNKNO  initial fluxes.
+* NADI    number of inner ADI iterations.
+*
+*Parameters: output
+* FUNKNO  fluxes.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IELEM,ICOL,NBLOS,NLF,NVD,L4,IPERT(NBLOS),
+     1 KN(NBLOS,4+6*IELEM*(IELEM+1)),MU(L4),IIMAX,LC,NADI
+      REAL SIDE,QFR(NBLOS,6),V(LC,LC-1),SYS(IIMAX),SUNKNO(L4*NLF/2),
+     1 FUNKNO(L4*NLF/2)
+*
+      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
+      MUMAX=MU(L4)
+      NELEM=IELEM*(IELEM+1)
+      DO 170 IADI=1,MAX(1,NADI)
+      DO 160 IL=0,NLF-1
+      FACT=REAL(2*IL+1)
+      IF(MOD(IL,2).EQ.0) THEN
+         DO 10 I=1,L4
+         FUNKNO((IL/2)*L4+I)=SUNKNO((IL/2)*L4+I)
+   10    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE SOLUTION AT ORDER IL.
+*----
+      NUM=0
+      DO 150 KEL=1,NBLOS
+      IF(IPERT(KEL).EQ.0) GO TO 150
+      NUM=NUM+1
+      IF(MOD(IL,2).EQ.0) THEN
+*        EVEN PARITY EQUATION
+         IF(IL.GE.2) THEN
+            DO 33 K4=0,1
+            DO 32 K3=0,IELEM-1
+            DO 31 K2=1,IELEM+1
+            KNW1=KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K2)
+            KNX1=KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K2)
+            KNY1=KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K2)
+            INW1=((IL-2)/2)*L4+ABS(KNW1)
+            INX1=((IL-2)/2)*L4+ABS(KNX1)
+            INY1=((IL-2)/2)*L4+ABS(KNY1)
+            DO 30 K1=0,IELEM-1
+            IF(V(K2,K1+1).EQ.0.0) GO TO 30
+            IF(K4.EQ.0) THEN
+               SSS=(-1.0)**K1
+               JND1=(IL/2)*L4+KN(NUM,1)+K3*IELEM+K1
+               JND2=(IL/2)*L4+KN(NUM,2)+K3*IELEM+K1
+               JND3=(IL/2)*L4+KN(NUM,3)+K3*IELEM+K1
+            ELSE
+               SSS=1.0
+               JND1=(IL/2)*L4+KN(NUM,2)+K1*IELEM+K3
+               JND2=(IL/2)*L4+KN(NUM,3)+K1*IELEM+K3
+               JND3=(IL/2)*L4+KN(NUM,4)+K1*IELEM+K3
+            ENDIF
+            IF(KNW1.NE.0) THEN
+               SG=REAL(SIGN(1,KNW1))
+               FUNKNO(JND1)=FUNKNO(JND1)-SG*SSS*REAL(IL)*SIDE*
+     1         V(K2,K1+1)*FUNKNO(INW1)
+            ENDIF
+            IF(KNX1.NE.0) THEN
+               SG=REAL(SIGN(1,KNX1))
+               FUNKNO(JND2)=FUNKNO(JND2)-SG*SSS*REAL(IL)*SIDE*
+     1         V(K2,K1+1)*FUNKNO(INX1)
+            ENDIF
+            IF(KNY1.NE.0) THEN
+               SG=REAL(SIGN(1,KNY1))
+               FUNKNO(JND3)=FUNKNO(JND3)-SG*SSS*REAL(IL)*SIDE*
+     1         V(K2,K1+1)*FUNKNO(INY1)
+            ENDIF
+   30       CONTINUE
+   31       CONTINUE
+   32       CONTINUE
+   33       CONTINUE
+         ENDIF
+      ELSE
+*        ODD PARITY EQUATION
+         DO 142 K4=0,1
+         DO 141 K3=0,IELEM-1
+         DO 140 K2=1,IELEM+1
+         KNW1=KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K2)
+         KNX1=KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K2)
+         KNY1=KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K2)
+         INW1=(IL/2)*L4+ABS(KNW1)
+         INX1=(IL/2)*L4+ABS(KNX1)
+         INY1=(IL/2)*L4+ABS(KNY1)
+         IF(KNW1.NE.0) THEN
+            DO 90 IL2=1,NLF-1,2
+            IF(IL2.EQ.IL) GO TO 90
+            ZMARS=PNMAR2(NZMAR,IL2,IL)
+            INW2=(IL2/2)*L4+ABS(KNW1)
+            IF((K2.EQ.1).AND.(K4.EQ.0)) THEN
+               FUNKNO(INW1)=FUNKNO(INW1)+0.5*FACT*QFR(NUM,1)*ZMARS*
+     1         FUNKNO(INW2)
+            ELSE IF((K2.EQ.IELEM+1).AND.(K4.EQ.1)) THEN
+               FUNKNO(INW1)=FUNKNO(INW1)+0.5*FACT*QFR(NUM,2)*ZMARS*
+     1         FUNKNO(INW2)
+            ENDIF
+   90       CONTINUE
+         ENDIF
+         IF(KNX1.NE.0) THEN
+            DO 100 IL2=1,NLF-1,2
+            IF(IL2.EQ.IL) GO TO 100
+            ZMARS=PNMAR2(NZMAR,IL2,IL)
+            INX2=(IL2/2)*L4+ABS(KNX1)
+            IF((K2.EQ.1).AND.(K4.EQ.0)) THEN
+               FUNKNO(INX1)=FUNKNO(INX1)+0.5*FACT*QFR(NUM,3)*ZMARS*
+     1         FUNKNO(INX2)
+            ELSE IF((K2.EQ.IELEM+1).AND.(K4.EQ.1)) THEN
+               FUNKNO(INX1)=FUNKNO(INX1)+0.5*FACT*QFR(NUM,4)*ZMARS*
+     1         FUNKNO(INX2)
+            ENDIF
+  100       CONTINUE
+         ENDIF
+         IF(KNY1.NE.0) THEN
+            DO 110 IL2=1,NLF-1,2
+            IF(IL2.EQ.IL) GO TO 110
+            ZMARS=PNMAR2(NZMAR,IL2,IL)
+            INY2=(IL2/2)*L4+ABS(KNY1)
+            IF((K2.EQ.1).AND.(K4.EQ.0)) THEN
+               FUNKNO(INY1)=FUNKNO(INY1)+0.5*FACT*QFR(NUM,5)*ZMARS*
+     1         FUNKNO(INY2)
+            ELSE IF((K2.EQ.IELEM+1).AND.(K4.EQ.1)) THEN
+               FUNKNO(INY1)=FUNKNO(INY1)+0.5*FACT*QFR(NUM,6)*ZMARS*
+     1         FUNKNO(INY2)
+            ENDIF
+  110       CONTINUE
+         ENDIF
+         IF(IL.LE.NLF-3) THEN
+            DO 130 K1=0,IELEM-1
+            IF(V(K2,K1+1).EQ.0.0) GO TO 130
+            IF(K4.EQ.0) THEN
+               SSS=(-1.0)**K1
+               JND1=((IL+2)/2)*L4+KN(NUM,1)+K3*IELEM+K1
+               JND2=((IL+2)/2)*L4+KN(NUM,2)+K3*IELEM+K1
+               JND3=((IL+2)/2)*L4+KN(NUM,3)+K3*IELEM+K1
+            ELSE
+               SSS=1.0
+               JND1=((IL+2)/2)*L4+KN(NUM,2)+K1*IELEM+K3
+               JND2=((IL+2)/2)*L4+KN(NUM,3)+K1*IELEM+K3
+               JND3=((IL+2)/2)*L4+KN(NUM,4)+K1*IELEM+K3
+            ENDIF
+            IF(KNW1.NE.0) THEN
+               SG=REAL(SIGN(1,KNW1))
+               FUNKNO(INW1)=FUNKNO(INW1)-SG*SSS*REAL(IL+1)*SIDE*
+     1         V(K2,K1+1)*FUNKNO(JND1)
+            ENDIF
+            IF(KNX1.NE.0) THEN
+               SG=REAL(SIGN(1,KNX1))
+               FUNKNO(INX1)=FUNKNO(INX1)-SG*SSS*REAL(IL+1)*SIDE*
+     1         V(K2,K1+1)*FUNKNO(JND2)
+            ENDIF
+            IF(KNY1.NE.0) THEN
+               SG=REAL(SIGN(1,KNY1))
+               FUNKNO(INY1)=FUNKNO(INY1)-SG*SSS*REAL(IL+1)*SIDE*
+     1         V(K2,K1+1)*FUNKNO(JND3)
+            ENDIF
+  130       CONTINUE
+         ENDIF
+  140    CONTINUE
+  141    CONTINUE
+  142    CONTINUE
+      ENDIF
+  150 CONTINUE
+      IF(MOD(IL,2).EQ.1) THEN
+         CALL ALLDLS(L4,MU,SYS((IL/2)*MUMAX+1),FUNKNO((IL/2)*L4+1))
+      ENDIF
+  160 CONTINUE
+  170 CONTINUE
+      RETURN
+      END