diff options
Diffstat (limited to 'Dragon/src/SNFBC1.f')
| -rw-r--r-- | Dragon/src/SNFBC1.f | 304 |
1 files changed, 304 insertions, 0 deletions
diff --git a/Dragon/src/SNFBC1.f b/Dragon/src/SNFBC1.f new file mode 100644 index 0000000..fa89ac3 --- /dev/null +++ b/Dragon/src/SNFBC1.f @@ -0,0 +1,304 @@ +*DECK SNFBC1 + SUBROUTINE SNFBC1(LX,NMAT,IELEM,NLF,NSCT,U,MAT,VOL,TOTAL, + 1 NCODE,ZCODE,QEXT,LFIXUP,LSHOOT,ISBS,NBS,ISBSM,BS,WX,CST, + 2 ISADPTX,NUN,FUNKNO,MN,DN) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Perform one inner iteration for solving SN equations in 1D slab +* geometry. Albedo boundary conditions. +* +*Copyright: +* Copyright (C) 2020 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, A. A. Calloo and C. Bienvenue +* +*Parameters: input +* LX number of regions. +* NMAT number of material mixtures. +* IELEM measure of order of the spatial approximation polynomial: +* =1 constant - default for HODD; +* =2 linear - default for DG; +* >3 higher orders. +* NLF number of $\\mu$ levels. +* NSCT number of Legendre components in the flux: +* =1: isotropic sources; +* =2: linearly anisotropic sources. +* U base points in $\\mu$ of the SN quadrature. +* W weights of the SN quadrature. +* MN moment-to-discrete matrix. +* DN discrete-to-moment matrix. +* MAT material mixture index in each region. +* VOL volumes of each region. +* TOTAL macroscopic total cross sections. +* NCODE boundary condition indices. +* ZCODE albedos. +* QEXT Legendre components of the fixed source. +* LFIXUP flag to enable negative flux fixup. +* LSHOOT flag to enable/disable shooting method. +* ISBS flag to indicate the presence or not of boundary fixed +* sources. +* NBS number of boundary fixed sources. +* ISBSM flag array to indicate the presence or not of boundary fixed +* source in each unit surface. +* BS boundary source array with their intensities. +* WX spatial closure relation weighting factors. +* CST constants for the polynomial approximations. +* ISADPTX flag to enable/disable spatial adaptive flux calculations. +* NUN total number of unknowns in vector FUNKNO +* +*Parameters: input/output +* FUNKNO Legendre components of the flux and boundary fluxes. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER LX,NMAT,IELEM,NLF,NSCT,MAT(LX),NCODE(2),ISBS,NBS, + 1 ISBSM(2*ISBS,NLF*ISBS),NUN + LOGICAL LFIXUP,LSHOOT,ISADPTX + REAL U(NLF),VOL(LX),TOTAL(0:NMAT),ZCODE(2),QEXT(IELEM,NSCT,LX), + 1 FUNKNO(NUN),BS(NBS*ISBS),WX(IELEM+1),CST(IELEM),MN(NLF,NSCT), + 2 DN(NSCT,NLF) +*---- +* LOCAL VARIABLES +*---- + REAL WX0(IELEM+1) + DOUBLE PRECISION XNI,XNI1,XNI2,XNIA,XNIB,XNIA1,XNIA2,XNIB1,XNIB2 + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: Q + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: Q2 + PARAMETER(RLOG=1.0E-8) + LOGICAL ISSHOOT,ISFIX +*---- +* ALLOCATABLE ARRAYS +*---- + ALLOCATE(Q(IELEM),Q2(IELEM,IELEM+1)) +*---- +* LENGTH OF FUNKNO COMPONENTS (IN ORDER) +*---- + LFLX=IELEM*LX*NSCT + LXNI=NLF +*---- +* INNER ITERATION +*---- + + FUNKNO(1:LFLX)=0.0 + XNI=0.0D0 + XNI1=0.0D0 + XNI2=0.0D0 + XNIA=0.0D0 + XNIA2=0.0D0 + WX0=WX + + ! SHOOTING METHOD (ONLY IF THERE IS A NON-VACUUM RIGHT + ! BOUNDARY CONDITION. + ISSHOOT=(ZCODE(2).NE.0.0).AND.LSHOOT + IF(ISSHOOT) THEN + NS=6 + ELSE + NS=2 + ENDIF + + ! LOOP OVER ALL DIRECTIONS + DO 200 M0=1,NLF/2 + + ! LOOP FOR SHOOTING METHOD + DO 500 IS=1,NS + + ! CHOOSE DIRECTION + IF(MOD(IS,2).EQ.0) THEN + M=NLF-M0+1 ! FORWARD + ELSE + M=M0 ! BACKWARD + ENDIF + + ! SHOOTING METHOD BOUNDARY CONDITIONS. + IF(ISSHOOT) THEN + ! 1ST BACKWARD SWEEP + IF(IS.EQ.1) THEN + XNI=0.0D0 + XNI1=0.0D0 + XNI2=0.0D0 + ! 1ST FORWARD SWEEP + ELSEIF(IS.EQ.2) THEN + XNIA1=0.0D0 + IF(NCODE(1).EQ.4) THEN + XNIA1=REAL(XNI) + XNI=0.0D0 + ELSE + XNI=ZCODE(1)*REAL(XNI) + ENDIF + ! 2ND BACKWARD SWEEP + ELSEIF(IS.EQ.3) THEN + XNIA2=0.0D0 + XNIA=0.0D0 + IF(NCODE(1).EQ.4) THEN + XNIA2=REAL(XNI) + ELSE + XNIA=REAL(XNI) + ENDIF + XNI=1.0D0 + ! 2ND FORWARD SWEEP + ELSEIF(IS.EQ.4) THEN + IF(NCODE(1).EQ.4) THEN + XNIB1=REAL(XNI) + XNI1=XNIA1/(1.0D0+XNIA1-XNIB1) + XNI=1.0D0 + ELSE + XNI=ZCODE(1)*REAL(XNI) + ENDIF + ! 3RD BACKWARD SWEEP + ELSEIF(IS.EQ.5) THEN + IF(NCODE(1).EQ.4) THEN + XNIB2=REAL(XNI) + XNI2=XNIA2/(1.0D0+XNIA2-XNIB2) + XNI=XNI1 + ELSE + XNIB=REAL(XNI) + XNI=ZCODE(2)*XNIA/(1.0D0+ZCODE(2)*(XNIA-XNIB)) + ENDIF + ! 3RD FORWARD SWEEP + ELSEIF(IS.EQ.6) THEN + XNI=ZCODE(1)*XNI + IF(NCODE(1).EQ.4) XNI=XNI2 + ENDIF + ! NO SHOOTING METHOD BOUNDARY CONDITIONS + ELSE + IF(.NOT.LSHOOT) THEN + IF(U(M).GT.0.0) THEN + IF(NCODE(1).NE.4) FUNKNO(LFLX+M)=FUNKNO(LFLX+NLF-M+1) + ELSE + IF(NCODE(2).NE.4) FUNKNO(LFLX+M)=FUNKNO(LFLX+NLF-M+1) + ENDIF + XNI=0.0D0 + ELSE + IF(IS.EQ.1) THEN + XNI=0.0D0 + ELSE + XNI=ZCODE(1)*XNI + ENDIF + ENDIF + ENDIF + + ! X-BOUNDARIES CONDITIONS (NO SHOOTING) + IF(.NOT.LSHOOT) THEN + IF(U(M).GT.0.0) THEN + XNI=FUNKNO(LFLX+M)*ZCODE(1) + ELSE + XNI=FUNKNO(LFLX+M)*ZCODE(2) + ENDIF + ENDIF + + ! BOUNDARY FIXED SOURCES + IF(U(M).GT.0.0) THEN + IF(ISBS.EQ.1.AND.ISBSM(1,M).NE.0) XNI=XNI+BS(ISBSM(1,M)) + ELSE + IF(ISBS.EQ.1.AND.ISBSM(2,M).NE.0) XNI=XNI+BS(ISBSM(2,M)) + ENDIF + + ! SWEEPING OVER ALL VOXELS + DO 30 I0=1,LX + I=I0 + IF(U(M).LT.0.0) I=LX+1-I0 + + ! DATA + IBM=MAT(I) + SIGMA=TOTAL(IBM) + + ! SOURCE DENSITY TERM + DO IEL=1,IELEM + Q(IEL)=0.0 + DO L=1,NSCT + Q(IEL)=Q(IEL)+QEXT(IEL,L,I)*MN(M,L) + ENDDO + ENDDO + + ISFIX=.FALSE. + DO WHILE (.NOT.ISFIX) ! LOOP FOR ADAPTIVE CALCULATION + + ! FLUX MOMENTS CALCULATIONS + Q2(:IELEM,:IELEM+1)=0.0D0 + DO II=1,IELEM + DO JJ=1,IELEM + + ! MOMENT COEFFICIENTS + IF(II.EQ.JJ) THEN + Q2(II,JJ)=SIGMA*VOL(I)+CST(II)**2*WX(JJ+1)*ABS(U(M)) + ELSEIF(II.LT.JJ) THEN + IF(MOD(II+JJ,2).EQ.1) THEN + Q2(II,JJ)=CST(II)*CST(JJ)*WX(JJ+1)*U(M) + ELSE + Q2(II,JJ)=CST(II)*CST(JJ)*WX(JJ+1)*ABS(U(M)) + ENDIF + ELSE + IF(MOD(II+JJ,2).EQ.1) THEN + Q2(II,JJ)=CST(II)*CST(JJ)*(WX(JJ+1)-2.0D0)*U(M) + ELSE + Q2(II,JJ)=CST(II)*CST(JJ)*WX(JJ+1)*ABS(U(M)) + ENDIF + ENDIF + ENDDO + ENDDO + + ! SOURCE TERMS + DO II=1,IELEM + IF(MOD(II,2).EQ.1) THEN + Q2(II,IELEM+1)=Q(II)*VOL(I)+CST(II)*(1-WX(1))*ABS(U(M))*XNI + ELSE + Q2(II,IELEM+1)=Q(II)*VOL(I)-CST(II)*(1+WX(1))*U(M)*XNI + ENDIF + ENDDO + + CALL ALSBD(IELEM,1,Q2,IER,IELEM) + IF(IER.NE.0) CALL XABORT('SNFBC1: SINGULAR MATRIX.') + + ! ADAPTIVE CORRECTION OF WEIGHTING PARAMETERS + IF(ISADPTX) THEN + CALL SNADPT(IELEM,IELEM,1,Q2(:IELEM,IELEM+1),XNI, + 1 1.0,WX,ISFIX) + ELSE + ISFIX=.TRUE. + ENDIF + + END DO ! END OF ADAPTIVE LOOP + + ! CLOSURE RELATIONS + IF(IELEM.EQ.1.AND.LFIXUP.AND.(Q2(1,2).LE.RLOG)) Q2(1,2)=0.0D0 + XNI=WX(1)*XNI + DO II=1,IELEM + IF(MOD(II,2).EQ.1) THEN + XNI=XNI+CST(II)*WX(II+1)*Q2(II,IELEM+1) + ELSE + XNI=XNI+CST(II)*WX(II+1)*Q2(II,IELEM+1)*SIGN(1.0,U(M)) + ENDIF + ENDDO + IF(IELEM.EQ.1.AND.LFIXUP.AND.(XNI.LE.RLOG)) XNI=0.0D0 + WX=WX0 + + IF(ISSHOOT.AND.IS.LT.5) GO TO 30 + + ! SAVE LEGENDRE MOMENT OF THE FLUX + DO L=1,NSCT + DO IEL=1,IELEM + IOF=(I-1)*NSCT*IELEM+(L-1)*IELEM+IEL + FUNKNO(IOF)=FUNKNO(IOF)+REAL(Q2(IEL,IELEM+1))*DN(L,M) + ENDDO + ENDDO + + 30 CONTINUE ! END OF X-LOOP + + ! SAVE BOUNDARIES FLUX + IF(.NOT.LSHOOT) FUNKNO(LFLX+M)=REAL(XNI) + + 500 CONTINUE ! END OF SHOOTING METHOD LOOP + 200 CONTINUE ! END OF DIRECTION LOOP + + DEALLOCATE(Q,Q2) + RETURN + END |