diff options
Diffstat (limited to 'Dragon/src/SNQU01.f')
| -rw-r--r-- | Dragon/src/SNQU01.f | 143 |
1 files changed, 143 insertions, 0 deletions
diff --git a/Dragon/src/SNQU01.f b/Dragon/src/SNQU01.f new file mode 100644 index 0000000..0959d2f --- /dev/null +++ b/Dragon/src/SNQU01.f @@ -0,0 +1,143 @@ +*DECK SNQU01 + SUBROUTINE SNQU01(NLF,JOP,U,W,TPQ,UPQ,VPQ,WPQ) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Set the level-symmetric (type 1) quadratures. +* +*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 +* NLF order of the SN approximation (even number). +* +*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 axial quadrature in $\\xi$. +* 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. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NLF,JOP(NLF/2) + REAL U(NLF/2),W(NLF/2),TPQ(NLF*(NLF/2+1)/4),UPQ(NLF*(NLF/2+1)/4), + 1 VPQ(NLF*(NLF/2+1)/4),WPQ(NLF*(NLF/2+1)/4) +*---- +* LOCAL VARIABLES +*---- + PARAMETER(PI=3.141592654,MAXNLF=20,MAXNBA=55,MAXW=12) + INTEGER INWEI(MAXNBA) + DOUBLE PRECISION ZMU1,ZMU2,WSUM2,WEI(MAXW) +*---- +* SET THE UNIQUE QUADRATURE VALUES. +*---- + IF(NLF.GT.MAXNLF) CALL XABORT('SNQU01: MAXNLF OVERFLOW.') + M2=NLF/2 + NPQ=M2*(M2+1)/2 + NW=1+(NLF*(NLF+8)-1)/48 + IF(NW.GT.MAXW) CALL XABORT('SNQU01: MAXW OVERFLOW.') + ZMU1=0.0D0 + WEI(:NW)=0.0D0 + IF(NLF.EQ.2) THEN + ZMU1=1.0D0/3.0D0 + WEI(1)=1.0D0 + ELSE IF(NLF.EQ.4) THEN + ZMU1=0.3500212**2 + WEI(1)=1.0D0/3.0D0 + ELSE IF(NLF.EQ.6) THEN + ZMU1=0.2666355**2 + WEI(1)=0.1761263 + WEI(2)=0.1572071 + ELSE IF(NLF.EQ.8) THEN + ZMU1=1.0D0/21.0D0 + WEI(1)=0.1209877 + WEI(2)=0.0907407 + WEI(3)=0.0925926 + ELSE IF(NLF.EQ.12) THEN + ZMU1=0.1672126**2 + WEI(1)=0.0707626 + WEI(2)=0.0558811 + WEI(3)=0.0373377 + WEI(4)=0.0502819 + WEI(5)=0.0258513 + ELSE IF(NLF.EQ.16) THEN + ZMU1=0.1389568**2 + WEI(1)=0.0489872 + WEI(2)=0.0413296 + WEI(3)=0.0212326 + WEI(4)=0.0256207 + WEI(5)=0.0360486 + WEI(6)=0.0144586 + WEI(7)=0.0344958 + WEI(8)=0.0085179 + ELSE + CALL XABORT('SNQU01: ORDER NOT AVAILABLE.') + ENDIF + U(1)=REAL(SQRT(ZMU1)) + DO I=2,M2 + ZMU2=ZMU1+2.0D0*DBLE(I-1)*(1.0D0-3.0D0*ZMU1)/DBLE(NLF-2) + U(I)=REAL(SQRT(ZMU2)) + ENDDO +*---- +* COMPUTE THE POSITION OF WEIGHTS. +*---- + IPR=0 + INMAX=0 + DO IP=1,M2 + JOP(IP)=M2-IP+1 + DO IQ=1,JOP(IP) + IPR=IPR+1 + IF(IPR.GT.MAXNBA) CALL XABORT('SNQU01: MAXNBA OVERFLOW.') + TPQ(IPR)=U(IP) + UPQ(IPR)=U(M2+2-IP-IQ) + VPQ(IPR)=U(IQ) + IS=MIN(IP,IQ,M2+2-IP-IQ) + NW0=0 + DO II=1,IS-1 + NW0=NW0+(M2-3*(II-1)+1)/2 + ENDDO + KK=IP-IS+1 + LL=IQ-IS+1 + IF(KK.EQ.1)THEN + INWEI(IPR)=NW0+MIN(LL,M2-3*(IS-1)+1-LL) + ELSEIF(LL.EQ.1)THEN + INWEI(IPR)=NW0+MIN(KK,M2-3*(IS-1)+1-KK) + ELSE + INWEI(IPR)=NW0+MIN(KK,LL) + ENDIF + INMAX=MAX(INMAX,INWEI(IPR)) + ENDDO + ENDDO + IF(INMAX.NE.NW) CALL XABORT('SNQU01: INVALID VALUE OD NW.') + IF(IPR.NE.NPQ) CALL XABORT('SNQU01: BAD VALUE ON NPQ.') +*---- +* SET THE LEVEL-SYMMETRIC QUADRATURES. +*---- + IPQ=0 + WSUM=0.0 + DO IP=1,M2 + WSUM2=0.0D0 + DO IQ=1,JOP(IP) + IPQ=IPQ+1 + WPQ(IPQ)=REAL(WEI(INWEI(IPQ))*PI/2.0) + WSUM2=WSUM2+WEI(INWEI(IPQ)) + ENDDO + W(IP)=REAL(WSUM2) + WSUM=WSUM+REAL(WSUM2*PI/2.0) + ENDDO + RETURN + END |