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/SHISN2.f | 435 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 435 insertions(+) create mode 100644 Dragon/src/SHISN2.f (limited to 'Dragon/src/SHISN2.f') diff --git a/Dragon/src/SHISN2.f b/Dragon/src/SHISN2.f new file mode 100644 index 0000000..edd4fee --- /dev/null +++ b/Dragon/src/SHISN2.f @@ -0,0 +1,435 @@ +*DECK SHISN2 + SUBROUTINE SHISN2 (IPLIB,IPTRK,IFTRAK,NGRO,NBISO,NBM,NREG,NUN, + 1 CDOOR,NRES,NBM2,IMPX,ISONAM,MIX,DEN,SN,SB,LSHI,IPHASE,MAT,VOL, + 2 KEYFLX,LEAKSW,TITR,START,SIGT,SIGT3,NOCONV,BIEFF,LGC,SIGE) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Perform one multidimensional self-shielding iteration using the +* generalized Stamm'ler algorithm without Nordheim (PIC) approximation. +* +*Copyright: +* Copyright (C) 2002 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 +* IPLIB pointer to the internal microscopic cross section library +* (L_LIBRARY signature). +* IPTRK pointer to the tracking (L_TRACK signature). +* IFTRAK unit number of the sequential binary tracking file. +* NGRO number of energy groups. +* NBISO number of isotopes present in the calculation domain. +* NBM number of mixtures in the macrolib. +* NREG number of volumes. +* NUN number of unknowns in the flux or source vector in one +* energy group. +* CDOOR name of the geometry/solution module. +* NRES number of resonant mixtures. +* NBM2 number of resonant isotopes. +* IMPX print flag. +* ISONAM alias name of isotopes. +* MIX mix number of each isotope (can be zero). +* DEN density of each isotope. +* LSHI resonant region number associated with each isotope. +* Infinite dilution will be assumed if LSHI(i)=0. +* IPHASE type of flux solution (=1 use a native flux solution door; +* =2 use collision probabilities). +* MAT index-number of the mixture type assigned to each volume. +* VOL volumes. +* KEYFLX pointers of fluxes in unknown vector. +* LEAKSW leakage flag (.TRUE. only if leakage is present on the outer +* surface). +* TITR title. +* START beginning-of-iteration flag (.TRUE. if SHISN2 is called +* for the first time). +* SIGT3 transport correction. +* NOCONV mixture convergence flag (.TRUE. if mixture IBM +* is not converged in group L). +* BIEFF Livolant-Jeanpierre normalization flag (.TRUE. to +* activate). +* LGC Goldstein-Cohen approximation flag (.TRUE. to activate). +* +*Parameters: output +* SN on input, estimate of the dilution cross section in each +* energy group of each isotope. A value of 1.0e10 is used +* for infinite dilution. +* On output, computed dilution cross section in each energy +* group of each isotope. +* SIGT total macroscopic cross sections as modified by Shiba. +* +*Parameters: output +* SB dilution cross section as used in Livolant-Jeanpierre +* normalization. +* SIGE computed macroscopic dilution cross section in each resonant +* mixture and each energy group. +* +*Reference: +* A. Hebert and G. Marleau, Generalization of the Stamm'ler Method +* for the Self-Shielding of Resonant isotopes in Arbitrary Geometries, +* Nucl. Sci. Eng. 108, 230 (1991). +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + PARAMETER (NALPHA=5) + TYPE(C_PTR) IPLIB,IPTRK + INTEGER IFTRAK,NGRO,NBISO,NBM,NREG,NUN,NRES,NBM2,IMPX, + 1 ISONAM(3,NBISO),MIX(NBISO),LSHI(NBISO),IPHASE,MAT(NREG), + 2 KEYFLX(NREG) + REAL DEN(NBISO),SN(NGRO,NBISO),SB(NGRO,NBISO),VOL(NREG), + 1 SIGT(NBM,NGRO),SIGT3(NBM,NGRO),SIGE(NRES,NGRO) + CHARACTER CDOOR*12,TITR*72,CGRPNM*12 + LOGICAL LEAKSW,START,NOCONV(NBM,NGRO),BIEFF,LGC +*---- +* LOCAL VARIABLES +*---- + CHARACTER HSMG*131 + LOGICAL LOGDO + COMPLEX COEF(3),DENOM(3),EAV + PARAMETER (NRAT=(NALPHA+1)/2) + TYPE(C_PTR) KPLIB + REAL FACT(NALPHA),SIGX(NALPHA) +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: IRES,MIX2,IRNBM,NPSYS + REAL, ALLOCATABLE, DIMENSION(:) :: GAR,SIGRES,DILAV,FUN + REAL, ALLOCATABLE, DIMENSION(:,:) :: SIG0,SIG1,SIG3,TOTAL,SIGOLD, + 1 DILUT + LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKI + TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO +*---- +* DATA STATEMENTS +*---- + DATA FACT/0.01,0.1,1.0,10.0,100.0/ +*---- +* SCRATCH STORAGE ALLOCATION +* SIG0 macroscopic xs of the resonant isotopes as interpolated. +* SIG1 macroscopic xs of the resonant isotopes at various SIGX. +* SIG3 macroscopic transport correction. +*---- + ALLOCATE(IRES(NBM),MIX2(NBISO),IRNBM(NBM),NPSYS(NGRO)) + ALLOCATE(SIG0(NBM,NGRO),SIG1(NBM,NGRO),SIG3(NBM,NGRO), + 1 TOTAL(NGRO,NBM2),SIGOLD(NGRO,NBM2),GAR(NGRO),SIGRES(NBM), + 2 DILAV(NGRO),DILUT(NALPHA,NGRO)) + ALLOCATE(MASKI(NBISO)) + ALLOCATE(IPISO(NBISO)) +*---- +* SET THE LCM MICROLIB ISOTOPEWISE DIRECTORIES. +*---- + CALL LIBIPS(IPLIB,NBISO,IPISO) +*---- +* UNLOAD MICROSCOPIC X-S FROM LCM TO SCRATCH STORAGE +*---- + IBM=0 + DO 20 ISO=1,NBISO + MIX2(ISO)=0 + IF(LSHI(ISO).GT.0) THEN + IBM=IBM+1 + MIX2(ISO)=IBM + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMGET(KPLIB,'NTOT0',TOTAL(1,IBM)) + CALL LCMLEN(KPLIB,'NGOLD',LENGT,ITYLCM) + IF((LENGT.EQ.NGRO).AND.(.NOT.START).AND.LGC) THEN + IF(IMPX.GE.5) WRITE (6,390) (ISONAM(I0,ISO),I0=1,3) + CALL LCMGET(KPLIB,'SIGS00',SIGOLD(1,IBM)) + CALL LCMGET(KPLIB,'NGOLD',GAR) + DO 10 LLL=1,NGRO + SIGOLD(LLL,IBM)=(1.0-GAR(LLL))*SIGOLD(LLL,IBM) + 10 CONTINUE + ELSE + SIGOLD(:NGRO,IBM)=0.0 + ENDIF + ENDIF + 20 CONTINUE +*---- +* LOOP OVER RESONANT REGIONS. THE CP ARE STORED ON DIRECTORY SHIBA +*---- + CALL LCMSIX(IPLIB,'SHIBA',1) + DO 260 INRS=1,NRES +*---- +* FIND THE RESONANT MIXTURE NUMBERS (IRNBM) ASSOCIATED WITH REGION INRS +*---- + NBNRS=0 + DO 50 IBM=1,NBM + IRES(IBM)=0 + DO 40 ISO=1,NBISO + IF((MIX(ISO).EQ.IBM).AND.(LSHI(ISO).EQ.INRS)) THEN + NBNRS=NBNRS+1 + IRNBM(NBNRS)=IBM + IRES(IBM)=1 + GO TO 50 + ENDIF + 40 CONTINUE + 50 CONTINUE + IF(NBNRS.EQ.0) THEN + IF(START.AND.(IMPX.GE.1)) WRITE(6,385) 'SHISN2',INRS + GO TO 260 + ELSE IF(START.AND.(NBNRS.GT.1).AND.(IMPX.GE.5)) THEN + WRITE (6,380) NBNRS,INRS + ENDIF +* + NPSYS(:NGRO)=0 + DO 120 LLL=1,NGRO + LOGDO=.FALSE. + DO 60 I=1,NBNRS + LOGDO=LOGDO.OR.NOCONV(IRNBM(I),LLL) + 60 CONTINUE + IF(LOGDO) THEN + NPSYS(LLL)=LLL +* +* COMPUTE THE LIGHT AND RESONANT COMPONENTS OF THE MACROSCOPIC +* CROSS SECTIONS IN EACH RESONANT MIXTURE. + DO 80 I=1,NBNRS + SIGRES(I)=0.0 + DO 70 ISO=1,NBISO + IF((MIX(ISO).EQ.IRNBM(I)).AND.(LSHI(ISO).EQ.INRS)) THEN + SIGRES(I)=SIGRES(I)+TOTAL(LLL,MIX2(ISO))*DEN(ISO) + ENDIF + 70 CONTINUE + SIGT(IRNBM(I),LLL)=SIGT(IRNBM(I),LLL)-SIGRES(I) + 80 CONTINUE + DO 90 IBM=1,NBM + SIG0(IBM,LLL)=0.0 + SIG1(IBM,LLL)=0.0 + SIG3(IBM,LLL)=SIGT3(IBM,LLL) + 90 CONTINUE + DO 110 I=1,NBNRS + SIG0(IRNBM(I),LLL)=SIGRES(I) + SIG3(IRNBM(I),LLL)=0.0 + 110 CONTINUE + IF(IMPX.GE.10) THEN + WRITE (6,400) LLL,(SIG0(I,LLL),I=1,NBM) + WRITE (6,410) LLL,(SIGT(I,LLL),I=1,NBM) + WRITE (6,420) LLL,(SIGT3(I,LLL),I=1,NBM) + ENDIF + ENDIF + 120 CONTINUE +*---- +* SET UP VECTORS DILUT AND SIGX. +*---- + DILAV(:NGRO)=0.0 + IF(START) THEN +* USE A VERY CHEAP APPROXIMATION TO START ITERATIONS. + ALLOCATE(FUN(NUN*NGRO)) + CALL LCMSIX(IPLIB,'--AVERAGE--',1) + CALL SHIDST (IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG, + 1 NUN,NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG0,SIGT, + 2 SIGT3(1,1),TITR,FUN,DILAV) + CALL LCMSIX(IPLIB,' ',2) + DEALLOCATE(FUN) + DO 135 LLL=1,NGRO + DO 130 IALP=1,NALPHA + DILUT(IALP,LLL)=DILAV(LLL) + 130 CONTINUE + 135 CONTINUE + ELSE + DO 165 IALP=1,NALPHA + DO 150 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + DO 140 I=1,NBNRS + SIG1(IRNBM(I),LLL)=FACT(IALP)*SIGE(INRS,LLL) + 140 CONTINUE + ENDIF + 150 CONTINUE + ALLOCATE(FUN(NUN*NGRO)) + WRITE(CGRPNM,'(8H--BAND--,I4.4)') IALP + CALL LCMSIX(IPLIB,CGRPNM,1) + CALL SHIDST (IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG, + 1 NUN,NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG1,SIGT, + 2 SIG3(1,1),TITR,FUN,DILAV) + CALL LCMSIX(IPLIB,' ',2) + DEALLOCATE(FUN) + DO 160 LLL=1,NGRO + DILUT(IALP,LLL)=DILAV(LLL) + 160 CONTINUE + 165 CONTINUE + ENDIF +*---- +* COMPUTE AVERAGE MACROSCOPIC DILUTION X-S (SIGE) USING A THREE-TERM +* RATIONAL APPROXIMATION. +*---- + DO 200 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + DO 170 IALP=1,NALPHA + SIGX(IALP)=FACT(IALP)*SIGE(INRS,LLL) + 170 CONTINUE + IMPX2=IMPX + IF(START) IMPX2=MAX(0,IMPX-10) +* ********************************************************** + CALL SHIRAT(IMPX2,NRAT,SIGX,DILUT(1,LLL),LLL,A,COEF,DENOM) +* ********************************************************** + EAV=(COEF(1)*SQRT(DENOM(1))+COEF(2)*SQRT(DENOM(2))+ + 1 COEF(3)*SQRT(DENOM(3)))**2 + SIGE(INRS,LLL)=REAL(EAV) + IF((.NOT.START).AND.(BIEFF).AND.(NBNRS.EQ.1)) THEN +* COMPUTE DILAV FOR THE L-J NORMALIZATION. + SIGXX=SIG0(IRNBM(1),LLL) + PXX=REAL(COEF(1)/(SIGXX+DENOM(1))+COEF(2)/(SIGXX+DENOM(2)) + 1 +COEF(3)/(SIGXX+DENOM(3))) + DILAV(LLL)=1.0/PXX-SIGXX + ENDIF +*---- +* COMPUTE THE ISOTOPE DILUTION MICROSCOPIC CROSS SECTIONS (SN) USED +* FOR LIBRARY INTERPOLATION. +*---- + DO 190 ISO=1,NBISO + IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND. + 1 (DEN(ISO).NE.0.)) THEN + SUM=0.0 + DO 180 JSO=1,NBISO + IBM=MIX(JSO) + IF((LSHI(JSO).EQ.INRS).AND.(IBM.EQ.MIX(ISO)).AND. + 1 (ISO.NE.JSO)) SUM=SUM+(TOTAL(LLL,MIX2(JSO))- + 2 SIGOLD(LLL,MIX2(JSO)))*DEN(JSO) + 180 CONTINUE + SN(LLL,ISO)=REAL(COEF(1)*SQRT(DENOM(1)+SUM)+COEF(2)* + 1 SQRT(DENOM(2)+SUM)+COEF(3)*SQRT(DENOM(3)+SUM))**2/DEN(ISO) + IF(SN(LLL,ISO).LE.0.0) THEN + WRITE (HSMG,510) (ISONAM(I0,ISO),I0=1,3),SN(LLL,ISO),LLL + CALL XABORT(HSMG) + ENDIF + ELSE IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND. + 1 (DEN(ISO).EQ.0.)) THEN + SN(LLL,ISO)=1.0E10 + ENDIF + 190 CONTINUE + IF((.NOT.START).AND.(IMPX.GE.10)) THEN + DO 195 I=1,NBNRS + PP=A-SIGT(IRNBM(I),LLL) + QQ=SIGE(INRS,LLL)-SIGT(IRNBM(I),LLL) + IF(ABS(PP).GT.1.0E-4*SIGT(IRNBM(I),LLL)) THEN + BEL=QQ/PP + ELSE + BEL=0.0 + ENDIF + WRITE (6,610) I,SIGE(INRS,LLL),BEL + 195 CONTINUE + ENDIF + ENDIF + 200 CONTINUE +*---- +* COMPUTE THE ISOTOPE DILUTION MICROSCOPIC CROSS SECTIONS (SB) USED +* IN L-J NORMALIZATION. +*---- + IF((.NOT.START).AND.(BIEFF).AND.(NBNRS.GT.1)) THEN +* COMPUTE DILAV FOR THE L-J NORMALIZATION. + ALLOCATE(FUN(NUN*NGRO)) + CALL LCMSIX(IPLIB,'--AVERAGE--',1) + CALL SHIDST (IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG, + 1 NUN,NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG0,SIGT, + 2 SIGT3(1,1),TITR,FUN,DILAV) + CALL LCMSIX(IPLIB,' ',2) + DEALLOCATE(FUN) + ENDIF + DO 250 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + DO 220 ISO=1,NBISO + IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND. + 1 (DEN(ISO).NE.0.)) THEN + SUM=0.0 + DO 210 JSO=1,NBISO + IBM=MIX(JSO) + IF((LSHI(JSO).EQ.INRS).AND.(IBM.EQ.MIX(ISO)).AND. + 1 (ISO.NE.JSO)) SUM=SUM+TOTAL(LLL,MIX2(JSO))*DEN(JSO) + 210 CONTINUE + IF(START.OR.(.NOT.BIEFF)) THEN + SB(LLL,ISO)=SN(LLL,ISO) + ELSE + SB(LLL,ISO)=(DILAV(LLL)+SUM)/DEN(ISO) + IF(SB(LLL,ISO).LT.0.0) THEN + WRITE (HSMG,515) (ISONAM(I0,ISO),I0=1,3),SB(LLL,ISO), + 1 LLL + CALL XABORT(HSMG) + ELSE IF(SB(LLL,ISO).LT.SN(LLL,ISO)) THEN + IF(SB(LLL,ISO).LT.0.99*SN(LLL,ISO)) WRITE (6,520) + 1 (ISONAM(I0,ISO),I0=1,3),SB(LLL,ISO)/SN(LLL,ISO),LLL + SB(LLL,ISO)=SN(LLL,ISO) + ENDIF + ENDIF + ELSE IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND. + 1 (DEN(ISO).EQ.0.)) THEN + SB(LLL,ISO)=1.0E10 + ENDIF + 220 CONTINUE +* +* RESTORE SIGT ARRAY. + DO 240 I=1,NBNRS + SIGRES(I)=0.0 + DO 230 ISO=1,NBISO + IF((MIX(ISO).EQ.IRNBM(I)).AND.(LSHI(ISO).EQ.INRS)) THEN + SIGRES(I)=SIGRES(I)+TOTAL(LLL,MIX2(ISO))*DEN(ISO) + ENDIF + 230 CONTINUE + SIGT(IRNBM(I),LLL)=SIGT(IRNBM(I),LLL)+SIGRES(I) + 240 CONTINUE + ENDIF + 250 CONTINUE + 260 CONTINUE + CALL LCMSIX(IPLIB,' ',2) +*---- +* SAVE THE GROUP- AND ISOTOPE-DEPENDENT DILUTIONS +*---- + CALL LCMPUT(IPLIB,'ISOTOPESDSB',NBISO*NGRO,2,SB) + CALL LCMPUT(IPLIB,'ISOTOPESDSN',NBISO*NGRO,2,SN) +*---- +* COMPUTE THE SELF-SHIELDED MICROSCOPIC CROSS SECTIONS AND UPDATE +* VECTOR SIGT +*---- + DO 290 ISO=1,NBISO + LOGDO=START.OR.(DEN(ISO).NE.0.) + MASKI(ISO)=(LSHI(ISO).GT.0).AND.LOGDO + 290 CONTINUE + IMPX2=MAX(0,IMPX-1) + CALL LIBLIB (IPLIB,NBISO,MASKI,IMPX2) + DO 320 ISO=1,NBISO + IBM=MIX(ISO) + IF((LSHI(ISO).GT.0).AND.(IBM.GT.0).AND.(DEN(ISO).NE.0.)) THEN + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMGET(KPLIB,'NTOT0',GAR) + DO 300 LLL=1,NGRO + TOTAL(LLL,MIX2(ISO))=TOTAL(LLL,MIX2(ISO))-GAR(LLL) + 300 CONTINUE + DO 310 LLL=1,NGRO + IF(NOCONV(IBM,LLL)) SIGT(IBM,LLL)=SIGT(IBM,LLL)-DEN(ISO)* + 1 TOTAL(LLL,MIX2(ISO)) + 310 CONTINUE + ENDIF + 320 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(IPISO) + DEALLOCATE(MASKI) + DEALLOCATE(DILUT,DILAV,SIGRES,GAR,SIGOLD,TOTAL,SIG3,SIG1,SIG0) + DEALLOCATE(NPSYS,IRNBM,MIX2,IRES) + RETURN +* + 380 FORMAT(/16H SHISN2: MERGING,I3,30H RESONANT MIXTURES IN RESONANT, + 1 14H REGION NUMBER,I3,1H.) + 385 FORMAT(A6,1X,': RESONANT REGION =',I10,1X,'NOT USED.') + 390 FORMAT(/53H SHISN2: GOLDSTEIN AND COHEN APPROXIMATION USED FOR I, + 1 8HSOTOPE ',3A4,2H'.) + 400 FORMAT(1X,'TOTAL MACROSCOPIC CROSS SECTIONS OF THE RESONANT ', + 1'MATERIALS IN EACH MIXTURE (GROUP',I5,'):'/(1X,1P,11E11.3)) + 410 FORMAT(1X,'TOTAL MACROSCOPIC CROSS SECTIONS OF THE OTHER ', + 1'MATERIALS IN EACH MIXTURE (GROUP',I5,'):'/(1X,1P,11E11.3)) + 420 FORMAT(//1X,'TRANSPORT CORRECTION CROSS SECTIONS OF THE OTHER ', + 1'MATERIALS IN EACH MIXTURE (GROUP',I5,'):'/(1X,1P,11E11.3)) + 510 FORMAT(30HSHISN2: THE RESONANT ISOTOPE ',3A4,14H' HAS A NEGATI, + 1 27HVE DILUTION CROSS-SECTION (,1P,E14.4,0P,10H) IN GROUP,I4,1H.) + 515 FORMAT(30HSHISN2: THE RESONANT ISOTOPE ',3A4,14H' HAS A NEGATI, + 1 22HVE L-J CROSS-SECTION (,1P,E14.4,0P,10H) IN GROUP,I4,1H.) + 520 FORMAT(54H SHISN2: THE L-J EQUIVALENCE FACTOR OF RESONANT ISOTOP, + 1 3HE ',3A4,18H' WAS CHANGED FROM,F6.3,16H TO 1.0 IN GROUP,I4,1H.) + 610 FORMAT(8X,8HAVERAGE(,I2,1H),1P,E13.5/8X,11HBELL FACTOR,E13.5) + END -- cgit v1.2.3