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/AUTSPH.f | 329 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 329 insertions(+) create mode 100644 Dragon/src/AUTSPH.f (limited to 'Dragon/src/AUTSPH.f') diff --git a/Dragon/src/AUTSPH.f b/Dragon/src/AUTSPH.f new file mode 100644 index 0000000..65bfe26 --- /dev/null +++ b/Dragon/src/AUTSPH.f @@ -0,0 +1,329 @@ +*DECK AUTSPH + SUBROUTINE AUTSPH(IPLI0,IPTRK,IFTRAK,NREG,NUN,NBMIX,NBISO,NIRES, + 1 NL,NED,NDEL,HCAL,MAT,VOL,KEYFLX,CDOOR,LEAKSW,IMPX,DEN,MIX,IAPT, + 2 ITRANC,IPHASE,NGRP,MASKG,IREX,TITR,SIGGAR,UNGAR,PHGAR,STGAR, + 3 SFGAR,SSGAR,S0GAR,SAGAR,SDGAR,DELTAU,SPH) +* +*----------------------------------------------------------------------- +* +*Purpose: +* SPH equivalence procedure over the self-shielded cross sections. Use +* all the standard solution doors of Dragon. Autosecol specific version. +* +*Copyright: +* Copyright (C) 2023 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 +* IPLI0 pointer to the LCM object containing subgroup-related +* information. +* 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. +* NBMIX number of mixtures in the internal library. +* NBISO number of isotopes. +* NIRES number of correlated resonant isotopes. +* NL number of Legendre orders required in the calculation +* (NL=1 or higher). +* NED number of extra vector edits. +* NDEL number of delayed neutron precursor groups. +* HCAL name of the self-shielding calculation. +* 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. +* 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). +* IAPT resonant isotope index associated with isotope I. Mixed +* moderator if IAPT(I)=NIRES+1. Out-of-fuel isotope if +* IAPT(I)=0. +* ITRANC type of transport correction. +* IPHASE type of flux solution (=1 use a native flux solution door; +* =2 use collision probabilities). +* NGRP number of energy groups. +* MASKG energy group mask pointing on self-shielded groups. +* IREX fuel region index assigned to each mixture. Equal to zero +* in non-resonant mixtures or in mixtures not used. +* TITR title. +* SIGGAR macroscopic x-s of the non-resonant isotopes in each mixture: +* (*,*,*,1) total; (*,*,*,2) transport correction; +* (*,*,*,3) P0 scattering. +* UNGAR averaged fluxes. +* STGAR microscopic self-shielded total x-s. +* +*Parameters: input/output +* PHGAR uncorrected and SPH-corrected averaged fluxes. +* SFGAR uncorrected and SPH-corrected microscopic self-shielded +* fission x-s. +* SSGAR uncorrected and SPH-corrected microscopic +* self-shielded scattering x-s. +* S0GAR uncorrected and SPH-corrected microscopic +* transfer scattering x-s +* (isotope,secondary,primary). +* SAGAR uncorrected and SPH-corrected microscopic +* additional x-s. +* SDGAR uncorrected and SPH-corrected microscopic +* self-shielded delayed nu-sigf x-s. +* DELTAU lethargy width of each energy group. +* +*Parameters: output +* SPH SPH factors. +* +*----------------------------------------------------------------------- +* + USE GANLIB + USE DOORS_MOD +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPLI0,IPTRK + INTEGER IFTRAK,NREG,NUN,NBMIX,NBISO,NIRES,NL,NED,NDEL,MAT(NREG), + 1 KEYFLX(NREG),IMPX,MIX(NBISO),IAPT(NBISO),ITRANC,IPHASE,NGRP, + 2 IREX(NBMIX) + REAL VOL(NREG),DEN(NBISO),SIGGAR(NBMIX,0:NIRES,NGRP,3), + 1 UNGAR(NUN,NGRP),PHGAR(NIRES,NGRP),STGAR(NIRES,NGRP), + 2 SFGAR(NIRES,NGRP),SSGAR(NIRES,NL,NGRP),S0GAR(NIRES,NL,NGRP,NGRP), + 3 SAGAR(NIRES,NED,NGRP),SDGAR(NIRES,NDEL,NGRP),DELTAU(NGRP), + 4 SPH(NIRES,NGRP) + LOGICAL LEAKSW,MASKG(NGRP) + CHARACTER CDOOR*12,HCAL*12,TITR*72 +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) JPLI0,KPLI0,IPMACR,IPSOU + LOGICAL LHOMOG,LPROB,LTIT,LEXAC,REBFLG + INTEGER NALBP +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: NPSYS + REAL, ALLOCATABLE, DIMENSION(:) :: SIGTXS,SIGS0X,SIGG + REAL, ALLOCATABLE, DIMENSION(:,:) :: SUNKNO,FUNKNO,SIGTI + LOGICAL, ALLOCATABLE, DIMENSION(:) :: LVOL +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(NPSYS(NGRP)) + ALLOCATE(SIGTI(NBMIX,5),SIGTXS(0:NBMIX),SIGS0X(0:NBMIX), + 1 SIGG(0:NBMIX),SUNKNO(NUN,NGRP),FUNKNO(NUN,NGRP),LVOL(NREG)) +*---- +* SET LHOMOG. +*---- + NALBP=0 + LHOMOG=.TRUE. + DO 10 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 10 + IF(IREX(IBM).EQ.0) LHOMOG=.FALSE. + 10 CONTINUE + IF(LHOMOG) GO TO 260 +*---- +* SET MACRO CALCULATION +*---- + ICPIJ=0 + CALL KDRCPU(TK1) + CALL LCMSIX(IPLI0,'SHIBA_SG',1) + CALL LCMSIX(IPLI0,HCAL,1) + LTIT=.TRUE. + JPLI0=LCMLID(IPLI0,'GROUP',NGRP) +*---- +* LOOP OVER SELF-SHIELDED ENERGY GROUPS. +*---- + FUNKNO(:NUN,:NGRP)=0.0 + SUNKNO(:NUN,:NGRP)=0.0 + NPSYS(:NGRP)=0 + DO 80 IGRP=1,NGRP + IF(.NOT.MASKG(IGRP)) GO TO 80 + NPSYS(IGRP)=IGRP +*---- +* SET THE MIXTURE-DEPENDENT MACROSCOPIC XS. +*---- + SIGTI(:NBMIX,:5)=0.0 + DO 50 IBM=1,NBMIX + DO 40 IRES=0,NIRES + IF(IRES.EQ.0) THEN + SIGTI(IBM,1)=SIGTI(IBM,1)+SIGGAR(IBM,0,IGRP,1) + SIGTI(IBM,3)=SIGTI(IBM,3)+SIGGAR(IBM,0,IGRP,3) + IF(ITRANC.NE.0) SIGTI(IBM,2)=SIGTI(IBM,2)+ + 1 SIGGAR(IBM,0,IGRP,2) + ELSE IF((IRES.GT.0).AND.(IREX(IBM).NE.0)) THEN + DENN=0.0 + DO 20 ISO=1,NBISO + IF((IAPT(ISO).EQ.IRES).AND.(MIX(ISO).EQ.IBM)) DENN=DEN(ISO) + 20 CONTINUE + SIGTI(IBM,5)=SIGTI(IBM,5)+DENN*STGAR(IRES,IGRP) + DO 30 JGRP=1,NGRP + SIGTI(IBM,4)=SIGTI(IBM,4)+PHGAR(IRES,JGRP)*DENN* + 1 S0GAR(IRES,1,IGRP,JGRP)*DELTAU(JGRP)/DELTAU(IGRP) + 30 CONTINUE + ENDIF + 40 CONTINUE + 50 CONTINUE +*---- +* COMPUTE THE SOURCES. +*---- + SIGG(0)=0.0 + DO 55 IBM=1,NBMIX + SIGG(IBM)=SIGTI(IBM,3) + IF(IREX(IBM).GT.0) SIGG(IBM)=SIGG(IBM)+SIGTI(IBM,4) + 55 CONTINUE + CALL DOORS(CDOOR,IPTRK,NBMIX,0,NUN,SIGG,SUNKNO(1,IGRP)) + SIGG(0:NBMIX)=0.0 + DO 60 IBM=1,NBMIX + IF(IREX(IBM).GT.0) THEN + SIGG(IBM)=SIGG(IBM)-SIGTI(IBM,5) + IF(.NOT.LHOMOG) SIGG(IBM)=SIGG(IBM)-SIGTI(IBM,1) + ENDIF + 60 CONTINUE + CALL DOORS(CDOOR,IPTRK,NBMIX,0,NUN,SIGG,SUNKNO(1,IGRP), + 1 UNGAR(1,IGRP)) +* + IF(NPSYS(IGRP).NE.0) THEN + ICPIJ=ICPIJ+1 + SIGTXS(0:)=0.0 + SIGS0X(0:)=0.0 + DO 70 IBM=1,NBMIX + SIGTXS(IBM)=SIGTI(IBM,1)-SIGTI(IBM,2) + IND=IREX(IBM) + IF(IND.EQ.0) THEN +* REMOVE TRANSPORT CORRECTION. + SIGS0X(IBM)=-SIGTI(IBM,2) + ELSE IF(IND.GT.0) THEN +* BELL ACCELERATION. + SIGTXS(IBM)=SIGTXS(IBM)+SIGTI(IBM,5) + SIGS0X(IBM)=SIGTXS(IBM) + IF(LHOMOG) SIGS0X(IBM)=SIGS0X(IBM)-SIGTI(IBM,1) + ENDIF + 70 CONTINUE + KPLI0=LCMDIL(JPLI0,IGRP) + CALL LCMPUT(KPLI0,'DRAGON-TXSC',NBMIX+1,2,SIGTXS) + CALL LCMPUT(KPLI0,'DRAGON-S0XSC',NBMIX+1,2,SIGS0X) + ENDIF + 80 CONTINUE +*---- +* SOLVE FOR THE FLUX USING DIRECT SELF-SHIELDED CROSS SECTIONS +*---- + ISTRM=1 + NANI=1 + KNORM=1 + IMPY=MAX(0,IMPX-3) + IF(IPHASE.EQ.1) THEN +* USE A NATIVE DOOR. + NW=0 + CALL DOORAV(CDOOR,JPLI0,NPSYS,IPTRK,IFTRAK,IMPY,NGRP,NREG, + 1 NBMIX,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,JPLI0,NPSYS,IPTRK,IFTRAK,IMPY,NGRP,NREG, + 1 NBMIX,NANI,MAT,VOL,KNORM,IPIJK,LEAKSW,ITPIJ,.FALSE.,TITR, + 2 NALBP) + ENDIF + IDIR=0 + LEXAC=.FALSE. + IPMACR=C_NULL_PTR + IPSOU=C_NULL_PTR + REBFLG=.FALSE. + CALL DOORFV(CDOOR,JPLI0,NPSYS,IPTRK,IFTRAK,IMPX,NGRP,NBMIX,IDIR, + 1 NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUNKNO,FUNKNO,IPMACR, + 2 IPSOU,REBFLG) +*---- +* LOOP OVER THE RESONANT ISOTOPES. +*---- + LVOL(:NREG)=.FALSE. + SPH(:NGRP,:NIRES)=1.0 + DO 240 IRES=1,NIRES +*---- +* HOMOGENIZE THE FLUX +*---- + VOLMER=0.0 + DO 100 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 100 + DO 90 ISO=1,NBISO + IF((IAPT(ISO).EQ.IRES).AND.(MIX(ISO).EQ.IBM)) LVOL(I)=.TRUE. + 90 CONTINUE + IF(LVOL(I)) VOLMER=VOLMER+VOL(I) + 100 CONTINUE + DO 230 IGRP=1,NGRP + IF(NPSYS(IGRP).NE.0) THEN + FLNEW=0.0 + DO 110 I=1,NREG + IF(LVOL(I)) FLNEW=FLNEW+FUNKNO(KEYFLX(I),IGRP)*VOL(I) + 110 CONTINUE + FLNEW=FLNEW/VOLMER +*---- +* SPH FACTOR CONTROL. +*---- + SPHNEW=PHGAR(IRES,IGRP)/FLNEW + LPROB=(SPHNEW.LE.0.0).OR.(SPHNEW.GT.1.0).OR.(FLNEW.LT.0.05) + IF(LPROB) SPHNEW=1.0 + SPH(IRES,IGRP)=SPHNEW + ENDIF + IF(MASKG(IGRP)) THEN + SPHNEW=SPH(IRES,IGRP) + PHGAR(IRES,IGRP)=PHGAR(IRES,IGRP)/SPHNEW + SFGAR(IRES,IGRP)=SFGAR(IRES,IGRP)*SPHNEW + DO 170 IL=1,NL + IF(MOD(IL-1,2).EQ.0) THEN + SSGAR(IRES,IL,IGRP)=SSGAR(IRES,IL,IGRP)*SPHNEW+ + 1 STGAR(IRES,IGRP)*(1.0-SPHNEW) + ELSE + SSGAR(IRES,IL,IGRP)=0.0 + ENDIF + DO 160 JGRP=1,NGRP + IF(MOD(IL-1,2).EQ.0) THEN + IF(IGRP.EQ.JGRP) THEN + S0GAR(IRES,IL,IGRP,IGRP)=S0GAR(IRES,IL,IGRP,IGRP)* + 1 SPHNEW+STGAR(IRES,IGRP)*(1.0-SPHNEW) + ELSE + S0GAR(IRES,IL,JGRP,IGRP)=S0GAR(IRES,IL,JGRP,IGRP)*SPHNEW + ENDIF + ELSE + IF(IGRP.EQ.JGRP) THEN + S0GAR(IRES,IL,IGRP,IGRP)=S0GAR(IRES,IL,IGRP,IGRP)/ + 1 SPHNEW+STGAR(IRES,IGRP)*(1.0-1.0/SPHNEW) + ELSE + S0GAR(IRES,IL,JGRP,IGRP)=S0GAR(IRES,IL,JGRP,IGRP)/ + 1 SPH(IRES,JGRP) + ENDIF + ENDIF + IF(MOD(IL-1,2).EQ.1) THEN + SSGAR(IRES,IL,IGRP)=SSGAR(IRES,IL,IGRP)+ + 1 S0GAR(IRES,IL,JGRP,IGRP) + ENDIF + 160 CONTINUE + 170 CONTINUE + DO 180 IED=1,NED + SAGAR(IRES,IED,IGRP)=SAGAR(IRES,IED,IGRP)*SPHNEW + 180 CONTINUE + DO 190 IDEL=1,NDEL + SDGAR(IRES,IDEL,IGRP)=SDGAR(IRES,IDEL,IGRP)*SPHNEW + 190 CONTINUE + ENDIF + 230 CONTINUE + 240 CONTINUE +* *************************************************************** + CALL LCMSIX(IPLI0,' ',2) + CALL LCMSIX(IPLI0,' ',2) + CALL KDRCPU(TK2) + IF(IMPX.GT.1) WRITE(6,'(/34H AUTSPH: CPU TIME SPENT TO COMPUTE, + 1 18H THE SPH FACTORS =,F8.1,8H SECOND./9X,17HNUMBER OF ASSEMBL, + 2 15HY DOORS CALLS =,I5,1H.)') TK2-TK1,ICPIJ +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + 260 DEALLOCATE(LVOL,FUNKNO,SUNKNO,SIGG,SIGS0X,SIGTXS,SIGTI) + DEALLOCATE(NPSYS) + RETURN + END -- cgit v1.2.3