*DECK SNTRK SUBROUTINE SNTRK(MAXPTS,IPTRK,IPGEOM,IMPX,ISCHM,IELEM,ISPLH,INSB, 1 NLF,MAXIT,EPSI,ISCAT,IQUAD,LFIXUP,LIVO,ICL1,ICL2,LDSA,NSTART, 2 NSDSA,IELEMSA,ISOLVSA,LBIHET,LSHOOT,IBFP,MCELL,NMPI,NFOU, 3 EELEM,ESCHM,IGLK) * *----------------------------------------------------------------------- * *Purpose: * Recover of the geometry and tracking for SN methods. * *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): A. Hebert and N. Martin * *Parameters: input * MAXPTS allocated storage for arrays of dimension NEL. * IPTRK L_TRACK pointer to the tracking information. * IPGEOM L_GEOM pointer to the geometry. * IMPX print flag. * ISCHM method of spatial discretisation: * =1 High-Order Diamond Differencing (HODD) - default; * =2 Discontinuous Galerkin finite element method (DG); * =3 Adaptive weighted method (AWD). * IELEM measure of order of the spatial approximation polynomial: * =1 constant - default for HODD; * =2 linear - default for DG; * >3 higher orders. * ISPLH mesh-splitting index for hexagons into lozenges. * INSB group vectorization flag (=0/1 group vectorization off/on). * NLF SN order for the flux (even number). * MAXIT maximum number of inner iterations (default=100). * EPSI convergence criterion on inner iterations (default=5.E-5). * ISCAT anisotropy of one-speed sources: * =1 isotropic sources; * =2 linearly anisotropic sources. * IQUAD type of SN quadrature (=1 Level symmetric, type IQUAD; * =4 Gauss-Legendre and Gauss-Chebyshev; =10 product). * LFIXUP flag to enable negative flux fixup. * LIVO flag to enable Livolant acceleration. * ICL1 Number of free iterations with Livolant acceleration. * ICL2 Number of accelerated iterations with Livolant acceleration. * LDSA flag to enable diffusion synthetic acceleration. * NSTART restarts the GMRES method every NSTART iterations. * NSDSA number of inner flux iterations before enabling SA. * IELEMSA degree of the Lagrangian finite elements for the SA: * <0 order -IELEMSA primal finite elements; * >0 order IELEMSA dual finite elements. * ISOLVSA type of solver to be used for the SA: * 1 - BIVAC ; * 2 - TRIVAC. * LBIHET flag to enable the double-heterogeneity model. * LSHOOT enablig flag for the shooting method. * IBFP type of energy proparation relation: * =0 no Fokker-Planck term; * =1 Galerkin type; * =2 heuristic Przybylski and Ligou type. * MCELL number of macrocells along each axis (in Cartesian geometry) * for the parallelisation using the OpenMP paradigm; OR * number of macrocells along the z-axis (in hexagonal geometry) * for the parallelisation using the OpenMP paradigm. * NMPI number of macrocells along each axis (in Cartesian geometry) * or along the z-axis for the hexagonal geometry for the * parallelisation using the MPI paradigm, when using WYVERN. * NFOU number of Fourier frequencies in the range (2*pi/L) to be * investigated. * ESCHM method of energy discretisation: * =1 High-Order Diamond Differencing (HODD) - default; * =2 Discontinuous Galerkin finite element method (DG); * =3 Adaptive weighted method (AWD). * EELEM measure of order of the energy approximation polynomial: * =1 constant - default for HODD; * =2 linear - default for DG; * >3 higher orders. * IGLK angular interpolation type: * =0 classical SN method. * =1 Galerkin quadrature method (M = inv(D)) * =2 Galerkin quadrature method (D = inv(M)) * *----------------------------------------------------------------------- * USE GANLIB *---- * SUBROUTINE ARGUMENTS *---- TYPE(C_PTR) IPTRK,IPGEOM INTEGER MAXPTS,IMPX,ISCHM,IELEM,ISPLH,INSB,NLF,ISCAT,IQUAD, 1 MAXIT,ICL1,ICL2,NSTART,NSDSA,IELEMSA,ISOLVSA,MCELL,NMPI,NFOU, 2 EELEM,ESCHM,IGLK REAL EPSI LOGICAL LFIXUP,LDSA,LBIHET,LIVO,LSHOOT *---- * LOCAL VARIABLES *---- PARAMETER(NSTATE=40) LOGICAL ILK CHARACTER HSMG*131 INTEGER ISTATE(NSTATE),IGP(NSTATE),NCODE(6),ICODE(6),LOZSWP(3,6),P REAL ZCODE(6) INTEGER, ALLOCATABLE, DIMENSION(:) :: MAT,IDL,ISPLX,ISPLY,ISPLZ, 1 JOP,MRMX,MRMY,MRMZ,IL,IM INTEGER, ALLOCATABLE, DIMENSION(:,:) :: COORDMAP INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: KEYANI REAL, ALLOCATABLE, DIMENSION(:) :: VOL,XXX,YYY,ZZZ,UU,WW,PL,TPQ, 1 UPQ,VPQ,WPQ,ALPHA,PLZ,SURF,DU,DE,DZ,DC,DB,DA,DAL,WX,WE,CST,MN,DN *---- * SCRATCH STORAGE ALLOCATION *---- ALLOCATE(MAT(MAXPTS),VOL(MAXPTS),IDL(MAXPTS)) * CALL LCMGET(IPGEOM,'STATE-VECTOR',ISTATE) ITYPE=ISTATE(1) ISUB2=ISTATE(9) *---- * PARAMETER VALIDATION. *---- IF(ISUB2.NE.0) CALL XABORT('SNTRK: DISCRETIZATION NOT AVAILABLE 1 .(1)') IF((ISCHM.NE.1).AND.(ISCHM.NE.2).AND.(ISCHM.NE.3)) THEN CALL XABORT('SNTRK: SPATIAL DISCRETIZATION SCHEME NOT AVAILABLE 1. ONLY VALUES OF 1 (=DIAMOND-DIFFERENCE) OR 2 (=DISCONTINUOUS GALE 2RKIN) OR 3 (=ADAPTIVE WEIGHTED DIFFERENCE) ARE ALLOWED.') ENDIF IF((ESCHM.NE.1).AND.(ESCHM.NE.2).AND.(ESCHM.NE.3)) THEN CALL XABORT('SNTRK: ENERGY DISCRETIZATION SCHEME NOT AVAILABLE 1. ONLY VALUES OF 1 (=DIAMOND-DIFFERENCE) OR 2 (=DISCONTINUOUS GALE 2RKIN) OR 3 (=ADAPTIVE WEIGHTED DIFFERENCE) ARE ALLOWED.') ENDIF IF(ISCHM.EQ.3.OR.ESCHM.EQ.3) THEN IF(.NOT.(ITYPE.EQ.2.OR.ITYPE.EQ.5.OR.ITYPE.EQ.7)) THEN CALL XABORT('SNTRK: ADAPTIVE SCHEME ONLY AVAILABLE IN CARTESIAN' 1 //' GEOMETRIES') ELSEIF(IELEM.GT.1.AND.ISCHM.EQ.3.OR. 1 EELEM.GT.1.AND.ESCHM.EQ.3) THEN CALL XABORT('SNTRK: ADAPTIVE SCHEMES ONLY AVAILABLE FOR CONSTANT' 1 //' ORDER CLOSURE RELATION.') ENDIF ENDIF IF((ITYPE.NE.2).AND.(ITYPE.NE.3).AND.(ITYPE.NE.4).AND. 1 (ITYPE.NE.5).AND.(ITYPE.NE.6).AND.(ITYPE.NE.7).AND. 2 (ITYPE.NE.8).AND.(ITYPE.NE.9)) THEN CALL XABORT('SNTRK: DISCRETIZATION NOT AVAILABLE.(2)') ENDIF IF((LDSA).AND.(ISOLVSA.EQ.1))THEN IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) > CALL XABORT('SNTRK: SYNTHETIC ACCELERATION WITH BIVAC NOT AV 1AILABLE IN 3D.') ENDIF * ALLOCATE(XXX(MAXPTS+1),YYY(MAXPTS+1),ZZZ(MAXPTS+1)) ALLOCATE(ISPLX(MAXPTS),ISPLY(MAXPTS),ISPLZ(MAXPTS)) CALL READ3D(MAXPTS,MAXPTS,MAXPTS,MAXPTS,IPGEOM,IHEX,IR,ILK,SIDE, 1 XXX,YYY,ZZZ,IMPX,LX,LY,LZ,MAT,NEL,NCODE,ICODE,ZCODE,ISPLX,ISPLY, 2 ISPLZ,ISPLTH,ISPLTL) DEALLOCATE(ISPLZ,ISPLY,ISPLX) * IF(LX*LY*LZ.GT.MAXPTS) THEN WRITE (HSMG,'(38HSNTRK: MAXPTS SHOULD BE INCREASED FROM,I7, 1 3H TO,I7)') MAXPTS,LX*LY*LZ CALL XABORT(HSMG) ENDIF IF((ITYPE.EQ.2).OR.(ITYPE.EQ.3)) THEN ! 1-D AND 2-D CYLINDRICAL CASES. NCODE(3)=2 NCODE(4)=5 ZCODE(3)=1.0 ZCODE(4)=1.0 YYY(1)=0.0 YYY(2)=2.0 ELSE IF(ITYPE.EQ.6) THEN LY=LZ DO I=1,LZ+1 YYY(I)=ZZZ(I) ENDDO NCODE(3)=NCODE(5) NCODE(4)=NCODE(6) ZCODE(3)=ZCODE(5) ZCODE(4)=ZCODE(6) ICODE(3)=ICODE(5) ICODE(4)=ICODE(6) ENDIF IF(IBFP.EQ.0) EELEM=1 *---- * UNFOLD THE DOMAIN IN DIAGONAL SYMMETRY CASES. *---- IF(ITYPE.EQ.7) THEN ! CARTESIAN 3D CASE IDIAG=0 IF((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) THEN IDIAG=1 NCODE(3)=NCODE(1) NCODE(2)=NCODE(4) ZCODE(3)=ZCODE(1) ZCODE(2)=ZCODE(4) ICODE(3)=ICODE(1) ICODE(2)=ICODE(4) K=NEL DO IZ=LZ,1,-1 IOFF=(IZ-1)*LX*LY DO IY=LY,1,-1 DO IX=LX,IY+1,-1 MAT(IOFF+(IY-1)*LX+IX)=MAT(IOFF+(IX-1)*LY+IY) ENDDO DO IX=IY,1,-1 MAT(IOFF+(IY-1)*LX+IX)=MAT(K) K=K-1 ENDDO ENDDO ENDDO NEL=LX*LY*LZ IF(K.NE.0) THEN CALL XABORT('SNTRK: UNABLE TO UNFOLD THE DOMAIN.') ENDIF ELSE IF((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) THEN IDIAG=1 NCODE(1)=NCODE(3) NCODE(4)=NCODE(2) ZCODE(1)=ZCODE(3) ZCODE(4)=ZCODE(2) ICODE(1)=ICODE(3) ICODE(4)=ICODE(2) K=NEL DO IZ=LZ,1,-1 IOFF=(IZ-1)*LX*LY DO IY=LY,1,-1 DO IX=LX,IY,-1 MAT(IOFF+(IY-1)*LX+IX)=MAT(K) K=K-1 ENDDO ENDDO ENDDO DO IZ=1,LZ IOFF=(IZ-1)*LX*LY DO IY=1,LY DO IX=1,IY-1 MAT(IOFF+(IY-1)*LX+IX)=MAT(IOFF+(IX-1)*LY+IY) ENDDO ENDDO ENDDO NEL=LX*LY*LZ IF(K.NE.0) THEN CALL XABORT('SNTRK: UNABLE TO UNFOLD THE DOMAIN.') ENDIF ENDIF ELSE ! OTHER CASES IF((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) THEN NCODE(3)=NCODE(1) NCODE(2)=NCODE(4) ZCODE(3)=ZCODE(1) ZCODE(2)=ZCODE(4) ICODE(3)=ICODE(1) ICODE(2)=ICODE(4) K=LX*(LX+1)/2 DO IY=LY,1,-1 DO IX=LX,IY+1,-1 MAT((IY-1)*LX+IX)=MAT((IX-1)*LY+IY) ENDDO DO IX=IY,1,-1 MAT((IY-1)*LX+IX)=MAT(K) K=K-1 ENDDO ENDDO NEL=LX*LY IF(K.NE.0) THEN CALL XABORT('SNTRK: UNABLE TO UNFOLD THE DOMAIN.') ENDIF ELSE IF((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) THEN NCODE(1)=NCODE(3) NCODE(4)=NCODE(2) ZCODE(1)=ZCODE(3) ZCODE(4)=ZCODE(2) ICODE(1)=ICODE(3) ICODE(4)=ICODE(2) K=LX*(LX+1)/2 DO IY=LY,1,-1 DO IX=LX,IY,-1 MAT((IY-1)*LX+IX)=MAT(K) K=K-1 ENDDO ENDDO DO IY=1,LY DO IX=1,IY-1 MAT((IY-1)*LX+IX)=MAT((IX-1)*LY+IY) ENDDO ENDDO NEL=LX*LY IF(K.NE.0) THEN CALL XABORT('SNTRK: UNABLE TO UNFOLD THE DOMAIN.') ENDIF ENDIF ENDIF IF(IMPX.GT.5) THEN WRITE(6,600) 'NCODE',(NCODE(I),I=1,6) WRITE(6,600) 'ICODE',(ICODE(I),I=1,6) WRITE(6,600) 'MAT',(MAT(I),I=1,LX*LY*LZ) ENDIF *--- * CALL TO THE SN TRACKING MODULE RELEVANT TO EACH GEOMETRY *--- NDIM=0 IF(ITYPE.EQ.2) THEN * 1D SLAB GEOMETRY. NDIM=1 IF(ISCAT.EQ.0) CALL XABORT('SNTRK: SCAT NOT DEFINED.') IF(IGLK.NE.0) THEN NSCT=NLF ELSE NSCT=ISCAT ENDIF ALLOCATE(UU(NLF),WW(NLF),WX(IELEM+1),WE(EELEM+1), 1 CST(MAX(IELEM,EELEM)),MN(NSCT*NLF),DN(NLF*NSCT),IL(NSCT), 2 IM(NSCT),PL(NSCT*NLF)) CALL SNT1DP(IMPX,LX,IELEM,NCODE,ZCODE,XXX,NLF,NSCT,UU,WW,PL, 1 VOL,IDL,LL4,NUN,LSHOOT,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM, 2 IGLK,MN,DN,IL,IM,IQUAD) CALL LCMPUT(IPTRK,'U',NLF,2,UU) CALL LCMPUT(IPTRK,'W',NLF,2,WW) CALL LCMPUT(IPTRK,'WX',IELEM+1,2,WX) IF(IBFP.NE.0) CALL LCMPUT(IPTRK,'WE',EELEM+1,2,WE) CALL LCMPUT(IPTRK,'CST',MAX(IELEM,EELEM),2,CST) CALL LCMPUT(IPTRK,'MN',NSCT*NLF,2,MN) CALL LCMPUT(IPTRK,'DN',NLF*NSCT,2,DN) CALL LCMPUT(IPTRK,'IL',NSCT,1,IL) CALL LCMPUT(IPTRK,'IM',NSCT,1,IM) CALL LCMPUT(IPTRK,'PL',NSCT*NLF,2,PL) DEALLOCATE(WW,UU,WX,WE,CST,MN,DN,IL,IM,PL) ! For Fourier Analysis IF(NFOU.GT.0)THEN XLEN = XXX(LX+1) CALL LCMPUT(IPTRK,'XXX',LX+1,2,XXX) CALL LCMPUT(IPTRK,'XLEN',1,2,XLEN) ENDIF ELSE IF(ITYPE.EQ.3) THEN * 1D CYLINDRICAL GEOMETRY. NDIM=1 IF(ISCAT.EQ.0) CALL XABORT('SNTRK: SCAT NOT DEFINED.') NSCT=(ISCAT/2)*(ISCAT/2+1)+(ISCAT+1)*MOD(ISCAT,2)/2 IF(IQUAD.GE.10) THEN * PRODUCT QUADRATURE. NPQ=(NLF**2)/2 ELSE NPQ=NLF*(1+NLF/2)/2 ENDIF ALLOCATE(JOP(NLF/2),UU(NLF/2),WW(NLF/2),TPQ(NPQ),UPQ(NPQ), 1 VPQ(NPQ),WPQ(NPQ),ALPHA(NPQ),PLZ(NSCT*NLF/2),PL(NSCT*NPQ), 2 SURF(LX+1),IL(NSCT),IM(NSCT)) CALL SNT1DC(IMPX,LX,NCODE,ZCODE,XXX,NLF,NPQ,NSCT,IQUAD,JOP, 1 UU,WW,TPQ,UPQ,VPQ,WPQ,ALPHA,PLZ,PL,VOL,IDL,SURF,IL,IM) DEALLOCATE(VPQ,TPQ,WW) CALL LCMPUT(IPTRK,'JOP',NLF/2,1,JOP) CALL LCMPUT(IPTRK,'U',NLF/2,2,UU) CALL LCMPUT(IPTRK,'UPQ',NPQ,2,UPQ) CALL LCMPUT(IPTRK,'WPQ',NPQ,2,WPQ) CALL LCMPUT(IPTRK,'ALPHA',NPQ,2,ALPHA) CALL LCMPUT(IPTRK,'PLZ',NSCT*NLF/2,2,PLZ) CALL LCMPUT(IPTRK,'PL',NSCT*NPQ,2,PL) CALL LCMPUT(IPTRK,'SURF',LX+1,2,SURF) CALL LCMPUT(IPTRK,'IL',NSCT,1,IL) CALL LCMPUT(IPTRK,'IM',NSCT,1,IM) DEALLOCATE(SURF,PL,PLZ,ALPHA,WPQ,UPQ,UU,JOP,IL,IM) LL4=LX*NSCT NUN=LL4 ELSE IF(ITYPE.EQ.4) THEN * 1D SPHERICAL GEOMETRY. NDIM=1 IF(ISCAT.EQ.0) CALL XABORT('SNTRK: SCAT NOT DEFINED.') NSCT=ISCAT ALLOCATE(UU(NLF),WW(NLF),ALPHA(NLF),PLZ(NSCT),PL(NSCT*NLF), 1 SURF(LX+1),IL(NSCT),IM(NSCT)) CALL SNT1DS(IMPX,LX,NCODE,ZCODE,XXX,NLF,NSCT,UU,WW,ALPHA, 1 PLZ,PL,VOL,IDL,SURF,IQUAD,IL,IM) CALL LCMPUT(IPTRK,'U',NLF,2,UU) CALL LCMPUT(IPTRK,'W',NLF,2,WW) CALL LCMPUT(IPTRK,'ALPHA',NLF,2,ALPHA) CALL LCMPUT(IPTRK,'PLZ',NSCT,2,PLZ) CALL LCMPUT(IPTRK,'PL',NSCT*NLF,2,PL) CALL LCMPUT(IPTRK,'SURF',LX+1,2,SURF) CALL LCMPUT(IPTRK,'XXX',LX+1,2,XXX) CALL LCMPUT(IPTRK,'IL',NSCT,1,IL) CALL LCMPUT(IPTRK,'IM',NSCT,1,IM) DEALLOCATE(SURF,PL,PLZ,ALPHA,WW,UU,IL,IM) LL4=LX*NSCT NUN=LL4 ELSE IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) THEN * 2D GEOMETRIES: CARTESIAN; TUBE; HEXAGONAL NDIM=2 IF(ISCAT.EQ.0) CALL XABORT('SNTRK: SCAT NOT DEFINED.') IF(IQUAD.GE.10) THEN NPQ=NLF**2 ELSE NPQ=(NLF+4)*NLF/2 ENDIF IF(IGLK.NE.0) THEN NPQ=(NLF+2)*NLF/2 NSCT=NPQ ELSE NSCT=ISCAT*(ISCAT+1)/2 ENDIF IGE=0 * IF(ITYPE.EQ.5) THEN ! 2D Cartesian IF(NFOU.GT.0)THEN XLEN = XXX(LX+1) YLEN = YYY(LY+1) CALL LCMPUT(IPTRK,'XXX',LX+1,2,XXX) CALL LCMPUT(IPTRK,'XLEN',1,2,XLEN) CALL LCMPUT(IPTRK,'YYY',LY+1,2,YYY) CALL LCMPUT(IPTRK,'YLEN',1,2,YLEN) ENDIF ELSEIF(ITYPE.EQ.6) THEN ! 2D Tube IGE=1 ELSEIF(ITYPE.EQ.8) THEN ! 2D Hexagonal IGE=2 NHEX=LX/(3*ISPLH**2) ALLOCATE(COORDMAP(3,NHEX)) COORDMAP(:,:)=0 CALL SNTSFH(IMPX,ITYPE,NHEX,LZ,MCELL,ISPLH,MAT,LOZSWP, > COORDMAP) CALL LCMPUT(IPTRK,'LOZSWP',3*6,1,LOZSWP) CALL LCMPUT(IPTRK,'COORDMAP',3*NHEX,1,COORDMAP) CALL LCMPUT(IPTRK,'SIDE',1,2,SIDE) DEALLOCATE(COORDMAP) ENDIF * ALLOCATE(MRMX(NPQ),MRMY(NPQ)) ALLOCATE(DU(NPQ),DE(NPQ),WW(NPQ),DB(LX*NPQ),DA(LX*LY*NPQ), 1 DAL(LX*LY*NPQ),WX(IELEM+1), 2 WE(EELEM+1),CST(MAX(IELEM,EELEM)),MN(NSCT*NPQ),DN(NPQ*NSCT), 3 IL(NSCT),IM(NSCT),PL(NSCT*NPQ)) CALL SNTT2D(IGE,IMPX,LX,LY,SIDE,IELEM,NLF,NPQ,NSCT,IQUAD, 1 NCODE,ZCODE,MAT,XXX,YYY,VOL,IDL,DU,DE,WW,MRMX,MRMY,DB,DA,DAL, 2 PL,LL4,NUN,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM,IGLK,MN,DN,IL,IM, 3 ISCAT) CALL LCMPUT(IPTRK,'DU',NPQ,2,DU) CALL LCMPUT(IPTRK,'DE',NPQ,2,DE) CALL LCMPUT(IPTRK,'W',NPQ,2,WW) CALL LCMPUT(IPTRK,'MRM',NPQ,1,MRMX) CALL LCMPUT(IPTRK,'MRMY',NPQ,1,MRMY) CALL LCMPUT(IPTRK,'DB',LX*NPQ,2,DB) CALL LCMPUT(IPTRK,'DA',LX*LY*NPQ,2,DA) IF(IGE.EQ.1) CALL LCMPUT(IPTRK,'DAL',LX*LY*NPQ,2,DAL) CALL LCMPUT(IPTRK,'PL',NSCT*NPQ,2,PL) CALL LCMPUT(IPTRK,'WX',IELEM+1,2,WX) IF(IBFP.NE.0) CALL LCMPUT(IPTRK,'WE',EELEM+1,2,WE) CALL LCMPUT(IPTRK,'CST',MAX(IELEM,EELEM),2,CST) CALL LCMPUT(IPTRK,'MN',NSCT*NPQ,2,MN) CALL LCMPUT(IPTRK,'DN',NPQ*NSCT,2,DN) CALL LCMPUT(IPTRK,'IL',NSCT,1,IL) CALL LCMPUT(IPTRK,'IM',NSCT,1,IM) CALL LCMPUT(IPTRK,'PL',NSCT*NPQ,2,PL) DEALLOCATE(DAL,DA,DB,WW,DE,DU,WX,WE,CST,MN,DN,IL,IM,PL) DEALLOCATE(MRMY,MRMX) ELSE IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) THEN * 3D GEOMETRIES: CARTESIAN; HEXAGONAL NDIM=3 IF(ISCAT.EQ.0) CALL XABORT('SNTRK: SCAT NOT DEFINED.') IF(IQUAD.GE.10) THEN NPQ=2*NLF**2 ELSE NPQ=(NLF+2)*NLF ENDIF IF(IGLK.NE.0) THEN NSCT=NPQ ELSE NSCT=ISCAT*(ISCAT+1) ENDIF IGE=0 * IF(ITYPE.EQ.9) THEN ! 3D Hexagonal IGE=2 NHEX =LX/(3*ISPLH**2) ALLOCATE(COORDMAP(3,NHEX)) COORDMAP(:,:)=0 CALL SNTSFH(IMPX,ITYPE,NHEX,LZ,MCELL,ISPLH,MAT,LOZSWP, > COORDMAP) CALL LCMPUT(IPTRK,'LOZSWP',3*6,1,LOZSWP) CALL LCMPUT(IPTRK,'COORDMAP',3*NHEX,1,COORDMAP) CALL LCMPUT(IPTRK,'SIDE',1,2,SIDE) DEALLOCATE(COORDMAP) ENDIF * ALLOCATE(MRMX(NPQ),MRMY(NPQ),MRMZ(NPQ)) ALLOCATE(DU(NPQ),DE(NPQ),DZ(NPQ),WW(NPQ),DC(LX*LY*NPQ), 1 DB(LX*LZ*NPQ),DA(LY*LZ*NPQ),WX(IELEM+1), 2 WE(EELEM+1),CST(MAX(IELEM,EELEM)),MN(NSCT*NPQ),DN(NPQ*NSCT), 3 IL(NSCT),IM(NSCT),PL(NSCT*NPQ)) CALL SNTT3D(IGE,IMPX,LX,LY,LZ,SIDE,IELEM,NLF,NPQ,NSCT,IQUAD, 1 NCODE,ZCODE,MAT,XXX,YYY,ZZZ,VOL,IDL,DU,DE,DZ,WW,MRMX,MRMY, 2 MRMZ,DC,DB,DA,PL,LL4,NUN,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM, 3 IGLK,MN,DN,IL,IM,ISCAT) * CALL LCMPUT(IPTRK,'DU',NPQ,2,DU) CALL LCMPUT(IPTRK,'DE',NPQ,2,DE) CALL LCMPUT(IPTRK,'DZ',NPQ,2,DZ) CALL LCMPUT(IPTRK,'W',NPQ,2,WW) CALL LCMPUT(IPTRK,'MRMX',NPQ,1,MRMX) CALL LCMPUT(IPTRK,'MRMY',NPQ,1,MRMY) CALL LCMPUT(IPTRK,'MRMZ',NPQ,1,MRMZ) CALL LCMPUT(IPTRK,'DC',LX*LY*NPQ,2,DC) CALL LCMPUT(IPTRK,'DB',LX*LZ*NPQ,2,DB) CALL LCMPUT(IPTRK,'DA',LY*LZ*NPQ,2,DA) CALL LCMPUT(IPTRK,'PL',NSCT*NPQ,2,PL) CALL LCMPUT(IPTRK,'WX',IELEM+1,2,WX) IF(IBFP.NE.0) CALL LCMPUT(IPTRK,'WE',EELEM+1,2,WE) CALL LCMPUT(IPTRK,'CST',MAX(IELEM,EELEM),2,CST) CALL LCMPUT(IPTRK,'MN',NSCT*NPQ,2,MN) CALL LCMPUT(IPTRK,'DN',NPQ*NSCT,2,DN) CALL LCMPUT(IPTRK,'IL',NSCT,1,IL) CALL LCMPUT(IPTRK,'IM',NSCT,1,IM) CALL LCMPUT(IPTRK,'PL',NSCT*NPQ,2,PL) DEALLOCATE(DA,DB,DC,WW,DZ,DE,DU,WX,WE,CST,MN,DN,IL,IM,PL) DEALLOCATE(MRMZ,MRMY,MRMX) ELSE CALL XABORT('SNTRK: UNKNOWN GEOMETRY.') ENDIF DEALLOCATE(YYY,ZZZ) *---- * THE NUMBER OF UNKNOWNS OF A BOLTZMANN-FOKKER-PLANCK DISCRETIZATION IS * INCREASED TO HOLD SLOWING-DOWN ANGULAR FLUXES. *---- IF(IBFP.GT.0) THEN IF(NDIM.EQ.1) THEN NUN=NUN+IELEM*NLF*NEL ELSE IF(NDIM.EQ.2) THEN NUN=NUN+IELEM**2*NPQ*NEL ELSE IF(NDIM.EQ.3) THEN NUN=NUN+IELEM**3*NPQ*NEL ELSE CALL XABORT('SNTRK: FOKKER-PLANCK NOT IMPLEMENTED.') ENDIF ENDIF IF(IMPX.GT.0) WRITE (6,'(/33H SNTRK: ORDER OF LINEAR SYSTEMS =, 1 I10/8X,37HNUMBER OF UNKNOWNS PER ENERGY GROUP =,I10)') LL4,NUN * IF(IMPX.GT.5) THEN I1=1 DO I=1,(NEL-1)/8+1 I2=I1+7 IF(I2.GT.NEL) I2=NEL WRITE (6,620) (J,J=I1,I2) WRITE (6,630) (MAT(J),J=I1,I2) WRITE (6,640) (IDL(J),J=I1,I2) WRITE (6,650) (VOL(J),J=I1,I2) I1=I1+8 ENDDO ENDIF *---- * SYNTHETIC ACCELERATION TRACKING INFORMATION. *---- IF(LDSA) THEN IF(IMPX.GT.0) WRITE (6,'(/32H SNTRK: SYNTHETIC ACCELERATION I, 1 19HNFORMATION FOLLOWS:)') CALL LCMSIX(IPTRK,'DSA',1) ICOL=3 ! Gauss-Legendre quadrature NLFSA=2 ! P1 method ISPN=1 ! simplified PN method ISCSA=1 ! isotropic scattering NVD=1 ! SN-type VOID boundary conditions NADI=2 ! ADI iteration ICHX=2 ! Raviart-Thomas finite elements ISEG=0 ! scalar algorithm IMPV=0 ! print parameter for vector operations IF(MAXPTS.EQ.0) CALL XABORT('SNTRK: MAXPTS NOT DEFINED.') IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) THEN CALL TRITRK(MAXPTS,IPTRK,IPGEOM,IMPX,IELEMSA,ICOL,ICHX,ISEG, 1 IMPV,NLFSA,NVD,ISPN,ISCSA,NADI) ELSE IF(ISOLVSA.EQ.1)THEN CALL BIVTRK(MAXPTS,IPTRK,IPGEOM,IMPX,IELEMSA,ICOL, 1 NLFSA,NVD,ISPN,ISCSA) ELSEIF(ISOLVSA.EQ.2)THEN CALL TRITRK(MAXPTS,IPTRK,IPGEOM,IMPX,IELEMSA,ICOL,ICHX, 1 ISEG,IMPV,NLFSA,NVD,ISPN,ISCSA,NADI) ELSE CALL XABORT('SNTRK: UNDEFINED SOLVER OPTION FOR ' 1 //'SYNTHETIC ACCELERATION.') ENDIF ENDIF CALL LCMSIX(IPTRK,' ',2) ENDIF *---- * SAVE GENERAL AND SN-SPECIFIC TRACKING INFORMATION. *---- IGP(:NSTATE)=0 IGP(1)=NEL IGP(2)=NUN IF(ILK) THEN IGP(3)=0 ELSE IGP(3)=1 ENDIF IGP(4)=ISTATE(7) IGP(5)=1 IGP(6)=ITYPE IGP(7)=NSCT IGP(8)=IELEM IGP(9)=NDIM IGP(10)=ISCHM IGP(11)=LL4 IGP(12)=LX IGP(13)=LY IGP(14)=LZ IGP(15)=NLF IGP(16)=ISCAT IGP(17)=IQUAD IGP(18)=0 IF(LFIXUP) IGP(18)=1 IGP(19)=0 IF(LDSA) IGP(19)=1 IGP(20)=NSTART IGP(21)=NSDSA IGP(22)=MAXIT IF(LIVO) IGP(23)=1 IGP(24)=ICL1 IGP(25)=ICL2 IGP(26)=ISPLH IGP(27)=INSB IGP(28)=MCELL IF((ITYPE.EQ.3).OR.(ITYPE.GE.4)) IGP(29)=1 IGP(30)=0 IF(LSHOOT) IGP(30)=1 IGP(31)=IBFP IGP(32)=NMPI IGP(33)=ISOLVSA IGP(34)=NFOU IGP(35)=EELEM IGP(36)=ESCHM IGP(37)=IGLK IF(LBIHET) IGP(40)=1 CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,IGP) CALL LCMPUT(IPTRK,'MATCOD',NEL,1,MAT) CALL LCMPUT(IPTRK,'VOLUME',NEL,2,VOL) CALL LCMPUT(IPTRK,'KEYFLX',NEL,1,IDL) CALL LCMPUT(IPTRK,'NCODE',6,1,NCODE) CALL LCMPUT(IPTRK,'ICODE',6,1,ICODE) CALL LCMPUT(IPTRK,'ZCODE',6,2,ZCODE) CALL LCMPUT(IPTRK,'EPSI',1,2,EPSI) IF(ITYPE.EQ.4) CALL LCMPUT(IPTRK,'XXX',LX+1,2,XXX) DEALLOCATE(XXX) *---- * SET KEYFLX$ANIS *---- NLIN=IELEM**NDIM*EELEM ALLOCATE(KEYANI(NEL,NLIN,NSCT)) DO IR=1,NEL IND=IDL(IR) DO IE=1,NLIN DO P=1,NSCT IF(IND.EQ.0) THEN KEYANI(IR,IE,P)=0 ELSE KEYANI(IR,IE,P)=IND+(P-1)*NLIN+IE-1 ENDIF ENDDO ENDDO ENDDO CALL LCMPUT(IPTRK,'KEYFLX$ANIS',NEL*NLIN*NSCT,1,KEYANI) DEALLOCATE(KEYANI) *---- * SCRATCH STORAGE DEALLOCATION *---- DEALLOCATE(IDL,VOL,MAT) RETURN * 600 FORMAT(/25H SNTRK: VALUES OF VECTOR ,A6,4H ARE/(1X,1P,20I6)) 620 FORMAT (///11H REGION ,8(I8,6X,1HI)) 630 FORMAT ( 11H MIXTURE ,8(I8,6X,1HI)) 640 FORMAT ( 11H POINTER ,8(I8,6X,1HI)) 650 FORMAT ( 11H VOLUME ,8(1P,E13.6,2H I)) END