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/BIVA02.f | 210 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 210 insertions(+) create mode 100755 Trivac/src/BIVA02.f (limited to 'Trivac/src/BIVA02.f') diff --git a/Trivac/src/BIVA02.f b/Trivac/src/BIVA02.f new file mode 100755 index 0000000..da999ea --- /dev/null +++ b/Trivac/src/BIVA02.f @@ -0,0 +1,210 @@ +*DECK BIVA02 + SUBROUTINE BIVA02(ITY,SGD,CYLIND,IELEM,ICOL,NREG,LL4,NBMIX,IIMAX, + 1 XX,YY,DD,MAT,KN,QFR,VOL,MU,LC,R,V,SYS) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Assembly of a within-group (leakage and removal) or out-of-group +* system matrix in mixed-dual finite element diffusion approximation +* (Cartesian geometry). +* +*Copyright: +* Copyright (C) 2002 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. +* SGD nuclear properties. SGD(:,1) and SGD(:,2) are diffusion +* coefficients. SGD(:,3) are removal macroscopic cross sections. +* CYLIND cylinderization flag (=.true. for cylindrical geometry). +* 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). +* NREG number of elements in BIVAC. +* LL4 number of unknowns per group in BIVAC. +* NBMIX number of macro-mixtures. +* IIMAX allocated dimension of array SYS. +* XX X-directed mesh spacings. +* YY Y-directed mesh spacings. +* DD value used used with a cylindrical geometry. +* MAT mixture index per region. +* KN element-ordered unknown list. +* QFR element-ordered boundary conditions. +* VOL volume of regions. +* MU indices used with compressed diagonal storage mode matrix SYS. +* LC number of polynomials in a complete 1-D basis. +* R cartesian mass matrix. +* V nodal coupling matrix. +* +*Parameters: output +* SYS system matrix. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ITY,IELEM,ICOL,NREG,LL4,NBMIX,IIMAX,MAT(NREG),KN(5*NREG), + 1 MU(LL4),LC + REAL SGD(NBMIX,3),XX(NREG),YY(NREG),DD(NREG),QFR(4*NREG), + 1 VOL(NREG),R(LC,LC),V(LC,LC-1),SYS(IIMAX) + LOGICAL CYLIND +*---- +* LOCAL VARIABLES +*---- + REAL QQ(5,5) +* + IF((CYLIND).AND.((IELEM.GT.1).OR.(ICOL.NE.2))) + 1 CALL XABORT('BIVA02: TYPE OF DISCRETIZATION NOT IMPLEMENTED.') +*---- +* ASSEMBLY OF A SYSTEM MATRIX. +*---- + IF(ITY.EQ.0) THEN +* LEAKAGE-REMOVAL SYSTEM MATRIX ASSEMBLY. + DO 12 I0=1,IELEM + DO 11 J0=1,IELEM + QQ(I0,J0)=0.0 + DO 10 K0=2,IELEM + QQ(I0,J0)=QQ(I0,J0)+V(K0,I0)*V(K0,J0)/R(K0,K0) + 10 CONTINUE + 11 CONTINUE + 12 CONTINUE + NUM1=0 + NUM2=0 + DO 80 K=1,NREG + L=MAT(K) + IF(L.EQ.0) GO TO 80 + VOL0=VOL(K) + IF(VOL0.EQ.0.0) GO TO 70 + DX=XX(K) + DY=YY(K) + IF(CYLIND) THEN + DIN=1.0-0.5*DX/DD(K) + DOT=1.0+0.5*DX/DD(K) + ELSE + DIN=1.0 + DOT=1.0 + ENDIF +* + DO 60 I0=1,IELEM + INX1=ABS(KN(NUM1+2))+I0-1 + INX2=ABS(KN(NUM1+3))+I0-1 + INY1=ABS(KN(NUM1+4))+I0-1 + INY2=ABS(KN(NUM1+5))+I0-1 + DO 50 J0=1,IELEM + JND1=KN(NUM1+1)+(I0-1)*IELEM+J0-1 + KEY=MU(JND1) + SYS(KEY)=SYS(KEY)+VOL0*SGD(L,3) + DO 20 K0=1,J0 + IF(QQ(J0,K0).EQ.0.0) GO TO 20 + KND1=KN(NUM1+1)+(I0-1)*IELEM+K0-1 + KEY=MU(JND1)-JND1+KND1 + SYS(KEY)=SYS(KEY)+VOL0*QQ(J0,K0)*SGD(L,1)/(DX*DX) + 20 CONTINUE + IF(KN(NUM1+2).NE.0) THEN + IF(JND1.GT.INX1) KEY=MU(JND1)-JND1+INX1 + IF(JND1.LT.INX1) KEY=MU(INX1)-INX1+JND1 + SG=REAL(SIGN(1,KN(NUM1+2))) + SYS(KEY)=SYS(KEY)+SG*(VOL0/DX)*DIN*V(1,J0) + ENDIF + IF(KN(NUM1+3).NE.0) THEN + IF(INX2.GT.JND1) KEY=MU(INX2)-INX2+JND1 + IF(INX2.LT.JND1) KEY=MU(JND1)-JND1+INX2 + SG=REAL(SIGN(1,KN(NUM1+3))) + SYS(KEY)=SYS(KEY)+SG*(VOL0/DX)*DOT*V(IELEM+1,J0) + ENDIF + JND1=KN(NUM1+1)+(J0-1)*IELEM+I0-1 + DO 30 K0=1,J0 + IF(QQ(J0,K0).EQ.0.0) GO TO 30 + KND1=KN(NUM1+1)+(K0-1)*IELEM+I0-1 + KEY=MU(JND1)-JND1+KND1 + SYS(KEY)=SYS(KEY)+VOL0*QQ(J0,K0)*SGD(L,2)/(DY*DY) + 30 CONTINUE + IF(KN(NUM1+4).NE.0) THEN + IF(JND1.GT.INY1) KEY=MU(JND1)-JND1+INY1 + IF(JND1.LT.INY1) KEY=MU(INY1)-INY1+JND1 + SG=REAL(SIGN(1,KN(NUM1+4))) + SYS(KEY)=SYS(KEY)+SG*(VOL0/DY)*V(1,J0) + ENDIF + IF(KN(NUM1+5).NE.0) THEN + IF(INY2.GT.JND1) KEY=MU(INY2)-INY2+JND1 + IF(INY2.LT.JND1) KEY=MU(JND1)-JND1+INY2 + SG=REAL(SIGN(1,KN(NUM1+5))) + SYS(KEY)=SYS(KEY)+SG*(VOL0/DY)*V(IELEM+1,J0) + ENDIF + 50 CONTINUE + IF(KN(NUM1+2).NE.0) THEN + KEY=MU(INX1) + SYS(KEY)=SYS(KEY)-DIN*(VOL0*R(1,1)/SGD(L,1)+QFR(NUM2+1)) + ENDIF + IF(KN(NUM1+3).NE.0) THEN + KEY=MU(INX2) + SYS(KEY)=SYS(KEY)-DOT*(VOL0*R(IELEM+1,IELEM+1)/SGD(L,1) + 1 +QFR(NUM2+2)) + ENDIF + IF(KN(NUM1+4).NE.0) THEN + KEY=MU(INY1) + SYS(KEY)=SYS(KEY)-VOL0*R(1,1)/SGD(L,2)-QFR(NUM2+3) + ENDIF + IF(KN(NUM1+5).NE.0) THEN + KEY=MU(INY2) + SYS(KEY)=SYS(KEY)-VOL0*R(IELEM+1,IELEM+1)/SGD(L,2) + 1 -QFR(NUM2+4) + ENDIF + IF(ICOL.NE.2) THEN + IF((KN(NUM1+2).NE.0).AND.(KN(NUM1+3).NE.0)) THEN + IF(INX2.GT.INX1) KEY=MU(INX2)-INX2+INX1 + IF(INX2.LE.INX1) KEY=MU(INX1)-INX1+INX2 + SG=REAL(SIGN(1,KN(NUM1+2))*SIGN(1,KN(NUM1+3))) + IF(INX1.EQ.INX2) SG=2.0*SG + SYS(KEY)=SYS(KEY)-SG*VOL0*R(IELEM+1,1)/SGD(L,1) + ENDIF + IF((KN(NUM1+4).NE.0).AND.(KN(NUM1+5).NE.0)) THEN + IF(INY2.GT.INY1) KEY=MU(INY2)-INY2+INY1 + IF(INY2.LE.INY1) KEY=MU(INY1)-INY1+INY2 + SG=REAL(SIGN(1,KN(NUM1+4))*SIGN(1,KN(NUM1+5))) + IF(INY1.EQ.INY2) SG=2.0*SG + SYS(KEY)=SYS(KEY)-SG*VOL0*R(IELEM+1,1)/SGD(L,2) + ENDIF + ENDIF + 60 CONTINUE + 70 NUM1=NUM1+5 + NUM2=NUM2+4 + 80 CONTINUE + ELSE +* CROSS SECTION SYSTEM MATRIX ASSEMBLY. COMPONENTS WITH 1E-10 +* FACTORS ARE INTRODUCED TO MAKE THE MATRIX INVERTIBLE. + NUM1=0 + DO 110 K=1,NREG + L=MAT(K) + IF(L.EQ.0) GO TO 110 + VOL0=VOL(K) + IF(VOL0.EQ.0.0) GO TO 100 + DO 95 I0=1,IELEM + INX1=ABS(KN(NUM1+2))+I0-1 + INX2=ABS(KN(NUM1+3))+I0-1 + INY1=ABS(KN(NUM1+4))+I0-1 + INY2=ABS(KN(NUM1+5))+I0-1 + IF(KN(NUM1+2).NE.0) SYS(MU(INX1))=SYS(MU(INX1))+1.0E-30 + IF(KN(NUM1+3).NE.0) SYS(MU(INX2))=SYS(MU(INX2))+1.0E-30 + IF(KN(NUM1+4).NE.0) SYS(MU(INY1))=SYS(MU(INY1))+1.0E-30 + IF(KN(NUM1+5).NE.0) SYS(MU(INY2))=SYS(MU(INY2))+1.0E-30 + DO 90 J0=1,IELEM + JND1=KN(NUM1+1)+(I0-1)*IELEM+J0-1 + KEY=MU(JND1) + SYS(KEY)=SYS(KEY)+VOL0*SGD(L,1) + 90 CONTINUE + 95 CONTINUE + 100 NUM1=NUM1+5 + 110 CONTINUE + ENDIF + RETURN + END -- cgit v1.2.3