Initial commit from Polytechnique Montreal

1 files changed, 357 insertions, 0 deletions

diff --git a/Trivac/src/NSSFL4.f b/Trivac/src/NSSFL4.f
new file mode 100755
index 0000000..bdf4e73
--- /dev/null
+++ b/Trivac/src/NSSFL4.f
@@ -0,0 +1,357 @@
+*DECK NSSFL4
+      SUBROUTINE NSSFL4(IPFLX,NUN,NG,NX,NY,LL4F,LL4X,LL4Y,NMIX,NALB,
+     > ICL1,ICL2,NADI,EPSNOD,MAXNOD,EPSTHR,MAXTHR,EPSOUT,MAXOUT,MAT,
+     > XX,YY,XXX,YYY,IDL,VOL,KN,IQFR,QFR,DIFF,SIGR,CHI,SIGF,SCAT,BETA,
+     > FD,BNDTL,NPASS,MUX,MUY,IMAX,IMAY,IPY,IMPX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Flux calculation for the analytic nodal method in Cartesian 2D
+* 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.
+* NG      number of energy groups.
+* NX      number of nodes in the X direction.
+* NY      number of nodes in the Y direction.
+* LL4F    number of nodal flux unknowns.
+* LL4X    number of nodal X-directed net currents unknowns.
+* LL4Y    number of nodal Y-directed net currents unknowns.
+* NMIX    number of mixtures in the nodal calculation.
+* NALB    number of physical albedos.
+* 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.
+* MAT     material mixtures.
+* XX      mesh spacings in the X direction.
+* YY      mesh spacings in the Y direction.
+* XXX     Cartesian coordinates along the X axis.
+* YYY     Cartesian coordinates along the Y axis.
+* IDL     position of averaged fluxes in unknown vector.
+* VOL     node volumes.
+* KN      node-ordered interface net current unknown list.
+* 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.
+* BNDTL   set to 'flat' or 'quadratic'.
+* NPASS   number of transverse current iterations.
+* MUX     X-oriented compressed storage mode indices.
+* MUY     Y-oriented compressed storage mode indices.
+* IMAX    X-oriented position of each first non-zero column element.
+* IMAY    Y-oriented position of each first non-zero column element.
+* IPY     Y-oriented permutation matrices.
+* IMPX    edition flag.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPFLX
+      INTEGER NUN,NG,NX,NY,LL4F,LL4X,LL4Y,NMIX,NALB,ICL1,ICL2,NADI,
+     1 MAXNOD,MAXTHR,MAXOUT,IMPX,MAT(NX*NY),IDL(NX*NY),KN(6,NX,NY),
+     2 NPASS,IQFR(6,NX,NY),MUX(LL4F),MUY(LL4F),IMAX(LL4F),IMAY(LL4F),
+     3 IPY(LL4F)
+      REAL EPSNOD,EPSTHR,EPSOUT,XX(NX*NY),YY(NX*NY),XXX(NX+1),YYY(NY+1),
+     1 VOL(NX*NY),QFR(6,NX*NY),DIFF(NMIX,NG),SIGR(NMIX,NG),CHI(NMIX,NG),
+     2 SIGF(NMIX,NG),SCAT(NMIX,NG,NG),BETA(NALB,NG,NG),FD(NMIX,4,NG,NG)
+      CHARACTER(LEN=12) :: BNDTL
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPFLX
+      INTEGER, PARAMETER :: NZ=1,NDIM=2
+      INTEGER :: MUZ(1),IMAZ(1),IPZ(1)
+      REAL :: COEF(6),KEFF,KEFF_OLD
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: ZZ,EVECT
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: A11X,A11Y,A11Z,SAVG
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: QFR2,DRIFT
+*
+      ALB(X)=0.5*(1.0-X)/(1.0+X)
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      NEL=NX*NY
+      N=LL4F*NG
+      ALLOCATE(QFR2(6,NEL,NG),ZZ(NEL),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,4
+          DO I=1,NX
+            DO J=1,NY
+              IEL=(J-1)*NX+I
+              IALB=IQFR(IQW,I,J)
+              IF(IALB > 0) THEN
+                IF(IALB.GT.NALB) CALL XABORT('NSSFL4: 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
+      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('NSSFL4: 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('NSSFL4: INVALID DRIFT.')
+          CALL LCMGDL(JPFLX,IG,DRIFT(1,1,IG))
+        ENDDO
+      ENDIF
+      DO IEL=1,LL4F
+        DO IG=1,NG
+          EVECT((IG-1)*LL4F+IEL)=SAVG(IEL,IG)
+        ENDDO
+      ENDDO
+*----
+*  NODAL CORRECTION LOOP
+*----
+      NMAX=IMAX(LL4F)
+      NMAY=IMAY(LL4F)
+      NMAZ=1
+      ALLOCATE(A11X(NMAX,NG),A11Y(NMAY,NG),A11Z(NMAZ,NG))
+      ZZ(:NEL)=1.0
+      MUZ(1)=0
+      IMAZ(1)=0
+      IPZ(1)=0
+      JTER=0
+      SAVG(:NUN,:NG)=0.0
+      IOFY=5*LL4F+LL4X
+      DO WHILE((ABS(KEFF_OLD-KEFF) >= EPSNOD).OR.(JTER==0))
+        JTER=JTER+1
+        IF(IMPX.GT.0) THEN
+          WRITE(6,'(36H NSSFL4: Nodal correction iteration=,I5)')
+     >    JTER
+        ENDIF
+        IF(JTER > MAXNOD) THEN
+          WRITE(6,'(/22H ACCURACY AT ITERATION,I4,2H =,1P,E12.5)')
+     >    JTER,ABS(KEFF_OLD-KEFF)
+          CALL XABORT('NSSFL4: NODAL ITERATION FAILURE')
+        ENDIF
+        !
+        ! set coarse mesh finite difference matrix
+        IOF=0
+        DO IG=1,NG
+          CALL NSSMXYZ(LL4F,NDIM,NX,NY,NZ,NMIX,MAT,XX,YY,ZZ,IDL,VOL,
+     >    IQFR,QFR2(1,1,IG),DIFF(1,IG),DRIFT(1,1,IG),SIGR(1,IG),
+     >    MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,A11X(1,IG),A11Y(1,IG),
+     >    A11Z(1,IG))
+        ENDDO
+        !
+        ! CMFD power iteration
+        DELTA=ABS(KEFF_OLD-KEFF)
+        KEFF_OLD=KEFF
+        CALL NSSEIG(NMAX,NMAY,NMAZ,LL4F,NDIM,NEL,NMIX,NG,MAT,IDL,VOL,
+     >  MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ,CHI,SIGF,SCAT,A11X,A11Y,A11Z,
+     >  EPSTHR,MAXTHR,NADI,EPSOUT,MAXOUT,ICL1,ICL2,ITER,EVECT,KEFF,IMPX)
+        IF(IMPX > 0) WRITE(6,10) JTER,KEFF,ITER,DELTA
+        IF(IMPX > 2) THEN
+          WRITE(6,'(1X,A)') 'NSSFL4: EVECT='
+          IOF=0
+          DO IG=1,NG
+            WRITE(6,'(1X,1P,14E12.4)') EVECT(IOF+1:IOF+LL4F)
+            IOF=IOF+LL4F
+          ENDDO
+        ENDIF
+        !
+        ! begin construct SAVG
+        IF(NUN /= IOFY+LL4Y) CALL XABORT('NSSFL4: INVALID NUN.')
+        DO IND1=1,LL4F
+          DO IG=1,NG
+            SAVG(IND1,IG)=EVECT((IG-1)*LL4F+IND1)
+          ENDDO
+        ENDDO
+        !
+        ! one- and two-node anm relations
+        CALL NSSANM2(NUN,NX,NY,LL4F,LL4X,NG,BNDTL,NPASS,NMIX,IDL,KN,
+     >  IQFR,QFR2,MAT,XXX,YYY,KEFF,DIFF,SIGR,CHI,SIGF,SCAT,FD,SAVG)
+        !
+        ! compute new drift coefficients
+        DRIFT(:6,:NEL,:NG)=0.0
+        DO IG=1,NG
+          DO J=1,NY
+            DO I=1,NX
+              IEL=(J-1)*NX+I
+              IND1=IDL(IEL)
+              IF(IND1 == 0) CYCLE
+              KK1=IQFR(1,I,J)
+              KK2=IQFR(2,I,J)
+              JXM=KN(1,I,J) ; JXP=KN(2,I,J)
+              JYM=KN(3,I,J) ; JYP=KN(4,I,J)
+              CALL NSSCO(NX,NY,NZ,NMIX,I,J,1,MAT,XX,YY,ZZ,DIFF(1,IG),
+     >        IQFR(1,I,J),QFR2(1,IEL,IG),COEF)
+              IF((KK1 < 0).AND.(KK2 < 0)) THEN
+                DRIFT(1,IEL,IG)=-(SAVG(5*LL4F+JXM,IG)+COEF(1)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+                DRIFT(2,IEL,IG)=-(SAVG(5*LL4F+JXP,IG)-COEF(2)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+              ELSE IF(KK1 < 0) THEN
+                DRIFT(1,IEL,IG)=-(SAVG(5*LL4F+JXM,IG)+COEF(1)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+                IND3=IDL((J-1)*NX+I+1)
+                IF(IND3 /= 0) DRIFT(2,IEL,IG)=-(SAVG(5*LL4F+JXP,IG)+
+     >          COEF(2)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+
+     >          SAVG(IND1,IG))
+              ELSE IF(KK2 < 0) THEN
+                IND2=IDL((J-1)*NX+I-1)
+                IF(IND2 /= 0) DRIFT(1,IEL,IG)=-(SAVG(5*LL4F+JXM,IG)+
+     >          COEF(1)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+
+     >          SAVG(IND2,IG))
+                DRIFT(2,IEL,IG)=-(SAVG(5*LL4F+JXP,IG)-COEF(2)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+              ELSE
+                IND2=IDL((J-1)*NX+I-1)
+                IND3=IDL((J-1)*NX+I+1)
+                IF(IND2 /= 0) DRIFT(1,IEL,IG)=-(SAVG(5*LL4F+JXM,IG)+
+     >          COEF(1)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+
+     >          SAVG(IND2,IG))
+                IF(IND3 /= 0) DRIFT(2,IEL,IG)=-(SAVG(5*LL4F+JXP,IG)+
+     >          COEF(2)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+
+     >          SAVG(IND1,IG))
+              ENDIF
+              KK3=IQFR(3,I,J)
+              KK4=IQFR(4,I,J)
+              IF((KK3 < 0).AND.(KK4 < 0)) THEN
+                DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+COEF(3)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+                DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)-COEF(4)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+              ELSE IF(KK3 < 0) THEN
+                DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+COEF(3)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+                IND3=IDL(J*NX+I)
+                IF(IND3 /= 0) DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)+
+     >          COEF(4)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+
+     >          SAVG(IND1,IG))
+              ELSE IF(KK4 < 0) THEN
+                IND2=IDL((J-2)*NX+I)
+                IF(IND2 /= 0) DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+
+     >          COEF(3)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+
+     >          SAVG(IND2,IG))
+                DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)-COEF(4)*
+     >          SAVG(IND1,IG))/SAVG(IND1,IG)
+              ELSE
+                IND2=IDL((J-2)*NX+I)
+                IND3=IDL(J*NX+I)
+                IF(IND2 /= 0) DRIFT(3,IEL,IG)=-(SAVG(IOFY+JYM,IG)+
+     >          COEF(3)*(SAVG(IND1,IG)-SAVG(IND2,IG)))/(SAVG(IND1,IG)+
+     >          SAVG(IND2,IG))
+                IF(IND3 /= 0) DRIFT(4,IEL,IG)=-(SAVG(IOFY+JYP,IG)+
+     >          COEF(4)*(SAVG(IND3,IG)-SAVG(IND1,IG)))/(SAVG(IND3,IG)+
+     >          SAVG(IND1,IG))
+              ENDIF
+            ENDDO
+          ENDDO
+        ENDDO
+        IF(IMPX > 5) THEN
+          DO IG=1,NG
+            WRITE(6,'(28H NSSFL4: DRIFT COEFFICIENTS(,I5,2H):)') IG
+            DO IEL=1,NX*NY
+              WRITE(6,'(1P,I7,4E12.4)') IEL,DRIFT(:4,IEL,IG)
+            ENDDO
+          ENDDO
+        ENDIF
+      ENDDO
+      DEALLOCATE(A11Z,A11Y,A11X)
+*----
+*  END OF NODAL CORRECTION LOOP
+*----
+      IF(IMPX.GT.0) WRITE(6,20) KEFF,JTER
+      IF(IMPX > 2) THEN
+        WRITE(6,'(/21H NSSFL4: UNKNOWNS----)')
+        DO IG=1,NG
+          WRITE(6,'(14H NSSFL4: SAVG(,I4,2H)=)') IG
+          WRITE(6,'(1P,12E12.4)') SAVG(:LL4F,IG)
+          WRITE(6,'(19H X-BOUNDARY FLUXES:)')
+          WRITE(6,'(1P,12E12.4)') SAVG(LL4F+1:2*LL4F,IG)
+          WRITE(6,'(1P,12E12.4)') SAVG(2*LL4F+1:3*LL4F,IG)
+          WRITE(6,'(19H Y-BOUNDARY FLUXES:)')
+          WRITE(6,'(1P,12E12.4)') SAVG(3*LL4F+1:4*LL4F,IG)
+          WRITE(6,'(1P,12E12.4)') SAVG(4*LL4F+1:5*LL4F,IG)
+          WRITE(6,'(12H X-CURRENTS:)')
+          WRITE(6,'(1P,12E12.4)') SAVG(5*LL4F+1:IOFY,IG)
+          WRITE(6,'(12H Y-CURRENTS:)')
+          WRITE(6,'(1P,12E12.4)') SAVG(IOFY+1:NUN,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,ZZ,QFR2)
+      RETURN
+*
+   10 FORMAT(14H NSSFL4: 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 NSSFL4: ANM KEFF=,F11.8,12H OBTAINED IN,I5,
+     1 11H ITERATIONS)
+      END