From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/SNTT3D.f | 609 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 609 insertions(+) create mode 100644 Dragon/src/SNTT3D.f (limited to 'Dragon/src/SNTT3D.f') diff --git a/Dragon/src/SNTT3D.f b/Dragon/src/SNTT3D.f new file mode 100644 index 0000000..3162d44 --- /dev/null +++ b/Dragon/src/SNTT3D.f @@ -0,0 +1,609 @@ +*DECK SNTT3D + SUBROUTINE SNTT3D(IGE,IMPX,LX,LY,LZ,SIDE,IELEM,NLF,NPQ,NSCT,IQUAD, + 1 NCODE,ZCODE,MAT,XXX,YYY,ZZZ,VOL,IDL,DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC, + 2 DB,DA,PL,LL4,NUN,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM,IGLK,MN,DN,IL, + 3 IM,ISCAT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Numbering corresponding to a 3-D Cartesian with discrete ordinates +* approximation of the flux. +* +*Copyright: +* Copyright (C) 2007 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): N. Martin and C. Bienvenue +* +*Parameters: input +* IGE geometry type with (=0 for Cartesian, =2 for hexagonal). +* IMPX print parameter. +* LX number of elements along the X axis; +* OR, number of elements in X-Y plane in hexagonal geometry, +* including lozenges and further submeshing of lozenges +* e.g. domain of 7 hex. by 5 levels with a submeshing of 2 +* will have LX=7*3*2*2=84 +* LY number of elements along the Y axis. +* LZ number of elements along the Z axis. +* SIDE side of hexagon. +* IELEM measure of order of the spatial approximation polynomial: +* =1 constant - only for HODD, classical diamond scheme +* (default for HODD); +* =2 linear - default for DG; +* =3 parabolic; +* =4 cubic - only for DG. +* NLF SN order for the flux (even number). +* NPQ number of SN directions in eight octants. +* NSCT maximum number of spherical harmonics moments of the flux. +* IQUAD type of SN quadrature (1 Level symmetric, type IQUAD; +* 4 Legendre-Chebyshev; 5 symmetric Legendre-Chebyshev; +* 6 quadruple range). +* 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+): +* =1 VOID; =2 REFL; =4 TRAN. +* ZCODE ZCODE(I) is the 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. +* EELEM measure of order of the energy approximation polynomial: +* =1 constant - default for HODD; +* =2 linear - default for DG; +* >3 higher orders. +* IBFP type of energy proparation relation: +* =0 no Fokker-Planck term; +* =1 Galerkin type; +* =2 heuristic Przybylski and Ligou type. +* ISCHM method of spatial discretisation: +* =1 High-Order Diamond Differencing (HODD) - default; +* =2 Discontinuous Galerkin finite element method (DG); +* =3 Adaptive weighted method (AWD). +* ESCHM method of energy discretisation: +* =1 High-Order Diamond Differencing (HODD) - default; +* =2 Discontinuous Galerkin finite element method (DG); +* =3 Adaptive weighted method (AWD). +* IGLK angular interpolation type: +* =0 classical SN method. +* =1 Galerkin quadrature method (M = inv(D)) +* =2 Galerkin quadrature method (D = inv(M)) +* ISCAT maximum number of spherical harmonics moments of the flux. +* +*Parameters: output +* VOL volume of each element. +* IDL isotropic flux indices. +* DU first direction cosines ($\\mu$). +* DE second direction cosines ($\\eta$). +* DZ third direction cosines ($\\xi$). +* W weights. +* MRMX quadrature index. +* MRMY quadrature index. +* MRMZ quadrature index. +* DC diamond-scheme parameter. +* DB diamond-scheme parameter. +* DA diamond-scheme parameter. +* PL discrete values of the spherical harmonics corresponding +* to the 3D SN quadrature. +* LL4 number of unknowns being solved for, over the domain. This +* includes the various moments of the isotropic (and if present, +* anisotropic) flux. +* NUN total number of unknowns stored in the FLUX vector per group. +* This includes LL4 (see above) as well as any surface boundary +* fluxes, if present. +* WX spatial closure relation weighting factors. +* WE energy closure relation weighting factors. +* CST constants for the polynomial approximations. +* MN moment-to-discrete matrix. +* DN discrete-to-moment matrix. +* IL indexes (l) of each spherical harmonics in the +* interpolation basis. +* IM indexes (m) of each spherical harmonics in the +* interpolation basis. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IMPX,LX,LY,LZ,IELEM,NLF,NPQ,NSCT,IQUAD,NCODE(6), + 1 MAT(LX,LY,LZ),IDL(LX*LY*LZ),MRMX(NPQ),MRMY(NPQ),MRMZ(NPQ),EELEM, + 2 IBFP,ISCHM,ESCHM,IL(NSCT),IM(NSCT),ISCAT,IGLK + REAL ZCODE(6),VOL(LX,LY,LZ),XXX(LX+1),YYY(LY+1),ZZZ(LZ+1), + 1 DU(NPQ),DE(NPQ),DZ(NPQ),W(NPQ),DC(LX,LY,NPQ),DB(LX,LZ,NPQ), + 2 DA(LY,LZ,NPQ),PL(NSCT,NPQ),WX(IELEM+1), + 3 WE(EELEM+1),CST(MAX(IELEM,EELEM)),MN(NPQ,NSCT),DN(NSCT,NPQ) +*---- +* LOCAL VARIABLES +*---- + CHARACTER HSMG*131 + LOGICAL L1,L2,L3,L4,L5,L6 + PARAMETER(RLOG=1.0E-8,PI=3.141592654) + REAL PX,PE + DOUBLE PRECISION MND(NPQ,NPQ),NORM,IPROD + INTEGER, ALLOCATABLE, DIMENSION(:) :: JOP + REAL, ALLOCATABLE, DIMENSION(:) :: XX,YY,ZZ,UU,WW,TPQ,UPQ,VPQ,WPQ + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: V,V2 + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: U + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: RLM +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(XX(LX),YY(LY),ZZ(LZ)) +*---- +* UNFOLD HEIGHT-OCTANT QUADRATURES. +*---- + IF(MOD(NLF,2).EQ.1) CALL XABORT('SNTT3D: EVEN NLF EXPECTED.') + IF(IQUAD.EQ.10) THEN + NPQ0=NLF**2/4 + ELSE + NPQ0=NLF*(NLF/2+1)/4 + ENDIF + ALLOCATE(JOP(NLF/2),UU(NLF/2),WW(NLF/2),TPQ(NPQ0),UPQ(NPQ0), + 1 VPQ(NPQ0),WPQ(NPQ0)) + IF(IQUAD.EQ.1) THEN +!Level-symmetric quadrature of type 1 + CALL SNQU01(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ) + ELSE IF(IQUAD.EQ.2) THEN +!Level-symmetric quadrature of type 2 + CALL SNQU02(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ) + ELSE IF(IQUAD.EQ.3) THEN +!Snow Level-symmetric type quadrature + CALL SNQU03(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ) + ELSE IF(IQUAD.EQ.4) THEN +!Legendre-Chebyshev quadrature + CALL SNQU04(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ) + ELSE IF(IQUAD.EQ.5) THEN +!Symmetric Legendre-Chebyshev quadrature + CALL SNQU05(NLF,TPQ,UPQ,VPQ,WPQ) + ELSE IF(IQUAD.EQ.6) THEN +!Quadruple Range quadrature + CALL SNQU06(NLF,TPQ,UPQ,VPQ,WPQ) + ELSE IF(IQUAD.EQ.10) THEN + CALL SNQU10(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ) + ELSE + CALL XABORT('SNTT3D: UNKNOWN QUADRATURE TYPE.') + ENDIF + N=0 + IOF=0 + DO 320 I=1,NLF/2 + JOF=IOF+NLF-2*I+2 + DO 330 J=0,NLF/2-I + IOF=IOF+1 + KOF=IOF+(NLF+2)*NLF/4 + LOF=IOF+(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=-UPQ(N+J+1) + DE(IOF)=-VPQ(N+J+1) + DZ(IOF)=-TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 330 CONTINUE + DO 340 J=NLF/2-I,0,-1 + IOF=IOF+1 + KOF=IOF+(NLF+2)*NLF/4 + LOF=IOF+(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=UPQ(N+J+1) + DE(IOF)=-VPQ(N+J+1) + DZ(IOF)=-TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 340 CONTINUE + N=N+NLF/2-I+1 + 320 CONTINUE + N=0 + DO 350 I=1,NLF/2 + JOF=IOF+NLF-2*I+2 + DO 360 J=0,NLF/2-I + IOF=IOF+1 + KOF=IOF-(NLF+2)*NLF/4 + LOF=IOF+(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=-UPQ(N+J+1) + DE(IOF)=VPQ(N+J+1) + DZ(IOF)=-TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 360 CONTINUE + DO 370 J=NLF/2-I,0,-1 + IOF=IOF+1 + KOF=IOF-(NLF+2)*NLF/4 + LOF=IOF+(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=UPQ(N+J+1) + DE(IOF)=VPQ(N+J+1) + DZ(IOF)=-TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 370 CONTINUE + N=N+NLF/2-I+1 + 350 CONTINUE + N=0 + DO 380 I=1,NLF/2 + JOF=IOF+NLF-2*I+2 + DO 390 J=0,NLF/2-I + IOF=IOF+1 + KOF=IOF+(NLF+2)*NLF/4 + LOF=IOF-(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=-UPQ(N+J+1) + DE(IOF)=-VPQ(N+J+1) + DZ(IOF)=TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 390 CONTINUE + DO 400 J=NLF/2-I,0,-1 + IOF=IOF+1 + KOF=IOF+(NLF+2)*NLF/4 + LOF=IOF-(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=UPQ(N+J+1) + DE(IOF)=-VPQ(N+J+1) + DZ(IOF)=TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 400 CONTINUE + N=N+NLF/2-I+1 + 380 CONTINUE + N=0 + DO 410 I=1,NLF/2 + JOF=IOF+NLF-2*I+2 + DO 420 J=0,NLF/2-I + IOF=IOF+1 + KOF=IOF-(NLF+2)*NLF/4 + LOF=IOF-(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=-UPQ(N+J+1) + DE(IOF)=VPQ(N+J+1) + DZ(IOF)=TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 420 CONTINUE + DO 430 J=NLF/2-I,0,-1 + IOF=IOF+1 + KOF=IOF-(NLF+2)*NLF/4 + LOF=IOF-(NLF+2)*NLF/2 + MRMX(IOF)=JOF + MRMY(IOF)=KOF + MRMZ(IOF)=LOF + DU(IOF)=UPQ(N+J+1) + DE(IOF)=VPQ(N+J+1) + DZ(IOF)=TPQ(N+J+1) + W(IOF)=WPQ(N+J+1) + JOF=JOF-1 + 430 CONTINUE + N=N+NLF/2-I+1 + 410 CONTINUE + DEALLOCATE(WPQ,VPQ,UPQ,TPQ,WW,UU,JOP) + IF(IMPX.GE.4) THEN + WRITE(6,'(/41H SNTT3D:HEIGHT-OCTANT ANGULAR QUADRATURES:/26X, + 1 2HMU,9X,3HETA,10X,2HXI,6X,6HWEIGHT)') + SUM=0.0 + DO 70 N=1,NPQ + SUM=SUM+W(N) + WRITE(6,'(1X,4I5,1P,4E12.4)') N,MRMX(N),MRMY(N),MRMZ(N),DU(N), + 1 DE(N),DZ(N),W(N) + 70 CONTINUE + WRITE(6,'(54X,10(1H-)/52X,1P,E12.4)') SUM + ENDIF +*---- +* IDENTIFICATION OF THE GEOMETRY. +*---- + IF(IGE.EQ.0) THEN +* ---------- +* 3D CARTESIAN +* ---------- + DO 83 N=1,NPQ + VU=DU(N) + VE=DE(N) + VZ=DZ(N) + DO 82 I=1,LX + DO 81 J=1,LY + DO 80 K=1,LZ + XX(I)=XXX(I+1)-XXX(I) + YY(J)=YYY(J+1)-YYY(J) + ZZ(K)=ZZZ(K+1)-ZZZ(K) + DA(J,K,N)=VU*YY(J)*ZZ(K) + DB(I,K,N)=VE*XX(I)*ZZ(K) + DC(I,J,N)=VZ*XX(I)*YY(J) + VOL(I,J,K)=XX(I)*YY(J)*ZZ(K) + 80 CONTINUE + 81 CONTINUE + 82 CONTINUE + 83 CONTINUE + ELSEIF(IGE.EQ.2) THEN +* ---------- +* 3D HEXAGONAL +* ---------- + DET = SQRT(3.0)*(SIDE**2)/2.0 + DO 93 N=1,NPQ + VU=DU(N) + VE=DE(N) + VZ=DZ(N) + DO 92 K=1,LZ + DO 91 J=1,LY + DO 90 I=1,LX + ZZ(K)=ZZZ(K+1)-ZZZ(K) + DA(J,K,N)=VU*ZZ(K) + DB(I,K,N)=VE*ZZ(K) + DC(I,J,N)=VZ*DET + VOL(I,J,K)=DET*ZZ(K) + 90 CONTINUE + 91 CONTINUE + 92 CONTINUE + 93 CONTINUE + ENDIF +*---- +* GENERATE SPHERICAL HARMONICS FOR SCATTERING SOURCE. +*---- + IOF=0 + DO 211 L=0,ISCAT-1 + DO 210 M=-L,L + IOF=IOF+1 + IF(IOF.GT.ISCAT) GO TO 211 + DO 200 N=1,NPQ + PL(IOF,N)=PNSH(L,M,DU(N),DE(N),DZ(N)) + 200 CONTINUE + 210 CONTINUE + 211 CONTINUE +*---- +* GENERATE MAPPING MATRIX FOR GALERKIN QUADRATURE METHOD +*---- + MN(:NPQ,:NSCT)=0.0 + DN(:NSCT,:NPQ)=0.0 + IL(:NSCT)=0 + IM(:NSCT)=0 + IF(IGLK.NE.0) THEN + ALLOCATE(U(NPQ,NPQ),RLM(NPQ,ISCAT,2*ISCAT-1),V(NPQ),V2(NPQ)) + RLM(:NPQ,:ISCAT-1,:2*ISCAT-1)=0.0 + DO L=0,ISCAT-1 + DO M=-L,L + DO N=1,NPQ + RLM(N,L+1,M+L+1)=PNSH(L,M,DU(N),DE(N),DZ(N)) + ENDDO + ENDDO + ENDDO + ! GRAM-SCHMIDT PROCEDURE TO FIND INDEPENDANT SET + ! OF SPHERICAL HARMONICS WITH ANY QUADRATURE + U(:NPQ,:NPQ)=0.0D0 + NORM=0.0D0 + DO N=1,NPQ + NORM=NORM+RLM(N,1,1)**2 + ENDDO + NORM=SQRT(NORM) + DO N=1,NPQ + IF(IGLK.EQ.1) THEN + MND(1,N)=W(N)*RLM(N,1,1) + ELSEIF(IGLK.EQ.2) THEN + MND(N,1)=(2.0*L+1.0)/(4.0*PI)*RLM(N,1,1) + ELSE + CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.') + ENDIF + U(N,1)=RLM(N,1,1)/NORM + ENDDO + IND=1 + ! ITERATE OVER THE SPHERICAL HARMONICS + DO 212 L=1,ISCAT-1 + DO 213 M=-L,L + V2(:NPQ)=0.0D0 + DO N=1,IND + IPROD=0.0D0 + DO N2=1,NPQ + IPROD=IPROD+U(N2,N)*RLM(N2,L+1,M+L+1) + ENDDO + DO N2=1,NPQ + V2(N2)=V2(N2)+IPROD*U(N2,N) + ENDDO + ENDDO + V(:NPQ)=0.0D0 + DO N=1,NPQ + V(N)=RLM(N,L+1,M+L+1)-V2(N) + ENDDO + NORM=0.0D0 + DO N=1,NPQ + NORM=NORM+V(N)**2 + ENDDO + NORM=SQRT(NORM) + ! KEEP THE SPHERICAL HARMONICS IF IT IS INDEPENDANT + IF(NORM.GE.1.0E-5) THEN + IND=IND+1 + DO N=1,NPQ + U(N,IND)=V(N)/NORM + IF(IGLK.EQ.1) THEN + MND(IND,N)=W(N)*RLM(N,L+1,M+L+1) + ELSEIF(IGLK.EQ.2) THEN + MND(N,IND)=(2.0*L+1.0)/(4.0*PI)*RLM(N,L+1,M+L+1) + ELSE + CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.') + ENDIF + ENDDO + IL(IND)=L + IM(IND)=M + ENDIF + IF(IND.EQ.NPQ) GOTO 217 + 213 ENDDO + 212 ENDDO + CALL XABORT('SNTT3D: THE'// + 1 ' GRAM-SCHMIDTH PROCEDURE TO FIND A SUITABLE INTERPOLATION'// + 2 ' BASIS REQUIRE HIGHER LEGENDRE ORDER.') + ! FIND INVERSE MATRIX + 217 IF(IGLK.EQ.1) THEN + DN=REAL(MND) + CALL ALINVD(NPQ,MND,NPQ,IER) + IF(IER.NE.0) CALL XABORT('SNTT3D: SINGULAR MATRIX.') + MN=REAL(MND) + ELSEIF(IGLK.EQ.2) THEN + MN=REAL(MND) + CALL ALINVD(NPQ,MND,NPQ,IER) + IF(IER.NE.0) CALL XABORT('SNTT3D: SINGULAR MATRIX.') + DN=REAL(MND) + ELSE + CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.') + ENDIF + DEALLOCATE(U,RLM,V,V2) + ELSE + IND=1 + DO L=0,ISCAT-1 + DO 218 M=-L,L + IL(IND)=L + IM(IND)=M + DO N=1,NPQ + DN(IND,N)=W(N)*PNSH(L,M,DU(N),DE(N),DZ(N)) + MN(N,IND)=(2.0*L+1.0)/(4.0*PI) + 1 *PNSH(L,M,DU(N),DE(N),DZ(N)) + ENDDO + IND=IND+1 + 218 ENDDO + ENDDO + ENDIF +*---- +* GENERATE THE WEIGHTING PARAMETERS OF THE CLOSURE RELATION. +*---- + PX=1 + PE=1 + IF(ISCHM.EQ.1.OR.ISCHM.EQ.3) THEN + PX=1 + ELSEIF(ISCHM.EQ.2) THEN + PX=0 + ELSE + CALL XABORT('SNTT3D: UNKNOWN TYPE OF SPATIAL CLOSURE RELATION.') + ENDIF + IF(MOD(IELEM,2).EQ.1) THEN + WX(1)=-PX + WX(2:IELEM+1:2)=1+PX + IF(IELEM.GE.2) WX(3:IELEM+1:2)=1-PX + ELSE + WX(1)=PX + WX(2:IELEM+1:2)=1-PX + IF(IELEM.GE.2) WX(3:IELEM+1:2)=1+PX + ENDIF + IF(IBFP.NE.0) THEN + IF(ESCHM.EQ.1.OR.ESCHM.EQ.3) THEN + PE=1 + ELSEIF(ESCHM.EQ.2) THEN + PE=0 + ELSE + CALL XABORT('SNTT3D: UNKNOWN TYPE OF ENERGY CLOSURE RELATION.') + ENDIF + IF(MOD(EELEM,2).EQ.1) THEN + WE(1)=-PE + WE(2:EELEM+1:2)=1+PE + IF(EELEM.GE.2) WE(3:EELEM+1:2)=1-PE + ELSE + WE(1)=PE + WE(2:EELEM+1:2)=1-PE + IF(EELEM.GE.2) WE(3:EELEM+1:2)=1+PE + ENDIF + ENDIF + ! NORMALIZED LEGENDRE POLYNOMIAL CONSTANTS + DO IEL=1,MAX(IELEM,EELEM) + CST(IEL)=SQRT(2.0*IEL-1.0) + ENDDO +*---- +* COMPUTE ISOTROPIC FLUX INDICES. +*---- + NM=IELEM**3*EELEM + NMX=IELEM**2*EELEM + NMY=NMX + NMZ=NMX + NME=IELEM**3 + LL4=LX*LY*LZ*NSCT*NM + IF(IGE.EQ.0)THEN + NUN=LL4+(NMX*LY*LZ+NMY*LX*LZ+NMZ*LX*LY)*NPQ + ELSEIF(IGE.EQ.2)THEN + IF((NCODE(5)==1).and.(NCODE(5)==NCODE(6)))THEN + NUN=LL4 + ELSE + NUN=LL4+(LX*LY*LZ)*NMZ*NPQ + ENDIF + ELSE + CALL XABORT('SNTT3D: CHECK SPATIAL SCHEME DISCRETISATION '// + 1 'PARAMETER.') + ENDIF + DO I=1,LX*LY*LZ + IDL(I)=(I-1)*NSCT*NM+1 + ENDDO +*---- +* SET BOUNDARY CONDITIONS. +*---- + DO 240 I=1,6 + IF(NCODE(I).NE.1) ZCODE(I)=1.0 + IF(NCODE(I).EQ.5) CALL XABORT('SNTT3D: SYME BC NOT ALLOWED.') + IF(NCODE(I).EQ.7) CALL XABORT('SNTT3D: ZERO FLUX BC NOT ALLOWED.') + 240 CONTINUE +*---- +* CHECK FOR INVALID VIRTUAL ELEMENTS. +*---- + DO 292 I=2,LX-1 + DO 291 J=2,LY-1 + DO 290 K=2,LZ-1 + IF(MAT(I,J,K).EQ.0) THEN + L1=(NCODE(1).NE.1) + DO 251 J1=1,J-1 + DO 250 K1=1,K-1 + L1=L1.OR.(MAT(I,J1,K1).NE.0) + 250 CONTINUE + 251 CONTINUE + L2=(NCODE(2).NE.1) + DO 261 J1=J+1,LY + DO 260 K1=K+1,LZ + L2=L2.OR.(MAT(I,J1,K1).NE.0) + 260 CONTINUE + 261 CONTINUE + L3=(NCODE(3).NE.1) + DO 271 I1=1,I-1 + DO 270 K1=1,K-1 + L3=L3.OR.(MAT(I1,J,K1).NE.0) + 270 CONTINUE + 271 CONTINUE + L4=(NCODE(4).NE.1) + DO 281 I1=I+1,LX + DO 280 K1=K+1,LZ + L4=L4.OR.(MAT(I1,J,K1).NE.0) + 280 CONTINUE + 281 CONTINUE + L5=(NCODE(5).NE.1) + DO 301 I1=1,I-1 + DO 300 J1=1,J-1 + L5=L5.OR.(MAT(I1,J1,K).NE.0) + 300 CONTINUE + 301 CONTINUE + L6=(NCODE(6).NE.0) + DO 311 I1=I+1,LX + DO 310 J1=I+1,LY + L6=L6.OR.(MAT(I1,J1,K).NE.0) + 310 CONTINUE + 311 CONTINUE + IF(L1.AND.L2.AND.L3.AND.L4.AND.L5.AND.L6) THEN + WRITE(HSMG,'(17HSNTT3D: ELEMENT (,I3,1H,,I3,11H) CANNOT BE, + 1 9H VIRTUAL.)') I,J,K + CALL XABORT(HSMG) + ENDIF + ENDIF + 290 CONTINUE + 291 CONTINUE + 292 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(ZZ,YY,XX) + RETURN + END -- cgit v1.2.3