Initial commit from Polytechnique Montreal

1 files changed, 343 insertions, 0 deletions

diff --git a/Trivac/src/NSSDRV.f b/Trivac/src/NSSDRV.f
new file mode 100755
index 0000000..4504051
--- /dev/null
+++ b/Trivac/src/NSSDRV.f
@@ -0,0 +1,343 @@
+*DECK NSSDRV
+      SUBROUTINE NSSDRV(IPTRK,IPMAC,IPFLX,ICHX,IDIM,NUN,NG,NEL,NMIX,
+     1 ITRIAL,ICL1,ICL2,NADI,EPSNOD,MAXNOD,EPSTHR,MAXTHR,EPSOUT,MAXOUT,
+     2 LNODF,BNDTL,NPASS,BB2,IPRINT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Driver for the flux calculation with the nodal expansion method.
+*
+*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
+* IPTRK   nodal tracking.
+* IPMAC   nodal macrolib.
+* IPFLX   nodal flux.
+* ICHX    solution flag (=4.:CMFD; =5: NEM; =6: ANM).
+* IDIM    number of dimensions (1, 2 or 3).
+* NUN     number of unknowns per energy group.
+* NG      number of energy groups.
+* NEL     number of nodes in the nodal calculation.
+* NMIX    number of mixtures in the nodal calculation.
+* ITRIAL  type of expansion functions in the nodal calculation
+*         (=1: polynomial; =2: hyperbolic).
+* ICL1    number of free iterations in one cycle of the inverse power
+*         method (used for thermal iterations).
+* ICL2    number of accelerated iterations in one cycle.
+* NADI    number of inner ADI iterations.
+* EPSNOD  nodal correction epsilon.
+* MAXNOD  maximum number of nodal correction iterations.
+* EPSTHR  thermal iteration epsilon.
+* MAXTHR  maximum number of thermal iterations.
+* EPSOUT  convergence epsilon for the power method.
+* MAXOUT  maximum number of iterations for the power method.
+* LNODF   flag set to .true. to force discontinuity factors to one.
+* BNDTL   set to 'flat', 'linear' or 'quadratic' in 2D cases.
+* BB2     imposed leakage used in non-regression tests.
+* NPASS   number of transverse current iterations.
+* IPRINT  edition flag.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPMAC,IPFLX
+      INTEGER ICHX,IDIM,NUN,NG,NEL,NMIX,ITRIAL(NMIX,NG),ICL1,ICL2,
+     > MAXNOD,NADI,MAXTHR,MAXOUT,NPASS,IPRINT
+      REAL EPSNOD,EPSTHR,EPSOUT,BB2
+      LOGICAL LNODF
+      CHARACTER*12 BNDTL
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (NSTATE=40)
+      INTEGER ISTATE(NSTATE),ICODE(6)
+      TYPE(C_PTR) JPMAC,KPMAC
+      CHARACTER HSMG*131
+      CHARACTER(LEN=8) HADF(6)
+      CHARACTER(LEN=72) TITLE
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: MAT,IJJ,NJJ,IPOS,IDL,MUX,
+     1 MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KN,IQFR
+      REAL, ALLOCATABLE, DIMENSION(:) :: XX,YY,ZZ,XXX,YYY,ZZZ,WORK,VOL
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: DIFF,SIGR,CHI,SIGF,QFR,ALBP,
+     1 GAR2,GAR3
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: BETA,SCAT,FDXM,FDXP,GAR4
+      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: FD
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(MAT(NEL),IDL(NEL),KN(6,NEL),IQFR(6,NEL))
+      ALLOCATE(XX(NEL),YY(NEL),ZZ(NEL),VOL(NEL),DIFF(NMIX,NG),
+     1 SIGR(NMIX,NG),CHI(NMIX,NG),SIGF(NMIX,NG),SCAT(NMIX,NG,NG),
+     2 QFR(6,NEL),FD(NMIX,2*IDIM,NG,NG))
+*----
+*  RECOVER TRACKING INFORMATION
+*----
+      TITLE=' '
+      CALL LCMLEN(IPTRK,'TITLE',LENGT,ITYLCM)
+      IF(LENGT.GT.0) THEN
+        CALL LCMGTC(IPTRK,'TITLE',72,TITLE)
+        IF(IPRINT.GT.0) WRITE(6,'(/9H NSSDRV: ,A72)') TITLE
+      ENDIF
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      NX=ISTATE(14)
+      NY=ISTATE(15)
+      NZ=ISTATE(16)
+      LL4F=ISTATE(25)
+      LL4X=ISTATE(27)
+      LL4Y=ISTATE(28)
+      LL4Z=ISTATE(29)
+      ALLOCATE(MUX(LL4F),MUY(LL4F),MUZ(LL4F),IMAX(LL4F),IMAY(LL4F),
+     1 IMAZ(LL4F),IPY(LL4F),IPZ(LL4F))
+      CALL LCMGET(IPTRK,'ICODE',ICODE)
+      CALL LCMGET(IPTRK,'MATCOD',MAT)
+      CALL LCMGET(IPTRK,'KEYFLX',IDL)
+      CALL LCMGET(IPTRK,'VOLUME',VOL)
+      CALL LCMGET(IPTRK,'XX',XX)
+      CALL LCMGET(IPTRK,'KN',KN)
+      IF(IDIM.GE.2) CALL LCMGET(IPTRK,'YY',YY)
+      IF(IDIM.EQ.3) CALL LCMGET(IPTRK,'ZZ',ZZ)
+      CALL LCMGET(IPTRK,'QFR',QFR)
+      CALL LCMGET(IPTRK,'IQFR',IQFR)
+      ALLOCATE(XXX(NX+1),YYY(NY+1),ZZZ(NZ+1))
+      CALL LCMGET(IPTRK,'XXX',XXX)
+      CALL LCMGET(IPTRK,'MUX',MUX)
+      CALL LCMGET(IPTRK,'IMAX',IMAX)
+      IF(IDIM.GE.2) THEN
+        CALL LCMGET(IPTRK,'YYY',YYY)
+        CALL LCMGET(IPTRK,'MUY',MUY)
+        CALL LCMGET(IPTRK,'IMAY',IMAY)
+        CALL LCMGET(IPTRK,'IPY',IPY)
+      ENDIF
+      IF(IDIM.EQ.3) THEN
+        CALL LCMGET(IPTRK,'ZZZ',ZZZ)
+        CALL LCMGET(IPTRK,'MUZ',MUZ)
+        CALL LCMGET(IPTRK,'IMAZ',IMAZ)
+        CALL LCMGET(IPTRK,'IPZ',IPZ)
+      ENDIF
+*----
+*  RECOVER MACROLIB INFORMATION
+*----
+      IF(BB2.NE.0.0) THEN
+        IF(IPRINT.GT.0) WRITE(6,'(/32H NSSDRV: INCLUDE LEAKAGE IN THE ,
+     >  13HMACROLIB (B2=,1P,E12.5,2H).)') BB2
+      ENDIF
+      CALL LCMGET(IPMAC,'STATE-VECTOR',ISTATE)
+      NALB=ISTATE(8) ! number of physical albedos
+      JPMAC=LCMGID(IPMAC,'GROUP')
+      ALLOCATE(WORK(NMIX*NG),IJJ(NMIX),NJJ(NMIX),IPOS(NMIX))
+      DO IGR=1,NG
+        KPMAC=LCMGIL(JPMAC,IGR)
+        CALL LCMGET(KPMAC,'NTOT0',SIGR(1,IGR))
+        CALL LCMGET(KPMAC,'DIFF',DIFF(1,IGR))
+        CALL LCMGET(KPMAC,'CHI',CHI(1,IGR))
+        CALL LCMGET(KPMAC,'NUSIGF',SIGF(1,IGR))
+        CALL LCMGET(KPMAC,'IJJS00',IJJ)
+        CALL LCMGET(KPMAC,'NJJS00',NJJ)
+        CALL LCMGET(KPMAC,'IPOS00',IPOS)
+        CALL LCMGET(KPMAC,'SCAT00',WORK)
+        DO IBM=1,NMIX
+          SCAT(IBM,IGR,:)=0.0
+          IPOSDE=IPOS(IBM)-1
+          DO JGR=IJJ(IBM),IJJ(IBM)-NJJ(IBM)+1,-1
+            IPOSDE=IPOSDE+1
+            IF(IPOSDE.GT.NMIX*NG) CALL XABORT('NSSDRV: SCAT OVERFLOW.')
+            SCAT(IBM,IGR,JGR)=WORK(IPOSDE) ! IGR <-- JGR
+          ENDDO
+          SIGR(IBM,IGR)=SIGR(IBM,IGR)-SCAT(IBM,IGR,IGR)
+        ENDDO
+        IF(BB2.NE.0.0) THEN
+          DO IBM=1,NMIX
+            SIGR(IBM,IGR)=SIGR(IBM,IGR)+BB2*DIFF(IBM,IGR)
+          ENDDO
+        ENDIF
+        DO IBM=1,NMIX
+          IF(SIGR(IBM,IGR).LE.0.0) CALL XABORT('NSSDRV: SIGR<=0.')
+        ENDDO
+      ENDDO
+      DEALLOCATE(IPOS,NJJ,IJJ,WORK)
+      ALLOCATE(FDXM(NMIX,NG,NG),FDXP(NMIX,NG,NG),BETA(NALB,NG,NG),
+     > GAR2(NMIX,NG),GAR3(NMIX,NG),GAR4(NMIX,NG,NG))
+      IF(NALB.GT.0) THEN
+        CALL LCMLEN(IPMAC,'ALBEDO',ILONG,ITYLCM)
+        IF(ILONG.EQ.NALB*NG) THEN
+          ALLOCATE(ALBP(NALB,NG))
+          CALL LCMGET(IPMAC,'ALBEDO',ALBP)
+          BETA(:,:,:)=1.0
+          DO IGR=1,NG
+            BETA(:NALB,IGR,IGR)=ALBP(:NALB,IGR)
+          ENDDO
+          DEALLOCATE(ALBP)
+        ELSE IF(ILONG.EQ.NALB*NG*NG) THEN
+          CALL LCMGET(IPMAC,'ALBEDO',BETA)
+        ELSE
+          CALL XABORT('NSSDRV: INVALID ALBEDO LENGTH.')
+        ENDIF
+        IF(IPRINT.GT.1) THEN
+          DO IALB=1,NALB
+            WRITE(6,'(/35H NSSDRV: PHYSICAL ALBEDO MATRIX ID=,I4)') IALB
+            DO IGR=1,NG
+              WRITE(6,'(5X,1P,10E12.4)') BETA(IALB,IGR,:)
+            ENDDO
+          ENDDO
+        ENDIF
+      ENDIF
+      FD(:,:,:,:)=0.0
+      IF(LNODF.OR.ISTATE(12).EQ.0) THEN
+        DO IGR=1,NG
+          FD(:NMIX,:2*IDIM,IGR,IGR)=1.0
+        ENDDO
+      ELSE IF(ISTATE(12).EQ.2) then
+        CALL LCMSIX(IPMAC,'ADF',1)
+        CALL LCMGET(IPMAC,'NTYPE',NSURFD)
+        CALL LCMGET(IPMAC,'AVG_FLUX',GAR3)
+        CALL LCMGTC(IPMAC,'HADF',8,NSURFD,HADF)
+        IF(NSURFD.EQ.1) THEN
+          CALL LCMGET(IPMAC,HADF(1),GAR2)
+          DO IBM=1,NMIX
+            DO IGR=1,NG
+              FD(IBM,:2*IDIM,IGR,IGR)=GAR2(IBM,IGR)/GAR3(IBM,IGR)
+            ENDDO
+          ENDDO
+        ELSE IF(NSURFD.EQ.2*IDIM) THEN
+          DO I=1,NSURFD
+            IF(HADF(I)(1:3).NE.'FD_') THEN
+              WRITE(HSMG,'(7HNSSDRV:,A,28H FOUND; FD_ PREFIX EXPECTED.)
+     1        ') TRIM(HADF(I))
+              CALL XABORT(HSMG)
+            ENDIF
+            CALL LCMGET(IPMAC,HADF(I),GAR2)
+            DO IGR=1,NG
+              FD(:NMIX,I,IGR,IGR)=GAR2(:NMIX,IGR)/GAR3(:NMIX,IGR)
+            ENDDO
+          ENDDO
+        ELSE
+          WRITE(HSMG,'(12HNSSDRV: 1 OR,I3,25HDISCONTINUITY FACTORS EXP,
+     1    6HECTED.)') 2*IDIM
+          CALL XABORT(HSMG)
+        ENDIF
+        CALL LCMSIX(IPMAC,' ',2)
+      ELSE IF(ISTATE(12).EQ.3) THEN
+        CALL LCMSIX(IPMAC,'ADF',1)
+        CALL LCMGET(IPMAC,'NTYPE',NSURFD)
+        CALL LCMGTC(IPMAC,'HADF',8,NSURFD,HADF)
+        IF(NSURFD.EQ.1) THEN
+          CALL LCMGET(IPMAC,HADF(1),GAR2)
+          DO IBM=1,NMIX
+            DO IGR=1,NG
+              FD(IBM,:2*IDIM,IGR,IGR)=GAR2(IBM,IGR)
+            ENDDO
+          ENDDO
+        ELSE IF(NSURFD.EQ.2*IDIM) THEN
+          DO I=1,NSURFD
+            IF(HADF(I)(1:3).NE.'FD_') THEN
+              WRITE(HSMG,'(7HNSSDRV:,A,28H FOUND; FD_ PREFIX EXPECTED.)
+     1        ') TRIM(HADF(I))
+              CALL XABORT(HSMG)
+            ENDIF
+            CALL LCMGET(IPMAC,HADF(I),GAR2)
+            DO IGR=1,NG
+              FD(:NMIX,I,IGR,IGR)=GAR2(:NMIX,IGR)
+            ENDDO
+          ENDDO
+        ELSE
+          WRITE(HSMG,'(12HNSSDRV: 1 OR,I3,25HDISCONTINUITY FACTORS EXP,
+     1    6HECTED.)') 2*IDIM
+          CALL XABORT(HSMG)
+        ENDIF
+        CALL LCMSIX(IPMAC,' ',2)
+      ELSE IF(ISTATE(12).EQ.4) THEN
+        ! matrix discontinuity factors
+        CALL LCMSIX(IPMAC,'ADF',1)
+        CALL LCMGET(IPMAC,'NTYPE',NSURFD)
+        IF(NSURFD.NE.2*IDIM) THEN
+          WRITE(HSMG,'(7HNSSDRV:,I3,30HDISCONTINUITY FACTORS EXPECTED)'
+     1    ) 2*IDIM
+          CALL XABORT(HSMG)
+        ENDIF
+        CALL LCMGTC(IPMAC,'HADF',8,NSURFD,HADF)
+        DO I=1,NSURFD
+          IF((HADF(I)(1:4).NE.'ERM_').AND.(HADF(I)(1:3).NE.'FD_')) THEN
+            WRITE(HSMG,'(7HNSSDRV:,A,30H FOUND; ERM_ OR FD_ PREFIX EXP,
+     1      6HECTED.)') TRIM(HADF(I))
+            CALL XABORT(HSMG)
+          ENDIF
+          CALL LCMGET(IPMAC,HADF(I),GAR4)
+          DO JGR=1,NG
+            DO IGR=1,NG
+              FD(:NMIX,I,IGR,JGR)=GAR4(:NMIX,IGR,JGR)
+            ENDDO
+          ENDDO
+        ENDDO
+        CALL LCMSIX(IPMAC,' ',2)
+      ELSE
+        WRITE(6,'(13H NSSDRV: IDF=,I3)') ISTATE(12)
+        CALL XABORT('NSSDRV: FLUX/CURRENT INFORMATION NOT SUPPORTED.')
+      ENDIF
+      IF(IPRINT.GT.3) THEN
+        DO I=1,NSURFD
+          WRITE(6,'(/31H NSSDRV: discontinuity factors ,A8)') HADF(I)
+          DO IBM=1,NMIX
+            DO IGR=1,NG
+              WRITE(6,'(4H FD(,2I4,2H)=,1p,12E12.4/(8X,12E12.4))')
+     1        IBM,IGR,FD(IBM,:,IGR,IGR)
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDIF
+      DEALLOCATE(GAR4,GAR3,GAR2,FDXP,FDXM)
+*----
+*  COMPUTE THE FLUX AND STORE NODAL SOLUTION IN IPFLX
+*----
+      IF(ICHX.EQ.5) THEN ! NEM
+        CALL NSSFL1(IPFLX,NUN,NG,NEL,NMIX,NALB,ITRIAL,EPSOUT,MAXOUT,
+     1  MAT,XX,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,FD,IPRINT)
+      ELSE IF(ICHX.EQ.4) THEN ! CMFD
+        CALL NSSFL2(IPFLX,NUN,NG,NEL,NMIX,NALB,EPSOUT,MAXOUT,MAT,XX,
+     1  IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,FD,IPRINT)
+      ELSE IF((ICHX.EQ.6).AND.(IDIM.EQ.1)) THEN ! ANM-1D
+        CALL NSSFL3(IPFLX,NUN,NG,NEL,NMIX,NALB,EPSNOD,MAXNOD,EPSOUT,
+     1  MAXOUT,MAT,XX,XXX,IDL,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,
+     2  FD,IPRINT)
+      ELSE IF((ICHX.EQ.6).AND.(IDIM.EQ.2)) THEN ! ANM-2D
+        CALL NSSFL4(IPFLX,NUN,NG,NX,NY,LL4F,LL4X,LL4Y,NMIX,NALB,ICL1,
+     1  ICL2,NADI,EPSNOD,MAXNOD,EPSTHR,MAXTHR,EPSOUT,MAXOUT,MAT,XX,YY,
+     2  XXX,YYY,IDL,VOL,KN,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,FD,
+     3  BNDTL,NPASS,MUX,MUY,IMAX,IMAY,IPY,IPRINT)
+      ELSE IF((ICHX.EQ.6).AND.(IDIM.EQ.3)) THEN ! ANM-3D
+        CALL NSSFL5(IPFLX,NUN,NG,NX,NY,NZ,LL4F,LL4X,LL4Y,LL4Z,NMIX,NALB,
+     1  ICL1,ICL2,NADI,EPSNOD,MAXNOD,EPSTHR,MAXTHR,EPSOUT,MAXOUT,MAT,XX,
+     2  YY,ZZ,XXX,YYY,ZZZ,IDL,VOL,KN,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,
+     3  BETA,FD,BNDTL,NPASS,MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,IPRINT)
+      ELSE
+        CALL XABORT('NSSDRV: OPTION NOT AVAILABLE.')
+      ENDIF
+      ISTATE(:)=0
+      ISTATE(1)=NG
+      ISTATE(2)=NUN
+      ISTATE(6)=2
+      CALL LCMPUT(IPFLX,'STATE-VECTOR',NSTATE,1,ISTATE)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IPZ,IPY,IMAZ,IMAY,IMAX,MUZ,MUY,MUX)
+      DEALLOCATE(BETA)
+      DEALLOCATE(ZZZ,YYY,XXX)
+      DEALLOCATE(FD,QFR,SCAT,SIGF,CHI,SIGR,DIFF,VOL,ZZ,YY,XX)
+      DEALLOCATE(IQFR,KN,IDL,MAT)
+      RETURN
+      END