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 --- Donjon/src/IDET01.f | 440 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 440 insertions(+) create mode 100644 Donjon/src/IDET01.f (limited to 'Donjon/src/IDET01.f') diff --git a/Donjon/src/IDET01.f b/Donjon/src/IDET01.f new file mode 100644 index 0000000..6cb4e56 --- /dev/null +++ b/Donjon/src/IDET01.f @@ -0,0 +1,440 @@ +*DECK IDET01 + SUBROUTINE IDET01(IPTRK,IPFLU,IPLIB,IPMAP,IMPX,NDETC,MAXNI,NINX, + > NINY,NINZ,COORD1,COORD2,COORD3,DETNAM,REANAM,ICORN,DETECT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute detector integrated response on Cartesian geometry +* +*Copyright: +* Copyright (C) 2019 Ecole Polytechnique de Montreal +* +*Author(s): +* A. Hebert +* +*Parameters: input +* IPTRK pointer to the tracking. +* IPFLU pointer to the finite-element flux. +* IPLIB pointer to the interpolated microlib. +* IPMAP pointer to the fuelmap. +* IMPX print parameter. +* NDETC number of detectors +* MAXNI first dimension of matrices NIN and COORD. +* NINX number of interpolation points per detector along x axis. +* NINY number of interpolation points per detector along y axis. +* NINZ number of interpolation points per detector along z axis. +* COORD1 interpolation points per detector along x axis. +* COORD2 interpolation points per detector along y axis. +* COORD3 interpolation points per detector along z axis. +* DETNAM character*12 alias name of the isotope used as detector. +* REANAM character*12 name of the nuclear reaction used as detector. +* ICORN flag to activate corner flux correction (0/1: ON/OFF). +* +*Parameters: output +* DETECT detector response. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK,IPFLU,IPLIB,IPMAP + INTEGER IMPX,NDETC,MAXNI,NINX(MAXNI),NINY(MAXNI),NINZ(MAXNI),ICORN + REAL COORD1(MAXNI,NDETC),COORD2(MAXNI,NDETC),COORD3(MAXNI,NDETC), + > DETECT(NDETC) + CHARACTER DETNAM*12,REANAM*12 +*---- +* LOCAL VARIABLES +*---- + INTEGER NSTATE + PARAMETER (NSTATE=40) + INTEGER ISTATE(NSTATE) + LOGICAL L3D + CHARACTER HSMG*131,CMODUL*12 + TYPE(C_PTR) JPFLU,JPLIB,KPLIB,IPMAC,JPMAC,KPMAC +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, DIMENSION(:), ALLOCATABLE :: MAT,KFLX,KN,IMIX + INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: MIXT + REAL, DIMENSION(:), ALLOCATABLE :: XX,YY,ZZ,XXX,YYY,ZZZ,MXD,MYD, + > MZD,FLXD,GAR,TERPX,TERPY,TERPZ + REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: TFLUX,AFLUX,DFLUX,SIGF + CHARACTER(LEN=12), DIMENSION(:), ALLOCATABLE :: HNAMIS + TYPE(C_PTR), DIMENSION(:), ALLOCATABLE :: IPISO +*---- +* RECOVER GENERAL TRACKING INFORMATION +*---- + IF(.NOT.C_ASSOCIATED(IPTRK)) CALL XABORT('IDET01: IPTRK not set.') + CALL LCMGTC(IPTRK,'TRACK-TYPE',12,CMODUL) + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + NREG=ISTATE(1) + NUN=ISTATE(2) + ITYPE=ISTATE(6) + NLF=0 + ICHX=0 + IDIM=1 + IF(ITYPE.EQ.5) THEN + IDIM=2 + ELSE IF(ITYPE.EQ.7) THEN + IDIM=3 + ELSE + CALL XABORT('IDET01: Cartesian geometry expected.') + ENDIF + IELEM=ISTATE(9) + L4=ISTATE(11) + ICHX=ISTATE(12) + NLF=ISTATE(30) + NXD=ISTATE(14) + NYD=ISTATE(15) + NZD=ISTATE(16) + LL4F=0 + LL4X=0 + LL4Y=0 + IF(CMODUL.EQ.'TRIVAC') THEN + LL4F=ISTATE(25) + LL4X=ISTATE(27) + LL4Y=ISTATE(28) + ENDIF + L3D=(NZD.GT.0) + IF(.NOT.L3D) CALL XABORT('IDET01: 3D geometry expected.') + NZD=MAX(1,NZD) + ALLOCATE(MAT(NREG),KFLX(NREG)) + CALL LCMGET(IPTRK,'MATCOD',MAT) + CALL LCMGET(IPTRK,'KEYFLX',KFLX) + ALLOCATE(MXD(NXD+1),MYD(NYD+1),MZD(NZD+1)) + ALLOCATE(XX(NREG),YY(NREG),ZZ(NREG)) + CALL LCMGET(IPTRK,'XX',XX) + CALL LCMGET(IPTRK,'YY',YY) + CALL LCMGET(IPTRK,'ZZ',ZZ) +*---- +* RECOVER FINITE-ELEMENT FLUX INFORMATION +*---- + IF(.NOT.C_ASSOCIATED(IPFLU)) CALL XABORT('IDET01: IPFLU not set.') + CALL LCMGET(IPFLU,'STATE-VECTOR',ISTATE) + NG=ISTATE(1) +*---- +* RECOVER RENUMBERED MIXTURE INDICES FROM FUELMAP +*---- + IF(C_ASSOCIATED(IPMAP)) THEN + CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE) + IF(ISTATE(4).NE.NG) CALL XABORT('IDET01: invalid group nb(1).') + CALL LCMSIX(IPMAP,'GEOMAP',1) + CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE) + NEL=NXD*NYD*NZD + IF(ISTATE(3).NE.NXD) CALL XABORT('IDET01: invalid NXD.') + IF(ISTATE(4).NE.NYD) CALL XABORT('IDET01: invalid NYD.') + IF(ISTATE(5).NE.NZD) CALL XABORT('IDET01: invalid NZD.') + IF(ISTATE(6).NE.NEL) CALL XABORT('IDET01: invalid NEL.') + IF(NREG.NE.NEL) CALL XABORT('IDET01: invalid NREG.') + CALL LCMSIX(IPMAP,' ',2) + CALL LCMGET(IPMAP,'BMIX',MAT) + ENDIF +*---- +* RECOVER MICROLIB INFORMATION +*---- + IF(.NOT.C_ASSOCIATED(IPLIB)) CALL XABORT('IDET01: IPLIB not set.') + CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE) + MAXMIX=ISTATE(1) + NBISO=ISTATE(2) + IF(ISTATE(3).NE.NG) CALL XABORT('IDET01: invalid group nb(2).') + NMIX=ISTATE(14) + ALLOCATE(IMIX(NBISO),HNAMIS(NBISO),IPISO(NBISO),GAR(NG)) + CALL LCMGET(IPLIB,'ISOTOPESMIX',IMIX) + CALL LCMGTC(IPLIB,'ISOTOPESUSED',12,NBISO,HNAMIS) + DO ISO=1,NBISO + IF(HNAMIS(ISO).EQ.DETNAM) GO TO 10 + ENDDO + DO ISO=1,NBISO + WRITE(6,'(5X,3H-->,A12)') HNAMIS(ISO) + ENDDO + WRITE(HSMG,'(48HIDET01: NO DETECTOR ISOTOPE FOUND IN MICROLIB WI, + > 8HTH NAME=,A12)') DETNAM + CALL XABORT(HSMG) + 10 CALL LIBIPS(IPLIB,NBISO,IPISO) + JPLIB=LCMGID(IPLIB,'ISOTOPESLIST') +*---- +* COMPUTE MESH FROM L_TRACK +*---- + ALLOCATE(XXX(NXD),YYY(NYD)) + XXX(:NXD)=0.0 + YYY(:NYD)=0.0 + IREG=0 + IF(L3D) THEN + ALLOCATE(ZZZ(NZD)) + ZZZ(:NZD)=0.0 + DO K=1,NZD + DO J=1,NYD + DO I=1,NXD + IREG=IREG+1 + IF(XX(IREG).NE.0.0) THEN + IF(XXX(I).EQ.0.0) THEN + XXX(I)=XX(IREG) + ELSE IF(ABS(XXX(I)-XX(IREG)).GT.1.0E-6) THEN + CALL XABORT('IDET01: inconsistent tracking in X') + ENDIF + ENDIF + IF(YY(IREG).NE.0.0) THEN + IF(YYY(J).EQ.0.0) THEN + YYY(J)=YY(IREG) + ELSE IF(ABS(YYY(J)-YY(IREG)).GT.1.0E-6) THEN + CALL XABORT('IDET01: inconsistent tracking in Y') + ENDIF + ENDIF + IF(ZZ(IREG).NE.0.0) THEN + IF(ZZZ(K).EQ.0.0) THEN + ZZZ(K)=ZZ(IREG) + ELSE IF(ABS(ZZZ(K)-ZZ(IREG)).GT.1.0E-6) THEN + CALL XABORT('IDET01: inconsistent tracking in Z') + ENDIF + ENDIF + ENDDO + ENDDO + ENDDO + ELSE + DO J=1,NYD + DO I=1,NXD + IREG=IREG+1 + IF(XX(IREG).NE.0.0) THEN + IF(XXX(I).EQ.0.0) THEN + XXX(I)=XX(IREG) + ELSE IF(ABS(XXX(I)-XX(IREG)).GT.1.0E-6) THEN + CALL XABORT('IDET01: inconsistent tracking in X') + ENDIF + ENDIF + IF(YY(IREG).NE.0.0) THEN + IF(YYY(J).EQ.0.0) THEN + YYY(J)=YY(IREG) + ELSE IF(ABS(YYY(J)-YY(IREG)).GT.1.0E-6) THEN + CALL XABORT('IDET01: inconsistent tracking in Y') + ENDIF + ENDIF + ENDDO + ENDDO + ENDIF + IF(IREG.NE.NREG) CALL XABORT('IDET01: invalid tracking') + MXD(1)=0.0 + MYD(1)=0.0 + MZD(1)=0.0 + DO I=1,NXD + MXD(I+1)=MXD(I)+XXX(I) + ENDDO + MYD(1)=0.0 + DO I=1,NYD + MYD(I+1)=MYD(I)+YYY(I) + ENDDO + MZD(1)=0.0 + IF(L3D) THEN + DO I=1,NZD + MZD(I+1)=MZD(I)+ZZZ(I) + ENDDO + DEALLOCATE(ZZZ) + ELSE + MZD(2)=0.0 + ENDIF + DEALLOCATE(YYY,XXX) +*---- +* PERFORM FLUX INTERPOLATION OVER DETECTOR LOCATIONS +*---- + IF(IMPX.GT.1) THEN + WRITE(6,'(/29H IDET01: DETECTOR INFORMATION)') + WRITE(6,'(5X,12HENERGY GROUP,1X,8HDETECTOR,2X,7HMIXTURE,5X, + 1 13HDETECTOR FLUX,3X,11HDONJON FLUX,5X,9HFLUX RATO,7X, + 2 11HDRAGON FLUX,5X,10HFISSION XS)') + ENDIF + ALLOCATE(FLXD(NUN)) + JPFLU=LCMGID(IPFLU,'FLUX') + DO I=1,NDETC + ININX=NINX(I) + ININY=NINY(I) + ININZ=NINZ(I) + ALLOCATE(TFLUX(ININX,ININY,ININZ,NG)) + DO IG=1,NG + CALL LCMGDL(JPFLU,IG,FLXD) + IF(ICHX.EQ.1) THEN +* Variational collocation method + CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM) + MKN=MAXKN/(NXD*NYD*NZD) + ALLOCATE(KN(MAXKN)) + CALL LCMGET(IPTRK,'KN',KN) + CALL LCMSIX(IPTRK,'BIVCOL',1) + CALL LCMLEN(IPTRK,'T',LC,ITYLCM) + CALL LCMGET(IPTRK,'E',E) + CALL LCMSIX(IPTRK,' ',2) + CALL VALUE2(LC,MKN,NXD,NYD,NZD,L4,COORD1(1,I),COORD2(1,I), + 1 COORD3(1,I),MXD,MYD,MZD,FLXD,MAT,KN,ININX,ININY,ININZ,E, + 2 TFLUX(1,1,1,IG)) + DEALLOCATE(KN) + ELSE IF(ICHX.EQ.2) THEN +* Raviart-Thomas finite element method + CALL VALUE4(IELEM,NUN,NXD,NYD,NZD,COORD1(1,I),COORD2(1,I), + 1 COORD3(1,I),MXD,MYD,MZD,FLXD,MAT,KFLX,ININX,ININY,ININZ, + 2 TFLUX(1,1,1,IG)) + ELSE IF(ICHX.EQ.3) THEN +* Nodal collocation method (MCFD) + CALL VALUE1(IDIM,NXD,NYD,NZD,L4,COORD1(1,I),COORD2(1,I), + 1 COORD3(1,I),MXD,MYD,MZD,FLXD,MAT,IELEM,ININX,ININY,ININZ, + 2 TFLUX(1,1,1,IG)) + ELSE IF(ICHX.EQ.6) THEN +* Analytic nodal method (ANM) + IF(IMPX.GT.0) WRITE(6,100) ICORN + IPMAC=LCMGID(IPLIB,'MACROLIB') + JPMAC=LCMGID(IPMAC,'GROUP') + CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM) + ALLOCATE(KN(MAXKN)) + CALL LCMGET(IPTRK,'KN',KN) + KPMAC=LCMGIL(JPMAC,IG) + CALL VALU5(KPMAC,NXD,NYD,NZD,LL4F,LL4X,LL4Y,NUN,NMIX, + 1 COORD1(1,I),COORD2(1,I),COORD3(1,I),MXD,MYD,MZD,FLXD,MAT, + 2 KFLX,KN,ININX,ININY,ININZ,ICORN,TFLUX(1,1,1,IG)) + DEALLOCATE(KN) + ELSE + CALL XABORT('IDET01: interpolation not implemented.') + ENDIF + ENDDO +*---- +* RECOVER AVERAGED FLUX FROM FINITE-ELEMENT CALCULATION +*---- + ALLOCATE(AFLUX(ININX,ININY,ININZ,NG),MIXT(ININX,ININY,ININZ)) + DO INX=1,ININX + NX=0 + DO IX=1,NXD + IF(COORD1(INX,I).LE.MXD(IX)) EXIT + NX=NX+1 + ENDDO + DO INY=1,ININY + NY=0 + DO IY=1,NYD + IF(COORD2(INY,I).LE.MXD(IY)) EXIT + NY=NY+1 + ENDDO + DO INZ=1,ININZ + NZ=0 + DO IZ=1,NZD + IF(COORD3(INZ,I).LE.MZD(IZ)) EXIT + NZ=NZ+1 + ENDDO + IF(NX*NY*NZ.EQ.0) THEN + WRITE(HSMG,'(38HIDET01: element not found for detector, + 1 I5,7h(1). x=,1p,e12.4,3H y=,e12.4,3H z=,e12.4)') I, + 2 COORD1(INX,I),COORD2(INY,I),COORD3(INZ,I) + CALL XABORT(HSMG) + ENDIF + IEL=(NZ-1)*NXD*NYD+(NY-1)*NXD+NX + IF(MAT(IEL).EQ.0) THEN + WRITE(HSMG,'(38HIDET01: element not found for detector, + 1 I5,7h(1). x=,1p,e12.4,3H y=,e12.4,3H z=,e12.4)') I, + 2 COORD1(INX,I),COORD2(INY,I),COORD3(INZ,I) + CALL XABORT(HSMG) + ENDIF + MIXT(INX,INY,INZ)=MAT(IEL) + IUN=KFLX(IEL) + IF(IUN.EQ.0) CALL XABORT('IDET01: flux not defined.') + DO IG=1,NG + CALL LCMGDL(JPFLU,IG,FLXD) + AFLUX(INX,INY,INZ,IG)=FLXD(IUN) + ENDDO + ENDDO + ENDDO + ENDDO +*---- +* RECOVER FLUX AND FISSION CROSS SECTION FROM MICROLIB +*---- + ALLOCATE(DFLUX(ININX,ININY,ININZ,NG),SIGF(ININX,ININY,ININZ,NG)) + DO INX=1,ININX + DO INY=1,ININY + DO INZ=1,ININZ + IBM=MIXT(INX,INY,INZ) + DFLUX(INX,INY,INZ,:NG)=0.0 + SIGF(INX,INY,INZ,:NG)=0.0 + DO ISO=1,NBISO + IF((HNAMIS(ISO).EQ.DETNAM).AND.(IMIX(ISO).EQ.IBM)) THEN + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMLEN(KPLIB,REANAM,LENGT,ITYLCM) + IF(LENGT.NE.NG) THEN + CALL LCMLIB(KPLIB) + WRITE(HSMG,'(23HIDET01: unable to find ,A,6H for i, + > 7Hsotope ,A,11H in mixture,I6)') REANAM,DETNAM,IBM + CALL XABORT(HSMG) + ENDIF + CALL LCMGET(KPLIB,'NWT0',GAR) + DFLUX(INX,INY,INZ,:NG)=GAR(:NG) + CALL LCMGET(KPLIB,REANAM,GAR) + SIGF(INX,INY,INZ,:NG)=GAR(:NG) + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO +*---- +* PRINT DETECTOR-DEPENDENT VALUES +*---- + IF(IMPX.GT.1) THEN + DO IG=1,NG + DO INZ=1,ININZ + DO INY=1,ININY + DO INX=1,ININX + IBM=MIXT(INX,INY,INZ) + TFLUX_I=TFLUX(INX,INY,INZ,IG) + AFLUX_I=AFLUX(INX,INY,INZ,IG) + DFLUX_I=DFLUX(INX,INY,INZ,IG) + SIGF_I=SIGF(INX,INY,INZ,IG) + WRITE(6,'(8X,3I9,1P,5E16.5)') IG,I,IBM,TFLUX_I, + > AFLUX_I,TFLUX_I/AFLUX_I,DFLUX_I,SIGF_I + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF +*---- +* COMPUTE DETECTOR RESPONSE +*---- + ALLOCATE(TERPX(ININX),TERPY(ININY),TERPZ(ININZ)) + IF(ININX.EQ.1) THEN + TERPX(1)=1.0 + ELSE + CALL ALTERI(.TRUE.,ININX,COORD1(1,I),COORD1(1,I), + > COORD1(ININX,I),TERPX) + ENDIF + IF(ININY.EQ.1) THEN + TERPY(1)=1.0 + ELSE + CALL ALTERI(.TRUE.,ININY,COORD2(1,I),COORD2(1,I), + > COORD2(ININY,I),TERPY) + ENDIF + IF(ININZ.EQ.1) THEN + TERPZ(1)=1.0 + ELSE + CALL ALTERI(.TRUE.,ININZ,COORD3(1,I),COORD3(1,I), + > COORD3(ININZ,I),TERPZ) + ENDIF +* integrate along axial direction + DETECT(I)=0.0 + DO IG=1,NG + ZNUM=0.0 + ZDEN=0.0 + DO INX=1,ININX + DO INY=1,ININY + DO INZ=1,ININZ + TRP=TERPX(INX)*TERPY(INY)*TERPZ(INZ) + ZNUM=ZNUM+TRP*TFLUX(INX,INY,INZ,IG)*SIGF(INX,INY,INZ,IG) + ZDEN=ZDEN+TRP + ENDDO + ENDDO + ENDDO + DETECT(I)=DETECT(I)+ZNUM/ZDEN + ENDDO + DEALLOCATE(TERPZ,TERPY,TERPX) + DEALLOCATE(SIGF,DFLUX,MIXT,AFLUX,TFLUX) + ENDDO +*---- +* RELEASE MEMORY +*---- + DEALLOCATE(FLXD) + DEALLOCATE(GAR,IPISO,HNAMIS,IMIX) + DEALLOCATE(XX,YY,ZZ,KFLX,MAT) + RETURN + 100 FORMAT(/46H IDET01: CORNER FLUX CORRECTION (0/1: OFF/ON)=,I3) + END -- cgit v1.2.3