summaryrefslogtreecommitdiff
path: root/Dragon/src/SYBJJ1.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SYBJJ1.f')
-rw-r--r--Dragon/src/SYBJJ1.f370
1 files changed, 370 insertions, 0 deletions
diff --git a/Dragon/src/SYBJJ1.f b/Dragon/src/SYBJJ1.f
new file mode 100644
index 0000000..a342976
--- /dev/null
+++ b/Dragon/src/SYBJJ1.f
@@ -0,0 +1,370 @@
+*DECK SYBJJ1
+ SUBROUTINE SYBJJ1 (IPAS,NMCEL,NMERGE,NGEN,NPIJ,NPIS,EPSJ,NUNKNO,
+ 1 FUNKNO,SUNKNO,IMPX,NCOUR,XX,YY,NMC,IFR,ALB,INUM,IGEN,PIJW,PISW,
+ 2 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
+* 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. Equal to the total number
+* of distinct out-currents (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).
+* 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).
+* XX X-thickness of the generating cells.
+* YY Y-thickness of the generating cells.
+* NMC offset of the first volume in each generating cell.
+* IFR index-number of in-currents.
+* ALB transmission/albedo associated with each in-current.
+* Note: IFR and ALB contains information to rebuild the
+* geometrical 'A' matrix.
+* INUM index-number of the merged cell associated to each cell.
+* 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,NPIJ,NPIS,NUNKNO,IMPX,NCOUR,
+ 1 NMC(NGEN+1),IFR(NCOUR*NMCEL),INUM(NMCEL),IGEN(NMERGE)
+ REAL EPSJ,FUNKNO(NUNKNO),SUNKNO(NUNKNO),XX(NGEN),YY(NGEN),
+ 1 ALB(NCOUR*NMCEL),PIJW(NPIJ),PISW(NPIS),PSJW(NPIS),PSSW(NGEN)
+*----
+* LOCAL VARIABLES
+*----
+ REAL PIJ,PIS
+ DOUBLE PRECISION PIBB(6)
+ 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,NMERGE),CIT0(NMERGE),AITR(2,NMERGE))
+ ALLOCATE(WCURR(NMERGE))
+*
+ 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,NMERGE
+ WCURR(I)=1.0D0
+ CIT0(I)=0.0D0
+ CITR(1,I)=FUNKNO(IPAS+I)
+ 30 CONTINUE
+*----
+* COMPUTE PSJW * Q(*) CONTRIBUTION
+*----
+ DO 45 IKK=1,NMERGE
+ IKG=IGEN(IKK)
+ I1P=NMC(IKG)
+ I2=NMC(IKG+1)-I1P
+ KNMC=INDNMC(IKK)
+ DO 40 I=1,I2
+ CIT0(IKK)=CIT0(IKK)+PSJW(I1P+I)*SUNKNO(KNMC+I)
+ 40 CONTINUE
+ 45 CONTINUE
+*----
+* COMPUTE NORMALIZATION VECTOR WCURR
+*----
+ DO 65 ICEL=1,NMCEL
+ IKK=INUM(ICEL)
+ IS=NCOUR*(ICEL-1)
+ IKG=IGEN(IKK)
+ IF(NCOUR.EQ.4) THEN
+ A=XX(IKG)
+ B=YY(IKG)
+ DEN1=2.0D0*(A+B)
+ PIBB(1)=B/DEN1
+ PIBB(2)=B/DEN1
+ PIBB(3)=A/DEN1
+ PIBB(4)=A/DEN1
+ ELSE
+ DO 50 JC=1,NCOUR
+ PIBB(JC)=1.0D0/6.0D0
+ 50 CONTINUE
+ ENDIF
+ DO 60 JC=1,NCOUR
+ J1=IFR(IS+JC)
+ WCURR(J1)=WCURR(J1)-PSSW(IKG)*PIBB(JC)*ALB(IS+JC)
+ 60 CONTINUE
+ 65 CONTINUE
+*
+ ISTART=1
+ TEST=0.0D0
+ ITER=0
+ 70 ITER=ITER+1
+ IF(ITER.GT.MAXIT) THEN
+ WRITE(6,'(/47H SYBJJ1: *** 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,NMERGE
+ CITR(IT3,I)=CIT0(I)
+ 80 CONTINUE
+*----
+* COMPUTE PSSW * J(-) CONTRIBUTION
+*----
+ DO 95 ICEL=1,NMCEL
+ IKK=INUM(ICEL)
+ IS=NCOUR*(ICEL-1)
+ IKG=IGEN(IKK)
+ IF(NCOUR.EQ.4) THEN
+ A=XX(IKG)
+ B=YY(IKG)
+ DEN1=2.0D0*(A+B)
+ PIBB(1)=B/DEN1
+ PIBB(2)=B/DEN1
+ PIBB(3)=A/DEN1
+ PIBB(4)=A/DEN1
+ ELSE
+ DO 85 JC=1,NCOUR
+ PIBB(JC)=1.0D0/6.0D0
+ 85 CONTINUE
+ ENDIF
+ DO 90 JC=1,NCOUR
+ J1=IFR(IS+JC)
+ PSS=PSSW(IKG)*PIBB(JC)
+ CITR(IT3,IKK)=CITR(IT3,IKK)+PSS*ALB(IS+JC)*CITR(IT2,J1)
+ 90 CONTINUE
+ 95 CONTINUE
+*----
+* NORMALIZATION
+*----
+ S1=0.0D0
+ S2=0.0D0
+ DO 100 I=1,NMERGE
+ 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,NMERGE
+ 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,NMERGE
+ AITR(1,I)=CITR(IT3,I)-CITR(IT2,I)
+ AITR(2,I)=CITR(IT2,I)-CITR(IT1,I)
+ 120 CONTINUE
+ DO 135 ICEL=1,NMCEL
+ IKK=INUM(ICEL)
+ IS=NCOUR*(ICEL-1)
+ IKG=IGEN(IKK)
+ IF(NCOUR.EQ.4) THEN
+ A=XX(IKG)
+ B=YY(IKG)
+ DEN1=2.0D0*(A+B)
+ PIBB(1)=B/DEN1
+ PIBB(2)=B/DEN1
+ PIBB(3)=A/DEN1
+ PIBB(4)=A/DEN1
+ ELSE
+ DO 125 JC=1,NCOUR
+ PIBB(JC)=1.0D0/6.0D0
+ 125 CONTINUE
+ ENDIF
+ DO 130 JC=1,NCOUR
+ J1=IFR(IS+JC)
+ PSS=PSSW(IKG)*PIBB(JC)*ALB(IS+JC)
+ AITR(1,IKK)=AITR(1,IKK)-PSS*(CITR(IT3,J1)-CITR(IT2,J1))
+ AITR(2,IKK)=AITR(2,IKK)-PSS*(CITR(IT2,J1)-CITR(IT1,J1))
+ 130 CONTINUE
+ 135 CONTINUE
+ IF((LACCFC.EQ.1).OR.(MOD(ITER-ISTART,ICL1+ICL2).EQ.ICL1)) THEN
+ S1=0.0D0
+ S2=0.0D0
+ DO 140 I=1,NMERGE
+ 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,NMERGE
+ 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,NMERGE
+ 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,NMERGE
+ 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,NMERGE
+ 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 SYBJJ1: 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)) CALL XABORT('SYBJJ1: '
+ 1 //'CONVERGENCE FAILURE.')
+ IF(LOGTES.OR.(ERR1.GT.EPSJ*ERR2)) GO TO 70
+ IF(IMPX.GT.2) WRITE(6,'(37H SYBJJ1: 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,NMERGE
+ FUNKNO(IPAS+I)=REAL(CITR(IT3,I))
+ 210 CONTINUE
+*----
+* COMPUTE PISW * J(-) CONTRIBUTION
+*----
+ DO 250 ICEL=1,NMCEL
+ IKK=INUM(ICEL)
+ IS=NCOUR*(ICEL-1)
+ IKG=IGEN(IKK)
+ IF(NCOUR.EQ.4) THEN
+ A=XX(IKG)
+ B=YY(IKG)
+ DEN1=2.0D0*(A+B)
+ PIBB(1)=B/DEN1
+ PIBB(2)=B/DEN1
+ PIBB(3)=A/DEN1
+ PIBB(4)=A/DEN1
+ ELSE
+ DO 220 JC=1,NCOUR
+ PIBB(JC)=1.0D0/6.0D0
+ 220 CONTINUE
+ ENDIF
+ I1P=NMC(IKG)
+ I2=NMC(IKG+1)-I1P
+ KNMC=INDNMC(IKK)
+ DO 240 J=1,I2
+ DO 230 JC=1,NCOUR
+ J1=IFR(IS+JC)
+ PIS=PISW(I1P+J)*REAL(PIBB(JC))
+ FUNKNO(KNMC+J)=FUNKNO(KNMC+J)+PIS*ALB(IS+JC)*FUNKNO(IPAS+J1)
+ 230 CONTINUE
+ 240 CONTINUE
+ 250 CONTINUE
+*----
+* COMPUTE PIJW * Q(*) CONTRIBUTION
+*----
+ DO 280 IKK=1,NMERGE
+ IKG=IGEN(IKK)
+ I2=NMC(IKG+1)-NMC(IKG)
+ KNMC=INDNMC(IKK)
+ DO 270 I=1,I2
+ DO 260 J=1,I2
+ PIJ=PIJW(INDPIJ(IKG)+(I-1)*I2+J)
+ FUNKNO(KNMC+J)=FUNKNO(KNMC+J)+PIJ*SUNKNO(KNMC+I)
+ 260 CONTINUE
+ 270 CONTINUE
+ 280 CONTINUE
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(WCURR)
+ DEALLOCATE(AITR,CIT0,CITR)
+ DEALLOCATE(INDNMC,INDPIJ)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:30:55 +0000