Initial commit from Polytechnique Montreal

1 files changed, 598 insertions, 0 deletions

diff --git a/Dragon/src/SNFBH2.F b/Dragon/src/SNFBH2.F
new file mode 100644
index 0000000..6ff2f01
--- /dev/null
+++ b/Dragon/src/SNFBH2.F
@@ -0,0 +1,598 @@
+*DECK SNFBH2
+      SUBROUTINE SNFBH2(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,ISPLH,SIDE,
+     1 IELEM,NM,NMX,NMY,NMAT,NPQ,NSCT,MAT,VOL,TOTAL,QEXT,LFIXUP,DU,DE,W,
+     2 DB,DA,MN,DN,WX,WY,CST,ISADPT,LOZSWP,COORDMAP,FUNKNO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform one inner iteration for solving SN equations in 2D hexagonal
+* geometry for the HODD/DG method. Energy-angle multithreading. VOID
+* boundary conditions. Boltzmann (BTE) discretization.
+
+*
+*Copyright:
+* Copyright (C) 2025 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
+* NUN      total number of unknowns in vector FUNKNO.
+* NGEFF    number of energy groups processed in parallel.
+* IMPX     print flag (equal to zero for no print).
+* INCONV   energy group convergence flag (set to .FALSE. if converged).
+* NGIND    energy group indices assign to the NGEFF set.
+* NHEX     number of hexagons in X-Y plane.
+* ISPLH    splitting option for hexagons.
+* SIDE     side of an hexagon.
+* IELEM    measure of order of the spatial approximation polynomial:
+*          =1 constant - default for HODD;
+*          =2 linear - default for DG;
+*          >3 higher orders.
+* NM       number of moments in space and energy for flux components
+* NMX      number of moments for X axis boundaries components
+* NMY      number of moments for Y axis boundaries components
+* 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.
+* 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.
+* DB       diamond-scheme parameter.
+* DA       diamond-scheme parameter.
+* MN       moment-to-discrete matrix.
+* DN       discrete-to-moment matrix.
+* WX       spatial X axis closure relation weighting factors.
+* WY       spatial Y axis closure relation weighting factors.
+* CST      constants for the polynomial approximations.
+* ISADPT   flag to enable/disable adaptive flux calculations.
+* LOZSWP   lozenge sweep order depending on direction.
+* COORDMAP coordinate map - mapping the hexagons from the indices 
+*          within the DRAGON geometry to a Cartesian axial coordinate
+*          array (see redblobgames.com website).
+*
+*Parameters: input/output
+* FUNKNO  Legendre components of the flux and boundary fluxes.
+*
+*Comments:
+* 1. The direction of the axes I, J and D for the surface boundary 
+*    fluxes are shown in the diagram below. This means that 
+*    i) lozenge A has I- and D-boundaries (instead of I and J)
+*    i) lozenge B has I- and J-boundaries
+*    i) lozenge C has D- and J-boundaries (instead of I and J)
+*
+*                                  ^
+*                         j-axis   |
+*                                  |          ^
+*                              _________     /    d-axis
+*                             /       / \   /
+*                            /   B   /   \
+*                           /       /     \
+*                          (-------(   A   )
+*                           \       \     /
+*                            \  C    \   / 
+*                             \_______\_/   \
+*                                            \   i-axis
+*                                             ^
+*
+*-----------------------------------------------------------------------
+#if defined(_OPENMP)
+      USE omp_lib
+#endif
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NUN,NGEFF,IMPX,NGIND(NGEFF),NHEX,ISPLH,IELEM,NM,NMX,
+     1 NMY,NMAT,NPQ,NSCT,MAT(ISPLH,ISPLH,3,NHEX),LOZSWP(3,6),
+     2 COORDMAP(3,NHEX)
+      LOGICAL INCONV(NGEFF)
+      REAL SIDE,VOL(ISPLH,ISPLH,3,NHEX),TOTAL(0:NMAT,NGEFF),
+     1 QEXT(NUN,NGEFF),DU(NPQ),DE(NPQ),W(NPQ),
+     2 DB(ISPLH,ISPLH,3,NHEX,NPQ),DA(ISPLH,ISPLH,3,NHEX,NPQ),
+     3 MN(NPQ,NSCT),DN(NSCT,NPQ),FUNKNO(NUN,NGEFF),WX(IELEM+1),
+     3 WY(IELEM+1),CST(IELEM)
+      LOGICAL LFIXUP,ISADPT(2)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER :: NPQD(6),IIND(6),P,DCOORD
+      REAL :: JAC(2,2,3), MUH, ETAH, AAA, BBB, CCC, DDD, MUHTEMP,
+     1 ETAHTEMP, WX0(IELEM+1),WY0(IELEM+1)
+      DOUBLE PRECISION :: Q(NM), Q2(NM,NM+1), V,THETA, XNI(NMX),
+     1 XNJ(NMY)
+      PARAMETER(IUNOUT=6,RLOG=1.0E-8,PI=3.141592654)
+      LOGICAL :: LHEX(NHEX)
+      LOGICAL ISFIX(2)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: INDANG
+      INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: TMPMAT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: FLUX
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: FLUX_G
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: TMPXNI,
+     > TMPXNJ, TMPXND
+*----
+*  MAP MATERIAL VALUES TO CARTESIAN AXIAL COORDINATE MAP
+*----
+      NRINGS=INT((SQRT(  REAL((4*NHEX-1)/3)  )+1.)/2.)
+      NCOL=2*NRINGS -1
+      ALLOCATE(TMPMAT(ISPLH,ISPLH,3,NCOL,NCOL))
+      TMPMAT(:,:,:,:,:) = -1
+      DO IHEX_DOM=1,NHEX
+         TMPMAT(:,:,:,COORDMAP(1,IHEX_DOM),COORDMAP(2,IHEX_DOM)) = 
+     >         MAT(:,:,:,IHEX_DOM)
+      ENDDO
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(INDANG(NPQ,6))
+      ALLOCATE(FLUX(NM,NSCT,3*ISPLH**2,NHEX))
+      ALLOCATE(FLUX_G(NM,NSCT,3*ISPLH**2,NHEX,NGEFF))
+      ALLOCATE(TMPXNI(IELEM,ISPLH,NCOL))
+      ALLOCATE(TMPXNJ(IELEM,ISPLH,NCOL))
+      ALLOCATE(TMPXND(IELEM,ISPLH,NCOL))
+*----
+*  CONSTRUCT JACOBIAN MATRIX FOR EACH LOZENGE
+*----
+      JAC = RESHAPE((/ 1., -SQRT(3.), 1., SQRT(3.), 2., 0., 1.,
+     >    SQRT(3.), 2., 0., -1., SQRT(3.) /), SHAPE(JAC))
+      JAC = (SIDE/2.)*JAC
+*----
+*  LENGTH OF FUNKNO COMPONENTS (IN ORDER)
+*----
+      LFLX=3*NM*(ISPLH**2)*NHEX*NSCT
+*----
+*  SET DODECANT SWAPPING ORDER
+*----
+      NPQD(:6)=0
+      INDANG(:NPQ,:6)=0
+      IIND(:6)=0
+      DO M=1,NPQ
+        VU=DU(M)
+        VE=DE(M)
+        IF(W(M).EQ.0) CYCLE
+        THETA=0.0D0
+        IF(VE.GT.0.0)THEN
+          IF(VU.EQ.0.0)THEN
+            THETA = PI/2
+          ELSEIF(VU.GT.0.0)THEN
+            THETA = ATAN(ABS(VE/VU))
+          ELSEIF(VU.LT.0.0)THEN
+            THETA = PI - ATAN(ABS(VE/VU))
+          ENDIF
+        ELSEIF(VE.LT.0.0)THEN
+          IF(VU.EQ.0.0)THEN
+            THETA = 3*PI/2
+          ELSEIF(VU.LT.0.0)THEN
+            THETA = PI + ATAN(ABS(VE/VU))
+          ELSEIF(VU.GT.0.0)THEN
+            THETA = 2.*PI - ATAN(ABS(VE/VU))
+          ENDIF
+        ENDIF
+        IND=0
+        IF((THETA.GT.0.0).AND.(THETA.LT.(PI/3.)))THEN
+          IND=1
+        ELSEIF((THETA.GT.(PI/3.)).AND.(THETA.LT.(2.*PI/3.)))THEN
+          IND=2
+        ELSEIF((THETA.GT.(2.*PI/3.)).AND.(THETA.LT.(PI)))THEN
+          IND=3
+        ELSEIF((THETA.GT.(PI)).AND.(THETA.LT.(4.*PI/3.)))THEN
+          IND=4
+        ELSEIF((THETA.GT.(4.*PI/3.)).AND.(THETA.LT.(5.*PI/3.)))THEN
+          IND=5
+        ELSEIF((THETA.GT.(5.*PI/3.)).AND.(THETA.LT.(2.*PI)))THEN
+          IND=6
+        ENDIF
+        ! Assume IIND(I)=I in hexagonal geometry
+        IIND(IND)=IND
+        NPQD(IND)=NPQD(IND)+1
+        INDANG(NPQD(IND),IND)=M
+      ENDDO
+*----
+*  MAIN LOOP OVER DODECANTS
+*----
+
+      FLUX_G(:NM,:NSCT,:3*ISPLH**2,:NHEX,:NGEFF)=0.0D0
+      WX0=WX
+      WY0=WY
+
+      DO JND=1,6
+      IND=IIND(JND)
+      ! Needed because of S2 LS (4 dir. for 6 sextants)
+      IF(IND.EQ.0) CYCLE
+
+*----
+*  MAIN SWAPPING LOOPS FOR SN FLUX CALCULATION
+*  LOOP OVER ENERGY AND ANGLES
+*----
+
+*$OMP  PARALLEL DO 
+*$OMP+ PRIVATE(ITID,FLUX,M,IG,XNI,XNJ,Q,Q2,IOF,IER,II,JJ,IEL,I,J,P)
+*$OMP+ PRIVATE(IPQD,IBM,SIGMA,V,ISFIX,IX,JX,IY,JY,AAA,BBB,CCC,DDD)
+*$OMP+ PRIVATE(LHEX,IHEX,IIHEX,DCOORD,ILOZLOOP,ILOZ,IL,I2,JL,J2)
+*$OMP+ PRIVATE(MUHTEMP,MUH,ETAHTEMP,ETAH,I3,I_FETCH,III,JJJ,IIM,JIM)
+*$OMP+ PRIVATE(TMPXNI,TMPXNJ,TMPXND)
+*$OMP+ FIRSTPRIVATE(WX,WY,WX0,WY0) SHARED(FUNKNO)
+*$OMP+ REDUCTION(+:FLUX_G) COLLAPSE(2)
+
+      ! LOOP FOR GROUPS TO EXECUTE IN PARALLEL
+      DO IG=1,NGEFF
+
+      ! LOOP OVER ALL DIRECTIONS
+      DO IPQD=1,NPQD(IND)
+      IF(.NOT.INCONV(IG)) CYCLE
+      M=INDANG(IPQD,IND)
+      IF(W(M).EQ.0.0) CYCLE
+
+      ! GET AND PRINT THREAD NUMBER
+#if defined(_OPENMP)
+      ITID=omp_get_thread_num()
+#else
+      ITID=0
+#endif
+      IF(IMPX.GT.5) WRITE(IUNOUT,400) ITID,NGIND(IG),IPQD
+
+      ! INITIALIZE FLUXES AND BOUNDARY FLUXES
+      FLUX(:NM,:NSCT,:3*ISPLH**2,:NHEX)=0.0D0
+      TMPXNI(:IELEM,:ISPLH,:NCOL)=0.0D0
+      TMPXNJ(:IELEM,:ISPLH,:NCOL)=0.0D0
+      TMPXND(:IELEM,:ISPLH,:NCOL)=0.0D0
+
+      ! LOOP OVER CARTESIAN MAP OF HEXAGONAL DOMAIN
+      DO JJJ=1,NCOL
+      JIM=JJJ
+      ! Account for different sweep direction depending on angle 
+      IF((IND.EQ.1).OR.(IND.EQ.2).OR.(IND.EQ.3)) JIM=NCOL+1-JIM
+        
+      DO III=1,NCOL
+      IIM=III
+      ! Account for different sweep direction depending on angle 
+      IF((IND.EQ.2).OR.(IND.EQ.3).OR.(IND.EQ.4)) IIM=NCOL+1-IIM
+
+      ! For IND 3 or 6, Cartesian axial coordinate map is swept 
+      ! vertically instead of horizontally. IM suffix is for 'IMmutable'
+      I=IIM
+      J=JIM
+      IF((IND.EQ.3).OR.(IND.EQ.6))THEN
+        I=JIM
+        J=IIM
+      ENDIF
+
+      ! If within corners of Cartesian axial coordinate map (where
+      ! there are no hexagons), skip loop
+      IF(TMPMAT(1,1,1,I,J).EQ.-1) CYCLE
+
+      ! Find DRAGON geometry hexagonal index using I and J
+      LHEX=(COORDMAP(1,:).EQ.I .AND. COORDMAP(2,:).EQ.J)
+      IHEX=0
+      DO IIHEX=1,NHEX
+        IF(LHEX(IIHEX)) THEN
+          IHEX=IIHEX
+          EXIT
+        ENDIF
+      ENDDO
+      IF(IHEX.EQ.0) CALL XABORT('SNFTH1: IHEX FAILURE.')
+      ! Find D coordinate
+      DCOORD = ABS(COORDMAP(3,IHEX))-NRINGS
+
+      ! LOOP OVER LOZENGES
+      DO ILOZLOOP=1,3
+      ILOZ=LOZSWP(ILOZLOOP,IND)
+
+      ! Get Jacobian elements values
+      AAA = JAC(1,1,ILOZ)
+      BBB = JAC(1,2,ILOZ)
+      CCC = JAC(2,1,ILOZ)
+      DDD = JAC(2,2,ILOZ)
+
+
+      ! LOOP OVER SUBMESH WITHIN EACH LOZENGE
+      DO IL=1,ISPLH
+      I2=IL
+      ! Account for different sweep direction depending on angle 
+      IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+        IF((IND.EQ.2).OR.(IND.EQ.3).OR.(IND.EQ.4)) I2=ISPLH+1-I2
+      ELSEIF(ILOZ.EQ.3)THEN
+        IF((IND.EQ.3).OR.(IND.EQ.4).OR.(IND.EQ.5)) I2=ISPLH+1-I2
+      ENDIF
+
+      DO JL=1,ISPLH
+      J2=JL
+      ! Account for different sweep direction depending on angle 
+      IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+        IF((IND.EQ.4).OR.(IND.EQ.5).OR.(IND.EQ.6)) J2=ISPLH+1-J2
+      ELSEIF(ILOZ.EQ.1)THEN
+        IF((IND.EQ.3).OR.(IND.EQ.4).OR.(IND.EQ.5)) J2=ISPLH+1-J2
+      ENDIF
+
+      ! READ IN XNI AND XNJ DEPENDING ON LOZENGE
+      I_FETCH=0
+      IF((ILOZ.EQ.1))THEN
+        ! Read boundary fluxes in reverse for lozenge A since affine
+        ! transformation of lozenges causes the D and I directions
+        ! of lozenges C and A respectively to be reversed
+        I_FETCH=ISPLH+1-I2
+        XNI(:) = TMPXNI(:,J2,J)
+        XNJ(:) = TMPXND(:,I_FETCH,DCOORD)
+      ELSEIF((ILOZ.EQ.2))THEN
+        XNI(:) = TMPXNI(:,J2,J)
+        XNJ(:) = TMPXNJ(:,I2,I)
+      ELSEIF((ILOZ.EQ.3))THEN
+        XNI(:) = TMPXND(:,J2,DCOORD)
+        XNJ(:) = TMPXNJ(:,I2,I)
+      ENDIF
+
+      ! DATA
+      IBM=MAT(I2,J2,ILOZ,IHEX)
+      ! Skip loop if virtual element 
+      IF(IBM.EQ.0) CYCLE
+      SIGMA=TOTAL(IBM,IG)
+      V=VOL(I2,J2,ILOZ,IHEX)
+
+      ! COMPUTE ADJUSTED DIRECTION COSINES
+      MUHTEMP  =  DA(I2,J2,ILOZ,IHEX,M)
+      ETAHTEMP =  DB(I2,J2,ILOZ,IHEX,M)
+      MUH = (MUHTEMP*DDD) - (ETAHTEMP*BBB)
+      ETAH = (-MUHTEMP*CCC) + (ETAHTEMP*AAA)
+
+      ! SOURCE DENSITY TERM
+      DO IEL=1,NM
+      Q(IEL)=0.0D0
+      DO P=1,NSCT
+      IOF=((((((IHEX-1)*3+(ILOZ-1))*ISPLH+(J2-1))*ISPLH+
+     1   (I2-1))*NSCT+(P-1))*NM)+IEL
+      Q(IEL)=Q(IEL)+QEXT(IOF,IG)*MN(M,P)
+      ENDDO
+      ENDDO
+
+      ISFIX=.FALSE.
+      DO WHILE (.NOT.ALL(ISFIX)) ! LOOP FOR ADAPTIVE CALCULATION
+
+      ! FLUX MOMENT COEFFICIENTS MATRIX
+      Q2(:NM,:NM+1)=0.0D0
+      
+      DO IY=1,IELEM
+      DO JY=1,IELEM
+      DO IX=1,IELEM
+      DO JX=1,IELEM
+        II=IELEM*(IY-1)+IX
+        JJ=IELEM*(JY-1)+JX
+
+        ! DIAGONAL TERMS
+        IF(II.EQ.JJ) THEN
+          Q2(II,JJ)=SIGMA*V
+     1              +CST(IX)**2*WX(JX+1)*ABS(MUH)
+     2              +CST(IY)**2*WY(JY+1)*ABS(ETAH)
+
+        ! UPPER DIAGONAL TERMS
+        ELSEIF(II.LT.JJ) THEN
+          ! X-SPACE TERMS
+          IF(IY.EQ.JY) THEN
+          IF(MOD(IX+JX,2).EQ.1) THEN
+            Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*MUH
+          ELSE
+            Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(MUH)
+          ENDIF
+          ! Y-SPACE TERMS
+          ELSEIF(IX.EQ.JX) THEN
+          IF(MOD(IY+JY,2).EQ.1) THEN
+            Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ETAH
+          ELSE
+            Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ABS(ETAH)
+          ENDIF
+          ENDIF
+
+        ! UNDER DIAGONAL TERMS
+        ELSE
+          ! X-SPACE TERMS
+          IF(IY.EQ.JY) THEN
+          IF(MOD(IX+JX,2).EQ.1) THEN
+            Q2(II,JJ)=CST(IX)*CST(JX)*(WX(JX+1)-2)*MUH
+          ELSE
+            Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(MUH)
+          ENDIF 
+          ! Y-SPACE TERMS
+          ELSEIF(IX.EQ.JX) THEN
+          IF(MOD(IY+JY,2).EQ.1) THEN
+            Q2(II,JJ)=CST(IY)*CST(JY)*(WY(JY+1)-2)*ETAH
+          ELSE
+            Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ABS(ETAH)
+          ENDIF
+          ENDIF
+
+        ENDIF
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+
+      ! FLUX SOURCE VECTOR
+      DO IY=1,IELEM
+      DO IX=1,IELEM
+        II=IELEM*(IY-1)+IX
+        Q2(II,NM+1)=Q(II)*V
+        ! X-SPACE TERMS
+        IF(MOD(IX,2).EQ.1) THEN
+          Q2(II,NM+1)=Q2(II,NM+1)+CST(IX)*(1-WX(1))
+     1                *XNI(IY)*ABS(MUH)
+        ELSE
+          Q2(II,NM+1)=Q2(II,NM+1)-CST(IX)*(1+WX(1))
+     1                *XNI(IY)*MUH
+        ENDIF
+        ! Y-SPACE TERMS
+        IF(MOD(IY,2).EQ.1) THEN
+          Q2(II,NM+1)=Q2(II,NM+1)+CST(IY)*(1-WY(1))
+     1                *XNJ(IX)*ABS(ETAH)
+        ELSE
+          Q2(II,NM+1)=Q2(II,NM+1)-CST(IY)*(1+WY(1))
+     1                *XNJ(IX)*ETAH
+        ENDIF
+      ENDDO
+      ENDDO
+
+      CALL ALSBD(NM,1,Q2,IER,NM)
+      IF(IER.NE.0) CALL XABORT('SNFBH2: SINGULAR MATRIX.')
+
+      ! ADAPTIVE CORRECTION OF WEIGHTING PARAMETERS
+      IF(ANY(ISADPT)) THEN
+        IF(ISADPT(1)) THEN
+          CALL SNADPT(IELEM,NM,IELEM,Q2(1:IELEM:1,NM+1),
+     1    XNI(:NMX),1.0,WX,ISFIX(1))
+        ELSE
+          ISFIX(1)=.TRUE.
+        ENDIF
+        IF(ISADPT(2)) THEN
+          CALL SNADPT(IELEM,NM,IELEM,Q2(1:NM:IELEM,NM+1),
+     1    XNJ(:NMY),1.0,WY,ISFIX(2))
+        ELSE
+          ISFIX(2)=.TRUE.
+        ENDIF
+      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.0
+      ! Read XNI/XNI into TMPXNI/J/D
+      IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+        TMPXNI(:NMX,J2,J)=WX(1)*XNI(:NMX)
+      ELSE
+        TMPXND(:NMX,J2,DCOORD)=WX(1)*XNI(:NMX)
+      ENDIF
+      IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+        TMPXNJ(:NMY,I2,I)=WY(1)*XNJ(:NMY)
+      ELSE
+        I3=I_FETCH
+        TMPXND(:NMY,I3,DCOORD)=WY(1)*XNJ(:NMY)
+      ENDIF
+      DO IY=1,IELEM
+      DO IX=1,IELEM
+        II=IELEM*(IY-1)+IX
+        ! X-SPACE
+        ! Assign I-boundary fluxes if lozenges A or B
+        IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+        IF(MOD(IX,2).EQ.1) THEN
+          TMPXNI(IY,J2,J)=TMPXNI(IY,J2,J)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXNI(IY,J2,J)=TMPXNI(IY,J2,J)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,MUH)
+        ENDIF
+        ENDIF
+        ! Y-SPACE
+        ! Assign J-boundary fluxes if lozenges B or C
+        IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+        IF(MOD(IY,2).EQ.1) THEN
+          TMPXNJ(IX,I2,I)=TMPXNJ(IX,I2,I)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXNJ(IX,I2,I)=TMPXNJ(IX,I2,I)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,ETAH)
+        ENDIF
+        ENDIF
+        ! D-SPACE
+        ! Assign D-boundary fluxes if lozenge A using XNJ
+        IF((ILOZ.EQ.1))THEN
+        I3=I_FETCH
+        IF(MOD(IY,2).EQ.1) THEN
+          TMPXND(IX,I3,DCOORD)=TMPXND(IX,I3,DCOORD)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXND(IX,I3,DCOORD)=TMPXND(IX,I3,DCOORD)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,ETAH)
+        ENDIF
+        ENDIF
+        ! Assign D-boundary fluxes if lozenge C using XNI
+        IF((ILOZ.EQ.3))THEN
+        IF(MOD(IX,2).EQ.1) THEN
+          TMPXND(IY,J2,DCOORD)=TMPXND(IY,J2,DCOORD)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXND(IY,J2,DCOORD)=TMPXND(IY,J2,DCOORD)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,MUH)
+        ENDIF
+        ENDIF
+      ENDDO
+      ENDDO
+      ! FLIP GRADIENTS IF NECESSARY
+      DO IY=1,IELEM
+        IF((MOD(IY,2).EQ.0).AND.(ILOZ.EQ.3).AND.(IL.EQ.ISPLH))
+     1    TMPXND(IY,J2,DCOORD)=TMPXND(IY,J2,DCOORD)*(-1)
+      ENDDO
+      I3=I_FETCH
+      DO IX=1,IELEM
+        IF((MOD(IX,2).EQ.0).AND.(ILOZ.EQ.1).AND.(JL.EQ.ISPLH))
+     1    TMPXND(IX,I3,DCOORD)=TMPXND(IX,I3,DCOORD)*(-1)
+      ENDDO
+      ! LFIXUP
+      IF(IELEM.EQ.1.AND.LFIXUP)THEN
+        IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+          IF(TMPXNI(1,J2,J).LE.RLOG) TMPXNI(1,J2,J)=0.0
+        ELSE
+          IF(TMPXND(1,J2,DCOORD).LE.RLOG) TMPXND(1,J2,DCOORD)=0.0
+        ENDIF
+        IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+          IF(TMPXNJ(1,I2,I).LE.RLOG) TMPXNJ(1,I2,I)=0.0
+        ELSE
+          I3=I_FETCH
+          IF(TMPXND(1,I3,DCOORD).LE.RLOG) TMPXND(1,I3,DCOORD)=0.0
+        ENDIF
+      ENDIF
+      WX=WX0
+      WY=WY0
+
+      ! SAVE LEGENDRE MOMENT OF THE FLUX
+      IOF=((ILOZ-1)*ISPLH+(J2-1))*ISPLH+I2
+      DO P=1,NSCT
+      DO IEL=1,NM
+      FLUX(IEL,P,IOF,IHEX)=FLUX(IEL,P,IOF,IHEX)+Q2(IEL,NM+1)*DN(P,M)
+      ENDDO
+      ENDDO
+
+      ENDDO ! END OF WITHIN LOZENGE J-LOOP
+      ENDDO ! END OF WITHIN LOZENGE I-LOOP
+
+      ENDDO ! END OF LOZENGE LOOP
+
+      ENDDO ! END OF I COLUMNS OF CARTESIAN MAP LOOP
+      ENDDO ! END OF J COLUMNS OF CARTESIAN MAP LOOP
+
+      ! SAVE FLUX INFORMATION
+      FLUX_G(:,:,:,:,IG)=FLUX_G(:,:,:,:,IG)+FLUX(:,:,:,:)
+
+      ENDDO ! END OF DIRECTION LOOP
+      ENDDO ! END OF ENERGY LOOP
+*$OMP END PARALLEL DO
+      ENDDO ! END OF OCTANT LOOP
+
+      ! SAVE FLUX INFORMATION
+      DO IG=1,NGEFF
+        IF(.NOT.INCONV(IG)) CYCLE
+        FUNKNO(:LFLX,IG)=
+     1  RESHAPE(REAL(FLUX_G(:IELEM**2,:NSCT,:3*ISPLH**2,:NHEX,IG)),
+     2   (/ LFLX /) )
+      ENDDO
+
+      ! CALL XABORT('SNFBH2: testing')
+
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(FLUX_G,FLUX,INDANG,TMPXNI,TMPXNJ,TMPXND,TMPMAT)
+      RETURN
+  400 FORMAT(16H SNFBH2: thread=,I8,12H --->(group=,I4,7H angle=,I4,1H))
+      END