diff options
Diffstat (limited to 'Dragon/src/SNFT1C.f')
| -rw-r--r-- | Dragon/src/SNFT1C.f | 211 |
1 files changed, 211 insertions, 0 deletions
diff --git a/Dragon/src/SNFT1C.f b/Dragon/src/SNFT1C.f new file mode 100644 index 0000000..efaa775 --- /dev/null +++ b/Dragon/src/SNFT1C.f @@ -0,0 +1,211 @@ +*DECK SNFT1C + SUBROUTINE SNFT1C(NREG,NMAT,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, + 1 PLZ,PL,MAT,VOL,SURF,TOTAL,IGAV,QEXT,LFIXUP,CURR,FLUX) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Perform one inner iteration for solving SN equations in 1D cylindrical +* geometry. +* +*Copyright: +* Copyright (C) 2005 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 +* NREG number of regions. +* NMAT number of material mixtures. +* M2 number of axial $\\xi$ levels. +* NPQ number of SN directions in one octant. +* ISCAT anisotropy of one-speed sources: +* =1 isotropic sources; +* =2 linearly anisotropic sources. +* NSCT number of spherical harmonics components in the flux. +* JOP number of base points per axial level in one octant. +* U base points in $\\mu$ of the 1D quadrature. Used with +* zero-weight points. +* UPQ base points in $\\mu$ of the 2D SN quadrature. +* WPQ weights of the 2D SN quadrature. +* ALPHA angular redistribution terms. +* PLZ discrete values of the spherical harmonics corresponding +* to the 1D quadrature. Used with zero-weight points. +* MN moment-to-discrete matrix. +* DN discrete-to-moment matrix. +* MAT material mixture index in each region. +* VOL volumes of each region. +* SURF surfaces surrounding each region. +* TOTAL macroscopic total cross sections. +* IGAV type of condition at axial axis (=1 specular reflection; +* =2 zero-weight reflection; =3 averaged reflection). +* QEXT spherical harmonics components of the fixed source. +* LFIXUP flag to enable negative flux fixup. +* +*Parameters: input/output +* CURR entering current at input and leaving current at output. +* +*Parameters: output +* FLUX spherical harmonics components of the flux. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NREG,NMAT,M2,NPQ,ISCAT,NSCT,JOP(M2),MAT(NREG),IGAV + REAL U(M2),UPQ(NPQ),WPQ(NPQ),ALPHA(NPQ),PLZ(NSCT,M2),PL(NSCT,NPQ), + 1 VOL(NREG),SURF(NREG+1),TOTAL(0:NMAT),QEXT(NSCT,NREG),CURR, + 2 FLUX(NSCT,NREG) + LOGICAL LFIXUP +*---- +* LOCAL VARIABLES +*---- + PARAMETER(RLOG=1.0E-8,PI=3.141592654) + DOUBLE PRECISION Q,E2,AFB,Q1,Q2,PSIA,WSIA,CURSUM + REAL, ALLOCATABLE, DIMENSION(:) :: FLXB,AFGL +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(FLXB(M2),AFGL(NREG)) +*---- +* COMPUTE A NORMALIZATION CONSTANT. +*---- + DENOM=0.0 + DO 10 I=1,NPQ + IF(UPQ(I).GT.0.0) DENOM=DENOM+4.0*WPQ(I)*UPQ(I) + 10 CONTINUE +*---- +* OUTER LOOP OVER AXIAL LEVELS. +*---- + FLUX(:NSCT,:NREG)=0.0 + CURSUM=0.0D0 + IPQ=0 + DO 200 IP=1,M2 +*---- +* INITIALIZATION SWEEP (USING ZERO-WEIGHT POINTS). AFB IS THE EDGE +* FLUX VALUE AND AFGL(I) IS THE CENTERED FLUX VALUE. +*---- + AFB=CURR/DENOM + DO 30 I=NREG,1,-1 ! Spatial sweep + Q=0.0D0 + IBM=MAT(I) + IOF=0 + DO 25 IL=0,ISCAT-1 + DO 20 IM=0,IL + IF(MOD(IL+IM,2).EQ.1) GO TO 20 + IOF=IOF+1 + Q=Q+QEXT(IOF,I)*PLZ(IOF,IP)*(2.0*IL+1.0)/(4.0*PI) + 20 CONTINUE + 25 CONTINUE + Q1=-4.0D0*PI*SQRT(1.0-U(IP)*U(IP))/(SURF(I+1)-SURF(I)) + E2=(Q-Q1*AFB)/(TOTAL(IBM)-Q1) + IF(LFIXUP.AND.(E2.LE.RLOG)) E2=0.0 + AFB=2.0*E2-AFB + IF(LFIXUP.AND.(AFB.LE.RLOG)) AFB=0.0 + AFGL(I)=REAL(E2) + IF(LFIXUP.AND.(AFGL(I).LE.RLOG)) AFGL(I)=0.0 + 30 CONTINUE + WSIA=0.0D0 + IF(IGAV.EQ.2) THEN + PSIA=0.0D0 + ELSE + PSIA=AFB + ENDIF +*---- +* BACKWARD SWEEP (FROM EXTERNAL SURFACE TOWARD CENTRAL AXIS). +*---- + ALPMAX=0.0 + DO 80 M=2*JOP(IP),JOP(IP)+1,-1 ! Angular sweep + AFB=CURR/DENOM + ALPMIN=ALPHA(IPQ+M) + DO 70 I=NREG,1,-1 ! Spatial sweep + Q=0.0 + IBM=MAT(I) + IOF=0 + DO 45 IL=0,ISCAT-1 + DO 40 IM=0,IL + IF(MOD(IL+IM,2).EQ.1) GO TO 40 + IOF=IOF+1 + Q=Q+QEXT(IOF,I)*PL(IOF,IPQ+M)*(2.0*IL+1.0)/(4.0*PI) + 40 CONTINUE + 45 CONTINUE + Q1=Q*VOL(I)-UPQ(IPQ+M)*(SURF(I)+SURF(I+1))*AFB+(SURF(I+1)- + 1 SURF(I))*(ALPMIN+ALPMAX)*AFGL(I)/WPQ(IPQ+M) + Q2=TOTAL(IBM)*VOL(I)-2.0D0*UPQ(IPQ+M)*SURF(I)+2.0D0* + 1 (SURF(I+1)-SURF(I))*ALPMIN/WPQ(IPQ+M) + E2=Q1/Q2 + IF(LFIXUP.AND.(E2.LE.RLOG)) E2=0.0 + AFB=2.0*E2-AFB ! IN SPACE + IF(LFIXUP.AND.(AFB.LE.RLOG)) AFB=0.0 + AFGL(I)=2.0*REAL(E2)-AFGL(I) ! IN ANGLE + IF(LFIXUP.AND.(AFGL(I).LE.RLOG)) AFGL(I)=0.0 + DO 60 K=1,NSCT + FLUX(K,I)=FLUX(K,I)+4.0*WPQ(IPQ+M)*REAL(E2)*PL(K,IPQ+M) + 60 CONTINUE + 70 CONTINUE + IF(IGAV.EQ.1) THEN + FLXB(2*JOP(IP)-M+1)=REAL(AFB) + ELSE IF(IGAV.EQ.2) THEN + PSIA=PSIA+WPQ(IPQ+M)*REAL(AFB) + WSIA=WSIA+WPQ(IPQ+M) + ENDIF + ALPMAX=ALPMIN + 80 CONTINUE + IF(IGAV.EQ.2) PSIA=PSIA/WSIA +*---- +* FORWARD SWEEP (FROM CENTRAL AXIS TOWARD EXTERNAL SURFACE). +*---- + DO 130 M=JOP(IP),1,-1 ! Angular sweep + ALPMIN=ALPHA(IPQ+M) + IF(IGAV.EQ.1) THEN + AFB=FLXB(M) + ELSE IF(IGAV.EQ.2) THEN + AFB=PSIA + ELSE IF(IGAV.EQ.3) THEN + AFB=PSIA + ENDIF + DO 120 I=1,NREG ! Spatial sweep + Q=0.0 + IBM=MAT(I) + IOF=0 + DO 100 IL=0,ISCAT-1 + DO 90 IM=0,IL + IF(MOD(IL+IM,2).EQ.1) GO TO 90 + IOF=IOF+1 + Q=Q+QEXT(IOF,I)*PL(IOF,IPQ+M)*(2.0*IL+1.0)/(4.0*PI) + 90 CONTINUE + 100 CONTINUE + Q1=Q*VOL(I)+UPQ(IPQ+M)*(SURF(I)+SURF(I+1))*AFB+(SURF(I+1)- + 1 SURF(I))*(ALPMIN+ALPMAX)*AFGL(I)/WPQ(IPQ+M) + Q2=TOTAL(IBM)*VOL(I)+2.0D0*UPQ(IPQ+M)*SURF(I+1)+2.0D0* + 1 (SURF(I+1)-SURF(I))*ALPMIN/WPQ(IPQ+M) + E2=Q1/Q2 + IF(LFIXUP.AND.(E2.LE.RLOG)) E2=0.0 + AFB=2.0*E2-AFB ! IN SPACE + IF(LFIXUP.AND.(AFB.LE.RLOG)) AFB=0.0 + AFGL(I)=2.0*REAL(E2)-AFGL(I) ! IN ANGLE + IF(LFIXUP.AND.(AFGL(I).LE.RLOG)) AFGL(I)=0.0 + DO 110 K=1,NSCT + FLUX(K,I)=FLUX(K,I)+4.0*WPQ(IPQ+M)*REAL(E2)*PL(K,IPQ+M) + 110 CONTINUE + 120 CONTINUE + CURSUM=CURSUM+4.0*WPQ(IPQ+M)*UPQ(IPQ+M)*AFB + ALPMAX=ALPMIN + 130 CONTINUE +*---- +* END OF OUTER LOOP OVER AXIAL LEVELS +*---- + IPQ=IPQ+2*JOP(IP) + 200 CONTINUE + IF(IPQ.NE.NPQ) CALL XABORT('SN1T1C: BUG.') + CURR=REAL(CURSUM) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(AFGL,FLXB) + RETURN + END |