Initial commit from Polytechnique Montreal

1 files changed, 305 insertions, 0 deletions

diff --git a/Trivac/src/NSSFL3.f b/Trivac/src/NSSFL3.f
new file mode 100755
index 0000000..61ed29c
--- /dev/null
+++ b/Trivac/src/NSSFL3.f
@@ -0,0 +1,305 @@
+*DECK NSSFL3
+      SUBROUTINE NSSFL3(IPFLX,NUN,NG,NEL,NMIX,NALB,EPSNOD,MAXNOD,
+     1 EPSOUT,MAXOUT,MAT,XX,XXX,IDL,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,
+     2 BETA,FD,IMPX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Flux calculation for the analytic nodal method in Cartesian 1D
+* geometry using the nodal correction iteration strategy.
+*
+*Copyright:
+* Copyright (C) 2022 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 per energy group (=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.
+* EPSNOD  nodal correction epsilon.
+* MAXNOD  maximum number of nodal correction iterations.
+* EPSOUT  convergence epsilon for the power method.
+* MAXOUT  maximum number of iterations for the power method.
+* MAT     material mixtures.
+* XX      mesh spacings.
+* XXX     Cartesian coordinates along the X axis.
+* IDL     position of averaged fluxes in unknown vector.
+* 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.
+* IMPX    edition flag.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPFLX
+      INTEGER NUN,NG,NEL,NMIX,NALB,MAXNOD,MAXOUT,IMPX,MAT(NEL),IDL(NEL),
+     1 IQFR(6,NEL)
+      REAL EPSNOD,EPSOUT,XX(NEL),XXX(NEL+1),QFR(6,NEL),DIFF(NMIX,NG),
+     1 SIGR(NMIX,NG),CHI(NMIX,NG),SIGF(NMIX,NG),SCAT(NMIX,NG,NG),
+     2 BETA(NALB,NG,NG),FD(NMIX,2,NG,NG)
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPFLX
+      INTEGER, PARAMETER :: NY=1,NZ=1,NDIM=1
+      REAL :: COEF(6),CODR(6),KEFF,KEFF_OLD
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: YY,ZZ,EVECT
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: A,SAVG
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: QFR2,DRIFT
+*
+      ALB(X)=0.5*(1.0-X)/(1.0+X)
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      N=NEL*NG
+      ALLOCATE(QFR2(6,NEL,NG),YY(NEL),ZZ(NEL),A(N,2*N),EVECT(N))
+      ALLOCATE(DRIFT(6,NEL,NG),SAVG(NUN,NG))
+*----
+*  ALBEDO PROCESSING
+*----
+      QFR2(:6,:NEL,:NG)=0.0
+      DO IG=1,NG
+        DO IQW=1,2
+          DO IEL=1,NEL
+            IALB=IQFR(IQW,IEL)
+            IF(IALB > 0) THEN
+              IF(IALB.GT.NALB) CALL XABORT('NSSFL3: BETA OVERFLOW.')
+              QFR2(IQW,IEL,IG)=QFR(IQW,IEL)*ALB(BETA(IALB,IG,IG))
+            ELSE
+              QFR2(IQW,IEL,IG)=QFR(IQW,IEL)
+            ENDIF
+          ENDDO
+        ENDDO
+      ENDDO
+*----
+*  INITIALIZATIONS
+*----
+      KEFF_OLD=0.0
+      KEFF=1.0
+      CALL LCMLEN(IPFLX,'FLUX',ILONG,ITYLCM)
+      IF(ILONG == 0) THEN
+        JPFLX=LCMLID(IPFLX,'FLUX',NG)
+        SAVG(:NUN,:NG)=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,SAVG(1,IG))
+        ENDDO
+        CALL LCMGET(IPFLX,'K-EFFECTIVE',KEFF)
+      ENDIF
+      CALL LCMLEN(IPFLX,'DRIFT',ILONG,ITYLCM)
+      IF(ILONG == 0) THEN
+        JPFLX=LCMLID(IPFLX,'DRIFT',6*NEL)
+        DRIFT(:6,:NEL,:NG)=0.0
+      ELSE
+        JPFLX=LCMGID(IPFLX,'DRIFT')
+        DO IG=1,NG
+          CALL LCMLEL(JPFLX,IG,ILONG,ITYLCM)
+          IF(ILONG /= 6*NEL) CALL XABORT('NSSFL3: INVALID DRIFT.')
+          CALL LCMGDL(JPFLX,IG,DRIFT(1,1,IG))
+        ENDDO
+      ENDIF
+      DO IEL=1,NEL
+        DO IG=1,NG
+          EVECT((IG-1)*NEL+IEL)=SAVG(IEL,IG)
+        ENDDO
+      ENDDO
+*----
+*  NODAL CORRECTION LOOP
+*----
+      YY(:NEL)=1.0
+      ZZ(:NEL)=1.0
+      JTER=0
+      DO WHILE((ABS(KEFF_OLD-KEFF) >= EPSNOD).OR.(JTER==0))
+        JTER=JTER+1
+        IF(IMPX > 0) THEN
+          WRITE(6,'(36H NSSFL3: Nodal correction iteration=,I5)')
+     1    JTER
+        ENDIF
+        IF(JTER > MAXNOD) THEN
+          WRITE(6,'(/22H ACCURACY AT ITERATION,I4,2H =,1P,E12.5)')
+     1    JTER,ABS(KEFF_OLD-KEFF)
+          CALL XABORT('NSSFL3: NODAL ITERATION FAILURE')
+        ENDIF
+        !
+        ! set CMFD matrix for x-directed couplings
+        A(:N,:2*N)=0.D0
+        IOF=0
+        DO IG=1,NG
+          DO IEL=1,NEL
+            IBM=MAT(IEL)
+            IF(IBM <= 0) CYCLE
+            KEL=IDL(IEL)
+            IF(KEL == 0) CYCLE
+            VOL0=XX(IEL)
+            CALL NSSCO(NX,NY,NZ,NMIX,IEL,1,1,MAT,XX,YY,ZZ,DIFF(1,IG),
+     >      IQFR(1,IEL),QFR2(1,IEL,IG),COEF)
+            COEF(1:2)=COEF(1:2)*VOL0/XX(IEL)
+            CODR(1:2)=DRIFT(1:2,IEL,IG)*VOL0/XX(IEL)
+            KEL2=0
+            KK1=IQFR(1,IEL)
+            IF(KK1 == -4) THEN
+              KEL2=IDL(NX)
+            ELSE IF(KK1 == 0) THEN
+              KEL2=IDL(IEL-1)
+            ENDIF
+            IF(KEL2 /= 0) THEN
+              A(IOF+KEL,IOF+KEL2)=A(IOF+KEL,IOF+KEL2)-COEF(1)+CODR(1)
+            ENDIF
+            KEL2=0
+            KK2=IQFR(2,IEL)
+            IF(KK2 == -4) THEN
+              KEL2=IDL(1)
+            ELSE IF(KK2 == 0) THEN
+              KEL2=IDL(IEL+1)
+            ENDIF
+            IF(KEL2 /= 0) THEN
+              A(IOF+KEL,IOF+KEL2)=A(IOF+KEL,IOF+KEL2)-COEF(2)-CODR(2)
+            ENDIF
+            A(IOF+KEL,IOF+KEL)=A(IOF+KEL,IOF+KEL)+COEF(1)+CODR(1)+
+     >      COEF(2)-CODR(2)
+            A(IOF+KEL,IOF+KEL)=A(IOF+KEL,IOF+KEL)+SIGR(IBM,IG)*VOL0
+          ENDDO
+          JOF=0
+          DO JG=1,NG ! IG <-- JG
+            DO IEL=1,NEL
+              IBM=MAT(IEL)
+              IF(IBM <= 0) CYCLE
+              KEL=IDL(IEL)
+              IF(KEL == 0) CYCLE
+              IF(IG /= JG) A(IOF+KEL,JOF+KEL)=-XX(IEL)*SCAT(IBM,IG,JG)
+              A(IOF+KEL,N+JOF+KEL)=XX(IEL)*CHI(IBM,IG)*SIGF(IBM,JG)
+            ENDDO
+            JOF=JOF+NEL
+          ENDDO
+          IOF=IOF+NEL
+        ENDDO
+        CALL ALSB(N,N,A,IER,N)
+        IF(IER /= 0) CALL XABORT('NSSFL3: SINGULAR MATRIX.')
+        !
+        ! CMFD power iteration (use double precision)
+        DELTA=ABS(KEFF_OLD-KEFF)
+        KEFF_OLD=KEFF
+        CALL AL1EIG(N,A(1,N+1),EPSOUT,MAXOUT,ITER,EVECT,KEFF,IMPX)
+*----
+*  FLUX NORMALIZATION
+*----
+        FMAX=MAXVAL(EVECT(:N))
+        EVECT(:N)=EVECT(:N)/FMAX
+        IF(IMPX > 0) WRITE(6,10) JTER,KEFF,ITER,DELTA
+        IF(IMPX > 2) THEN
+          WRITE(6,'(1X,A)') 'NSSFL3: EVECT='
+          IOF=0
+          DO IG=1,NG
+            WRITE(6,'(1X,1P,14E12.4)') EVECT(IOF+1:IOF+NEL)
+            IOF=IOF+NEL
+          ENDDO
+        ENDIF
+        !
+        ! begin construct SAVG
+        IF(NUN /= 4*NEL+1) CALL XABORT('NSSFL3: INVALID NUN.')
+        SAVG(:NUN,:NG)=0.0
+        DO IEL=1,NEL
+          DO IG=1,NG
+            SAVG(IEL,IG)=EVECT((IG-1)*NEL+IEL)
+          ENDDO
+        ENDDO
+        !
+        ! one- and two-node anm relations
+        CALL NSSANM1(NEL,NG,NMIX,IQFR,QFR2,MAT,XXX,KEFF,DIFF,SIGR,CHI,
+     1  SIGF,SCAT,FD,SAVG)
+        !
+        ! compute new drift coefficients
+        DO IG=1,NG
+          DO IEL=1,NEL
+            IBM=MAT(IEL)
+            IF(IBM == 0) CYCLE
+            CALL NSSCO(NX,NY,NZ,NMIX,IEL,1,1,MAT,XX,YY,ZZ,DIFF(1,IG),
+     1      IQFR(1,IEL),QFR2(1,IEL,IG),COEF)
+            IF(IEL == 1) THEN
+              DRIFT(1,IEL,IG)=-(SAVG(3*NEL+IEL,IG)+COEF(1)*SAVG(IEL,IG))
+     1        /SAVG(IEL,IG)
+              DRIFT(2,IEL,IG)=-(SAVG(3*NEL+IEL+1,IG)+COEF(2)*
+     1        (SAVG(IEL+1,IG)-SAVG(IEL,IG)))/(SAVG(IEL+1,IG)+
+     2        SAVG(IEL,IG))
+            ELSE IF(IEL < NEL) THEN
+              DRIFT(1,IEL,IG)=-(SAVG(3*NEL+IEL,IG)+COEF(1)*(SAVG(IEL,IG)
+     1        -SAVG(IEL-1,IG)))/(SAVG(IEL,IG)+SAVG(IEL-1,IG))
+              DRIFT(2,IEL,IG)=-(SAVG(3*NEL+IEL+1,IG)+COEF(2)*
+     1        (SAVG(IEL+1,IG)-SAVG(IEL,IG)))/(SAVG(IEL+1,IG)+
+     2        SAVG(IEL,IG))
+            ELSE
+              DRIFT(1,IEL,IG)=-(SAVG(3*NEL+IEL,IG)+COEF(1)*(SAVG(IEL,IG)
+     1        -SAVG(IEL-1,IG)))/(SAVG(IEL,IG)+SAVG(IEL-1,IG))
+              DRIFT(2,IEL,IG)=-(SAVG(3*NEL+IEL+1,IG)-COEF(2)*
+     1        SAVG(IEL,IG))/SAVG(IEL,IG)
+            ENDIF
+          ENDDO
+        ENDDO
+      ENDDO
+*----
+*  END OF NODAL CORRECTION LOOP
+*----
+      IF(IMPX.GT.0) WRITE(6,20) KEFF,JTER
+      IF(IMPX > 2) THEN
+        WRITE(6,'(/21H NSSFL3: UNKNOWNS----)')
+        DO IG=1,NG
+          WRITE(6,'(14H NSSFL3: SAVG(,I4,2H)=)') IG
+          WRITE(6,'(1P,12E12.4)') SAVG(:NEL,IG)
+          WRITE(6,'(19H X-BOUNDARY FLUXES:)')
+          WRITE(6,'(1P,12E12.4)') SAVG(NEL+1:2*NEL,IG)
+          WRITE(6,'(1P,12E12.4)') SAVG(2*NEL+1:3*NEL,IG)
+          WRITE(6,'(12H X-CURRENTS:)')
+          WRITE(6,'(1P,12E12.4)') SAVG(3*NEL+1:,IG)
+          WRITE(6,'(5H ----)')
+        ENDDO
+      ENDIF
+*----
+*  SAVE SOLUTION
+*----
+      JPFLX=LCMGID(IPFLX,'FLUX')
+      DO IG=1,NG
+        CALL LCMPDL(JPFLX,IG,NUN,2,SAVG(1,IG))
+      ENDDO
+      JPFLX=LCMGID(IPFLX,'DRIFT')
+      DO IG=1,NG
+        CALL LCMPDL(JPFLX,IG,6*NEL,2,DRIFT(1,1,IG))
+      ENDDO
+      CALL LCMPUT(IPFLX,'K-EFFECTIVE',1,2,KEFF)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(SAVG,DRIFT)
+      DEALLOCATE(EVECT,A,ZZ,YY,QFR2)
+      RETURN
+*
+   10 FORMAT(14H NSSFL3: JTER=,I4,11H CMFD KEFF=,1P E13.6,
+     1 12H OBTAINED IN,I4,28H CMFD ITERATIONS WITH ERROR=,
+     2 1P,E11.4,1H.)
+   20 FORMAT(18H NSSFL3: ANM KEFF=,F11.8,12H OBTAINED IN,I5,
+     1 28H NODAL CORRECTION ITERATIONS)
+      END