diff options
Diffstat (limited to 'Dragon/src/SNF.f')
| -rw-r--r-- | Dragon/src/SNF.f | 257 |
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 |