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/NSSFL1.f | 220 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 220 insertions(+) create mode 100755 Trivac/src/NSSFL1.f (limited to 'Trivac/src/NSSFL1.f') diff --git a/Trivac/src/NSSFL1.f b/Trivac/src/NSSFL1.f new file mode 100755 index 0000000..838d2a3 --- /dev/null +++ b/Trivac/src/NSSFL1.f @@ -0,0 +1,220 @@ +*DECK NSSFL1 + SUBROUTINE NSSFL1(IPFLX,NUN,NG,NEL,NMIX,NALB,ITRIAL,EPSOUT,MAXOUT, + 1 MAT,XX,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,FD,IPRINT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Flux calculation for the nodal expansion method in Cartesian 1D +* geometry. +* +*Copyright: +* Copyright (C) 2021 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 +* IPFLX nodal flux. +* NUN number of unknowns (=4*NEL+1). +* NG number of energy groups. +* NEL number of nodes in the nodal calculation. +* NMIX number of mixtures in the nodal calculation. +* NALB number of physical albedos. +* ITRIAL type of expansion functions in the nodal calculation +* (=1: polynomial; =2: hyperbolic). +* EPSOUT convergence epsilon for the power method. +* MAXOUT maximum number of iterations for the power method. +* MAT material mixtures. +* XX mesh spacings. +* IQFR boundary condition information. +* QFR albedo function information. +* DIFF diffusion coefficients +* SIGR removal cross sections. +* CHI fission spectra. +* SIGF nu times fission cross section. +* SCAT scattering cross section. +* BETA albedos. +* FD discontinuity factors. +* IPRINT edition flag. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPFLX + INTEGER NUN,NG,NEL,NMIX,NALB,ITRIAL(NMIX,NG),MAXOUT,IPRINT, + 1 MAT(NEL),IQFR(6,NEL) + REAL EPSOUT,XX(NEL),QFR(6,NEL),DIFF(NMIX,NG),SIGR(NMIX,NG), + 1 CHI(NMIX,NG),SIGF(NMIX,NG),SCAT(NMIX,NG,NG),BETA(NALB,NG,NG), + 2 FD(NMIX,2,NG,NG) +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) :: JPFLX + INTEGER :: DIM + REAL :: KEFF + REAL, ALLOCATABLE, DIMENSION(:) :: WORK + REAL, ALLOCATABLE, DIMENSION(:,:) :: EVECT,A,B,AI,A11,QFR2,FUNKN +* + ALB(X)=0.5*(1.0-X)/(1.0+X) +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(EVECT(NUN,NG)) + DIM=5*NEL + ALLOCATE(FUNKN(DIM,NG),A(DIM*NG,DIM*NG),B(DIM*NG,DIM*NG)) + ALLOCATE(WORK(NMIX),A11(DIM,DIM),QFR2(6,NEL)) +*---- +* INITIALIZATIONS +*---- + CALL LCMLEN(IPFLX,'FLUX',ILONG,ITYLCM) + IF(ILONG == 0) THEN + JPFLX=LCMLID(IPFLX,'FLUX',NG) + EVECT(:NUN,:NG)=1.0 + KEFF=1.0 + ELSE + JPFLX=LCMGID(IPFLX,'FLUX') + DO IG=1,NG + CALL LCMLEL(JPFLX,IG,ILONG,ITYLCM) + IF(ILONG /= NUN) CALL XABORT('NSSFL3: INVALID FLUX.') + CALL LCMGDL(JPFLX,IG,EVECT(1,IG)) + ENDDO + CALL LCMGET(IPFLX,'K-EFFECTIVE',KEFF) + ENDIF + FUNKN(:,:)=0.0 + DO IEL=1,NEL + IOF=(IEL-1)*5 + FUNKN(IOF+1,:)=EVECT(IEL,:) + ENDDO +*---- +* COMPUTE NODAL SOLUTION +*---- + DIM=5*NEL + A(:DIM*NG,:DIM*NG)=0.0 + B(:DIM*NG,:DIM*NG)=0.0 + QFR2(:6,:NEL)=0.0 + DO J=1,NG + IOF1=(J-1)*DIM + DO I=1,NG + DO IQW=1,2 + DO IEL=1,NEL + IALB=IQFR(IQW,IEL) + IF(IALB.GT.0) THEN + IF(IALB.GT.NALB) CALL XABORT('NSSFL1: BETA OVERFLOW.') + QFR2(IQW,IEL)=QFR(IQW,IEL)*ALB(BETA(IALB,I,J)) + ELSE IF(I == J) THEN + QFR2(IQW,IEL)=QFR(IQW,IEL) + ELSE + QFR2(IQW,IEL)=0.0 + ENDIF + ENDDO + ENDDO + DO IBM=1,NMIX + WORK(IBM)=CHI(IBM,I)*SIGF(IBM,J) + ENDDO + IOF2=(I-1)*DIM + IF(I == J) THEN + CALL NSS1TR(ITRIAL(1,J),NEL,NMIX,MAT,XX,IQFR,QFR2,DIFF(:,I), + 1 SIGR(:,I),FD(:,:,I,J),A11) + A(IOF1+1:IOF1+DIM,IOF1+1:IOF1+DIM)=A11(:,:) + ELSE + CALL NSS2TR(ITRIAL(1,J),NEL,NMIX,MAT,XX,IQFR,QFR2,DIFF(:,J), + 1 SIGR(:,J),SCAT(:,I,J),FD(:,:,I,J),A11) + A(IOF2+1:IOF2+DIM,IOF1+1:IOF1+DIM)=-A11(:,:) + ENDIF + CALL NSS3TR(ITRIAL(1,J),NEL,NMIX,MAT,XX,DIFF(:,J),SIGR(:,J), + 1 WORK(:),A11) + B(IOF2+1:IOF2+DIM,IOF1+1:IOF1+DIM)=A11(:,:) + ENDDO + ENDDO + DEALLOCATE(QFR2,A11,WORK) +*---- +* SOLVE EIGENVALUE MATRIX SYSTEM +*---- + CALL ALINV(DIM*NG,A,DIM*NG,IER) + IF(IER.NE.0) CALL XABORT('NSSFL1: SINGULAR MATRIX') + ALLOCATE(AI(DIM*NG,DIM*NG)) + AI(:DIM*NG,:DIM*NG)=MATMUL(A(:DIM*NG,:DIM*NG),B(:DIM*NG,:DIM*NG)) + CALL AL1EIG(DIM*NG,AI,EPSOUT,MAXOUT,ITER,FUNKN,KEFF,IPRINT) + IF(IPRINT.GT.0) WRITE(6,10) KEFF,ITER + DEALLOCATE(AI,B,A) +*---- +* NORMALIZE THE FLUX +*---- + FLMAX=0.0 + DO IG=1,NG + NUM1=0 + DO IEL=1,NEL + IF(ABS(FUNKN(NUM1+1,IG)).GT.ABS(FLMAX)) FLMAX=FUNKN(NUM1+1,IG) + NUM1=NUM1+5 + ENDDO + ENDDO + FUNKN(:,:)=FUNKN(:,:)/FLMAX +*---- +* COMPUTE INTERFACE FLUXES AND CURRENTS +*---- + IOF1=NEL + IOF2=2*NEL + IOF3=3*NEL + IF(IOF3+NEL+1.NE.NUN) CALL XABORT('NSSFL1: NUN ERROR.') + DO IG=1,NG + DO KEL=1,NEL + IBM=MAT(KEL) + IOF=(KEL-1)*5 + EVECT(KEL,IG)=FUNKN(IOF+1,IG) + EVECT(IOF1+KEL,IG)=FUNKN(IOF+1,IG)+0.5*(-FUNKN(IOF+2,IG)+ + 1 FUNKN(IOF+3,IG)) + EVECT(IOF2+KEL,IG)=FUNKN(IOF+1,IG)+0.5*(FUNKN(IOF+2,IG)+ + 1 FUNKN(IOF+3,IG)) + IF(ITRIAL(IBM,IG).EQ.1) THEN + EVECT(IOF3+KEL,IG)=-(DIFF(IBM,IG)/XX(KEL))*(FUNKN(IOF+2,IG)- + 1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)/2.0-FUNKN(IOF+5,IG)/5.0) + ELSE + ETA=XX(KEL)*SQRT(SIGR(IBM,IG)/DIFF(IBM,IG)) + ALP1=ETA*COSH(ETA/2.0)-2.0*SINH(ETA/2.0) + EVECT(IOF1+KEL,IG)=EVECT(IOF1+KEL,IG)-FUNKN(IOF+4,IG)* + 1 SINH(ETA/2.0)+FUNKN(IOF+5,IG)*ALP1/ETA + EVECT(IOF2+KEL,IG)=EVECT(IOF2+KEL,IG)+FUNKN(IOF+4,IG)* + 1 SINH(ETA/2.0)+FUNKN(IOF+5,IG)*ALP1/ETA + EVECT(IOF3+KEL,IG)=-(DIFF(IBM,IG)/XX(KEL))*(FUNKN(IOF+2,IG)- + 1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)*ETA*COSH(ETA/2.0)- + 2 FUNKN(IOF+5,IG)*ETA*SINH(ETA/2.0)) + ENDIF + ENDDO + IBM=MAT(NEL) + IOF=(NEL-1)*5 + IF(ITRIAL(IBM,IG).EQ.1) THEN + EVECT(IOF3+NEL+1,IG)=-(DIFF(IBM,IG)/XX(NEL))*(FUNKN(IOF+2,IG)+ + 1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)/2.0+FUNKN(IOF+5,IG)/5.0) + ELSE + ETA=XX(NEL)*SQRT(SIGR(IBM,IG)/DIFF(IBM,IG)) + EVECT(IOF3+NEL+1,IG)=-(DIFF(IBM,IG)/XX(NEL))*(FUNKN(IOF+2,IG)+ + 1 3.0*FUNKN(IOF+3,IG)+FUNKN(IOF+4,IG)*ETA*COSH(ETA/2.0)+ + 2 FUNKN(IOF+5,IG)*ETA*SINH(ETA/2.0)) + ENDIF + IF(IPRINT.GT.2) THEN + WRITE(6,'(/33H NSSFL1: AVERAGED FLUXES IN GROUP,I5)') IG + WRITE(6,'(1P,10e12.4)') (EVECT(I,IG),I=1,NEL) + WRITE(6,'(/39H NSSFL1: SURFACIC NET CURRENTS IN GROUP,I5)') IG + WRITE(6,'(1P,10e12.4)') (EVECT(IOF3+I,IG),I=1,NEL+1) + ENDIF + ENDDO +*---- +* SAVE SOLUTION +*---- + JPFLX=LCMGID(IPFLX,'FLUX') + DO IG=1,NG + CALL LCMPDL(JPFLX,IG,NUN,2,EVECT(1,IG)) + ENDDO + CALL LCMPUT(IPFLX,'K-EFFECTIVE',1,2,KEFF) + DEALLOCATE(FUNKN,EVECT) + RETURN +* + 10 FORMAT(14H NSSFL1: KEFF=,F11.8,12H OBTAINED IN,I5,11H ITERATIONS) + END -- cgit v1.2.3