diff options
Diffstat (limited to 'Trivac/src/NSST.f')
| -rwxr-xr-x | Trivac/src/NSST.f | 370 |
1 files changed, 370 insertions, 0 deletions
diff --git a/Trivac/src/NSST.f b/Trivac/src/NSST.f new file mode 100755 index 0000000..b7cb74d --- /dev/null +++ b/Trivac/src/NSST.f @@ -0,0 +1,370 @@ +*DECK NSST + SUBROUTINE NSST(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Nodal expansion method (NEM) tracking operator. +* +*Copyright: +* Copyright (C) 2020 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/output +* NENTRY number of LCM objects or files used by the operator. +* HENTRY name of each LCM object or file: +* HENTRY(1): create or modification type(L_TRACK); +* HENTRY(2): read-only type(L_GEOM). +* IENTRY type of each LCM object or file: +* =1 LCM memory object; =2 XSM file; =3 sequential binary file; +* =4 sequential ascii file. +* JENTRY access of each LCM object or file: +* =0 the LCM object or file is created; +* =1 the LCM object or file is open for modifications; +* =2 the LCM object or file is open in read-only mode. +* KENTRY LCM object address or file unit number. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY) + CHARACTER HENTRY(NENTRY)*12 + TYPE(C_PTR) KENTRY(NENTRY) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (NSTATE=40) + CHARACTER TEXT4*4,TEXT12*12,TITLE*72,HSIGN*12 + DOUBLE PRECISION DFLOTT + LOGICAL ILK + INTEGER ISTATE(NSTATE),NCODE(6),ICODE(6) + REAL ZCODE(6) + TYPE(C_PTR) IPGEO,IPTRK +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: MAT,IDL,ISPLX,ISPLY,ISPLZ, + 1 MUX,MUY,MUZ,IMAX,IMAY,IMAZ,IPY,IPZ + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KN,IQFR + REAL, ALLOCATABLE, DIMENSION(:) :: XXX,YYY,ZZZ,XX,YY,ZZ,VOL + REAL, ALLOCATABLE, DIMENSION(:,:) :: QFR +*---- +* PARAMETER VALIDATION. +*---- + IF(NENTRY.NE.2) CALL XABORT('NSST: TWO PARAMETERS EXPECTED.') + IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2)) CALL XABORT('NSST: L' + 1 //'CM OBJECT EXPECTED AT LHS.') + IF((JENTRY(1).NE.0).AND.(JENTRY(1).NE.1)) CALL XABORT('NSST: E' + 1 //'NTRY IN CREATE OR MODIFICATION MODE EXPECTED.') + IF((JENTRY(2).NE.2).OR.((IENTRY(2).NE.1).AND.(IENTRY(2).NE.2))) + 1 CALL XABORT('NSST: LCM OBJECT IN READ-ONLY MODE EXPECTED AT R' + 2 //'HS.') + CALL LCMGTC(KENTRY(2),'SIGNATURE',12,HSIGN) + IF(HSIGN.NE.'L_GEOM') THEN + TEXT12=HENTRY(2) + CALL XABORT('NSST: SIGNATURE OF '//TEXT12//' IS '//HSIGN// + 1 '. L_GEOM EXPECTED.') + ENDIF + IPTRK=KENTRY(1) + IPGEO=KENTRY(2) + HSIGN='L_TRACK' + CALL LCMPTC(IPTRK,'SIGNATURE',12,HSIGN) + HSIGN='TRIVAC' + CALL LCMPTC(IPTRK,'TRACK-TYPE',12,HSIGN) + CALL LCMGET(IPGEO,'STATE-VECTOR',ISTATE) + IDIM=0 + ITYPE=ISTATE(1) + IF(ITYPE.EQ.2) THEN + IDIM=1 + ELSE IF(ITYPE.EQ.5) THEN + IDIM=2 + ELSE IF(ITYPE.EQ.7) THEN + IDIM=3 + ELSE + CALL XABORT('NSST: 1D, 2D OR 3D CARTESIAN GEOMETRY EXPECTED.') + ENDIF + NX=ISTATE(3) + NY=ISTATE(4) + NZ=ISTATE(5) + CALL LCMLEN(IPGEO,'BIHET',ILONG,ITYLCM) + IF(ILONG.NE.0) CALL XABORT('NSST: DOUBLE-HETEROGENEITY NOT SUPPO' + 1 //'RTED.') +* + IMPX=1 + TITLE=' ' + IGMAX=0 + ICHX=5 + NADI=2 + IF(IDIM.EQ.1) THEN + MAXPTS=NX + ELSE IF(IDIM.EQ.2) THEN + MAXPTS=NX*NY + ELSE + MAXPTS=NX*NY*NZ + ENDIF + IF(JENTRY(1).EQ.1) THEN + CALL LCMGTC(IPTRK,'SIGNATURE',12,HSIGN) + IF(HSIGN.NE.'L_TRACK') THEN + TEXT12=HENTRY(1) + CALL XABORT('NSST: SIGNATURE OF '//TEXT12//' IS '//HSIGN// + 1 '. L_TRACK EXPECTED.') + ENDIF + CALL LCMGTC(IPTRK,'TRACK-TYPE',12,HSIGN) + IF(HSIGN.NE.'TRIVAC') THEN + TEXT12=HENTRY(3) + CALL XABORT('NSST: TRACK-TYPE OF '//TEXT12//' IS '//HSIGN + 1 //'. TRIVAC EXPECTED.') + ENDIF + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + ICHX=ISTATE(12) ! CMFD/NEM/ANM + IGMAX=ISTATE(39) + CALL LCMLEN(IPTRK,'TITLE',LENGT,ITYLCM) + IF(LENGT.GT.0) CALL LCMGTC(IPTRK,'TITLE',72,TITLE) + ENDIF + 10 CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) + IF(INDIC.EQ.10) GO TO 30 + IF(INDIC.NE.3) CALL XABORT('NSST: CHARACTER DATA EXPECTED.') + IF(TEXT4.EQ.'EDIT') THEN + CALL REDGET(INDIC,IMPX,FLOTT,TEXT4,DFLOTT) + IF(INDIC.NE.1) CALL XABORT('NSST: INTEGER DATA EXPECTED(1).') + ELSE IF(TEXT4.EQ.'TITL') THEN + CALL REDGET(INDIC,NITMA,FLOTT,TITLE,DFLOTT) + IF(INDIC.NE.3) CALL XABORT('NSST: TITLE EXPECTED.') + ELSE IF(TEXT4.EQ.'MAXR') THEN + CALL REDGET(INDIC,MAXPTS,FLOTT,TEXT4,DFLOTT) + IF(INDIC.NE.1) CALL XABORT('NSST: INTEGER DATA EXPECTED(2).') + ELSE IF(TEXT4.EQ.'CMFD') THEN + ICHX=4 + ELSE IF(TEXT4.EQ.'NEM') THEN + ICHX=5 + ELSE IF(TEXT4.EQ.'ANM') THEN + ICHX=6 + ELSE IF(TEXT4.EQ.'HYPE') THEN + CALL REDGET(INDIC,IGMAX,FLOTT,TEXT4,DFLOTT) + IF(INDIC.NE.1) CALL XABORT('NSST: INTEGER DATA EXPECTED(3).') + ELSE IF(TEXT4.EQ.'ADI') THEN + CALL REDGET(INDIC,NADI,FLOTT,TEXT4,DFLOTT) + IF(INDIC.NE.1) CALL XABORT('NSST: INTEGER DATA EXPECTED(4).') + ELSE IF(TEXT4.EQ.';') THEN + GO TO 30 + ELSE + CALL XABORT('NSST: '//TEXT4//' IS AN INVALID KEYWORD.') + ENDIF + GO TO 10 +*---- +* SCRATCH STORAGE ALLOCATION +*---- + 30 IF(IMPX.GT.1) WRITE(6,100) TITLE + ALLOCATE(XXX(MAXPTS+1),YYY(MAXPTS+1),ZZZ(MAXPTS+1),MAT(MAXPTS), + 1 IDL(MAXPTS),VOL(MAXPTS),XX(MAXPTS),YY(MAXPTS),ZZ(MAXPTS), + 2 KN(6,MAXPTS),QFR(6,MAXPTS),IQFR(6,MAXPTS)) +*---- +* RECOVER TRACKING INFORMATION +*---- + ALLOCATE(ISPLX(MAXPTS),ISPLY(MAXPTS),ISPLZ(MAXPTS)) + CALL READ3D(MAXPTS,MAXPTS,MAXPTS,MAXPTS,IPGEO,IHEX,IR,ILK,SIDE, + 1 XXX,YYY,ZZZ,IMPX,NX,NY,NZ,MAT,NEL,NCODE,ICODE,ZCODE,ISPLX,ISPLY, + 2 ISPLZ,ISPLH,ISPLL) + DEALLOCATE(ISPLX,ISPLY,ISPLZ) + IF(IDIM.EQ.1) THEN +* 1D GEOMETRY + NY=1 + NCODE(3)=2 + NCODE(4)=2 + ZCODE(3)=1.0 + ZCODE(4)=1.0 + YYY(1)=0.0 + YYY(2)=1.0 + ENDIF + IF(IDIM.LE.2) THEN +* 1D OR 2D GEOMETRY + NZ=1 + NCODE(5)=2 + NCODE(6)=2 + ZCODE(5)=1.0 + ZCODE(6)=1.0 + ZZZ(1)=0.0 + ZZZ(2)=1.0 + ENDIF +*---- +* UNFOLD THE DOMAIN IN DIAGONAL SYMMETRY CASES. +*---- + IDIAG=0 + IF((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) THEN + IDIAG=1 + NCODE(3)=NCODE(1) + NCODE(2)=NCODE(4) + ICODE(3)=ICODE(1) + ICODE(2)=ICODE(4) + ZCODE(3)=ZCODE(1) + ZCODE(2)=ZCODE(4) + K=NEL + DO IZ=NZ,1,-1 + IOFF=(IZ-1)*NX*NY + DO IY=NY,1,-1 + DO IX=NX,IY+1,-1 + MAT(IOFF+(IY-1)*NX+IX)=MAT(IOFF+(IX-1)*NY+IY) + ENDDO + DO IX=IY,1,-1 + MAT(IOFF+(IY-1)*NX+IX)=MAT(K) + K=K-1 + ENDDO + ENDDO + ENDDO + NEL=NX*NY*NZ + IF(K.NE.0) THEN + CALL XABORT('TRITRK: UNABLE TO UNFOLD THE DOMAIN(1).') + ENDIF + ELSE IF((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) THEN + IDIAG=1 + NCODE(1)=NCODE(3) + NCODE(4)=NCODE(2) + ICODE(1)=ICODE(3) + ICODE(4)=ICODE(2) + ZCODE(1)=ZCODE(3) + ZCODE(4)=ZCODE(2) + K=NEL + DO IZ=NZ,1,-1 + IOFF=(IZ-1)*NX*NY + DO IY=NY,1,-1 + DO IX=NX,IY,-1 + MAT(IOFF+(IY-1)*NX+IX)=MAT(K) + K=K-1 + ENDDO + ENDDO + ENDDO + DO IZ=1,NZ + IOFF=(IZ-1)*NX*NY + DO IY=1,NY + DO IX=1,IY-1 + MAT(IOFF+(IY-1)*NX+IX)=MAT(IOFF+(IX-1)*NY+IY) + ENDDO + ENDDO + ENDDO + NEL=NX*NY*NZ + IF(K.NE.0) THEN + CALL XABORT('TRITRK: UNABLE TO UNFOLD THE DOMAIN(2).') + ENDIF + ENDIF + IF(IMPX.GT.5) THEN + WRITE(6,120) 'NCODE',(NCODE(I),I=1,6) + WRITE(6,120) 'MAT',(MAT(I),I=1,NX*NY*NZ) + ENDIF +*---- +* SET TRACKING INFORMATION +*---- + LL4F=0 + DO KEL=1,NEL + IF(MAT(KEL).GT.0) LL4F=LL4F+1 + ENDDO + ALLOCATE(MUX(LL4F),MUY(LL4F),MUZ(LL4F),IMAX(LL4F),IMAY(LL4F), + 1 IMAZ(LL4F),IPY(LL4F),IPZ(LL4F)) + CALL NSSDFC(IMPX,IDIM,NX,NY,NZ,NCODE,ICODE,ZCODE,MAT,XXX,YYY,ZZZ, + 1 LL4F,LL4X,LL4Y,LL4Z,VOL,XX,YY,ZZ,IDL,KN,QFR,IQFR,MUX,MUY,MUZ, + 2 IMAX,IMAY,IMAZ,IPY,IPZ) + IF(IDIM.EQ.1) THEN + NUN=LL4F*3+LL4X + ELSE IF(IDIM.EQ.2) THEN + NUN=LL4F*5+LL4X+LL4Y + ELSE + NUN=LL4F*7+LL4X+LL4Y+LL4Z + ENDIF +*---- +* SAVE INFORMATION ON LCM +*---- + ISTATE(:)=0 + ISTATE(1)=NEL + ISTATE(2)=NUN + ISTATE(4)=MAXVAL(MAT(:NEL)) + ISTATE(6)=ITYPE ! Geometry type + ISTATE(12)=ICHX ! CMFD/NEM/ANM + ISTATE(14)=NX + IF(IDIM.GE.2) ISTATE(15)=NY + IF(IDIM.EQ.3) ISTATE(16)=NZ + ISTATE(25)=LL4F + ISTATE(27)=LL4X + ISTATE(28)=LL4Y + ISTATE(29)=LL4Z + ISTATE(33)=NADI + ISTATE(39)=IGMAX + CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,ISTATE) + CALL LCMPUT(IPTRK,'NCODE',6,1,NCODE) + CALL LCMPUT(IPTRK,'ZCODE',6,2,ZCODE) + CALL LCMPUT(IPTRK,'ICODE',6,1,ICODE) + CALL LCMPUT(IPTRK,'MATCOD',NEL,1,MAT) + CALL LCMPUT(IPTRK,'VOLUME',NEL,2,VOL) + CALL LCMPUT(IPTRK,'KEYFLX',NEL,1,IDL) + CALL LCMPUT(IPTRK,'XX',NEL,2,XX) + CALL LCMPUT(IPTRK,'XXX',NX+1,2,XXX) + IF(IDIM.GE.2) THEN + CALL LCMPUT(IPTRK,'YY',NEL,2,YY) + CALL LCMPUT(IPTRK,'YYY',NY+1,2,YYY) + ENDIF + IF(IDIM.EQ.3) THEN + CALL LCMPUT(IPTRK,'ZZ',NEL,2,ZZ) + CALL LCMPUT(IPTRK,'ZZZ',NZ+1,2,ZZZ) + ENDIF + CALL LCMPUT(IPTRK,'KN',6*NEL,1,KN) + CALL LCMPUT(IPTRK,'QFR',6*NEL,2,QFR) + CALL LCMPUT(IPTRK,'IQFR',6*NEL,1,IQFR) + CALL LCMPUT(IPTRK,'MUX',LL4F,1,MUX) + CALL LCMPUT(IPTRK,'IMAX',LL4F,1,IMAX) + IF(IDIM.GE.2) THEN + CALL LCMPUT(IPTRK,'MUY',LL4F,1,MUY) + CALL LCMPUT(IPTRK,'IMAY',LL4F,1,IMAY) + CALL LCMPUT(IPTRK,'IPY',LL4F,1,IPY) + ENDIF + IF(IDIM.EQ.3) THEN + CALL LCMPUT(IPTRK,'MUZ',LL4F,1,MUZ) + CALL LCMPUT(IPTRK,'IMAZ',LL4F,1,IMAZ) + CALL LCMPUT(IPTRK,'IPZ',LL4F,1,IPZ) + ENDIF + IF(TITLE.NE.' ') CALL LCMPTC(IPTRK,'TITLE',72,TITLE) + TEXT12=HENTRY(2) + CALL LCMPTC(IPTRK,'LINK.GEOM',12,TEXT12) + IF(IMPX.GT.1) THEN + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + WRITE(6,110) ISTATE(1:2),ISTATE(4),ISTATE(6),ISTATE(12), + 1 ISTATE(14:16),ISTATE(25),ISTATE(27:29),ISTATE(33),ISTATE(39) + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(IPZ,IPY,IMAZ,IMAY,IMAX,MUZ,MUY,MUX) + DEALLOCATE(XX,ZZ,YY,ZZZ,YYY,XXX) + DEALLOCATE(IQFR,QFR,KN,VOL,IDL,MAT) + RETURN +* + 100 FORMAT(1H1,24HNN NN SSSSS SSSSS, + 1 97(1H*)/26H NNN NN SSSSSSS SSSSSSS, + 2 58(1H*),38H MULTIGROUP VERSION. A. HEBERT (2021)/ + 3 26H NNNN NN SS SS SS SS/26H NN NN NN SSS SSS / + 4 26H NN NN NN SSS SSS /26H NN NNNN SS SS SS SS/ + 5 26H NN NNN SSSSSSS SSSSSSS/26H NN NN SSSSS SSSSS // + 6 1X,A72//) + 110 FORMAT(/14H STATE VECTOR:/ + 1 7H NREG ,I8,22H (NUMBER OF REGIONS)/ + 2 7H NUN ,I8,23H (NUMBER OF UNKNOWNS)/ + 3 7H NMIX ,I8,23H (NUMBER OF MIXTURES)/ + 4 7H ITYPE ,I8,41H (TYPE OF GEOMETRY -- 2:1D; 5:2D; 7:3D)/ + 5 7H ICHX ,I8,40H (TYPE OF SOLUTION 4/5/6:CMFD/NEM/ANM)/ + 6 7H NX ,I8,40H (NUMBER OF ELEMENTS ALONG THE X AXIS)/ + 7 7H NY ,I8,40H (NUMBER OF ELEMENTS ALONG THE Y AXIS)/ + 8 7H NZ ,I8,40H (NUMBER OF ELEMENTS ALONG THE Z AXIS)/ + 9 7H LL4F ,I8,29H (NUMBER OF AVERAGE FLUXES)/ + 1 7H LL4X ,I8,38H (NUMBER OF X-DIRECTED NET CURRENTS)/ + 2 7H LL4Y ,I8,38H (NUMBER OF Y-DIRECTED NET CURRENTS)/ + 3 7H LL4Z ,I8,38H (NUMBER OF Z-DIRECTED NET CURRENTS)/ + 4 7H NADI ,I8,29H (NUMBER OF ADI ITERATIONS)/ + 5 7H IGMAX ,I8,47H (ENERGY GROUP LIMIT WITH HYPERBOLIC TRIAL FU, + 6 8HNCTIONS)) + 120 FORMAT(/24H NSST: VALUES OF VECTOR ,A6,4H ARE/(1X,1P,20I6)) + END |