Initial commit from Polytechnique Montreal

1 files changed, 343 insertions, 0 deletions

diff --git a/Dragon/src/SYBJJ2.f b/Dragon/src/SYBJJ2.f
new file mode 100644
index 0000000..76dffa4
--- /dev/null
+++ b/Dragon/src/SYBJJ2.f
@@ -0,0 +1,343 @@
+*DECK SYBJJ2
+      SUBROUTINE SYBJJ2 (IPAS,NMCEL,NMERGE,NGEN,IJAT,NPIJ,NPIS,EPSJ,
+     1 NUNKNO,FUNKNO,SUNKNO,IMPX,NCOUR,NMC,IFR,ALB,INUM,MIX,DVX,IGEN,
+     2 PIJW,PISW,PSJW,PSSW)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the neutron flux and interface currents in a 2-D Cartesian
+* or hexagonal assembly using the current iteration method with
+* Roth X 4, DP0 or DP1 approximation.
+*
+*Copyright:
+* Copyright (C) 2002 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
+* IPAS    total number of regions.
+* NMCEL   total number of cells in the domain.
+* NMERGE  total number of merged cells for which specific values
+*         of the neutron flux and reactions rates are required.
+*         Many cells with different position in the domain can
+*         be merged before the neutron flux calculation if they
+*         own the same generating cell (NMERGE.le.NMCEL).
+* NGEN    total number of generating cells. A generating cell is
+*         defined by its material and dimensions, irrespective of
+*         its position in the domain (NGEN.le.NMERGE).
+* IJAT    total number of distinct out-currents.
+* NPIJ    size of cellwise scattering-reduced collision probability
+*         matrices.
+* NPIS    size of cellwise scattering-reduced escape probability
+*         matrices.
+* EPSJ    stopping criterion for flux-current iterations.
+* NUNKNO  total number of unknowns in vectors SUNKNO and FUNKNO.
+* SUNKNO  input source vector.
+* IMPX    print flag (equal to 0 for no print).
+* NCOUR   number of incoming currents (=4 Cartesian lattice;
+*         =6 hexagonal lattice).
+* NMC     offset of the first volume in each generating cell.
+* IFR     index-number of in-currents.
+* ALB     transmission/albedo associated with each in-current.
+* INUM    index-number of the merged cell associated to each cell.
+* MIX     index-number of out-currents.
+* DVX     weight associated with each out-current.
+*         Note: IFR, ALB, MIX and DVX contains information to rebuild
+*         the geometrical 'A' matrix.
+* IGEN    index-number of the generating cell associated with each
+*         merged cell.
+* PIJW    cellwise scattering-reduced collision probability matrices.
+* PISW    cellwise scattering-reduced escape probability matrices.
+* PSJW    cellwise scattering-reduced collision probability matrices
+*         for incoming neutrons.
+* PSSW    cellwise scattering-reduced transmission probability
+*         matrices.
+* 
+*Parameters: input/output
+* FUNKNO  unknown vector.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IPAS,NMCEL,NMERGE,NGEN,IJAT,NPIJ,NPIS,NUNKNO,IMPX,NCOUR,
+     1 NMC(NGEN+1),IFR(NCOUR*NMCEL),INUM(NMCEL),MIX(NCOUR*NMERGE),
+     2 IGEN(NMERGE)
+      REAL EPSJ,FUNKNO(NUNKNO),SUNKNO(NUNKNO),ALB(NCOUR*NMCEL),
+     1 DVX(NCOUR*NMERGE),PIJW(NPIJ),PISW(NCOUR*NPIS),PSJW(NCOUR*NPIS),
+     2 PSSW(NGEN*NCOUR*NCOUR)
+*----
+*  LOCAL VARIABLES
+*----
+      REAL PIJ,PIS
+      LOGICAL LOGTES
+      PARAMETER (MAXIT=400,LACCFC=2,ICL1=3,ICL2=3)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, DIMENSION(:), POINTER :: INDPIJ,INDNMC
+      DOUBLE PRECISION, DIMENSION(:), POINTER :: CIT0
+      DOUBLE PRECISION, DIMENSION(:,:), POINTER :: CITR,AITR
+      DOUBLE PRECISION, DIMENSION(:), POINTER :: WCURR
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(INDPIJ(NGEN),INDNMC(NMERGE))
+      ALLOCATE(CITR(3,IJAT),CIT0(IJAT),AITR(2,IJAT))
+      ALLOCATE(WCURR(IJAT))
+*
+      IPIJ=0
+      DO 10 JKG=1,NGEN
+      J2=NMC(JKG+1)-NMC(JKG)
+      INDPIJ(JKG)=IPIJ
+      IPIJ=IPIJ+J2*J2
+   10 CONTINUE
+      KNMC=0
+      DO 20 JKK=1,NMERGE
+      JKG=IGEN(JKK)
+      J2=NMC(JKG+1)-NMC(JKG)
+      INDNMC(JKK)=KNMC
+      KNMC=KNMC+J2
+   20 CONTINUE
+*
+      DO 30 I=1,IJAT
+      WCURR(I)=1.0D0
+      CIT0(I)=0.0D0
+      CITR(1,I)=FUNKNO(IPAS+I)
+   30 CONTINUE
+*----
+*  COMPUTE PSJW * Q(*) CONTRIBUTION
+*----
+      DO 42 IKK=1,NMERGE
+      IKG=IGEN(IKK)
+      I1P=NMC(IKG)
+      I2=NMC(IKG+1)-I1P
+      IT=NCOUR*(IKK-1)
+      KNMC=INDNMC(IKK)
+      DO 41 I=1,I2
+      DO 40 IC=1,NCOUR
+      JCC=MIX(IT+IC)
+      PBJ=PSJW(I1P*NCOUR+(I-1)*NCOUR+IC)
+      CIT0(JCC)=CIT0(JCC)+PBJ*DVX(IT+IC)*SUNKNO(KNMC+I)
+   40 CONTINUE
+   41 CONTINUE
+   42 CONTINUE
+*----
+*  COMPUTE NORMALIZATION VECTOR WCURR
+*----
+      DO 52 ICEL=1,NMCEL
+      IKK=INUM(ICEL)
+      IT=NCOUR*(IKK-1)
+      IS=NCOUR*(ICEL-1)
+      IKG=IGEN(IKK)
+      IPSS=(IKG-1)*NCOUR*NCOUR
+      DO 51 JC=1,NCOUR
+      J1=IFR(IS+JC)
+      DO 50 IC=1,NCOUR
+      PSS=PSSW(IPSS+(JC-1)*NCOUR+IC)
+      WCURR(J1)=WCURR(J1)-PSS*ALB(IS+JC)*DVX(IT+IC)
+   50 CONTINUE
+   51 CONTINUE
+   52 CONTINUE
+*
+      ISTART=1
+      TEST=0.0D0
+      ITER=0
+   70 ITER=ITER+1
+      IF(ITER.GT.MAXIT) THEN
+         WRITE(6,'(/47H SYBJJ2: *** WARNING *** MAXIMUM NUMBER OF ITER,
+     1   15HATIONS REACHED.)')
+         GO TO 190
+      ENDIF
+      IT3=MOD(ITER,3)+1
+      IT2=MOD(ITER-1,3)+1
+      IT1=MOD(ITER-2,3)+1
+      DO 80 I=1,IJAT
+      CITR(IT3,I)=CIT0(I)
+   80 CONTINUE
+*----
+*  COMPUTE PSSW * J(-) CONTRIBUTION
+*----
+      DO 92 ICEL=1,NMCEL
+      IKK=INUM(ICEL)
+      IT=NCOUR*(IKK-1)
+      IS=NCOUR*(ICEL-1)
+      IKG=IGEN(IKK)
+      IPSS=(IKG-1)*NCOUR*NCOUR
+      DO 91 JC=1,NCOUR
+      J1=IFR(IS+JC)
+      DO 90 IC=1,NCOUR
+      J2=MIX(IT+IC)
+      PSS=PSSW(IPSS+(JC-1)*NCOUR+IC)
+      CITR(IT3,J2)=CITR(IT3,J2)+PSS*ALB(IS+JC)*DVX(IT+IC)*CITR(IT2,J1)
+   90 CONTINUE
+   91 CONTINUE
+   92 CONTINUE
+*----
+*  NORMALIZATION
+*----
+      S1=0.0D0
+      S2=0.0D0
+      DO 100 I=1,IJAT
+      S1=S1+WCURR(I)*CITR(IT3,I)
+      S2=S2+CIT0(I)
+  100 CONTINUE
+      ZNORM=S2/S1
+      IF(ZNORM.LT.0.0D0) ZNORM=1.0D0
+      DO 110 I=1,IJAT
+      CITR(IT3,I)=CITR(IT3,I)*ZNORM
+  110 CONTINUE
+*----
+*  ONE/TWO PARAMETER ACCELERATION
+*----
+      ALP=1.0D0
+      BET=0.0D0
+      LOGTES=(1+MOD(ITER-ISTART,ICL1+ICL2).GT.ICL1)
+      IF(LOGTES) THEN
+         DO 120 I=1,IJAT
+         AITR(1,I)=CITR(IT3,I)-CITR(IT2,I)
+         AITR(2,I)=CITR(IT2,I)-CITR(IT1,I)
+  120    CONTINUE
+         DO 132 ICEL=1,NMCEL
+         IKK=INUM(ICEL)
+         IT=NCOUR*(IKK-1)
+         IS=NCOUR*(ICEL-1)
+         IKG=IGEN(IKK)
+         IPSS=(IKG-1)*NCOUR*NCOUR
+         DO 131 JC=1,NCOUR
+         J1=IFR(IS+JC)
+         DO 130 IC=1,NCOUR
+         J2=MIX(IT+IC)
+         PSS=PSSW(IPSS+(JC-1)*NCOUR+IC)*ALB(IS+JC)*DVX(IT+IC)
+         AITR(1,J2)=AITR(1,J2)-PSS*(CITR(IT3,J1)-CITR(IT2,J1))
+         AITR(2,J2)=AITR(2,J2)-PSS*(CITR(IT2,J1)-CITR(IT1,J1))
+  130    CONTINUE
+  131    CONTINUE
+  132    CONTINUE
+         IF((LACCFC.EQ.1).OR.(MOD(ITER-ISTART,ICL1+ICL2).EQ.ICL1)) THEN
+            S1=0.0D0
+            S2=0.0D0
+            DO 140 I=1,IJAT
+            S1=S1+(CITR(IT3,I)-CITR(IT2,I))*AITR(1,I)
+            S2=S2+AITR(1,I)*AITR(1,I)
+  140       CONTINUE
+            IF(S2.EQ.0.0D0) THEN
+               ISTART=ITER+1
+            ELSE
+               ALP=S1/S2
+               IF(ALP.LE.0.0D0) THEN
+                  ISTART=ITER+1
+                  ALP=1.0D0
+               ENDIF
+            ENDIF
+            DO 150 I=1,IJAT
+            CITR(IT3,I)=CITR(IT2,I)+ALP*(CITR(IT3,I)-CITR(IT2,I))
+  150       CONTINUE
+         ELSE IF(LACCFC.EQ.2) THEN
+            S1=0.0D0
+            S2=0.0D0
+            S3=0.0D0
+            S4=0.0D0
+            S5=0.0D0
+            DO 160 I=1,IJAT
+            S1=S1+(CITR(IT3,I)-CITR(IT2,I))*AITR(1,I)
+            S2=S2+AITR(1,I)*AITR(1,I)
+            S3=S3+(CITR(IT3,I)-CITR(IT2,I))*AITR(2,I)
+            S4=S4+AITR(1,I)*AITR(2,I)
+            S5=S5+AITR(2,I)*AITR(2,I)
+  160       CONTINUE
+            DET=S2*S5-S4*S4
+            IF(DET.EQ.0.0D0) THEN
+               ISTART=ITER+1
+            ELSE
+               ALP=(S5*S1-S4*S3)/DET
+               BET=(S2*S3-S4*S1)/DET
+               IF(ALP.LE.0.0D0) THEN
+                  ISTART=ITER+1
+                  ALP=1.0D0
+                  BET=0.0D0
+               ENDIF
+            ENDIF
+            DO 170 I=1,IJAT
+            CITR(IT3,I)=CITR(IT2,I)+ALP*(CITR(IT3,I)-CITR(IT2,I))+
+     1      BET*(CITR(IT2,I)-CITR(IT1,I))
+  170       CONTINUE
+         ENDIF
+      ENDIF
+*----
+*  CHECK THE CONVERGENCE ERROR
+*----
+      ERR1=0.0D0
+      ERR2=0.0D0
+      DO 180 I=1,IJAT
+      ERR1=MAX(ERR1,ABS(CITR(IT3,I)-CITR(IT2,I)))
+      ERR2=MAX(ERR2,ABS(CITR(IT3,I)))
+  180 CONTINUE
+      IF(IMPX.GT.3) WRITE(6,'(30H SYBJJ2: CURRENT ITERATION NB.,I4,
+     1 7H ERROR=,1P,E10.3,5H OVER,E10.3,15H NORMALIZATION=,E10.3,
+     2 14H ACCELERATION=,2E11.3,1H.)') ITER,ERR1,ERR2,ZNORM,ALP,
+     3 BET/ALP
+      IF(ITER.EQ.1) TEST=ERR1/ERR2
+      IF((ITER.GT.20).AND.(ERR1/ERR2.GT.TEST)) THEN
+         WRITE(6,'(/45H SYBJJ2: *** WARNING *** CONVERGENCE DIFFICUL,
+     1   5HTIES.)')
+         GO TO 190
+      ENDIF
+      IF(LOGTES.OR.(ERR1.GT.EPSJ*ERR2)) GO TO 70
+      IF(IMPX.GT.2) WRITE(6,'(37H SYBJJ2: CURRENT CONVERGENCE AT ITERA,
+     1 8HTION NB.,I4,7H ERROR=,1P,E10.3,5H OVER,E10.3,1H.)') ITER,ERR1,
+     2 ERR2
+*
+  190 DO 200 I=1,IPAS
+      FUNKNO(I)=0.0
+  200 CONTINUE
+      DO 210 I=1,IJAT
+      FUNKNO(IPAS+I)=REAL(CITR(IT3,I))
+  210 CONTINUE
+*----
+*  COMPUTE PISW * J(-) CONTRIBUTION
+*----
+      DO 240 ICEL=1,NMCEL
+      IKK=INUM(ICEL)
+      IS=NCOUR*(ICEL-1)
+      IKG=IGEN(IKK)
+      I1P=NMC(IKG)
+      I2=NMC(IKG+1)-I1P
+      KNMC=INDNMC(IKK)
+      DO 230 J=1,I2
+      DO 220 JC=1,NCOUR
+      J1=IFR(IS+JC)
+      PIS=PISW(I1P*NCOUR+(JC-1)*I2+J)
+      FUNKNO(KNMC+J)=FUNKNO(KNMC+J)+PIS*ALB(IS+JC)*FUNKNO(IPAS+J1)
+  220 CONTINUE
+  230 CONTINUE
+  240 CONTINUE
+*----
+*  COMPUTE PIJW * Q(*) CONTRIBUTION
+*----
+      DO 270 IKK=1,NMERGE
+      IKG=IGEN(IKK)
+      I2=NMC(IKG+1)-NMC(IKG)
+      KNMC=INDNMC(IKK)
+      DO 260 I=1,I2
+      DO 250 J=1,I2
+      PIJ=PIJW(INDPIJ(IKG)+(I-1)*I2+J)
+      FUNKNO(KNMC+J)=FUNKNO(KNMC+J)+PIJ*SUNKNO(KNMC+I)
+  250 CONTINUE
+  260 CONTINUE
+  270 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(WCURR)
+      DEALLOCATE(AITR,CIT0,CITR)
+      DEALLOCATE(INDNMC,INDPIJ)
+      RETURN
+      END