Initial commit from Polytechnique Montreal

1 files changed, 149 insertions, 0 deletions

diff --git a/Trivac/src/NSSCO.f b/Trivac/src/NSSCO.f
new file mode 100755
index 0000000..de8ff40
--- /dev/null
+++ b/Trivac/src/NSSCO.f
@@ -0,0 +1,149 @@
+*DECK NSSCO
+      SUBROUTINE NSSCO(NX,NY,NZ,NMIX,I,J,K,MAT,XX,YY,ZZ,DIFF,IQFR,QFR,
+     1 COEF)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the mesh centered finite difference coefficients.
+*
+*Copyright:
+* Copyright (C) 2023 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
+* NX      number of X-directed nodes.
+* NY      number of Y-directed nodes.
+* NZ      number of Z-directed nodes.
+* NMIX    number of material mixtures.
+* I       X-index of node under consideration.
+* J       Y-index of node under consideration.
+* K       Z-index of node under consideration.
+* MAT     mixture index assigned to each node.
+* DIF     diffusion coefficients.
+* XX      X-directed mesh spacings.
+* YY      Y-directed mesh spacings.
+* ZZ      Z-directed mesh spacings.
+* IQFR    node-ordered unknown list:
+*         =0: neighbour exists;
+*         =-1: void/albedo boundary condition;
+*         =-2: reflection boundary condition;
+*         =-3: ZERO flux boundary condition;
+*         =-4: SYME boundary condition (axial symmetry).
+* QFR     element-ordered boundary conditions.
+*
+*Parameters: output
+* COEF    mesh centered finite difference coefficients.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NX,NY,NZ,I,J,K,MAT(NX,NY,NZ),IQFR(6)
+      REAL DIFF(NMIX),XX(NX,NY,NZ),YY(NX,NY,NZ),ZZ(NX,NY,NZ),QFR(6),
+     1 COEF(6)
+*----
+*  LOCAL VARIABLES
+*----
+      DHARM(X1,X2,DIF1,DIF2)=2.0*DIF1*DIF2/(X1*DIF2+X2*DIF1)
+*
+      IBM=MAT(I,J,K)
+      DX=XX(I,J,K)
+      DY=YY(I,J,K)
+      DZ=ZZ(I,J,K)
+      KK1=IQFR(1)
+      KK2=IQFR(2)
+      KK3=IQFR(3)
+      KK4=IQFR(4)
+      KK5=IQFR(5)
+      KK6=IQFR(6)
+      COEF(:6)=0.
+      ! x- side:
+      IF(KK1 == 0) THEN
+        COEF(1)=DHARM(DX,XX(I-1,J,K),DIFF(IBM),DIFF(MAT(I-1,J,K)))
+      ELSE IF((KK1 > 0).OR.(KK1 == -1)) THEN
+        COEF(1)=DHARM(DX,DX,DIFF(IBM),DX*QFR(1)/2.0)
+      ELSE IF(KK1 == -2) THEN
+        COEF(1)=0.0
+      ELSE IF(KK1 == -3) THEN
+        COEF(1)=2.0*DHARM(DX,DX,DIFF(IBM),DIFF(IBM))
+      ENDIF
+      ! x+ side:
+      IF(KK2 == 0) THEN
+        COEF(2)=DHARM(DX,XX(I+1,J,K),DIFF(IBM),DIFF(MAT(I+1,J,K)))
+      ELSE IF((KK2 > 0).OR.(KK2 == -1)) THEN
+        COEF(2)=DHARM(DX,DX,DIFF(IBM),DX*QFR(2)/2.0)
+      ELSE IF(KK2 == -2) THEN
+        COEF(2)=0.0
+      ELSE IF(KK2 == -3) THEN
+        COEF(2)=2.0*DHARM(DX,DX,DIFF(IBM),DIFF(IBM))
+      ELSE IF(KK2 == -4) THEN
+        IF(KK1 == -4) CALL XABORT('NSSCO: INCONSISTENT SYME (1).')
+        COEF(2)=COEF(1)
+      ENDIF
+      IF(KK1 == -4) THEN
+        IF(KK2 == -4) CALL XABORT('NSSCO: INCONSISTENT SYME (2).')
+        COEF(1)=COEF(2)
+      ENDIF
+      ! y- side:
+      IF(KK3 == 0) THEN
+        COEF(3)=DHARM(DY,YY(I,J-1,K),DIFF(IBM),DIFF(MAT(I,J-1,K)))
+      ELSE IF((KK3 > 0).OR.(KK3 == -1)) THEN
+        COEF(3)=DHARM(DY,DY,DIFF(IBM),DY*QFR(3)/2.0)
+      ELSE IF(KK3 == -2) THEN
+        COEF(3)=0.0
+      ELSE IF(KK3 == -3) THEN
+        COEF(3)=2.0*DHARM(DY,DY,DIFF(IBM),DIFF(IBM))
+      ENDIF
+      ! y+ side:
+      IF(KK4 == 0) THEN
+        COEF(4)=DHARM(DY,YY(I,J+1,K),DIFF(IBM),DIFF(MAT(I,J+1,K)))
+      ELSE IF((KK4 > 0).OR.(KK4 == -1)) THEN
+        COEF(4)=DHARM(DY,DY,DIFF(IBM),DY*QFR(4)/2.0)
+      ELSE IF(KK4 == -2) THEN
+        COEF(4)=0.0
+      ELSE IF(KK4 == -3) THEN
+        COEF(4)=2.0*DHARM(DY,DY,DIFF(IBM),DIFF(IBM))
+      ELSE IF(KK4 == -4) THEN
+        IF(KK3 == -4) CALL XABORT('NSSCO: INCONSISTENT SYME (3).')
+        COEF(4)=COEF(3)
+      ENDIF
+      IF(KK3 == -4) THEN
+        IF(KK4 == -4) CALL XABORT('NSSCO: INCONSISTENT SYME (4).')
+        COEF(3)=COEF(4)
+      ENDIF
+      ! z- side:
+      IF(KK5 == 0) THEN
+        COEF(5)=DHARM(DZ,ZZ(I,J,K-1),DIFF(IBM),DIFF(MAT(I,J,K-1)))
+      ELSE IF((KK5 > 0).OR.(KK5 == -1)) THEN
+        COEF(5)=DHARM(DZ,DZ,DIFF(IBM),DZ*QFR(5)/2.0)
+      ELSE IF(KK5 == -2) THEN
+        COEF(5)=0.0
+      ELSE IF(KK5 == -3) THEN
+        COEF(5)=2.0*DHARM(DZ,DZ,DIFF(IBM),DIFF(IBM))
+      ENDIF
+      ! z+ side:
+      IF(KK6 == 0) THEN
+        COEF(6)=DHARM(DZ,ZZ(I,J,K+1),DIFF(IBM),DIFF(MAT(I,J,K+1)))
+      ELSE IF((KK6 > 0).OR.(KK6 == -1)) THEN
+        COEF(6)=DHARM(DZ,DZ,DIFF(IBM),DZ*QFR(6)/2.0)
+      ELSE IF(KK6 == -2) THEN
+        COEF(6)=0.0
+      ELSE IF(KK6 == -3) THEN
+        COEF(6)=2.0*DHARM(DZ,DZ,DIFF(IBM),DIFF(IBM))
+      ELSE IF(KK6 == -4) THEN
+        IF(KK5 == -4) CALL XABORT('NSSCO: INCONSISTENT SYME (5).')
+        COEF(6)=COEF(5)
+      ENDIF
+      IF(KK5 == -4) THEN
+        IF(KK6 == -4) CALL XABORT('NSSCO: INCONSISTENT SYME (6).')
+        COEF(5)=COEF(6)
+      ENDIF
+      RETURN
+      END