diff options
Diffstat (limited to 'Dragon/src/SNT1DC.f')
| -rw-r--r-- | Dragon/src/SNT1DC.f | 184 |
1 files changed, 184 insertions, 0 deletions
diff --git a/Dragon/src/SNT1DC.f b/Dragon/src/SNT1DC.f new file mode 100644 index 0000000..702c067 --- /dev/null +++ b/Dragon/src/SNT1DC.f @@ -0,0 +1,184 @@ +*DECK SNT1DC + SUBROUTINE SNT1DC (IMPX,LX,NCODE,ZCODE,XXX,NLF,NPQ,NSCT,IQUAD, + 1 JOP,U,W,TPQ,UPQ,VPQ,WPQ,ALPHA,PLZ,PL,VOL,IDL,SURF,IL,IM) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Numbering corresponding to a 1-D cylindrical geometry with discrete +* ordinates approximation of the flux. +* +*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 +* IMPX print parameter. +* LX number of elements along the R axis. +* NCODE type of boundary condition applied on each side +* (i=1 R-; i=2 R+): +* =1 VOID; =2 REFL; =7 ZERO. +* ZCODE ZCODE(I) is the albedo corresponding to boundary condition +* 'VOID' on each side (ZCODE(I)=0.0 by default). +* XXX coordinates along the R axis. +* NLF order of the SN approximation (even number). +* NPQ number of SN directions in two octants. +* NSCT maximum number of spherical harmonics moments of the flux. +* IQUAD type of SN quadrature (1 Level symmetric, type IQUAD; +* 4 Gauss-Legendre and Gauss-Chebyshev; 10 product). +* +*Parameters: output +* JOP number of base points per axial level in one octant. +* U base points in $\\xi$ of the axial quadrature. Used with +* zero-weight points. +* W weights for the quadrature in $\\mu$. +* TPQ base points in $\\xi$ of the 2D SN quadrature. +* UPQ base points in $\\mu$ of the 2D SN quadrature. +* VPQ base points in $\\eta$ 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. +* PL discrete values of the spherical harmonics corresponding +* to the 2D SN quadrature. +* VOL volume of each element. +* IDL position of integrated fluxes into unknown vector. +* SURF surfaces. +* IL indexes (l) of each spherical harmonics in the +* interpolation basis. +* IM indexes (m) of each spherical harmonics in the +* interpolation basis. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IMPX,LX,NCODE(2),NLF,NPQ,NSCT,IQUAD,JOP(NLF/2),IDL(LX), + 1 IL(NSCT),IM(NSCT) + REAL ZCODE(2),XXX(LX+1),U(NLF/2),W(NLF/2),PLZ(NSCT,NLF/2), + 1 PL(NSCT,NPQ),TPQ(NPQ),UPQ(NPQ),VPQ(NPQ),WPQ(NPQ),ALPHA(NPQ), + 2 VOL(LX),SURF(LX+1) +*---- +* LOCAL VARIABLES +*---- + PARAMETER(PI=3.141592654) + REAL, ALLOCATABLE, DIMENSION(:) :: TPQ2,UPQ2,VPQ2,WPQ2 +*---- +* GENERATE QUADRATURE BASE POINTS AND CORRESPONDING WEIGHTS. +*---- + IF(MOD(NLF,2).EQ.1) CALL XABORT('SNT1DC: EVEN NLF EXPECTED.') + M2=NLF/2 + ALLOCATE(TPQ2(NPQ/2),UPQ2(NPQ/2),VPQ2(NPQ/2),WPQ2(NPQ/2)) + IF(IQUAD.EQ.1) THEN + CALL SNQU01(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2) + ELSE IF(IQUAD.EQ.2) THEN + CALL SNQU02(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2) + ELSE IF(IQUAD.EQ.3) THEN + CALL SNQU03(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2) + ELSE IF(IQUAD.EQ.4) THEN + CALL SNQU04(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2) + ELSE IF(IQUAD.EQ.10) THEN + CALL SNQU10(NLF,JOP,U,W,TPQ2,UPQ2,VPQ2,WPQ2) + ELSE + CALL XABORT('SNT1DC: UNKNOWN QUADRATURE TYPE.') + ENDIF + IPQ=0 + IPR=0 + WSUM=0.0 + DO IP=1,M2 + DO IQ=1,JOP(IP) + IPQ=IPQ+1 + IPR=IPR+1 + TPQ(IPQ)=TPQ2(IPR) + UPQ(IPQ)=UPQ2(IPR) + VPQ(IPQ)=VPQ2(IPR) + WPQ(IPQ)=WPQ2(IPR) + WSUM=WSUM+WPQ(IPQ) + ENDDO + DO IQ=JOP(IP)+1,2*JOP(IP) + IPQ=IPQ+1 + JQ=IQ-(JOP(IP)+1)+1 + JPQ=IPQ-2*JQ+1 + TPQ(IPQ)=TPQ(JPQ) + UPQ(IPQ)=-UPQ(JPQ) + VPQ(IPQ)=VPQ(JPQ) + WPQ(IPQ)=WPQ(JPQ) + WSUM=WSUM+WPQ(IPQ) + ENDDO + ENDDO + DEALLOCATE(WPQ2,VPQ2,UPQ2,TPQ2) + IF(IPQ.NE.NPQ) CALL XABORT('SNT1DC: BAD VALUE ON NPQ.') +*---- +* COMPUTE ALPHA. +*---- + IPQ=0 + DO IP=1,M2 + SUMETA=0.0 + DO IQ=1,2*JOP(IP) + IPQ=IPQ+1 + ALPHA(IPQ)=SUMETA + SUMETA=SUMETA+WPQ(IPQ)*UPQ(IPQ) + ENDDO + ENDDO +*---- +* PRINT COSINES AND WEIGHTS. +*---- + IF(IMPX.GT.1) THEN + WRITE(6,'(/20H SNT1DC: WEIGHT SUM=,1P,E11.4)') WSUM + WRITE(6,60) (U(N),N=1,M2) + 60 FORMAT(//,1X,'THE POSITIVE QUADRATURE COSINES FOLLOW'// + 1 (1X,5E14.6)) + WRITE(6,70) (W(N),N=1,M2) + 70 FORMAT(//,1X,'THE CORRESPONDING QUADRATURE WEIGHTS FOLLOW'// + 1 (1X,5E14.6)) + WRITE(6,74) (UPQ(N),N=1,NPQ) + 74 FORMAT(//,1X,'THE BASE POINTS (MU) FOLLOW'//(1X,5E14.6)) + WRITE(6,75) (WPQ(N),N=1,NPQ) + 75 FORMAT(//,1X,'THE WEIGHTS FOLLOW'//(1X,5E14.6)) + WRITE(6,76) (ALPHA(N),N=1,NPQ) + 76 FORMAT(//,1X,'THE ALPHAS FOLLOW'//(1X,5E14.6)) + ENDIF +*---- +* GENERATE SPHERICAL HARMONICS FOR SCATTERING SOURCE. +*---- + IOF=0 + IL(:NSCT)=0 + IM(:NSCT)=0 + DO 130 L=0,NLF-1 + DO 120 M=0,L + IF(MOD(L+M,2).EQ.1) GO TO 120 + IOF=IOF+1 + IL(IOF)=L + IM(IOF)=M + IF(IOF.GT.NSCT) GO TO 140 + DO 100 N=1,M2 + PLZ(IOF,N)=PNSH(L,M,U(N),-SQRT(1.0-U(N)*U(N)),0.0) + 100 CONTINUE + DO 110 N=1,NPQ + PL(IOF,N)=PNSH(L,M,TPQ(N),UPQ(N),VPQ(N)) + 110 CONTINUE + 120 CONTINUE + 130 CONTINUE +*---- +* COMPUTE SURFACES AND VOLUMES. +*---- + 140 DO 150 I=1,LX+1 + SURF(I)=2.0*PI*XXX(I) + 150 CONTINUE + DO 160 I=1,LX + IDL(I)=(I-1)*NSCT+1 + VOL(I)=PI*(XXX(I+1)*XXX(I+1)-XXX(I)*XXX(I)) + 160 CONTINUE +*---- +* SET BOUNDARY CONDITIONS. +*---- + IF(NCODE(2).EQ.2) ZCODE(2)=1.0 +* + RETURN + END |