From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001
From: stainer_t <thomas.stainer@oecd-nea.org>
Date: Mon, 8 Sep 2025 13:48:49 +0200
Subject: Initial commit from Polytechnique Montreal

---
 Trivac/src/NSSDFC.f | 485 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 485 insertions(+)
 create mode 100755 Trivac/src/NSSDFC.f

(limited to 'Trivac/src/NSSDFC.f')

diff --git a/Trivac/src/NSSDFC.f b/Trivac/src/NSSDFC.f
new file mode 100755
index 0000000..1c910be
--- /dev/null
+++ b/Trivac/src/NSSDFC.f
@@ -0,0 +1,485 @@
+*DECK NSSDFC
+      SUBROUTINE NSSDFC(IMPX,IDIM,NX,NY,NZ,NCODE,ICODE,ZCODE,MAT,XXX,
+     1 YYY,ZZZ,LL4F,LL4X,LL4Y,LL4Z,VOL,XX,YY,ZZ,IDL,KN,QFR,IQFR,MUX,
+     2 MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a coarse mesh finite difference (NEM
+* type) in a 3-D geometry.
+*
+*Copyright:
+* Copyright (C) 2022 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.
+* IDIM    number of Cartesian dimensions.
+* NX      number of elements along the X axis.
+* NY      number of elements along the Y axis.
+* NZ      number of elements along the Z axis.
+* NCODE   type of boundary condition applied on each side:
+*         I=1: X-; I=2: X+; I=3: Y-; I=4: Y+; I=5: Z-; I=6: Z+;
+*         NCODE(I)=1: VOID;  NCODE(I)=2: REFL;  NCODE(I)=4: TRAN;
+*         NCODE(I)=7: ZERO.
+* ICODE   physical albedo index on each side of the domain.
+* ZCODE   albedo corresponding to boundary condition 'VOID' on each
+*         side (ZCODE(i)=0.0 by default).
+* MAT     mixture index assigned to each element.
+* XXX     Cartesian coordinates along the X axis.
+* YYY     Cartesian coordinates along the Y axis.
+* ZZZ     Cartesian coordinates along the Z axis.
+* LL4F    total number of averaged flux unknown per energy group.
+*
+*Parameters: output
+* LL4X    total number of X-direccted interface net currents.
+* LL4Y    total number of Y-direccted interface net currents.
+* LL4Z    total number of Z-direccted interface net currents.
+* VOL     volume of each element.
+* XX      X-directed mesh spacings.
+* YY      Y-directed mesh spacings.
+* ZZ      Z-directed mesh spacings.
+* IDL     position of averaged fluxes in unknown vector.
+* KN      element-ordered interface net current unknown list.
+* QFR     element-ordered boundary conditions.
+* IQFR    element-ordered physical albedo indices.
+* MUX     X-oriented compressed storage mode indices.
+* MUY     Y-oriented compressed storage mode indices.
+* MUZ     Z-oriented compressed storage mode indices.
+* IMAX    X-oriented position of each first non-zero column element.
+* IMAY    Y-oriented position of each first non-zero column element.
+* IMAZ    Z-oriented position of each first non-zero column element.
+* IPY     Y-oriented permutation matrices.
+* IPZ     Z-oriented permutation matrices.
+*
+*-----------------------------------------------------------------------
+*
+      INTEGER IMPX,IDIM,NX,NY,NZ,NCODE(6),ICODE(6),MAT(NX,NY,NZ),LL4F,
+     1 LL4X,LL4Y,LL4Z,IDL(NX,NY,NZ),KN(6,NX,NY,NZ),IQFR(6,NX,NY,NZ),
+     2 MUX(LL4F),MUY(LL4F),MUZ(LL4F),IMAX(LL4F),IMAY(LL4F),IMAZ(LL4F),
+     3 IPY(LL4F),IPZ(LL4F)
+      REAL ZCODE(6),XXX(NX+1),YYY(NY+1),ZZZ(NZ+1),VOL(NX,NY,NZ),
+     1 XX(NX,NY,NZ),YY(NX,NY,NZ),ZZ(NX,NY,NZ),QFR(6,NX,NY,NZ)
+*----
+*  LOCAL VARIABLES
+*----
+      LOGICAL LL1,LALB
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: JPX,JPY,JPZ
+*
+      ALB(X)=0.5*(1.0-X)/(1.0+X)
+*----
+*  IDENTIFICATION OF THE NON VIRTUAL NODES
+*----
+      IF(IMPX.GT.0) WRITE(6,700) NX,NY,NZ
+      ALLOCATE(JPX((NX+1)*NY*NZ),JPY((NY+1)*NX*NZ),JPZ((NZ+1)*NX*NY))
+      JPX(:)=0
+      JPY(:)=0
+      JPZ(:)=0
+      IND=0
+      DO K0=1,NZ
+        DO K1=1,NY
+          DO K2=1,NX
+            IDL(K2,K1,K0)=0
+            KN(:6,K2,K1,K0)=0
+            IF(MAT(K2,K1,K0).EQ.0) CYCLE
+            IND=IND+1
+            IDL(K2,K1,K0)=IND
+            KN(1,K2,K1,K0)=K2    +(NX+1)*(K1-1)+(NX+1)*NY*(K0-1)
+            KN(2,K2,K1,K0)=(K2+1)+(NX+1)*(K1-1)+(NX+1)*NY*(K0-1)
+            KN(3,K2,K1,K0)=K1    +(NY+1)*(K0-1)+(NY+1)*NZ*(K2-1)
+            KN(4,K2,K1,K0)=(K1+1)+(NY+1)*(K0-1)+(NY+1)*NZ*(K2-1)
+            KN(5,K2,K1,K0)=K0    +(NZ+1)*(K2-1)+(NZ+1)*NX*(K1-1)
+            KN(6,K2,K1,K0)=(K0+1)+(NZ+1)*(K2-1)+(NZ+1)*NX*(K1-1)
+            JPX(KN(1:2,K2,K1,K0))=1
+            JPY(KN(3:4,K2,K1,K0))=1
+            JPZ(KN(5:6,K2,K1,K0))=1
+          ENDDO
+        ENDDO
+      ENDDO
+      IF(IND.NE.LL4F) CALL XABORT('NSSDFC: WRONG VALUE OF LL4F.')
+      LL4X=0
+      DO I=1,(NX+1)*NY*NZ
+        IF(JPX(I).EQ.1) THEN
+          LL4X=LL4X+1
+          JPX(I)=LL4X
+        ENDIF
+      ENDDO
+      LL4Y=0
+      DO I=1,(NY+1)*NX*NZ
+        IF(JPY(I).EQ.1) THEN
+          LL4Y=LL4Y+1
+          JPY(I)=LL4Y
+        ENDIF
+      ENDDO
+      LL4Z=0
+      DO I=1,(NZ+1)*NX*NY
+        IF(JPZ(I).EQ.1) THEN
+          LL4Z=LL4Z+1
+          JPZ(I)=LL4Z
+        ENDIF
+      ENDDO
+      DO K0=1,NZ
+        DO K1=1,NY
+          DO K2=1,NX
+            IF(MAT(K2,K1,K0).EQ.0) CYCLE
+            KN(1:2,K2,K1,K0)=JPX(KN(1:2,K2,K1,K0))
+            KN(3:4,K2,K1,K0)=JPY(KN(3:4,K2,K1,K0))
+            KN(5:6,K2,K1,K0)=JPZ(KN(5:6,K2,K1,K0))
+          ENDDO
+        ENDDO
+      ENDDO
+      DEALLOCATE(JPZ,JPY,JPX)
+*----
+*  IDENTIFICATION OF THE GEOMETRY. MAIN LOOP OVER THE NODES
+*----
+      QFR(:6,:NX,:NY,:NZ)=0.0
+      IQFR(:6,:NX,:NY,:NZ)=-99
+      DO K0=1,NZ
+        DO K1=1,NY
+          DO K2=1,NX
+            XX(K2,K1,K0)=0.0
+            YY(K2,K1,K0)=0.0
+            ZZ(K2,K1,K0)=0.0
+            VOL(K2,K1,K0)=0.0
+            IF(MAT(K2,K1,K0).LE.0) CYCLE
+            XX(K2,K1,K0)=XXX(K2+1)-XXX(K2)
+            YY(K2,K1,K0)=YYY(K1+1)-YYY(K1)
+            ZZ(K2,K1,K0)=ZZZ(K0+1)-ZZZ(K0)
+*----
+*  VOID, REFL OR ZERO BOUNDARY CONTITION
+*----
+            IQFR(:2,K2,K1,K0)=0
+            IF(K2.EQ.1) THEN
+              LL1=.TRUE.
+            ELSE
+              LL1=(MAT(K2-1,K1,K0).EQ.0)
+            ENDIF
+            IF(LL1) THEN
+              LALB=(NCODE(1).EQ.1).OR.(NCODE(1).EQ.6)
+              IF(LALB.AND.(ICODE(1).EQ.0)) THEN
+                QFR(1,K2,K1,K0)=ALB(ZCODE(1))
+                IQFR(1,K2,K1,K0)=-1
+              ELSE IF(LALB) THEN
+                QFR(1,K2,K1,K0)=1.0
+                IQFR(1,K2,K1,K0)=ICODE(1)
+              ELSE IF(NCODE(1).EQ.2) THEN
+                IQFR(1,K2,K1,K0)=-2
+              ELSE IF(NCODE(1).EQ.7) THEN
+                IQFR(1,K2,K1,K0)=-3
+              ELSE IF(NCODE(1).EQ.5) THEN
+                CALL XABORT('NSSDFC: SYME NOT IMPLEMENTED(1).')
+              ENDIF
+            ENDIF
+*
+            IF(K2.EQ.NX) THEN
+              LL1=.TRUE.
+            ELSE
+              LL1=(MAT(K2+1,K1,K0).EQ.0)
+            ENDIF
+            IF(LL1) THEN
+              LALB=(NCODE(2).EQ.1).OR.(NCODE(2).EQ.6)
+              IF(LALB.AND.(ICODE(2).EQ.0)) THEN
+                QFR(2,K2,K1,K0)=ALB(ZCODE(2))
+                IQFR(2,K2,K1,K0)=-1
+              ELSE IF(LALB) THEN
+                QFR(2,K2,K1,K0)=1.0
+                IQFR(2,K2,K1,K0)=ICODE(2)
+              ELSE IF(NCODE(2).EQ.2) THEN
+                IQFR(2,K2,K1,K0)=-2
+              ELSE IF(NCODE(2).EQ.7) THEN
+                IQFR(2,K2,K1,K0)=-3
+              ELSE IF(NCODE(1).EQ.5) THEN
+                CALL XABORT('NSSDFC: SYME NOT IMPLEMENTED(2).')
+              ENDIF
+            ENDIF
+*
+            IF(IDIM == 1) GO TO 100
+            IQFR(3:4,K2,K1,K0)=0
+            IF(K1.EQ.1) THEN
+              LL1=.TRUE.
+            ELSE
+              LL1=(MAT(K2,K1-1,K0).EQ.0)
+            ENDIF
+            IF(LL1) THEN
+              LALB=(NCODE(3).EQ.1).OR.(NCODE(3).EQ.6)
+              IF(LALB.AND.(ICODE(3).EQ.0)) THEN
+                QFR(3,K2,K1,K0)=ALB(ZCODE(3))
+                IQFR(3,K2,K1,K0)=-1
+              ELSE IF(LALB) THEN
+                QFR(3,K2,K1,K0)=1.0
+                IQFR(3,K2,K1,K0)=ICODE(3)
+              ELSE IF(NCODE(3).EQ.2) THEN
+                IQFR(3,K2,K1,K0)=-2
+              ELSE IF(NCODE(3).EQ.7) THEN
+                IQFR(3,K2,K1,K0)=-3
+              ELSE IF(NCODE(1).EQ.5) THEN
+                CALL XABORT('NSSDFC: SYME NOT IMPLEMENTED(3).')
+              ENDIF
+            ENDIF
+*
+            IF(K1.EQ.NY) THEN
+              LL1=.TRUE.
+            ELSE
+              LL1=(MAT(K2,K1+1,K0).EQ.0)
+            ENDIF
+            IF(LL1) THEN
+              LALB=(NCODE(4).EQ.1).OR.(NCODE(4).EQ.6)
+              IF(LALB.AND.(ICODE(4).EQ.0)) THEN
+                QFR(4,K2,K1,K0)=ALB(ZCODE(4))
+                IQFR(4,K2,K1,K0)=-1
+              ELSE IF(LALB) THEN
+                QFR(4,K2,K1,K0)=1.0
+                IQFR(4,K2,K1,K0)=ICODE(4)
+              ELSE IF(NCODE(4).EQ.2) THEN
+                IQFR(4,K2,K1,K0)=-2
+              ELSE IF(NCODE(4).EQ.7) THEN
+                IQFR(4,K2,K1,K0)=-3
+              ELSE IF(NCODE(1).EQ.5) THEN
+                CALL XABORT('NSSDFC: SYME NOT IMPLEMENTED(4).')
+              ENDIF
+            ENDIF
+*
+            IF(IDIM == 2) GO TO 100
+            IQFR(5:6,K2,K1,K0)=0
+            IF(K0.EQ.1) THEN
+              LL1=.TRUE.
+            ELSE
+              LL1=(MAT(K2,K1,K0-1).EQ.0)
+            ENDIF
+            IF(LL1) THEN
+              LALB=(NCODE(5).EQ.1).OR.(NCODE(5).EQ.6)
+              IF(LALB.AND.(ICODE(5).EQ.0)) THEN
+                QFR(5,K2,K1,K0)=ALB(ZCODE(5))
+                IQFR(5,K2,K1,K0)=-1
+              ELSE IF(LALB) THEN
+                QFR(5,K2,K1,K0)=1.0
+                IQFR(5,K2,K1,K0)=ICODE(5)
+              ELSE IF(NCODE(5).EQ.2) THEN
+                IQFR(5,K2,K1,K0)=-2
+              ELSE IF(NCODE(5).EQ.7) THEN
+                IQFR(5,K2,K1,K0)=-3
+              ELSE IF(NCODE(1).EQ.5) THEN
+                CALL XABORT('NSSDFC: SYME NOT IMPLEMENTED(5).')
+              ENDIF
+            ENDIF
+*
+            IF(K0.EQ.NZ) THEN
+              LL1=.TRUE.
+            ELSE
+              LL1=(MAT(K2,K1,K0+1).EQ.0)
+            ENDIF
+            IF(LL1) THEN
+              LALB=(NCODE(6).EQ.1).OR.(NCODE(6).EQ.6)
+              IF(LALB.AND.(ICODE(6).EQ.0)) THEN
+                QFR(6,K2,K1,K0)=ALB(ZCODE(6))
+                IQFR(6,K2,K1,K0)=-1
+              ELSE IF(LALB) THEN
+                QFR(6,K2,K1,K0)=1.0
+                IQFR(6,K2,K1,K0)=ICODE(6)
+              ELSE IF(NCODE(6).EQ.2) THEN
+                IQFR(6,K2,K1,K0)=-2
+              ELSE IF(NCODE(6).EQ.7) THEN
+                IQFR(6,K2,K1,K0)=-3
+              ELSE IF(NCODE(1).EQ.5) THEN
+                CALL XABORT('NSSDFC: SYME NOT IMPLEMENTED(6).')
+              ENDIF
+            ENDIF
+*----
+*  TRAN BOUNDARY CONDITION
+*----
+  100       IF((K2.EQ.1).AND.(NCODE(1).EQ.4)) THEN
+              KN(1,K2,K1,K0)=KN(2,NX,K1,K0)
+            ENDIF
+            IF((K2.EQ.NX).AND.(NCODE(2).EQ.4)) THEN
+              KN(2,K2,K1,K0)=KN(1,1,K1,K0)
+            ENDIF
+            IF((K1.EQ.1).AND.(NCODE(3).EQ.4)) THEN
+              KN(3,K2,K1,K0)=KN(2,K2,NY,K0)
+            ENDIF
+            IF((K1.EQ.NY).AND.(NCODE(4).EQ.4)) THEN
+              KN(4,K2,K1,K0)=KN(1,K2,1,K0)
+            ENDIF
+            IF((K0.EQ.1).AND.(NCODE(5).EQ.4)) THEN
+              KN(5,K2,K1,K0)=KN(6,K2,K1,NZ)
+            ENDIF
+            IF((K0.EQ.NZ).AND.(NCODE(6).EQ.4)) THEN
+              KN(6,K2,K1,K0)=KN(5,K2,K1,1)
+            ENDIF
+*
+            VOL(K2,K1,K0)=XX(K2,K1,K0)*YY(K2,K1,K0)*ZZ(K2,K1,K0)
+          ENDDO
+        ENDDO
+      ENDDO
+* END OF THE MAIN LOOP OVER NODES.
+*
+      IF(IMPX.GE.2) THEN
+         WRITE(6,720) VOL(:NX,:NY,:NZ)
+         WRITE(6,750)
+         DO K0=1,NZ
+           DO K1=1,NY
+             DO K2=1,NX
+               IF(MAT(K2,K1,K0).LE.0) CYCLE
+               KEL=(K0-1)*NX*NY+(K1-1)*NX+K2
+               WRITE (6,760) KEL,(KN(I,K2,K1,K0),I=1,6),
+     1         (QFR(I,K2,K1,K0),I=1,6),(IQFR(I,K2,K1,K0),I=1,6)
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDIF
+*----
+*  COMPUTE THE PERMUTATION VECTORS IPY AND IPZ
+*----
+      IF(IDIM.GE.2) THEN
+         INX1=0
+         DO K2=1,NX
+           DO K0=1,NZ
+             DO K1=1,NY
+               INX2=IDL(K2,K1,K0)
+               IF(INX2.LE.0) CYCLE
+               INX1=INX1+1
+               IPY(INX2)=INX1
+             ENDDO
+           ENDDO
+         ENDDO
+         IF(INX1.NE.IND) CALL XABORT('NSSDFC: FAILURE OF THE RENUMBERI'
+     1   //'NG ALGORITHM(1)')
+         IF(IDIM.EQ.3) THEN
+            INX1=0
+            DO K1=1,NY
+              DO K2=1,NX
+                DO K0=1,NZ
+                  INX2=IDL(K2,K1,K0)
+                  IF(INX2.LE.0) CYCLE
+                  INX1=INX1+1
+                  IPZ(INX2)=INX1
+                ENDDO
+              ENDDO
+            ENDDO
+            IF(INX1.NE.IND) CALL XABORT('NSSDFC: FAILURE OF THE RENUMB'
+     1      //'ERING ALGORITHM(2)')
+         ENDIF
+      ENDIF
+*----
+*  COMPUTE VECTOR MUX
+*----
+      MUX(:LL4F)=1
+      DO K0=1,NZ
+        DO K1=1,NY
+*         X- SIDE:
+          DO K2=2,NX
+            KEL=IDL(K2,K1,K0)
+            IF(KEL.EQ.0) CYCLE
+            KK1=IDL(K2-1,K1,K0)
+            IF(KK1.GT.0) MUX(KEL)=MAX0(MUX(KEL),KEL-KK1+1)
+          ENDDO
+*         X+ SIDE:
+          DO K2=1,NX-1
+            KEL=IDL(K2,K1,K0)
+            IF(KEL.EQ.0) CYCLE
+            KK2=IDL(K2+1,K1,K0)
+            IF(KK2.GT.0) MUX(KEL)=MAX0(MUX(KEL),KEL-KK2+1)
+          ENDDO
+        ENDDO
+      ENDDO
+*----
+*  COMPUTE VECTOR MUY
+*----
+      IF(IDIM.GE.2) THEN
+        MUY(:LL4F)=1
+        DO K2=1,NX
+          DO K0=1,NZ
+*           Y- SIDE:
+            DO K1=2,NY
+              KEL=IDL(K2,K1,K0)
+              IF(KEL.EQ.0) CYCLE
+              INY1=IPY(KEL)
+              KK3=IDL(K2,K1-1,K0)
+              IF(KK3.GT.0) MUY(INY1)=MAX0(MUY(INY1),INY1-IPY(KK3)+1)
+            ENDDO
+*           Y- SIDE:
+            DO K1=1,NY-1
+              KEL=IDL(K2,K1,K0)
+              IF(KEL.EQ.0) CYCLE
+              INY1=IPY(KEL)
+              KK4=IDL(K2,K1+1,K0)
+              IF(KK4.GT.0) MUY(INY1)=MAX0(MUY(INY1),INY1-IPY(KK4)+1)
+            ENDDO
+          ENDDO
+        ENDDO
+      ELSE
+        MUY(:LL4F)=0
+      ENDIF
+*----
+*  COMPUTE VECTOR MUZ
+*----
+      IF(IDIM.EQ.3) THEN
+        MUZ(:LL4F)=1
+        DO K1=1,NY
+          DO K2=1,NX
+*           Z- SIDE:
+            DO K0=2,NZ
+              KEL=IDL(K2,K1,K0)
+              IF(KEL.EQ.0) CYCLE
+              INZ1=IPZ(KEL)
+              KK5=IDL(K2,K1,K0-1)
+              IF(KK5.GT.0) MUZ(INZ1)=MAX0(MUZ(INZ1),INZ1-IPZ(KK5)+1)
+            ENDDO
+*           Z+ SIDE:
+            DO K0=1,NZ-1
+              KEL=IDL(K2,K1,K0)
+              IF(KEL.EQ.0) CYCLE
+              INZ1=IPZ(KEL)
+              KK6=IDL(K2,K1,K0+1)
+              IF(KK6.GT.0) MUZ(INZ1)=MAX0(MUZ(INZ1),INZ1-IPZ(KK6)+1)
+            ENDDO
+          ENDDO
+        ENDDO
+      ELSE
+        MUZ(:LL4F)=0
+      ENDIF
+*
+      MUXMAX=0
+      MUYMAX=0
+      MUZMAX=0
+      IIMAXX=0
+      IIMAXY=0
+      IIMAXZ=0
+      DO I=1,LL4F
+        MUXMAX=MAX(MUXMAX,MUX(I))
+        MUYMAX=MAX(MUYMAX,MUY(I))
+        MUZMAX=MAX(MUZMAX,MUZ(I))
+        IBAND=MUX(I)
+        IIMAXX=IIMAXX+IBAND
+        MUX(I)=IIMAXX
+        IIMAXX=IIMAXX+IBAND-1
+        IMAX(I)=IIMAXX
+        IBAND=MUY(I)
+        IIMAXY=IIMAXY+IBAND
+        MUY(I)=IIMAXY
+        IIMAXY=IIMAXY+IBAND-1
+        IMAY(I)=IIMAXY
+        IBAND=MUZ(I)
+        IIMAXZ=IIMAXZ+IBAND
+        MUZ(I)=IIMAXZ
+        IIMAXZ=IIMAXZ+IBAND-1
+        IMAZ(I)=IIMAXZ
+      ENDDO
+      IF(IMPX.GT.0) WRITE (6,770) MUXMAX,MUYMAX,MUZMAX
+      RETURN
+*
+  700 FORMAT(/46H NSSDFC: COARSE MESH FINITE DIFFERENCE METHOD.//3H NU,
+     1 28HMBER OF NODES ALONG X AXIS =,I3/17X,14HALONG Y AXIS =,I3/
+     2 17X,14HALONG Z AXIS =,I3)
+  720 FORMAT(/17H VOLUMES PER NODE/(1X,1P,10E13.4))
+  750 FORMAT(/22H NUMBERING OF UNKNOWNS/1X,21(1H-)//4X,4HNODE,5X,3HINT,
+     1 26HERFACE NET CURRENT INDICES,28X,23HVOID BOUNDARY CONDITION)
+  760 FORMAT(1X,I6,7X,6I8,6X,6F9.2/68X,6I9)
+  770 FORMAT(/41H NSSDFC: MAXIMUM BANDWIDTH ALONG X AXIS =,I5/
+     1 27X,14HALONG Y AXIS =,I5/27X,14HALONG Z AXIS =,I5)
+      END
-- 
cgit v1.2.3