diff options
Diffstat (limited to 'Dragon/src/SNT1DS.f')
| -rw-r--r-- | Dragon/src/SNT1DS.f | 135 |
1 files changed, 135 insertions, 0 deletions
diff --git a/Dragon/src/SNT1DS.f b/Dragon/src/SNT1DS.f new file mode 100644 index 0000000..e446c78 --- /dev/null +++ b/Dragon/src/SNT1DS.f @@ -0,0 +1,135 @@ +*DECK SNT1DS + SUBROUTINE SNT1DS (IMPX,LX,NCODE,ZCODE,XXX,NLF,NSCT,U,W,ALPHA, + 1 PLZ,PL,VOL,IDL,SURF,IQUAD,IL,IM) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Numbering corresponding to a 1-D spherical 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). +* NSCT maximum number of spherical harmonics moments of the flux. +* IQUAD quadrature type: +* =1 Gauss-Lobatto; +* =2 Gauss-Legendre. +* +*Parameters: output +* U base points in $\\xi$ of the axial quadrature. +* W weights for the quadrature in $\\mu$. +* 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,NSCT,IDL(LX),IQUAD,IL(NSCT),IM(NSCT) + REAL ZCODE(2),XXX(LX+1),U(NLF),W(NLF),PLZ(NSCT),PL(NSCT,NLF), + 1 ALPHA(NLF),VOL(LX),SURF(LX+1) +*---- +* LOCAL VARIABLES +*---- + PARAMETER(PI=3.141592654) +*---- +* GENERATE QUADRATURE BASE POINTS AND CORRESPONDING WEIGHTS. +*---- + IF(MOD(NLF,2).EQ.1) CALL XABORT('SNT1DS: EVEN NLF EXPECTED.') + IF(IQUAD.EQ.1) THEN + CALL SNQU07(NLF,U,W) ! GAUSS-LOBATTO + ELSEIF(IQUAD.EQ.2) THEN + CALL ALGPT(NLF,-1.0,1.0,U,W) ! GAUSS-LEGENDRE + ELSE + CALL XABORT('SNT1DS: UNKNOWN QUADRATURE TYPE.') + ENDIF +*---- +* COMPUTE ALPHA. +*---- + SUMETA=0.0 + DO IP=1,NLF + SUMETA=SUMETA-2.0*W(IP)*U(IP) + ALPHA(IP)=SUMETA + ENDDO +*---- +* PRINT COSINES AND WEIGHTS. +*---- + IF(IMPX.GT.1) THEN + WRITE(6,60) (U(M),M=1,NLF) + 60 FORMAT(//,1X,'THE POSITIVE QUADRATURE COSINES FOLLOW'// + 1 (1X,5E14.6)) + WRITE(6,70) (W(M),M=1,NLF) + 70 FORMAT(//,1X,'THE CORRESPONDING QUADRATURE WEIGHTS FOLLOW'// + 1 (1X,5E14.6)) + WRITE(6,76) (ALPHA(M),M=1,NLF) + 76 FORMAT(//,1X,'THE ALPHAS FOLLOW'//(1X,5E14.6)) + ENDIF +*---- +* GENERATE LEGENDRE POLYNOMIALS FOR SCATTERING SOURCE. +*---- + IL(:NSCT)=0 + IM(:NSCT)=0 + DO 150 M=1,NLF + PL(1,M)=1.0 + IL(1)=0 + IF(NSCT.GT.1) THEN + PL(2,M)=U(M) + IL(2)=1 + DO 140 L=2,NSCT-1 + PL(L+1,M)=((2.0*L-1.0)*U(M)*PL(L,M)-(L-1)*PL(L-1,M))/REAL(L) + IL(L+1)=L + 140 CONTINUE + ENDIF + 150 CONTINUE + PLZ(1)=1.0 + IF(NSCT.GT.1) THEN + PLZ(2)=-1.0 + DO 160 L=2,NSCT-1 + PLZ(L+1)=(-(2.0*L-1.0)*PLZ(L)-(L-1)*PLZ(L-1))/REAL(L) + 160 CONTINUE + ENDIF +*---- +* COMPUTE SURFACES AND VOLUMES. +*---- + DO 200 I=1,LX+1 + SURF(I)=4.0*PI*XXX(I)*XXX(I) + 200 CONTINUE + DO 210 I=1,LX + IDL(I)=(I-1)*NSCT+1 + VOL(I)=4.0*PI*(XXX(I+1)**3-XXX(I)**3)/3.0 + 210 CONTINUE +*---- +* SET BOUNDARY CONDITIONS. +*---- + IF(NCODE(2).EQ.2) ZCODE(2)=1.0 +* + RETURN + END |