Initial commit from Polytechnique Montreal

1 files changed, 178 insertions, 0 deletions

diff --git a/Dragon/src/SNFC12.f b/Dragon/src/SNFC12.f
new file mode 100644
index 0000000..f0a8df2
--- /dev/null
+++ b/Dragon/src/SNFC12.f
@@ -0,0 +1,178 @@
+*DECK SNFC12
+      SUBROUTINE SNFC12(LX,LY,NMAT,NPQ,NSCT,MAT,VOL,TOTAL,NCODE,ZCODE,
+     1 QEXT,LFIXUP,DU,DE,W,MRM,MRMY,DB,DA,DAL,FLUX,XNEI,XNEJ,MN,DN)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform one inner iteration for solving SN equations in 2D R-Z
+* geometry for the diamond differencing method. Albedo boundary
+* conditions.
+*
+*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
+* LX      number of meshes along X axis.
+* LY      number of meshes along Y axis.
+* NMAT    number of material mixtures.
+* NPQ     number of SN directions in four octants (including zero-weight
+*         directions).
+* NSCT    maximum number of spherical harmonics moments of the flux.
+* 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.
+* DU      first direction cosines ($\\mu$).
+* DE      second direction cosines ($\\eta$).
+* W       weights.
+* MRM     quadrature index.
+* MRMY    quadrature index.
+* DB      diamond-scheme parameter.
+* DA      diamond-scheme parameter.
+* DAL     diamond-scheme parameter.
+* MN      moment-to-discrete matrix.
+* DN      discrete-to-moment matrix.
+*
+*Parameters: input/output
+* XNEI    X-directed SN boundary fluxes.
+* XNEJ    Y-directed SN boundary fluxes.
+*
+*Parameters: output
+* FLUX    Legendre components of the flux.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER LX,LY,NMAT,NPQ,NSCT,MAT(LX,LY),NCODE(4),MRM(NPQ),MRMY(NPQ)
+      REAL VOL(LX,LY),TOTAL(0:NMAT),ZCODE(4),QEXT(NSCT,LX,LY),DU(NPQ),
+     1 DE(NPQ),W(NPQ),DB(LX,NPQ),DA(LX,LY,NPQ),DAL(LX,LY,NPQ),
+     2 FLUX(NSCT,LX,LY),XNEI(LY,NPQ),XNEJ(LX,NPQ),MN(NPQ,NSCT),
+     3 DN(NSCT,NPQ)
+      LOGICAL LFIXUP
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER P
+      DOUBLE PRECISION QQ,C1,XNM,XNJ,Q2(1,2)
+      PARAMETER(RLOG=1.0E-8,PI=3.141592654)
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: XARN
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: XNI
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(XARN(LX,LY),XNI(LY))
+*----
+*  MAIN LOOP OVER SN ANGLES.
+*----
+      FLUX(:NSCT,:LX,:LY)=0.0
+      XARN(:LX,:LY)=0.0
+      C1=0.0
+      XNM=0.0
+      DO 170 M=1,NPQ
+      WEIGHT=W(M)
+      VU=DU(M)
+      VE=DE(M)
+      IF(NCODE(1).NE.4) THEN
+         M1=MRM(M)
+         IF(WEIGHT.EQ.0.0) THEN
+            DO 10 J=1,LY
+            XNEI(J,M)=XNEI(J,M1)
+   10       CONTINUE
+         ELSE IF(VU.GT.0.0) THEN
+            DO 20 J=1,LY
+            E1=XNEI(J,M)
+            XNEI(J,M)=XNEI(J,M1)
+            XNEI(J,M1)=E1
+   20       CONTINUE
+         ENDIF
+      ENDIF
+      IF(NCODE(3).NE.4) THEN
+         M1=MRMY(M)
+         IF(VE.GT.0) THEN
+            DO 40 I=1,LX
+            E1=XNEJ(I,M)
+            XNEJ(I,M)=XNEJ(I,M1)
+            XNEJ(I,M1)=E1
+   40       CONTINUE
+         ENDIF
+      ENDIF
+      IF(VE.GT.0.0) GOTO 70
+      IF(VU.GT.0.0) GOTO 60
+      IND=3
+      GOTO 90
+   60 IND=4
+      GOTO 90
+   70 IF(VU.GT.0.0) GOTO 80
+      IND=2
+      GOTO 90
+   80 IND=1
+*----
+*  LOOP OVER X- AND Y-DIRECTED AXES.
+*----
+   90 DO 155 II=1,LX
+      I=II
+      IF((IND.EQ.2).OR.(IND.EQ.3)) I=LX+1-II
+      IF((IND.EQ.1).OR.(IND.EQ.2)) THEN
+         XNJ=XNEJ(I,M)*ZCODE(3)
+      ELSE
+         XNJ=XNEJ(I,M)*ZCODE(4)
+      ENDIF
+      DO 140 JJ=1,LY
+      J=JJ
+      IF((IND.EQ.3).OR.(IND.EQ.4)) J=LY+1-JJ
+      C1=DAL(I,J,M)
+      IF(II.EQ.1) THEN
+         IF((IND.EQ.1).OR.(IND.EQ.4)) THEN
+            XNI(J)=XNEI(J,M)
+         ELSE
+            XNI(J)=XNEI(J,M)*ZCODE(2)
+         ENDIF
+      ENDIF
+      IF(MAT(I,J).EQ.0) GO TO 140
+      QQ=0.0D0
+      DO 110 P=1,NSCT
+      QQ=QQ+QEXT(P,I,J)*MN(M,P)
+  110 CONTINUE
+      VT=VOL(I,J)*TOTAL(MAT(I,J))
+      XNM=XARN(I,J)
+      Q2(1,1)=C1+2.0D0*ABS(DA(I,J,M))+2.0D0*ABS(DB(I,M))+VT
+      Q2(1,2)=C1*XNM+2.0D0*ABS(DA(I,J,M))*XNI(J)+2.0D0*ABS(DB(I,M))
+     1          *XNJ+VOL(I,J)*QQ
+      IF(Q2(1,1).EQ.0.0D0) CALL XABORT('SNFC12: SINGULAR MATRIX.')
+      Q2(1,2)=Q2(1,2)/Q2(1,1)
+      IF(LFIXUP.AND.(Q2(1,2).LE.RLOG)) Q2(1,2)=0.0
+      XNI(J)=2.0D0*Q2(1,2)-XNI(J)
+      XNJ=2.0D0*Q2(1,2)-XNJ
+      XARN(I,J)=REAL(2.0D0*Q2(1,2)-XNM)
+      IF(LFIXUP.AND.(XARN(I,J).LE.RLOG)) XARN(I,J)=0.0
+      IF(W(M).LE.RLOG) XARN(I,J)=REAL(Q2(1,2))
+      IF(LFIXUP.AND.(XNI(J).LE.RLOG)) XNI(J)=0.0
+      IF(LFIXUP.AND.(XNJ.LE.RLOG)) XNJ=0.0
+      DO 135 P=1,NSCT
+      FLUX(P,I,J)=FLUX(P,I,J)+REAL(Q2(1,2))*DN(P,M)
+  135 CONTINUE
+  140 CONTINUE
+      XNEJ(I,M)=REAL(XNJ)
+  155 CONTINUE
+      DO 165 J=1,LY
+      XNEI(J,M)=REAL(XNI(J))
+  165 ENDDO
+  170 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(XNI,XARN)
+      RETURN
+      END