Initial commit from Polytechnique Montreal

1 files changed, 220 insertions, 0 deletions

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