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/TONSPH.f | 351 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 351 insertions(+) create mode 100644 Dragon/src/TONSPH.f (limited to 'Dragon/src/TONSPH.f') diff --git a/Dragon/src/TONSPH.f b/Dragon/src/TONSPH.f new file mode 100644 index 0000000..812ac62 --- /dev/null +++ b/Dragon/src/TONSPH.f @@ -0,0 +1,351 @@ +*DECK TONSPH + SUBROUTINE TONSPH(IPLIB,IPTRK,IFTRAK,NREG,NUN,NBM,NBISO,ISONAM, + 1 MAT,VOL,KEYFLX,CDOOR,INRS,LEAKSW,IMPX,DEN,MIX,LSHI,ITRANC, + 2 IPHASE,NGRO,IGRMIN,IGRMAX,NBNRS,TITR,SIGT2,SIGT3,SN,SPH,ICPIJ, + 3 TK3,TK4) +* +*----------------------------------------------------------------------- +* +*Purpose: +* SPH equivalence procedure over the self-shielded cross sections. Use +* all the standard solution doors of Dragon. +* +*Copyright: +* Copyright (C) 2017 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 file unit number used to store the tracks. +* NREG number of regions. +* NUN number of unknowns per energy group. +* NBM number of mixtures in the internal library. +* NBISO number of isotopes. +* ISONAM alias name of isotopes in IPLIB. +* MAT index-number of the mixture type assigned to each volume. +* VOL volumes. +* KEYFLX pointers of fluxes in unknown vector. +* CDOOR name of the geometry/solution operator. +* INRS index of the resonant isotope under consideration. +* LEAKSW leakage flag (LEAKSW=.TRUE. if neutron leakage through +* external boundary is present). +* IMPX print flag (equal to zero for no print). +* DEN density of each isotope. +* MIX mix number of each isotope (can be zero). +* LSHI resonant region number associated with each isotope. +* Infinite dilution will be assumed if LSHI(i)=0. +* ITRANC type of transport correction. +* IPHASE type of flux solution (=1 use a native flux solution door; +* =2 use collision probabilities). +* NGRO number of energy groups. +* IGRMIN first group where the self-shielding is applied. +* IGRMAX most thermal group where the self-shielding is applied. +* NBNRS number of totally correlated fuel regions. NBNRS=max(IRES). +* TITR title. +* SIGT2 total macroscopic cross sections. +* SIGT3 transport correction. +* SN computed dilution cross section in each energy group of +* each isotope. +* +*Parameters: output +* SPH SPH factors. +* ICPIJ number of flux solution door calls. +* +*Parameters: input/output +* TK3 cpu time to compute system matrices. +* TK4 cpu time to compute fluxes. +* +*----------------------------------------------------------------------- +* + USE GANLIB + USE DOORS_MOD +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPLIB,IPTRK + INTEGER IFTRAK,NREG,NUN,NBM,NBISO,ISONAM(3,NBISO),MAT(NREG), + 1 KEYFLX(NREG),INRS,IMPX,MIX(NBISO),LSHI(NBISO),ITRANC,IPHASE, + 2 NGRO,IGRMIN,IGRMAX,NBNRS,ICPIJ + REAL VOL(NREG),DEN(NBISO),SIGT2(0:NBM,NGRO),SIGT3(0:NBM,NGRO), + 1 SN(NGRO,NBISO),SPH(NBM,NGRO),TK3,TK4 + LOGICAL LEAKSW + CHARACTER CDOOR*12,TITR*72,HNAMIS*12 +*---- +* LOCAL VARIABLES +*---- + PARAMETER(NSTATE=40) + TYPE(C_PTR) JPLIB,KPLIB,IPMACR,IPSOU + LOGICAL LHOMOG,LPROB,LEXAC,LOGDO,REBFLG + CHARACTER TEXT12*12 + INTEGER NALBP,ISTATE(NSTATE) +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: IRES,ISONR,NPSYS + REAL, ALLOCATABLE, DIMENSION(:) :: VST,SIGTXS,SIGS0X,SIGG,FLNEW + REAL, ALLOCATABLE, DIMENSION(:,:) :: FLUX,SIG0,SIG1,SIG3,TOTAL, + 1 SIGS0,TRANC,PHGAR,SUNKNO,FUNKNO + LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKI + TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(IRES(NBM),ISONR(NBISO),NPSYS(NGRO)) + ALLOCATE(VST(NBNRS),SIGTXS(0:NBM),SIGS0X(0:NBM),SIGG(0:NBM), + 1 FLNEW(NBNRS),SUNKNO(NUN,NGRO),FUNKNO(NUN,NGRO),FLUX(NBM,NGRO), + 2 SIG0(NBM,NGRO),SIG1(NBM,NGRO),SIG3(NBM,NGRO),TOTAL(NGRO,NBNRS), + 3 SIGS0(NGRO,NBNRS),TRANC(NGRO,NBNRS),PHGAR(NGRO,NBNRS)) + ALLOCATE(MASKI(NBISO)) + ALLOCATE(IPISO(NBISO)) +*---- +* FIND THE RESONANT MIXTURE NUMBERS AND THE CORRELATED ISOTOPES +* ASSOCIATED WITH REGION INRS +*---- + IRES(:NBM)=0 + ISONR(:NBISO)=0 + IRS=0 + TEXT12=' ' + 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 + WRITE(HNAMIS,'(3A4)') (ISONAM(I0,ISO),I0=1,3) + IF(HNAMIS.NE.TEXT12) THEN + IRS=IRS+1 + TEXT12=HNAMIS + ENDIF + ISONR(ISO)=IRS + IRES(IBM)=IRS + ENDIF + 20 CONTINUE + 30 CONTINUE + IF(IRS.NE.NBNRS) CALL XABORT('TONSPH: INVALID VALUE OF NBNRS.') +*---- +* SET THE LCM MICROLIB ISOTOPEWISE DIRECTORIES. +*---- + CALL LIBIPS(IPLIB,NBISO,IPISO) +*---- +* UNLOAD MICROSCOPIC X-S FROM LCM TO SCRATCH STORAGE. +*---- + DO 40 ISO=1,NBISO + IRS=ISONR(ISO) + IF(IRS.GT.0) THEN + KPLIB=IPISO(ISO) ! set ISO-th isotope + CALL LCMGET(KPLIB,'NTOT0',TOTAL(1,IRS)) + CALL LCMGET(KPLIB,'SIGS00',SIGS0(1,IRS)) + DO IGRP=IGRMIN,IGRMAX +* Compute a ST flux for the homogeneous equivalent medium. + PHGAR(IGRP,IRS)=MAX(0.0,SN(IGRP,ISO)/(SN(IGRP,ISO)+ + 1 (TOTAL(IGRP,IRS)-SIGS0(IGRP,IRS)))) + ENDDO + IF(ITRANC.NE.0) CALL LCMGET(KPLIB,'TRANC',TRANC(1,IRS)) + ENDIF + 40 CONTINUE +*---- +* COMPUTE THE MERGED VOLUMES. +*---- + NALBP=0 + LHOMOG=.TRUE. + VST(:NBNRS)=0.0 + DO 50 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 50 + IND=IRES(IBM) + IF(IND.EQ.0) THEN + LHOMOG=.FALSE. + ELSE + VST(IND)=VST(IND)+VOL(I) + ENDIF + 50 CONTINUE + IF(LHOMOG.AND.(NBNRS.EQ.1)) GO TO 260 + IF(IMPX.GE.3) WRITE(6,'(37H TONSPH: SPH FACTOR CALCULATION (NBNR, + 1 2HS=,I5,1H)/)') NBNRS +*---- +* SET THE MIXTURE-DEPENDENT MACROSCOPIC XS. +*---- + FUNKNO(:NUN,:NGRO)=0.0 + SUNKNO(:NUN,:NGRO)=0.0 + NPSYS(:NGRO)=0 + CALL LCMSIX(IPLIB,'SHIBA',1) + JPLIB=LCMLID(IPLIB,'GROUP',NGRO) + DO 110 IGRP=IGRMIN,IGRMAX + NPSYS(IGRP)=IGRP +* +* COMPUTE THE LIGHT AND RESONANT COMPONENTS OF THE MACROSCOPIC +* CROSS SECTIONS IN EACH RESONANT MIXTURE. + DO 70 IBM=1,NBM + SIG0(IBM,IGRP)=0.0 + SIG1(IBM,IGRP)=0.0 + SIG3(IBM,IGRP)=SIGT3(IBM,IGRP) + 70 CONTINUE + DO 80 ISO=1,NBISO + IRS=ISONR(ISO) + IF(IRS.GT.0) THEN + IBM=MIX(ISO) + FLUX(IBM,IGRP)=PHGAR(IGRP,IRS) + SIGT2(IBM,IGRP)=SIGT2(IBM,IGRP)-TOTAL(IGRP,IRS)*DEN(ISO) + SIG0(IBM,IGRP)=TOTAL(IGRP,IRS)*DEN(ISO) + SIG1(IBM,IGRP)=SIGS0(IGRP,IRS)*DEN(ISO) + IF(ITRANC.NE.0) THEN + SIG3(IBM,IGRP)=SIGT3(IBM,IGRP)-TRANC(IGRP,IRS)*DEN(ISO) + ENDIF + ENDIF + 80 CONTINUE + IF(IMPX.GE.10) THEN + WRITE (6,400) IGRP,(SIG0(I,IGRP),I=1,NBM) + WRITE (6,410) IGRP,(SIG1(I,IGRP),I=1,NBM) + WRITE (6,420) IGRP,(SIGT2(I,IGRP),I=1,NBM) + WRITE (6,430) IGRP,(FLUX(I,IGRP),I=1,NBM) + ENDIF +*---- +* COMPUTE THE SOURCES. +*---- + SIGG(0)=0.0 + DO 90 IBM=1,NBM + SIGG(IBM)=SIGT2(IBM,IGRP) + IF(IRES(IBM).GT.0) THEN + SIGG(IBM)=SIGG(IBM)+FLUX(IBM,IGRP)*(SIG1(IBM,IGRP)- + > SIG0(IBM,IGRP)) + IF(.NOT.LHOMOG) SIGG(IBM)=SIGG(IBM)-FLUX(IBM,IGRP)* + > SIGT2(IBM,IGRP) + ENDIF + 90 CONTINUE + SUNKNO(:NUN,IGRP)=0.0 + CALL DOORS(CDOOR,IPTRK,NBM,0,NUN,SIGG,SUNKNO(1,IGRP)) +* + IF(NPSYS(IGRP).NE.0) THEN + ICPIJ=ICPIJ+1 + SIGTXS(0)=0.0 + SIGS0X(0)=0.0 + DO 100 IBM=1,NBM + IND=IRES(IBM) + IF((ITRANC.NE.0).AND.(IND.EQ.0)) THEN + SIGTXS(IBM)=SIGT2(IBM,IGRP)-SIG3(IBM,IGRP) + ELSE + SIGTXS(IBM)=SIGT2(IBM,IGRP) + ENDIF + IF(IND.EQ.0) THEN +* REMOVE TRANSPORT CORRECTION. + IF(ITRANC.NE.0) THEN + SIGS0X(IBM)=-SIG3(IBM,IGRP) + ELSE + SIGS0X(IBM)=0.0 + ENDIF + ELSE +* BELL ACCELERATION. + SIGTXS(IBM)=SIGTXS(IBM)+SIG0(IBM,IGRP) + SIGS0X(IBM)=SIGTXS(IBM) + IF(LHOMOG) SIGS0X(IBM)=SIGS0X(IBM)-SIGT2(IBM,IGRP) + ENDIF + 100 CONTINUE + KPLIB=LCMDIL(JPLIB,IGRP) + CALL LCMPUT(KPLIB,'DRAGON-TXSC',NBM+1,2,SIGTXS) + CALL LCMPUT(KPLIB,'DRAGON-S0XSC',NBM+1,2,SIGS0X) + ENDIF + 110 CONTINUE +*---- +* SOLVE FOR THE FLUX USING DIRECT SELF-SHIELDED CROSS SECTIONS +*---- + CALL KDRCPU(TKA) + ISTRM=1 + NANI=1 + NW=0 + KNORM=1 + IMPY=MAX(0,IMPX-3) + IF(IPHASE.EQ.1) THEN +* USE A NATIVE DOOR. + CALL DOORAV(CDOOR,JPLIB,NPSYS,IPTRK,IFTRAK,IMPY,NGRO,NREG, + 1 NBM,NANI,NW,MAT,VOL,KNORM,LEAKSW,TITR,NALBP,ISTRM) + ELSE IF(IPHASE.EQ.2) THEN +* USE A COLLISION PROBABILITY DOOR. + IPIJK=1 + ITPIJ=1 + CALL DOORPV(CDOOR,JPLIB,NPSYS,IPTRK,IFTRAK,IMPY,NGRO,NREG, + 1 NBM,NANI,MAT,VOL,KNORM,IPIJK,LEAKSW,ITPIJ,.FALSE.,TITR,NALBP) + ENDIF + CALL KDRCPU(TKB) + TK3=TK3+(TKB-TKA) + CALL KDRCPU(TKA) + IDIR=0 + LEXAC=.FALSE. + IPMACR=C_NULL_PTR + IPSOU=C_NULL_PTR + REBFLG=.FALSE. + CALL DOORFV(CDOOR,JPLIB,NPSYS,IPTRK,IFTRAK,IMPX,NGRO,NBM,IDIR, + 1 NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUNKNO,FUNKNO,IPMACR, + 2 IPSOU,REBFLG) + CALL LCMSIX(IPLIB,' ',2) + TK4=TK4+(TKB-TKA) +*---- +* HOMOGENIZE THE FLUX +*---- + DO 150 IGRP=IGRMIN,IGRMAX + IF(NPSYS(IGRP).NE.0) THEN + FLNEW(:NBNRS)=0.0 + DO 120 I=1,NREG + IF(MAT(I).EQ.0) GO TO 120 + IND=IRES(MAT(I)) + IF(IND.GT.0) FLNEW(IND)=FLNEW(IND)+FUNKNO(KEYFLX(I),IGRP)*VOL(I) + 120 CONTINUE + DO 130 IND=1,NBNRS + FLNEW(IND)=FLNEW(IND)/VST(IND) + 130 CONTINUE +*---- +* SPH FACTOR CONTROL +*---- + DO 140 IBM=1,NBM + IND=IRES(IBM) + IF(IND.GT.0) THEN + SPHNEW=PHGAR(IGRP,IND)/FLNEW(IND) + LPROB=(SPHNEW.LE.0.).OR.(SPHNEW.GT.1.).OR.(FLNEW(IND).LT.0.05) + IF(LPROB) SPHNEW=1.0 + SPH(IBM,IGRP)=SPHNEW + ENDIF + 140 CONTINUE + ENDIF + 150 CONTINUE +*---- +* SPH CORRECTION OF THE MICROLIB +*---- + CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE) + NL=ISTATE(4) + NED=ISTATE(13) + NDEL=ISTATE(19) + DO 160 ISO=1,NBISO + MASKI(ISO)=(ISONR(ISO).GT.0) + 160 CONTINUE + CALL TONCMI(IPLIB,IMPX,NBM,NBISO,NGRO,NL,NED,NDEL,MASKI,SPH) + IF(IMPX.GT.3) THEN + DO 170 IGRP=IGRMIN,IGRMAX + WRITE (6,440) IGRP,(SPH(IBM,IGRP),IBM=1,NBM) + 170 CONTINUE + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + 260 DEALLOCATE(IPISO) + DEALLOCATE(MASKI) + DEALLOCATE(PHGAR,TRANC,SIGS0,TOTAL,SIG3,SIG1,SIG0,FLUX,FUNKNO, + 1 SUNKNO,FLNEW,SIGG,SIGS0X,SIGTXS,VST) + DEALLOCATE(NPSYS,ISONR,IRES) + RETURN + 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 SCATTERING MACROSCOPIC CROSS SECTIONS OF THE OTHER, + 1 33H MATERIALS IN EACH MIXTURE (GROUP,I5,2H):/(1X,1P,11E11.3)) + 420 FORMAT(/51H TOTAL MACROSCOPIC CROSS SECTIONS OF THE OTHER MATE, + 1 28HRIALS IN EACH MIXTURE (GROUP,I5,2H):/(1X,1P,11E11.3)) + 430 FORMAT(/19H TABSN3 FLUX (GROUP,I5,2H):/(1X,1P,11E11.3)) + 440 FORMAT(/19H SPH FACTORS (GROUP,I5,2H):/(1X,1P,11E11.3)) + END -- cgit v1.2.3