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/SHISN3.f | 464 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 464 insertions(+) create mode 100644 Dragon/src/SHISN3.f (limited to 'Dragon/src/SHISN3.f') diff --git a/Dragon/src/SHISN3.f b/Dragon/src/SHISN3.f new file mode 100644 index 0000000..c82e61c --- /dev/null +++ b/Dragon/src/SHISN3.f @@ -0,0 +1,464 @@ +*DECK SHISN3 + SUBROUTINE SHISN3 (IPLIB,IPTRK,IFTRAK,LEVEL,NGRO,NBISO,NBM,NREG, + 1 NUN,CDOOR,INRS,NBNRS,IMPX,ISONAM,MIX,DEN,SN,SB,LSHI,IPHASE,MAT, + 2 VOL,KEYFLX,LEAKSW,TITR,START,SIGT,SIGT3,NOCONV,BIEFF,LGC) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Perform one multidimensional self-shielding iteration using the +* generalized Stamm'ler algorithm with Nordheim (PIC) approximation. +* +*Copyright: +* Copyright (C) 2004 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. +* LEVEL type of self-shielding model (=1 original Stamm'ler model +* with Nordheim approximation; =2 Stamm'ler model with Nordheim +* approximation and Riemann integration method). +* 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. +* INRS index of the resonant isotope under consideration. +* NBNRS number of totaly correlated resonant regions. +* 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 SHISN3 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: input/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 infinitedilution. +* 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. +* +*Reference: +* A. Hebert, Revisiting the Stamm'ler Self-Shielding Method, Presented +* at the 25th CNS annnual conference, June 6-9, Toronto, 2004. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + PARAMETER (NALPHA=9,NRAT=(NALPHA+1)/2) + TYPE(C_PTR) IPLIB,IPTRK + INTEGER IFTRAK,LEVEL,NGRO,NBISO,NBM,NREG,NUN,INRS,NBNRS,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) + CHARACTER CDOOR*12,TITR*72 + LOGICAL LEAKSW,START,NOCONV(NBM,NGRO),BIEFF,LGC +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) KPLIB + CHARACTER HSMG*131,CGRPNM*12 + LOGICAL LOGDO + REAL FACT(NALPHA),SIGX(NALPHA) + COMPLEX SUM +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: IRES,ISONR,NPSYS + REAL, ALLOCATABLE, DIMENSION(:) :: GAR,GAS,SIGE,VST,DIST,FUN,DILG + REAL, ALLOCATABLE, DIMENSION(:,:) :: SIG0,SIG3,TOTAL,DILUT,GC + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: PICX + LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKI + COMPLEX, ALLOCATABLE, DIMENSION(:,:) :: COEF,DENOM + COMPLEX*16, ALLOCATABLE, DIMENSION(:,:,:) :: XCOEF,XDENO + TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO +*---- +* DATA STATEMENTS +*---- + DATA FACT/0.01,0.03162278,0.1,0.3162278,1.0,3.162278,10.0, + 1 31.62278,100.0/ +*---- +* SCRATCH STORAGE ALLOCATION +* SIG0 macroscopic xs of the resonant isotopes as interpolated. +* SIG3 macroscopic transport correction. +*---- + ALLOCATE(IRES(NBM),ISONR(NBNRS),NPSYS(NGRO)) + ALLOCATE(SIG0(NBM,NGRO),SIG3(NBM,NGRO),TOTAL(NGRO,NBNRS), + 1 GAR(NGRO),GAS(NGRO),SIGE(NGRO),DILUT(NALPHA,NGRO), + 2 GC(NGRO,NBNRS),VST(NBNRS),DIST(NBNRS)) + ALLOCATE(PICX(NALPHA,NBNRS,NGRO)) + ALLOCATE(MASKI(NBISO)) + ALLOCATE(COEF(NRAT,NGRO),DENOM(NRAT,NGRO)) + ALLOCATE(XCOEF(NRAT,NBNRS,NGRO),XDENO(NRAT,NBNRS,NGRO)) + ALLOCATE(IPISO(NBISO)) +*---- +* FIND THE RESONANT MIXTURE NUMBERS AND THE CORRELATED ISOTOPES +* ASSOCIATED WITH REGION INRS +*---- + IRS=0 + DO 30 IBM=1,NBM + LOGDO=.FALSE. + DO 10 I=1,NREG + LOGDO=LOGDO.OR.(MAT(I).EQ.IBM) + 10 CONTINUE + IF(.NOT.LOGDO) GO TO 30 + DO 20 ISO=1,NBISO + IF((MIX(ISO).EQ.IBM).AND.(LSHI(ISO).EQ.INRS)) THEN + IRS=IRS+1 + ISONR(IRS)=ISO + GO TO 30 + ENDIF + 20 CONTINUE + 30 CONTINUE + IF(IRS.NE.NBNRS) CALL XABORT('SHISN3: INVALID VALUE OF NBNRS.') + IRES(:NBM)=0 + DO 40 IRS=1,NBNRS + ISO=ISONR(IRS) + IRES(MIX(ISO))=IRS + 40 CONTINUE +*---- +* SET THE LCM MICROLIB ISOTOPEWISE DIRECTORIES. +*---- + CALL LIBIPS(IPLIB,NBISO,IPISO) +*---- +* UNLOAD MICROSCOPIC X-S FROM LCM TO SCRATCH STORAGE. SET THE +* GOLDSTEIN-COHEN TO ONE IN LEVEL 2 CALCULATIONS. +*---- + DO 50 IRS=1,NBNRS + ISO=ISONR(IRS) + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMGET(KPLIB,'NTOT0',TOTAL(1,IRS)) + GC(:NGRO,IRS)=1.0 + 50 CONTINUE +* + VST(:NBNRS)=0.0 + DO 60 I=1,NREG + IF(MAT(I).EQ.0) GO TO 60 + IND=IRES(MAT(I)) + IF(IND.GT.0) VST(IND)=VST(IND)+VOL(I) + 60 CONTINUE +* + NPSYS(:NGRO)=0 + DO 110 LLL=1,NGRO + LOGDO=.FALSE. + DO 70 IRS=1,NBNRS + LOGDO=LOGDO.OR.NOCONV(MIX(ISONR(IRS)),LLL) + 70 CONTINUE + IF(LOGDO) THEN + NPSYS(LLL)=LLL +* +* COMPUTE THE LIGHT AND RESONANT COMPONENTS OF THE MACROSCOPIC +* CROSS SECTIONS IN EACH RESONANT MIXTURE. + DO 80 IRS=1,NBNRS + ISO=ISONR(IRS) + IBM=MIX(ISO) + SIGT(IBM,LLL)=SIGT(IBM,LLL)-TOTAL(LLL,IRS)*DEN(ISO) + 80 CONTINUE + DO 90 IBM=1,NBM + SIG0(IBM,LLL)=0.0 + SIG3(IBM,LLL)=SIGT3(IBM,LLL) + 90 CONTINUE + DO 100 IRS=1,NBNRS + ISO=ISONR(IRS) + SIG0(MIX(ISO),LLL)=TOTAL(LLL,IRS)*DEN(ISO) + SIG3(MIX(ISO),LLL)=0.0 + 100 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 + 110 CONTINUE +*---- +* SET UP VECTORS SIGE AND SB. +*---- + CALL LCMSIX(IPLIB,'SHIBA',1) +* + SIGE(:NGRO)=0.0 + ALLOCATE(FUN(NUN*NGRO)) + CALL LCMSIX(IPLIB,'--AVERAGE--',1) + CALL SHIDST(IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG,NUN, + 1 NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG0,SIGT,SIGT3(1,1), + 2 TITR,FUN,SIGE) + CALL LCMSIX(IPLIB,' ',2) + DO 130 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + DO 120 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 120 + IND=IRES(IBM) + IF(IND.GT.0) THEN + IOF=(LLL-1)*NUN+KEYFLX(I) + ISO=ISONR(IND) + IF(NOCONV(IBM,LLL)) SB(LLL,ISO)=FUN(IOF)/ + 1 SIG0(IBM,LLL) + ENDIF + 120 CONTINUE + ENDIF + 130 CONTINUE + DEALLOCATE(FUN) +*---- +* SET UP VECTORS DILUT AND SIGX. +*---- + IF(START) THEN +* USE A VERY CHEAP APPROXIMATION TO START ITERATIONS. + DO 145 LLL=1,NGRO + DO 140 IALP=1,NALPHA + SIGX(IALP)=0.0 + DILUT(IALP,LLL)=SIGE(LLL) + 140 CONTINUE + 145 CONTINUE + ELSE + AVDEN=0.0 + VOLTOT=0.0 + DO 150 IRS=1,NBNRS + AVDEN=AVDEN+DEN(ISONR(IRS))*VST(IRS) + VOLTOT=VOLTOT+VST(IRS) + 150 CONTINUE + AVDEN=AVDEN/VOLTOT + DO 160 IRS=1,NBNRS + ISO=ISONR(IRS) + DIST(IRS)=DEN(ISO)/AVDEN + 160 CONTINUE + DO 220 IALP=1,NALPHA + DO 190 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + SIGX(IALP)=FACT(IALP)*SIGE(LLL) + DO 170 IBM=1,NBM + SIG0(IBM,LLL)=0.0 + SIG3(IBM,LLL)=SIGT3(IBM,LLL) + 170 CONTINUE + DO 180 IRS=1,NBNRS + ISO=ISONR(IRS) + SIG0(MIX(ISO),LLL)=SIGX(IALP)*DIST(IRS) + SIG3(MIX(ISO),LLL)=0.0 + 180 CONTINUE + ENDIF + 190 CONTINUE + ALLOCATE(DILG(NGRO),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,NUN, + 1 NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG0,SIGT,SIG3(1,1), + 2 TITR,FUN,DILG) + CALL LCMSIX(IPLIB,' ',2) + DO 210 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + DILUT(IALP,LLL)=DILG(LLL) + DO 200 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 200 + IND=IRES(IBM) + IF(IND.GT.0) THEN + IOF=(LLL-1)*NUN+KEYFLX(I) + PICX(IALP,IND,LLL)=FUN(IOF)/SIG0(IBM,LLL) + ENDIF + 200 CONTINUE + ENDIF + 210 CONTINUE + DEALLOCATE(FUN,DILG) + 220 CONTINUE + ENDIF + CALL LCMSIX(IPLIB,' ',2) +*---- +* COMPUTE AVERAGE MACROSCOPIC DILUTION X-S USING A N-TERM RATIONAL +* APPROXIMATION. +*---- + DO 260 LLL=1,NGRO + IF(NPSYS(LLL).NE.0) THEN + DO 230 IALP=1,NALPHA + SIGX(IALP)=FACT(IALP)*SIGE(LLL) + 230 CONTINUE +* ********************************************************** + CALL SHIRAT(IMPX,NRAT,SIGX,DILUT(1,LLL),LLL,A,COEF(1,LLL), + 1 DENOM(1,LLL)) +* ********************************************************** +* +* COMPUTE THE PIC BASE POINTS FOR A N-TERM RATIONAL APPROXIMATION + IF(START) THEN + DO 245 IRS=1,NBNRS + DO 240 I=1,NRAT + XCOEF(I,IRS,LLL)=COEF(I,LLL) + XDENO(I,IRS,LLL)=DENOM(I,LLL) + 240 CONTINUE + 245 CONTINUE + ELSE + CALL SHIDIL(NRAT,NALPHA,NBNRS,COEF(1,LLL),DENOM(1,LLL), + 1 DILUT(1,LLL),PICX(1,1,LLL),SIGX,DIST,VST,IMPX,LLL, + 2 XCOEF(1,1,LLL),XDENO(1,1,LLL)) + ENDIF + IF(.NOT.START.AND.BIEFF) THEN + DO 255 IRS=1,NBNRS + ISO=ISONR(IRS) + IF(NOCONV(MIX(ISO),LLL)) THEN + SIGRES=TOTAL(LLL,IRS)*DEN(ISO) + IF(NBNRS.EQ.1) THEN + SUM=0.0 + DO 250 I=1,NRAT + SUM=SUM+COEF(I,LLL)/(SIGRES+DENOM(I,LLL)) + 250 CONTINUE + PX0=REAL(SUM) + ELSE + PX0=SB(LLL,ISO) + ENDIF + SB(LLL,ISO)=(1.0/PX0-SIGRES)/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) + ENDIF + ENDIF + 255 CONTINUE + ENDIF + ENDIF + 260 CONTINUE +*---- +* APPLY A GOLDSTEIN-COHEN CORRECTION SIMILAR TO THE CORRECTION USED +* IN SHISN2. +*---- + IF((.NOT.START).AND.LGC) THEN + DO 295 IRS=1,NBNRS + ISO=ISONR(IRS) + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMLEN(KPLIB,'NGOLD',LENGT,ITYLCM) + IF(LENGT.EQ.NGRO) THEN + IF(IMPX.GE.5) WRITE (6,390) (ISONAM(I0,ISO),I0=1,3) + CALL LCMGET(KPLIB,'NGOLD',GC(1,IRS)) + ENDIF + IF(LEVEL.EQ.2) GO TO 295 + DO 290 JSO=1,NBISO + IF((MIX(JSO).EQ.MIX(ISO)).AND.(JSO.NE.ISO).AND. + 1 (LSHI(JSO).NE.0)) THEN + KPLIB=IPISO(JSO) ! set JSO-th isotope + CALL LCMLEN(KPLIB,'NGOLD',LENGT,ITYLCM) + IF(LENGT.EQ.NGRO) THEN + CALL LCMGET(KPLIB,'SIGS00',GAS) + CALL LCMGET(KPLIB,'NGOLD',GAR) + DO 280 LLL=1,NGRO + IF((NOCONV(MIX(JSO),LLL)).AND.(GAR(LLL).NE.1.0)) THEN + DDD=(1.0-GAR(LLL))*GAS(LLL)*DEN(JSO) + DO 270 I=1,NRAT + XDENO(I,IRS,LLL)=XDENO(I,IRS,LLL)-DDD + 270 CONTINUE + ENDIF + 280 CONTINUE + ENDIF + ENDIF + 290 CONTINUE + 295 CONTINUE + ENDIF +*---- +* COMPUTE THE DILUTION PARAMETERS (BARN) FOR EACH RESONANT ISOTOPE IN +* RESONANT MIXTURE INRS +*---- + CALL SHIEQU(IPLIB,LEVEL,NGRO,NBISO,NBM,NBNRS,NRAT,MIX,ISONAM, + 1 NOCONV,ISONR,GC,COEF,DENOM,XCOEF,XDENO,DEN,IMPX,SN) +* + DO 320 LLL=1,NGRO + LOGDO=.FALSE. + DO 300 IRS=1,NBNRS + LOGDO=LOGDO.OR.NOCONV(MIX(ISONR(IRS)),LLL) + 300 CONTINUE + IF(LOGDO) THEN + DO 310 IRS=1,NBNRS + ISO=ISONR(IRS) + IBM=MIX(ISO) + IF(START.OR.(.NOT.BIEFF)) THEN + SB(LLL,ISO)=SN(LLL,ISO) + ELSE IF(SB(LLL,ISO).LT.0.97*SN(LLL,ISO)) THEN + WRITE (6,520)(ISONAM(I0,ISO),I0=1,3),SB(LLL,ISO)/ + 1 SN(LLL,ISO),0.97,LLL + SB(LLL,ISO)=0.97*SN(LLL,ISO) + ENDIF + SIGT(IBM,LLL)=SIGT(IBM,LLL)+TOTAL(LLL,IRS)*DEN(ISO) + 310 CONTINUE + ENDIF + 320 CONTINUE +*---- +* 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 330 ISO=1,NBISO + LOGDO=START.OR.(DEN(ISO).NE.0.) + MASKI(ISO)=(LSHI(ISO).EQ.INRS).AND.LOGDO + 330 CONTINUE + IMPX2=MAX(0,IMPX-1) + CALL LIBLIB (IPLIB,NBISO,MASKI,IMPX2) + DO 345 IRS=1,NBNRS + ISO=ISONR(IRS) + IBM=MIX(ISO) + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMGET(KPLIB,'NTOT0',GAR) + DO 340 LLL=1,NGRO + TOTAL(LLL,IRS)=TOTAL(LLL,IRS)-GAR(LLL) + SIGT(IBM,LLL)=SIGT(IBM,LLL)-DEN(ISO)*TOTAL(LLL,IRS) + 340 CONTINUE + 345 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(IPISO) + DEALLOCATE(XDENO,XCOEF) + DEALLOCATE(DENOM,COEF) + DEALLOCATE(MASKI) + DEALLOCATE(PICX) + DEALLOCATE(DIST,VST,GC,DILUT,SIGE,GAS,GAR,TOTAL,SIG3,SIG0) + DEALLOCATE(NPSYS,ISONR,IRES) + RETURN +* + 390 FORMAT(/53H SHISN3: GOLDSTEIN AND COHEN APPROXIMATION USED FOR I, + 1 8HSOTOPE ',3A4,2H'.) + 400 FORMAT(//51H TOTAL MACROSCOPIC CROSS SECTIONS OF THE RESONANT M, + 1 31HATERIALS IN EACH MIXTURE (GROUP,I5,2H):/(1X,1P,11E11.3)) + 410 FORMAT(//51H TOTAL MACROSCOPIC CROSS SECTIONS OF THE OTHER MATE, + 1 28HRIALS IN EACH MIXTURE (GROUP,I5,2H):/(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)) + 515 FORMAT(30HSHISN3: 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 SHISN3: THE L-J EQUIVALENCE FACTOR OF RESONANT ISOTOP, + 1 3HE ',3A4,18H' WAS CHANGED FROM,F6.3,3H TO,F5.2,9H IN GROUP,I4, + 2 1H.) + END -- cgit v1.2.3