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/PNDH2E.f | 300 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 300 insertions(+) create mode 100755 Trivac/src/PNDH2E.f (limited to 'Trivac/src/PNDH2E.f') diff --git a/Trivac/src/PNDH2E.f b/Trivac/src/PNDH2E.f new file mode 100755 index 0000000..9614391 --- /dev/null +++ b/Trivac/src/PNDH2E.f @@ -0,0 +1,300 @@ +*DECK PNDH2E + SUBROUTINE PNDH2E(ITY,IELEM,ICOL,NBLOS,L4,NBMIX,IIMAX,SIDE,MAT, + 1 IPERT,SIGT,KN,QFR,NLF,NVD,NAN,MU,LC,R,V,H,SYS) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of a within-group (leakage and removal) or out-of-group +* system matrix in a Thomas-Raviart-Schneider (dual) finite element +* simplified PN method approximation (2D 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 +* ITY type of assembly: +* =0: leakage-removal matrix assembly; =1: cross section matrix +* assembly. +* 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. +* L4 number of unknowns per energy group and per set of two +* Legendre orders. +* NBMIX number of mixtures. +* IIMAX allocated dimension of array SYS. +* SIDE side of the hexagons. +* MAT mixture index assigned to each element. +* SIGT total minus self-scattering macroscopic cross sections. +* SIGT(:,NAN) generally contains the total cross section only. +* KN element-ordered unknown list. +* QFR element-ordered boundary conditions. +* 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. +* MU indices used with compressed diagonal storage mode matrix SYS. +* LC order of the unit matrices. +* R Cartesian mass matrix. +* V nodal coupling matrix. +* H Piolat (hexagonal) coupling matrix. +* +*Parameters: output +* SYS system matrix. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ITY,IELEM,ICOL,NBLOS,L4,NBMIX,IIMAX,MAT(3,NBLOS), + 1 IPERT(NBLOS),KN(NBLOS,4+6*IELEM*(IELEM+1)),NLF,NVD,NAN,MU(L4),LC + REAL SIDE,SIGT(NBMIX,NAN),QFR(NBLOS,6),R(LC,LC),V(LC,LC-1), + 1 H(LC,LC-1),SYS(IIMAX) +*---- +* LOCAL VARIABLES +*---- + PARAMETER(MAXIEL=3) + DOUBLE PRECISION CTRAN(MAXIEL*(MAXIEL+1),MAXIEL*(MAXIEL+1)),VAR1 +* + TTTT=REAL(0.5D0*SQRT(3.D00)*SIDE*SIDE) + IF(IELEM.GT.MAXIEL) CALL XABORT('PNDH2E: MAXIEL OVERFLOW.') + NZMAR=65 + 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 + ENDIF + MUMAX=MU(L4) + NELEM=IELEM*(IELEM+1) + COEF=REAL(2.0D0*SIDE*SIDE/SQRT(3.D00)) +*---- +* COMPUTE THE TRANVERSE COUPLING PIOLAT UNIT MATRIX +*---- + CTRAN(:MAXIEL*(MAXIEL+1),:MAXIEL*(MAXIEL+1))=0.0D0 + CNORM=REAL(SIDE*SIDE/SQRT(3.D00)) + I=0 + DO 22 JS=1,IELEM + DO 21 JT=1,IELEM+1 + J=0 + I=I+1 + SSS=1.0 + DO 20 IT=1,IELEM + DO 10 IS=1,IELEM+1 + J=J+1 + CTRAN(I,J)=SSS*CNORM*H(IS,JS)*H(JT,IT) + 10 CONTINUE + SSS=-SSS + 20 CONTINUE + 21 CONTINUE + 22 CONTINUE +*---- +* ASSEMBLY OF THE MAIN COEFFICIENT MATRIX AT ORDER IL. +*---- + DO 100 IL=0,NLF-1 + ZMARS=0.0 + IF(MOD(IL,2).EQ.1) ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) + NUM=0 + KEY=0 + DO 90 KEL=1,NBLOS + IF(IPERT(KEL).EQ.0) GO TO 90 + IBM=MAT(1,IPERT(KEL)) + IF(IBM.EQ.0) GO TO 90 + NUM=NUM+1 + GARS=SIGT(IBM,MIN(IL+1,NAN)) + IF(MOD(IL,2).EQ.0) THEN +* EVEN PARITY EQUATION. + DO 35 K2=0,IELEM-1 + DO 30 K1=0,IELEM-1 + JND1=KN(NUM,1)+K2*IELEM+K1 + JND2=KN(NUM,2)+K2*IELEM+K1 + JND3=KN(NUM,3)+K2*IELEM+K1 + KEY=(IL/2)*MUMAX+MU(JND1) + SYS(KEY)=SYS(KEY)+FACT*TTTT*GARS + KEY=(IL/2)*MUMAX+MU(JND2) + SYS(KEY)=SYS(KEY)+FACT*TTTT*GARS + KEY=(IL/2)*MUMAX+MU(JND3) + SYS(KEY)=SYS(KEY)+FACT*TTTT*GARS + 30 CONTINUE + 35 CONTINUE + ELSE +* ODD PARITY EQUATION. + DO 52 K4=0,1 + DO 51 K3=0,IELEM-1 + DO 50 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=ABS(KNW1) + INX1=ABS(KNX1) + INY1=ABS(KNY1) + DO 40 K1=1,IELEM+1 + KNW2=KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K1) + KNX2=KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K1) + KNY2=KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K1) + INW2=ABS(KNW2) + INX2=ABS(KNX2) + INY2=ABS(KNY2) + IF((KNW2.NE.0).AND.(KNW1.NE.0).AND.(INW1.GE.INW2)) THEN + KEY=(IL/2)*MUMAX+MU(INW1)-INW1+INW2 + SG=REAL(SIGN(1,KNW1)*SIGN(1,KNW2)) + SYS(KEY)=SYS(KEY)-SG*FACT*COEF*GARS*R(K2,K1) + ENDIF + IF((KNX2.NE.0).AND.(KNX1.NE.0).AND.(INX1.GE.INX2)) THEN + KEY=(IL/2)*MUMAX+MU(INX1)-INX1+INX2 + SG=REAL(SIGN(1,KNX1)*SIGN(1,KNX2)) + SYS(KEY)=SYS(KEY)-SG*FACT*COEF*GARS*R(K2,K1) + ENDIF + IF((KNY2.NE.0).AND.(KNY1.NE.0).AND.(INY1.GE.INY2)) THEN + KEY=(IL/2)*MUMAX+MU(INY1)-INY1+INY2 + SG=REAL(SIGN(1,KNY1)*SIGN(1,KNY2)) + SYS(KEY)=SYS(KEY)-SG*FACT*COEF*GARS*R(K2,K1) + ENDIF + 40 CONTINUE + IF(ITY.EQ.0) THEN + IF(KNW1.NE.0) THEN + KEY=(IL/2)*MUMAX+MU(INW1) + IF((K2.EQ.1).AND.(K4.EQ.0)) THEN + SYS(KEY)=SYS(KEY)-0.5*FACT*QFR(NUM,1)*ZMARS + ELSE IF((K2.EQ.IELEM+1).AND.(K4.EQ.1)) THEN + SYS(KEY)=SYS(KEY)-0.5*FACT*QFR(NUM,2)*ZMARS + ENDIF + ENDIF + IF(KNX1.NE.0) THEN + KEY=(IL/2)*MUMAX+MU(INX1) + IF((K2.EQ.1).AND.(K4.EQ.0)) THEN + SYS(KEY)=SYS(KEY)-0.5*FACT*QFR(NUM,3)*ZMARS + ELSE IF((K2.EQ.IELEM+1).AND.(K4.EQ.1)) THEN + SYS(KEY)=SYS(KEY)-0.5*FACT*QFR(NUM,4)*ZMARS + ENDIF + ENDIF + IF(KNY1.NE.0) THEN + KEY=(IL/2)*MUMAX+MU(INY1) + IF((K2.EQ.1).AND.(K4.EQ.0)) THEN + SYS(KEY)=SYS(KEY)-0.5*FACT*QFR(NUM,5)*ZMARS + ELSE IF((K2.EQ.IELEM+1).AND.(K4.EQ.1)) THEN + SYS(KEY)=SYS(KEY)-0.5*FACT*QFR(NUM,6)*ZMARS + ENDIF + ENDIF + ENDIF + 50 CONTINUE + 51 CONTINUE + 52 CONTINUE +* + ITRS=0 + DO I=1,NBLOS + IF(KN(I,1).EQ.KN(NUM,4)) THEN + ITRS=I + GO TO 60 + ENDIF + ENDDO + CALL XABORT('PNDH2E: ITRS FAILURE.') + 60 DO 75 I=1,NELEM + KNW1=KN(ITRS,4+I) + KNX1=KN(NUM,4+2*NELEM+I) + KNY1=KN(NUM,4+4*NELEM+I) + INW1=ABS(KNW1) + INX1=ABS(KNX1) + INY1=ABS(KNY1) + DO 70 J=1,NELEM + KNW2=KN(NUM,4+NELEM+J) + KNX2=KN(NUM,4+3*NELEM+J) + KNY2=KN(NUM,4+5*NELEM+J) + INW2=ABS(KNW2) + INX2=ABS(KNX2) + INY2=ABS(KNY2) + VAR1=FACT*GARS*CTRAN(I,J) + IF((KNY2.NE.0).AND.(KNW1.NE.0).AND.(INW1.LT.INY2)) THEN + KEY=(IL/2)*MUMAX+MU(INY2)-INY2+INW1 + SG=REAL(SIGN(1,KNW1)*SIGN(1,KNY2)) + SYS(KEY)=SYS(KEY)-SG*REAL(VAR1) ! y w + ELSE IF((KNY2.NE.0).AND.(KNW1.NE.0).AND.(INW1.GT.INY2)) THEN + KEY=(IL/2)*MUMAX+MU(INW1)-INW1+INY2 + SG=REAL(SIGN(1,KNW1)*SIGN(1,KNY2)) + SYS(KEY)=SYS(KEY)-SG*REAL(VAR1) ! w y + ENDIF + IF((KNW2.NE.0).AND.(KNX1.NE.0).AND.(INW2.LT.INX1)) THEN + KEY=(IL/2)*MUMAX+MU(INX1)-INX1+INW2 + SG=REAL(SIGN(1,KNX1)*SIGN(1,KNW2)) + SYS(KEY)=SYS(KEY)-SG*REAL(VAR1) ! x w + ELSE IF((KNW2.NE.0).AND.(KNX1.NE.0).AND.(INW2.GT.INX1)) THEN + KEY=(IL/2)*MUMAX+MU(INW2)-INW2+INX1 + SG=REAL(SIGN(1,KNX1)*SIGN(1,KNW2)) + SYS(KEY)=SYS(KEY)-SG*REAL(VAR1) ! w x + ENDIF + IF((KNX2.NE.0).AND.(KNY1.NE.0).AND.(INX2.LT.INY1)) THEN + KEY=(IL/2)*MUMAX+MU(INY1)-INY1+INX2 + SG=REAL(SIGN(1,KNY1)*SIGN(1,KNX2)) + SYS(KEY)=SYS(KEY)-SG*REAL(VAR1) ! y x + ELSE IF((KNX2.NE.0).AND.(KNY1.NE.0).AND.(INX2.GT.INY1)) THEN + KEY=(IL/2)*MUMAX+MU(INX2)-INX2+INY1 + SG=REAL(SIGN(1,KNY1)*SIGN(1,KNX2)) + SYS(KEY)=SYS(KEY)-SG*REAL(VAR1) ! x y + ENDIF + 70 CONTINUE + 75 CONTINUE +* + IF(ITY.EQ.0) THEN + DO 83 K4=0,1 + DO 82 K3=0,IELEM-1 + DO 81 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=ABS(KNW1) + INX1=ABS(KNX1) + INY1=ABS(KNY1) + DO 80 K1=0,IELEM-1 + IF(V(K2,K1+1).EQ.0.0) GO TO 80 + IF(K4.EQ.0) THEN + SSS=(-1.0)**K1 + JND1=KN(NUM,1)+K3*IELEM+K1 + JND2=KN(NUM,2)+K3*IELEM+K1 + JND3=KN(NUM,3)+K3*IELEM+K1 + ELSE + SSS=1.0 + JND1=KN(NUM,2)+K1*IELEM+K3 + JND2=KN(NUM,3)+K1*IELEM+K3 + JND3=KN(NUM,4)+K1*IELEM+K3 + ENDIF + IF(KNW1.NE.0) THEN + IF(JND1.GT.INW1) KEY=(IL/2)*MUMAX+MU(JND1)-JND1+INW1 + IF(JND1.LT.INW1) KEY=(IL/2)*MUMAX+MU(INW1)-INW1+JND1 + SG=REAL(SIGN(1,KNW1)) + SYS(KEY)=SYS(KEY)+SG*SSS*REAL(IL)*SIDE*V(K2,K1+1) + ENDIF + IF(KNX1.NE.0) THEN + IF(JND2.GT.INX1) KEY=(IL/2)*MUMAX+MU(JND2)-JND2+INX1 + IF(JND2.LT.INX1) KEY=(IL/2)*MUMAX+MU(INX1)-INX1+JND2 + SG=REAL(SIGN(1,KNX1)) + SYS(KEY)=SYS(KEY)+SG*SSS*REAL(IL)*SIDE*V(K2,K1+1) + ENDIF + IF(KNY1.NE.0) THEN + IF(JND3.GT.INY1) KEY=(IL/2)*MUMAX+MU(JND3)-JND3+INY1 + IF(JND3.LT.INY1) KEY=(IL/2)*MUMAX+MU(INY1)-INY1+JND3 + SG=REAL(SIGN(1,KNY1)) + SYS(KEY)=SYS(KEY)+SG*SSS*REAL(IL)*SIDE*V(K2,K1+1) + ENDIF + 80 CONTINUE + 81 CONTINUE + 82 CONTINUE + 83 CONTINUE + ENDIF + ENDIF + 90 CONTINUE + 100 CONTINUE + RETURN + END -- cgit v1.2.3