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/PN3DXX.f | 455 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 455 insertions(+) create mode 100755 Trivac/src/PN3DXX.f (limited to 'Trivac/src/PN3DXX.f') diff --git a/Trivac/src/PN3DXX.f b/Trivac/src/PN3DXX.f new file mode 100755 index 0000000..5f51913 --- /dev/null +++ b/Trivac/src/PN3DXX.f @@ -0,0 +1,455 @@ +*DECK PN3DXX + SUBROUTINE PN3DXX(NBMIX,IELEM,ICOL,NEL,NLF,NVD,NAN,LL4F,LL4X, + 1 SIGT,SIGTI,MAT,VOL,XX,YY,ZZ,KN,QFR,MUX,IPBBX,LC,R,V,BBX,TTF, + 2 AX,C11X) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of system matrices for a Thomas-Raviart (dual) finite element +* method in 3-D simplified PN approximation. Note: system matrices +* should be initialized by the calling program. +* +*Copyright: +* Copyright (C) 2005 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 +* NBMIX number of mixtures. +* IELEM degree of the Lagrangian finite elements: =1 (linear); +* =2 (parabolic); =3 (cubic). +* ICOL type of quadrature: =1 (analytical integration); +* =2 (Gauss-Lobatto); =3 (Gauss-Legendre). +* NEL total number of finite elements. +* 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. +* LL4F number of flux components. +* LL4X number of X-directed currents. +* LL4Y number of Y-directed currents. +* LL4Z number of Z-directed currents. +* SIGT total minus self-scattering macroscopic cross sections. +* SIGT(:,NAN) generally contains the total cross section only. +* SIGTI inverse macroscopic cross sections ordered by mixture. +* SIGTI(:,NAN) generally contains the inverse total cross +* section only. +* MAT mixture index assigned to each element. +* VOL volume of each element. +* XX X-directed mesh spacings. +* YY Y-directed mesh spacings. +* ZZ Z-directed mesh spacings. +* KN element-ordered unknown list. +* QFR element-ordered boundary conditions. +* MUX X-directed compressed storage mode indices. +* MUY Y-directed compressed storage mode indices. +* MUZ Z-directed compressed storage mode indices. +* IPBBX X-directed perdue storage indices. +* IPBBY Y-directed perdue storage indices. +* IPBBZ Z-directed perdue storage indices. +* LC order of the unit matrices. +* R unit matrix. +* V unit matrix. +* BBX X-directed flux-current matrices. +* BBY Y-directed flux-current matrices. +* BBZ Z-directed flux-current matrices. +* +*Parameters: output +* TTF flux-flux matrices. +* AX X-directed main current-current matrices. Dimensionned to +* MUX(LL4X)*NLF/2. +* AY Y-directed main current-current matrices. Dimensionned to +* MUY(LL4Y)*NLF/2. +* AZ Z-directed main current-current matrices. Dimensionned to +* MUZ(LL4Z)*NLF/2. +* C11X X-directed main current-current matrices to be factorized. +* Dimensionned to MUX(LL4X)*NLF/2. +* C11Y Y-directed main current-current matrices to be factorized. +* Dimensionned to MUY(LL4Y)*NLF/2. +* C11Z Z-directed main current-current matrices to be factorized. +* Dimensionned to MUZ(LL4Z)*NLF/2. +* +*Reference(s): +* J.J. Lautard, D. Schneider, A.M. Baudron, "Mixed Dual Methods for +* Neutronic Reactor Core Calculations in the CRONOS System," Proc. +* Int. Conf. on Mathematics and Computation, Reactor Physics and +* Environmental Analysis in Nuclear Applications, Madrid, Spain, +* September 27-30, 1999. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,IELEM,ICOL,NEL,NLF,NVD,NAN,LL4F,LL4X,MAT(NEL), + 1 KN(NEL*(1+6*IELEM**2)),MUX(LL4X),IPBBX(2*IELEM,LL4X),LC + REAL SIGT(NBMIX,NAN),SIGTI(NBMIX,NAN),VOL(NEL),XX(NEL),YY(NEL), + 1 ZZ(NEL),QFR(6*NEL),R(LC,LC),V(LC,LC-1),BBX(2*IELEM,LL4X), + 2 TTF(LL4F*NLF/2),AX(*),C11X(*) +*---- +* LOCAL VARIABLES +*---- + REAL QQ(5,5) +*---- +* X-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + 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 25 I0=1,IELEM + DO 20 J0=1,IELEM + FFF=0.0 + 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.0 + QQ(I0,J0)=FFF + 20 CONTINUE + 25 CONTINUE + MUMAX=MUX(LL4X) + DO 170 IL=0,NLF-1 + IF(MOD(IL,2).EQ.1) ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE X-ORIENTED COEFFICIENT MATRICES AT ORDER IL. +*---- + NUM1=0 + NUM2=0 + DO 120 IE=1,NEL + IBM=MAT(IE) + IF(IBM.EQ.0) GO TO 120 + VOL0=VOL(IE) + IF(VOL0.EQ.0.0) GO TO 110 + DX=XX(IE) + DY=YY(IE) + DZ=ZZ(IE) + GARS=SIGT(IBM,MIN(IL+1,NAN)) + IF(MOD(IL,2).EQ.0) THEN +* EVEN PARITY EQUATION. + DO 32 K3=0,IELEM-1 + DO 31 K2=0,IELEM-1 + DO 30 K1=0,IELEM-1 + KEY=(IL/2)*LL4F+KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1 + TTF(KEY)=TTF(KEY)+FACT*VOL0*GARS + 30 CONTINUE + 31 CONTINUE + 32 CONTINUE + ELSE + GARSI=SIGTI(IBM,MIN(IL+1,NAN)) + DO 105 K3=0,IELEM-1 + DO 100 K2=0,IELEM-1 +* PARTIAL INVERSION OF THE ODD PARITY EQUATION. MODIFICATION OF +* THE EVEN PARITY EQUATION. + DO 40 K1=0,IELEM-1 + JND1=KN(NUM1+1)+(K3*IELEM+K2)*IELEM+K1 + KEY=((IL-1)/2)*LL4F+JND1 + TTF(KEY)=TTF(KEY)+(REAL(IL)**2)*VOL0*QQ(K1+1,K1+1)*GARSI/(FACT* + 1 DX*DX) + IF(IL.LE.NLF-3) THEN + KEY=((IL+2)/2)*LL4F+JND1 + TTF(KEY)=TTF(KEY)+(REAL(IL+1)**2)*VOL0*QQ(K1+1,K1+1)*GARSI/ + 1 (FACT*DX*DX) + ENDIF + KEY=((IL-1)/2)*LL4F+JND1 + TTF(KEY)=TTF(KEY)+(REAL(IL)**2)*VOL0*QQ(K2+1,K2+1)*GARSI/ + 1 (FACT*DY*DY) + IF(IL.LE.NLF-3) THEN + KEY=((IL+2)/2)*LL4F+JND1 + TTF(KEY)=TTF(KEY)+(REAL(IL+1)**2)*VOL0*QQ(K2+1,K2+1)*GARSI/ + 1 (FACT*DY*DY) + ENDIF + KEY=((IL-1)/2)*LL4F+JND1 + TTF(KEY)=TTF(KEY)+(REAL(IL)**2)*VOL0*QQ(K3+1,K3+1)*GARSI/ + 1 (FACT*DZ*DZ) + IF(IL.LE.NLF-3) THEN + KEY=((IL+2)/2)*LL4F+JND1 + TTF(KEY)=TTF(KEY)+(REAL(IL+1)**2)*VOL0*QQ(K3+1,K3+1)*GARSI/ + 1 (FACT*DZ*DZ) + ENDIF + 40 CONTINUE +* +* ODD PARITY EQUATION. + DO 55 IC=1,2 + IF(IC.EQ.1) IIC=1 + IF(IC.EQ.2) IIC=IELEM+1 + KN1=KN(NUM1+2+(IC-1)*IELEM**2+K3*IELEM+K2) + IND1=ABS(KN1)-LL4F + S1=REAL(SIGN(1,KN1)) + DO 50 JC=1,2 + IF(JC.EQ.1) JJC=1 + IF(JC.EQ.2) JJC=IELEM+1 + KN2=KN(NUM1+2+(JC-1)*IELEM**2+K3*IELEM+K2) + IND2=ABS(KN2)-LL4F + IF((KN1.NE.0).AND.(KN2.NE.0).AND.(IND1.GE.IND2)) THEN + S2=REAL(SIGN(1,KN2)) + KEY=((IL-1)/2)*MUMAX+MUX(IND1)-IND1+IND2 + AX(KEY)=AX(KEY)-S1*S2*FACT*R(IIC,JJC)*VOL0*GARS + ENDIF + 50 CONTINUE + 55 CONTINUE +* + KN1=KN(NUM1+2+K3*IELEM+K2) + KN2=KN(NUM1+2+IELEM**2+K3*IELEM+K2) + IND1=ABS(KN1)-LL4F + IND2=ABS(KN2)-LL4F + IF((QFR(NUM2+1).NE.0.0).AND.(KN1.NE.0)) THEN + KEY=((IL-1)/2)*MUMAX+MUX(IND1) + AX(KEY)=AX(KEY)-0.5*FACT*QFR(NUM2+1)*ZMARS + ENDIF + IF((QFR(NUM2+2).NE.0.0).AND.(KN2.NE.0)) THEN + KEY=((IL-1)/2)*MUMAX+MUX(IND2) + AX(KEY)=AX(KEY)-0.5*FACT*QFR(NUM2+2)*ZMARS + ENDIF + 100 CONTINUE + 105 CONTINUE + ENDIF + 110 NUM1=NUM1+1+6*IELEM**2 + NUM2=NUM2+6 + 120 CONTINUE +* + IF(MOD(IL,2).EQ.1) THEN + DO 130 I0=1,MUMAX + C11X(((IL-1)/2)*MUMAX+I0)=-AX(((IL-1)/2)*MUMAX+I0) + 130 CONTINUE + MUIM1=0 + DO 160 I=1,LL4X + MUI=MUX(I) + DO 150 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 145 I0=1,2*IELEM + II=IPBBX(I0,I) + IF(II.EQ.0) GO TO 150 + DO 140 J0=1,2*IELEM + JJ=IPBBX(J0,J) + IF(II.EQ.JJ) C11X(KEY)=C11X(KEY)+REAL(IL**2)*BBX(I0,I)* + 1 BBX(J0,J)/TTF(((IL-1)/2)*LL4F+II) + 140 CONTINUE + 145 CONTINUE + 150 CONTINUE + MUIM1=MUI + 160 CONTINUE + ENDIF + 170 CONTINUE + RETURN + END +* + SUBROUTINE PN3DXY(NBMIX,IELEM,ICOL,NEL,NLF,NVD,NAN,LL4F,LL4X,LL4Y, + 1 SIGT,MAT,VOL,YY,KN,QFR,MUY,IPBBY,LC,R,BBY,TTF,AY,C11Y) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,IELEM,ICOL,NEL,NLF,NVD,NAN,LL4F,LL4X,LL4Y,MAT(NEL), + 1 KN(NEL*(1+6*IELEM**2)),MUY(LL4Y),IPBBY(2*IELEM,LL4Y),LC + REAL SIGT(NBMIX,NAN),VOL(NEL),YY(NEL),QFR(6*NEL),R(LC,LC), + 1 BBY(2*IELEM,LL4Y),TTF(LL4F*NLF/2),AY(*),C11Y(*) +*---- +* Y-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + 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=MUY(LL4Y) + DO 320 IL=1,NLF-1,2 + ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE Y-ORIENTED COEFFICIENT MATRICES AT ODD ORDER IL. +*---- + NUM1=0 + NUM2=0 + DO 270 IE=1,NEL + IBM=MAT(IE) + IF(IBM.EQ.0) GO TO 270 + VOL0=VOL(IE) + IF(VOL0.EQ.0.0) GO TO 260 + DY=YY(IE) + GARS=SIGT(IBM,MIN(IL+1,NAN)) +* + DO 255 K3=0,IELEM-1 + DO 250 K1=0,IELEM-1 + DO 205 IC=3,4 + IF(IC.EQ.3) IIC=1 + IF(IC.EQ.4) IIC=IELEM+1 + KN1=KN(NUM1+2+(IC-1)*IELEM**2+K3*IELEM+K1) + IND1=ABS(KN1)-LL4F-LL4X + S1=REAL(SIGN(1,KN1)) + DO 200 JC=3,4 + IF(JC.EQ.3) JJC=1 + IF(JC.EQ.4) JJC=IELEM+1 + KN2=KN(NUM1+2+(JC-1)*IELEM**2+K3*IELEM+K1) + IND2=ABS(KN2)-LL4F-LL4X + IF((KN1.NE.0).AND.(KN2.NE.0).AND.(IND1.GE.IND2)) THEN + S2=REAL(SIGN(1,KN2)) + KEY=((IL-1)/2)*MUMAX+MUY(IND1)-IND1+IND2 + AY(KEY)=AY(KEY)-S1*S2*FACT*R(IIC,JJC)*VOL0*GARS + ENDIF + 200 CONTINUE + 205 CONTINUE +* + KN1=KN(NUM1+2+2*IELEM**2+K3*IELEM+K1) + KN2=KN(NUM1+2+3*IELEM**2+K3*IELEM+K1) + IND1=ABS(KN1)-LL4F-LL4X + IND2=ABS(KN2)-LL4F-LL4X + IF((QFR(NUM2+3).NE.0.0).AND.(KN1.NE.0)) THEN + KEY=((IL-1)/2)*MUMAX+MUY(IND1) + AY(KEY)=AY(KEY)-0.5*FACT*QFR(NUM2+3)*ZMARS + ENDIF + IF((QFR(NUM2+4).NE.0.0).AND.(KN2.NE.0)) THEN + KEY=((IL-1)/2)*MUMAX+MUY(IND2) + AY(KEY)=AY(KEY)-0.5*FACT*QFR(NUM2+4)*ZMARS + ENDIF + 250 CONTINUE + 255 CONTINUE + 260 NUM1=NUM1+1+6*IELEM**2 + NUM2=NUM2+6 + 270 CONTINUE +* + DO 280 I0=1,MUMAX + C11Y(((IL-1)/2)*MUMAX+I0)=-AY(((IL-1)/2)*MUMAX+I0) + 280 CONTINUE + MUIM1=0 + DO 310 I=1,LL4Y + MUI=MUY(I) + DO 300 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 295 I0=1,2*IELEM + II=IPBBY(I0,I) + IF(II.EQ.0) GO TO 300 + DO 290 J0=1,2*IELEM + JJ=IPBBY(J0,J) + IF(II.EQ.JJ) C11Y(KEY)=C11Y(KEY)+REAL(IL**2)*BBY(I0,I)*BBY(J0,J)/ + 1 TTF(((IL-1)/2)*LL4F+II) + 290 CONTINUE + 295 CONTINUE + 300 CONTINUE + MUIM1=MUI + 310 CONTINUE + 320 CONTINUE + RETURN + END +* + SUBROUTINE PN3DXZ(NBMIX,IELEM,ICOL,NEL,NLF,NVD,NAN,LL4F,LL4X,LL4Y, + 1 LL4Z,SIGT,MAT,VOL,ZZ,KN,QFR,MUZ,IPBBZ,LC,R,BBZ,TTF,AZ,C11Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NBMIX,IELEM,ICOL,NEL,NLF,NVD,NAN,LL4F,LL4X,LL4Y,LL4Z, + 1 MAT(NEL),KN(NEL*(1+6*IELEM**2)),MUZ(LL4Z),IPBBZ(2*IELEM,LL4Z),LC + REAL SIGT(NBMIX,NAN),VOL(NEL),ZZ(NEL),QFR(6*NEL),R(LC,LC), + 1 BBZ(2*IELEM,LL4Z),TTF(LL4F*NLF/2),AZ(*),C11Z(*) +*---- +* Z-ORIENTED COUPLINGS +*---- + ZMARS=0.0 + 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=MUZ(LL4Z) + DO 470 IL=1,NLF-1,2 + ZMARS=PNMAR2(NZMAR,IL,IL) + FACT=REAL(2*IL+1) +*---- +* ASSEMBLY OF THE Z-ORIENTED COEFFICIENT MATRICES AT ORDER IL. +*---- + NUM1=0 + NUM2=0 + DO 420 IE=1,NEL + IBM=MAT(IE) + IF(IBM.EQ.0) GO TO 420 + VOL0=VOL(IE) + IF(VOL0.EQ.0.0) GO TO 410 + DZ=ZZ(IE) + GARS=SIGT(IBM,MIN(IL+1,NAN)) +* + DO 405 K2=0,IELEM-1 + DO 400 K1=0,IELEM-1 + DO 355 IC=5,6 + IF(IC.EQ.5) IIC=1 + IF(IC.EQ.6) IIC=IELEM+1 + KN1=KN(NUM1+2+(IC-1)*IELEM**2+K2*IELEM+K1) + IND1=ABS(KN1)-LL4F-LL4X-LL4Y + S1=REAL(SIGN(1,KN1)) + DO 350 JC=5,6 + IF(JC.EQ.5) JJC=1 + IF(JC.EQ.6) JJC=IELEM+1 + KN2=KN(NUM1+2+(JC-1)*IELEM**2+K2*IELEM+K1) + IND2=ABS(KN2)-LL4F-LL4X-LL4Y + IF((KN1.NE.0).AND.(KN2.NE.0).AND.(IND1.GE.IND2)) THEN + S2=REAL(SIGN(1,KN2)) + KEY=((IL-1)/2)*MUMAX+MUZ(IND1)-IND1+IND2 + AZ(KEY)=AZ(KEY)-S1*S2*FACT*R(IIC,JJC)*VOL0*GARS + ENDIF + 350 CONTINUE + 355 CONTINUE +* + KN1=KN(NUM1+2+4*IELEM**2+K2*IELEM+K1) + KN2=KN(NUM1+2+5*IELEM**2+K2*IELEM+K1) + IND1=ABS(KN1)-LL4F-LL4X-LL4Y + IND2=ABS(KN2)-LL4F-LL4X-LL4Y + IF((QFR(NUM2+5).NE.0.0).AND.(KN1.NE.0)) THEN + KEY=((IL-1)/2)*MUMAX+MUZ(IND1) + AZ(KEY)=AZ(KEY)-0.5*FACT*QFR(NUM2+5)*ZMARS + ENDIF + IF((QFR(NUM2+6).NE.0.0).AND.(KN2.NE.0)) THEN + KEY=((IL-1)/2)*MUMAX+MUZ(IND2) + AZ(KEY)=AZ(KEY)-0.5*FACT*QFR(NUM2+6)*ZMARS + ENDIF + 400 CONTINUE + 405 CONTINUE + 410 NUM1=NUM1+1+6*IELEM**2 + NUM2=NUM2+6 + 420 CONTINUE +* + DO 430 I0=1,MUMAX + C11Z(((IL-1)/2)*MUMAX+I0)=-AZ(((IL-1)/2)*MUMAX+I0) + 430 CONTINUE + MUIM1=0 + DO 460 I=1,LL4Z + MUI=MUZ(I) + DO 450 J=I-(MUI-MUIM1)+1,I + KEY=((IL-1)/2)*MUMAX+(MUI-I+J) + DO 445 I0=1,2*IELEM + II=IPBBZ(I0,I) + IF(II.EQ.0) GO TO 450 + DO 440 J0=1,2*IELEM + JJ=IPBBZ(J0,J) + IF(II.EQ.JJ) C11Z(KEY)=C11Z(KEY)+REAL(IL**2)*BBZ(I0,I)*BBZ(J0,J)/ + 1 TTF(((IL-1)/2)*LL4F+II) + 440 CONTINUE + 445 CONTINUE + 450 CONTINUE + MUIM1=MUI + 460 CONTINUE + 470 CONTINUE + RETURN + END -- cgit v1.2.3