summaryrefslogtreecommitdiff
path: root/Donjon/src/IDET01.f
diff options
context:
space:
mode:
Diffstat (limited to 'Donjon/src/IDET01.f')
-rw-r--r--Donjon/src/IDET01.f440
1 files changed, 440 insertions, 0 deletions
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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 22:41:11 +0000