Initial commit from Polytechnique Montreal

1 files changed, 257 insertions, 0 deletions

diff --git a/Dragon/src/SNF.f b/Dragon/src/SNF.f
new file mode 100644
index 0000000..f8a1c90
--- /dev/null
+++ b/Dragon/src/SNF.f
@@ -0,0 +1,257 @@
+*DECK SNF
+      SUBROUTINE SNF(KPSYS,IPTRK,IFTRAK,IMPX,NGEFF,NGIND,IDIR,NREG,
+     1 NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,TITR,NBS,KPSOU1,KPSOU2,
+     2 FLUXC,EVALRHO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solve N-group transport equation for fluxes using the discrete
+* ordinates (SN) method.
+*
+*Copyright:
+* Copyright (C) 2007 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
+* KPSYS   pointer to the assembly LCM object (L_PIJ signature). KPSYS is
+*         an array of directories.
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IFTRAK  not used.
+* IMPX    print flag (equal to zero for no print).
+* NGEFF   number of energy groups processed in parallel.
+* NGIND   energy group indices assign to the NGEFF set.
+* IDIR    not used.
+* NREG    total number of regions for which specific values of the
+*         neutron flux and reactions rates are required.
+* NBMIX   number of mixtures.
+* NUN     total number of unknowns in vectors SUNKNO and FUNKNO.
+* MAT     index-number of the mixture type assigned to each volume.
+* VOL     volumes.
+* KEYFLX  position of averaged flux elements in FUNKNO vector.
+* SUNKNO  input source vector.
+* TITR    title.
+* NBS
+* KPSOU1
+* KPSOU2
+*
+*Parameters: input/output
+* FUNKNO  unknown vector.
+* FLUXC   flux at the cutoff energy.
+* EVALRHO dominance ratio.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      CHARACTER   TITR*72
+      TYPE(C_PTR) KPSYS(NGEFF),IPTRK,KPSOU1(NGEFF),KPSOU2(NGEFF)
+      INTEGER     NGEFF,NGIND(NGEFF),IFTRAK,IMPX,IDIR,NREG,NBMIX,NUN,
+     1            MAT(NREG),KEYFLX(NREG),NBS(NGEFF)
+      REAL        VOL(NREG),FUNKNO(NUN,NGEFF),SUNKNO(NUN,NGEFF)
+      REAL,OPTIONAL :: FLUXC(NREG)
+      REAL,OPTIONAL :: EVALRHO
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER  (IUNOUT=6,NSTATE=40)
+      INTEGER     IPAR(NSTATE)
+      LOGICAL     LIVO
+      DOUBLE PRECISION F1,F2,R1,R2,DMU
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: FGAR,TEST
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: OLD1,OLD2,OLD3
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: INCONV
+*----
+*  RECOVER SN SPECIFIC PARAMETERS
+*----
+      IF(IMPX.GT.2) THEN
+        WRITE(IUNOUT,'(//6H SNF: ,A72)') TITR
+        CALL KDRCPU(TK1)
+      ENDIF
+      IF(IDIR.NE.0) CALL XABORT('SNF: EXPECTING IDIR=0')
+      IF(IFTRAK.NE.0) CALL XABORT('SNF: EXPECTING IFTRAK=0')
+      CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR)
+      NSTART=IPAR(20)
+      MAXIT=IPAR(22)
+      LIVO=(IPAR(23).EQ.1)
+      ICL1=IPAR(24)
+      ICL2=IPAR(25)
+      IBFP=IPAR(31)
+      NFOU=IPAR(34)
+      CALL LCMGET(IPTRK,'EPSI',EPSINR)
+      IF(IMPX.GT.3) THEN
+        ALLOCATE(FGAR(NREG))
+        DO II=1,NGEFF
+          FGAR(:NREG)=0.0
+          DO I=1,NREG
+            IF(KEYFLX(I).NE.0) FGAR(I)=SUNKNO(KEYFLX(I),II)
+          ENDDO
+          WRITE(IUNOUT,'(/33H N E U T R O N    S O U R C E S (,I5,
+     1    3H ):)') NGIND(II)
+          WRITE(IUNOUT,'(1P,6(5X,E15.7))') (FGAR(I),I=1,NREG)
+        ENDDO
+        DEALLOCATE(FGAR)
+      ENDIF
+*
+      IF(NSTART.GT.0) THEN
+*----
+*  GMRES(M) INNER ITERATION LOOP FOR ONE-SPEED TRANSPORT EQUATION
+*----
+         CALL SNGMRE (KPSYS,NGIND,IPTRK,IMPX,NGEFF,NREG,NBMIX,NUN,
+     1   NSTART,MAXIT,EPSINR,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,NBS,KPSOU1,
+     2   KPSOU2,FLUXC)
+      ELSE
+*----
+*  LIVOLANT INNER ITERATION LOOP FOR ONE-SPEED TRANSPORT EQUATION
+*----
+         LNCONV=NGEFF
+         ALLOCATE(INCONV(NGEFF),OLD1(NUN,NGEFF),OLD2(NUN,NGEFF),
+     1   TEST(NGEFF),OLD3(NUN,NGEFF))
+*
+         INCONV(:NGEFF)=.TRUE.
+         TEST(:NGEFF)=0.0
+*
+         OLD1(:NUN,:NGEFF)=0.0
+         OLD2(:NUN,:NGEFF)=0.0
+         IF(NFOU.GT.0) OLD3(:NUN,:NGEFF)=0.0
+*
+         ITER=0
+   10    ITER=ITER+1
+         IF(ITER.GT.MAXIT) THEN
+            WRITE(IUNOUT,'(40H SNF: MAXIMUM NUMBER OF ONE-SPEED ITERAT,
+     1      12HION REACHED.)')
+            IF(NFOU.GT.0) WRITE(6,310) TRHO1
+            GO TO 70
+         ENDIF
+*
+         IF(NFOU.GT.0)THEN
+            OLD1(:NUN,:NGEFF)=  OLD2(:NUN,:NGEFF)
+            OLD2(:NUN,:NGEFF)=  OLD3(:NUN,:NGEFF)
+            OLD3(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF)
+         ELSE
+            OLD1(:NUN,:NGEFF)=  OLD2(:NUN,:NGEFF)
+            OLD2(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF)
+         ENDIF
+*----
+*  UPDATE THE FIXED SOURCE AND COMPUTE THE FLUX
+*----
+         CALL SNFLUX(KPSYS,INCONV,NGIND,IPTRK,IMPX,NGEFF,NREG,
+     1   NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,ITER,NBS,KPSOU1,
+     2   KPSOU2,FLUXC)
+*----
+*  LOOP OVER ENERGY GROUPS
+*----
+         DO 60 II=1,NGEFF
+         IF(INCONV(II)) THEN
+*----
+*  FOURIER ANALYSIS NUMERICAL EIGENVALUE CALCULATION
+*----
+           IF(NFOU.GT.0)THEN
+           TRHO1=0.0
+           IF(ITER.GT.3)THEN
+             TOT1=0.0
+             TOT2=0.0
+             DO IR=1,NREG
+               IND=KEYFLX(IR)
+               IF(IND.GT.0)THEN
+                 TOT1=TOT1+(FUNKNO(IND,II)-OLD3(IND,II))**2
+                 TOT2=TOT2+(  OLD2(IND,II)-OLD1(IND,II))**2
+               ENDIF
+             ENDDO
+             TOT1 = SQRT(TOT1)
+             TOT2 = SQRT(TOT2)
+             TRHO1 = SQRT(TOT1/TOT2)
+             EVALRHO=TRHO1
+           ENDIF
+           ENDIF
+*----
+*  VARIATIONAL ACCELERATION. LIVOLANT INNER ITERATION LOOP FOR ONE-GROUP
+*  TRANSPORT EQUATION.
+*----
+           DMU=1.0D0
+           IF(LIVO.AND.(MOD(ITER-1,ICL1+ICL2).GE.ICL1)) THEN
+              F1=0.0
+              F2=0.0
+              DO 30 I=1,NUN
+              R1=OLD2(I,II)-OLD1(I,II)
+              R2=FUNKNO(I,II)-OLD2(I,II)
+              F1=F1+R1*(R2-R1)
+              F2=F2+(R2-R1)*(R2-R1)
+   30         CONTINUE
+              DMU=-F1/F2
+              IF(DMU.GT.0.0) THEN
+                RDMU=REAL(DMU)
+                DO 40 I=1,NUN
+                FUNKNO(I,II)=OLD2(I,II)+RDMU*(FUNKNO(I,II)-OLD2(I,II))
+                OLD2(I,II)=OLD1(I,II)+RDMU*(OLD2(I,II)-OLD1(I,II))
+   40           CONTINUE
+              ENDIF
+           ENDIF
+*----
+*  CALCULATE ERROR AND TEST FOR CONVERGENCE
+*----
+           AAA=0.0
+           BBB=0.0
+           DO 50 I=1,NREG
+           IF(KEYFLX(I).EQ.0) GO TO 50
+           AAA=MAX(AAA,ABS(FUNKNO(KEYFLX(I),II)-OLD2(KEYFLX(I),II)))
+           BBB=MAX(BBB,ABS(FUNKNO(KEYFLX(I),II)))
+   50      CONTINUE
+           IF(IMPX.GT.2) WRITE(IUNOUT,300) NGIND(II),ITER,AAA,BBB,
+     1     AAA/BBB,DMU
+           IF(IMPX.GT.5) THEN
+             ALLOCATE(FGAR(NREG))
+             FGAR(:NREG)=0.0
+             DO I=1,NREG
+              IF(KEYFLX(I).NE.0) FGAR(I)=FUNKNO(KEYFLX(I),II)
+             ENDDO
+             WRITE(IUNOUT,'(//33H N E U T R O N    F L U X E S   :)')
+             WRITE(IUNOUT,'(1P,6(5X,E15.7))') (FGAR(I),I=1,NREG)
+             DEALLOCATE(FGAR)
+           ENDIF
+           IF(AAA.LE.0.1*EPSINR*BBB) THEN
+             LNCONV=LNCONV-1
+             INCONV(II)=.FALSE.
+           ENDIF
+           IF(ITER.EQ.1) TEST(II)=AAA
+!          Be careful if the value of ITER is changed below.
+           IF((ITER.GE.10).AND.(AAA.GT.TEST(II))) THEN
+             WRITE(IUNOUT,'(39H SNF: UNABLE TO CONVERGE ONE-SPEED ITER,
+     1       15HATIONS IN GROUP,I5,1H.)') NGIND(II)
+             LNCONV=LNCONV-1
+             INCONV(II)=.FALSE.
+           ENDIF
+         ENDIF
+   60    CONTINUE
+         IF((NFOU.GT.0).AND.(LNCONV.EQ.0)) WRITE(6,310) TRHO1
+         IF(LNCONV.EQ.0) GO TO 70
+         GO TO 10
+*----
+* CONVERGENCE OF ONE-SPEED ITERATIONS IN ALL NGEFF GROUPS
+*----
+   70    IF(IMPX.GT.1) WRITE(IUNOUT,'(29H SNF: NUMBER OF ONE-SPEED ITE,
+     1   8HRATIONS=,I5,1H.)') ITER
+         DEALLOCATE(OLD3,TEST,OLD2,OLD1,INCONV)
+      ENDIF
+      IF(IMPX.GT.2) THEN
+        CALL KDRCPU(TK2)
+        WRITE(IUNOUT,'(15H SNF: CPU TIME=,1P,E11.3,8H SECOND./)')
+     1  TK2-TK1
+      ENDIF
+      RETURN
+*
+  300 FORMAT(11H SNF: GROUP,I5,20H ONE-SPEED ITERATION,I4,8H  ERROR=,
+     1 1P,E11.4,5H OVER,E11.4,5H PREC,E12.4,22H  ACCELERATION FACTOR=,
+     2 0P,F7.3)
+  310 FORMAT (44H SNF: EIGENVALUE FOR FOURIER ANALYSIS, RHO= ,E13.6)
+      END