diff options
Diffstat (limited to 'Trivac/src/NSSCO.f')
| -rwxr-xr-x | Trivac/src/NSSCO.f | 149 |
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 |