diff options
Diffstat (limited to 'Dragon/src/USSEXD.f')
| -rw-r--r-- | Dragon/src/USSEXD.f | 377 |
1 files changed, 377 insertions, 0 deletions
diff --git a/Dragon/src/USSEXD.f b/Dragon/src/USSEXD.f new file mode 100644 index 0000000..31f7622 --- /dev/null +++ b/Dragon/src/USSEXD.f @@ -0,0 +1,377 @@ +*DECK USSEXD + SUBROUTINE USSEXD(MAXNOR,CDOOR,IPLI0,IPTRK,IFTRAK,IMPX,NGRP,IG, + 1 IASM,NBMIX,NREG,NUN,IPHASE,MAT,VOL,KEYFLX,IREX,SIGGAR,TITR,NIRES, + 2 IRES,NBNRS,MRANK,CONR,GOLD,IPPT1,IPPT2,VOLMER,XFLUX,UNGAR) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solution of the flux for the resonance spectrum expansion (RSE) method +* using the response matrix method. This is a non-iterative approach +* which is useful in exceptional cases where the fixed-point approach +* fails. +* +*Copyright: +* Copyright (C) 2024 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 +* MAXNOR maximum order of the probability tables (RSE). +* CDOOR name of the geometry/solution operator. +* IPLI0 pointer to the internal microscopic cross section library +* builded by the self-shielding module. +* IPTRK pointer to the tracking (L_TRACK signature). +* IFTRAK file unit number used to store the tracks. +* IMPX print flag (equal to zero for no print). +* NGRP number of energy groups. +* IG index of energy group being processed. +* IASM offset of information computed by DOORAV or DOORPV. +* NBMIX number of mixtures in the internal library. +* NREG number of regions. +* NUN number of unknowns in the flux or source vector in one +* energy group and one band. +* 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. +* IREX fuel region index assigned to each mixture. Equal to zero +* in non-resonant mixtures or in mixtures not used. +* SIGGAR macroscopic x-s of the non-resonant isotopes in each mixture: +* (*,*,*,1) total; (*,*,*,2) transport correction; +* (*,*,*,3) P0 scattering. +* TITR title. +* NIRES exact number of correlated resonant isotopes. +* IRES index of the resonant isotope being processed. +* NBNRS number of correlated fuel regions. +* MRANK exact order of the probability table. +* CONR number density of the resonant isotopes. +* GOLD type of self-shielding model (=1.0 physical probability +* tables; =-1001.0 resonance spectrum expansion method). +* IPPT1 pointer to LCM directory of each resonant isotope. +* IPPT2 information related to each resonant isotope: +* IPPT2(:,1) index of a resonant region (used with infinite +* dilution case); +* IPPT2(:,2:4) alias name of resonant isotope. +* VOLMER volumes of the resonant and non-resonant regions. +* +*Parameters: input/output +* XFLUX subgroup flux. +* +*Parameters: output +* UNGAR averaged fluxes per volume. +* +*----------------------------------------------------------------------- +* + USE GANLIB + USE DOORS_MOD +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPLI0,IPTRK,IPPT1(NIRES) + INTEGER MAXNOR,IFTRAK,IMPX,NGRP,IG,IASM,NBMIX,NREG,NUN,IPHASE, + 1 MAT(NREG),KEYFLX(NREG),IREX(NBMIX),NIRES,IRES,NBNRS,MRANK(NGRP), + 2 IPPT2(NIRES,4) + REAL VOL(NREG),SIGGAR(NBMIX,0:NIRES,NGRP,3),CONR(NBNRS,NIRES), + 1 GOLD(NIRES,NGRP),VOLMER(0:NBNRS), + 2 XFLUX(NBNRS,MAXNOR,NGRP),UNGAR(NREG,NIRES,NGRP) + CHARACTER CDOOR*12,TITR*72 +*---- +* LOCAL VARIABLES +*---- + TYPE(C_PTR) IPSYS,KPSYS,IPMACR,IPSOU,IPLIB,JPLIB1,KPLIB,IOFSET + DOUBLE PRECISION QQQ,SSS,T1 + LOGICAL LEXAC,REBFLG + CHARACTER CBDPNM*12,TEXT12*12 + TYPE VECTOR_ARRAY + DOUBLE PRECISION, POINTER, DIMENSION(:) :: VECTOR + END TYPE VECTOR_ARRAY + TYPE(VECTOR_ARRAY) :: WEIGHT_V,GAMMA_V +*---- +* ALLOCATABLE ARRAYS +*---- + TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: JPLIB2 + INTEGER, ALLOCATABLE, DIMENSION(:) :: NPSYS + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NJJ + REAL, ALLOCATABLE, DIMENSION(:) :: AWPHI,FUN,SUN,SIGG + REAL, ALLOCATABLE, DIMENSION(:,:) :: PAV + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: MATRIX + TYPE MATRIX_ARRAY + DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: MATRIX + END TYPE MATRIX_ARRAY + TYPE(MATRIX_ARRAY), ALLOCATABLE, DIMENSION(:,:) :: SCAT_M +*---- +* STATEMENT FUNCTIONS +*---- + INM(IND,IM,NBNRS)=(IM-1)*NBNRS+IND +*---- +* SCRATCH STORAGE ALLOCATION +*---- + MI=MRANK(IG) + ALLOCATE(NJJ(NGRP,NIRES)) + ALLOCATE(PAV(0:NBNRS,0:NBNRS),AWPHI(0:NBNRS)) + ALLOCATE(MATRIX(NBNRS*MI,NBNRS*MI+1)) + ALLOCATE(JPLIB2(NIRES),SCAT_M(NGRP,NIRES),SIGG(0:NBMIX)) +*---- +* RECOVER RSE INFORMATION FROM MICROLIB +*---- + IPLIB=IPPT1(IRES) + JPLIB1=LCMGID(IPLIB,'GROUP-RSE') + DO JRES=1,NIRES + WRITE(TEXT12,'(3A4)') (IPPT2(JRES,I),I=2,4) + CALL LCMSIX(IPLIB,TEXT12,1) + JPLIB2(JRES)=LCMGID(IPLIB,'SCAT_M') ! SCAT_M information + CALL LCMGET(IPLIB,'NJJS00',NJJ(:NGRP,JRES)) + CALL LCMSIX(IPLIB,' ',2) + ENDDO + KPLIB=LCMGIL(JPLIB1,IG) + CALL LCMGPD(KPLIB,'WEIGHT_V',IOFSET) + CALL C_F_POINTER(IOFSET,WEIGHT_V%VECTOR,(/MI/)) + CALL LCMGPD(KPLIB,'GAMMA_V',IOFSET) + CALL C_F_POINTER(IOFSET,GAMMA_V%VECTOR,(/MI/)) + DO JRES=1,NIRES + IPOS=1 + DO JG=1,IG-1 + IPOS=IPOS+NJJ(JG,JRES) + ENDDO + DO JG=IG-NJJ(IG,JRES)+1,IG + MJ=MRANK(JG) + CALL LCMGPL(JPLIB2(JRES),IPOS+IG-JG,IOFSET) + CALL C_F_POINTER(IOFSET,SCAT_M(JG,JRES)%MATRIX,(/MI,MJ/)) + ENDDO + ENDDO +*---- +* RECOVER THE SPECIFIC DIRECTORY FOR IRES-TH RESONANT ISOTOPE. +*---- + WRITE(CBDPNM,'(3HCOR,I4.4,1H/,I4.4)') IRES,NIRES + CALL LCMSIX(IPLI0,CBDPNM,1) + IPSYS=LCMGID(IPLI0,'ASSEMB-RSE') + CALL LCMSIX(IPLI0,' ',2) +*---- +* COMPUTE THE AVERAGED COLLISION PROBABILITY MATRIX. +*---- + ALLOCATE(NPSYS(MI*(NBNRS+1))) + ALLOCATE(FUN(NUN*MI*(NBNRS+1)),SUN(NUN*MI*(NBNRS+1))) + FUN(:NUN*MI*(NBNRS+1))=0.0 + SUN(:NUN*MI*(NBNRS+1))=0.0 + DO 50 IM=1,MI + DO 40 JNBN=0,NBNRS + NPSYS((IM-1)*(NBNRS+1)+JNBN+1)=IASM+IM + T1=0.0D0 + DO 10 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 10 + IND=IREX(IBM) + IF((JNBN.EQ.0).AND.(IND.EQ.0)) THEN + SSS=SIGGAR(IBM,0,IG,3)*GAMMA_V%VECTOR(IM) + T1=T1+SSS*VOL(I) + ELSE IF(IND.EQ.JNBN) THEN + T1=T1+VOL(I) + ENDIF + 10 CONTINUE + IOF=(IM-1)*NUN*(NBNRS+1)+JNBN*NUN + SIGG(0:NBMIX)=0.0 + DO 20 IBM=1,NBMIX + IND=IREX(IBM) + IF((JNBN.EQ.0).AND.(IND.EQ.0)) THEN + SSS=SIGGAR(IBM,0,IG,3)*GAMMA_V%VECTOR(IM) + SIGG(IBM)=REAL(SSS,4) + ELSE IF(IND.EQ.JNBN) THEN + SIGG(IBM)=1.0 + ENDIF + 20 CONTINUE + CALL DOORS(CDOOR,IPTRK,NBMIX,0,NUN,SIGG,SUN(IOF+1)) + DO 30 I=1,NUN + IF(T1.NE.0.0) SUN(IOF+I)=SUN(IOF+I)/REAL(T1,4) + 30 CONTINUE + 40 CONTINUE + 50 CONTINUE +*---- +* SOLVE FOR THE MULTIBAND FLUX. +*---- + IDIR=0 + NABS=MI*(NBNRS+1) + LEXAC=.FALSE. + IPMACR=C_NULL_PTR + IPSOU=C_NULL_PTR + REBFLG=.FALSE. + CALL DOORFV(CDOOR,IPSYS,NPSYS,IPTRK,IFTRAK,IMPX,NABS,NBMIX, + 1 IDIR,NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUN,FUN,IPMACR, + 2 IPSOU,REBFLG) +*---- +* HOMOGENIZE THE MULTIBAND FLUX. +*---- + DO 100 IM=1,MI + PAV(0:NBNRS,0:NBNRS)=0.0 + DO 70 JNBN=0,NBNRS + DO 60 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 60 + IOF=(IM-1)*NUN*(NBNRS+1)+JNBN*NUN+KEYFLX(I)-1 + PAV(IREX(IBM),JNBN)=PAV(IREX(IBM),JNBN)+FUN(IOF+1)*VOL(I) + 60 CONTINUE + 70 CONTINUE + DO 90 I=0,NBNRS + DO 80 J=0,NBNRS + IF(VOLMER(I).NE.0.0) PAV(I,J)=PAV(I,J)*VOLMER(J)/VOLMER(I) + 80 CONTINUE + 90 CONTINUE + KPSYS=LCMGIL(IPSYS,IASM+IM) + CALL LCMPUT(KPSYS,'DRAGON-PAV',(NBNRS+1)**2,2,PAV(0,0)) + 100 CONTINUE + DEALLOCATE(SUN,FUN,NPSYS) +*---- +* RESPONSE MATRIX APPROACH. LOOP OVER THE SECONDARY SUBGROUPS. +*---- + MATRIX(:NBNRS*MI,:NBNRS*MI+1)=0.0D0 + DO 200 IM=1,MI + KPSYS=LCMGIL(IPSYS,IASM+IM) + CALL LCMGET(KPSYS,'DRAGON-PAV',PAV(0,0)) +*---- +* LOOP OVER THE PRIMARY SUBGROUPS. MI+1 IS THE SOURCE. +*---- + DO 190 JM=1,MI+1 + IF(JM.LE.MI) THEN + JNBMAX=NBNRS + ELSE + JNBMAX=1 + ENDIF + DO 180 JNBN=1,JNBMAX + AWPHI(1:NBNRS)=0.0 + DO 160 I=1,NREG + IBM=MAT(I) + IF(IBM.EQ.0) GO TO 160 + JND=IREX(IBM) + QQQ=0.0D0 + IF(JM.EQ.MI+1) THEN + QQQ=QQQ+SIGGAR(IBM,0,IG,3)*GAMMA_V%VECTOR(IM) + IF(JND.NE.0) THEN + DO 130 JRES=1,NIRES + DENSIT=CONR(JND,JRES) + DO 120 JG=IG-NJJ(IG,JRES)+1,IG-1 + IF(GOLD(IRES,JG).NE.-1001.) CYCLE + DO 110 KM=1,MRANK(JG) + QQQ=QQQ+DENSIT*SCAT_M(JG,JRES)%MATRIX(IM,KM)* + 1 XFLUX(JND,KM,JG) + 110 CONTINUE + 120 CONTINUE + 130 CONTINUE + ENDIF + ELSE IF((JND.EQ.JNBN).AND.(JM.NE.IM)) THEN + DO 140 JRES=1,NIRES + DENSIT=CONR(JND,JRES) + IF(GOLD(IRES,IG).NE.-1001.) CYCLE + IF(JM.EQ.IM) CYCLE + QQQ=QQQ-DENSIT*SCAT_M(IG,JRES)%MATRIX(IM,JM) + 140 CONTINUE + ENDIF + DO 150 IND=1,NBNRS + AWPHI(IND)=AWPHI(IND)+PAV(IND,JND)*REAL(QQQ,4)*VOL(I)/VOLMER(JND) + 150 CONTINUE + 160 CONTINUE + DO 170 IND=1,NBNRS + MATRIX(INM(IND,IM,NBNRS),INM(JNBN,JM,NBNRS))=AWPHI(IND) + 170 CONTINUE + 180 CONTINUE + 190 CONTINUE + 200 CONTINUE +* + DO 210 I=1,NBNRS*MI + MATRIX(I,I)=MATRIX(I,I)+1.0D0 + 210 CONTINUE + CALL ALSBD(NBNRS*MI,1,MATRIX,IER,NBNRS*MI) + IF(IER.NE.0) CALL XABORT('USSEXD: SINGULAR MATRIX.') + XFLUX(:NBNRS,:MAXNOR,IG)=0.0 + DO 230 IND=1,NBNRS + DO 220 IM=1,MI + I1=INM(IND,IM,NBNRS) + XFLUX(IND,IM,IG)=REAL(MATRIX(I1,NBNRS*MI+1)) + 220 CONTINUE + 230 CONTINUE +* END OF RESPONSE MATRIX APPROACH. +*---- +* COMPUTE THE AVERAGED SOURCE. +*---- + ALLOCATE(FUN(NUN*MI),SUN(NUN*MI)) + SUN(:NUN*MI)=0.0 + ALLOCATE(NPSYS(MI)) + DO 250 IM=1,MI + NPSYS(IM)=IASM+IM + KPSYS=LCMGIL(IPSYS,IASM+IM) + CALL LCMLEN(KPSYS,'FUNKNO$USS',ILENG,ITYLCM) + IF(ILENG.EQ.NUN) THEN + CALL LCMGET(KPSYS,'FUNKNO$USS',FUN((IM-1)*NUN+1)) + ELSE + FUN((IM-1)*NUN+1:IM*NUN)=0.0 + ENDIF + SIGG(0:NBMIX)=0.0 + DO 240 IBM=1,NBMIX + IND=IREX(IBM) + QQQ=SIGGAR(IBM,0,IG,3)*GAMMA_V%VECTOR(IM) + IF(IND.GT.0) THEN + DO JG=1,IG + DO JRES=1,NIRES + IF(GOLD(IRES,JG).NE.-1001.) CYCLE + IF(JG.LT.IG-NJJ(IG,JRES)+1) CYCLE + DENSIT=CONR(IND,JRES) + DO JM=1,MRANK(JG) + IF((JG.EQ.IG).AND.(JM.EQ.IM)) CYCLE + QQQ=QQQ+DENSIT*SCAT_M(JG,JRES)%MATRIX(IM,JM)* + 1 XFLUX(IND,JM,JG) + ENDDO + ENDDO + ENDDO + ENDIF + SIGG(IBM)=REAL(QQQ,4) + 240 CONTINUE + IOF=(IM-1)*NUN + CALL DOORS(CDOOR,IPTRK,NBMIX,0,NUN,SIGG,SUN(IOF+1)) + 250 CONTINUE +* + IF(IMPX.GT.0) THEN + WRITE(TEXT12,'(3A4)') (IPPT2(IRES,I),I=2,4) + WRITE(6,'(15H USSEXD: GROUP=,I5,24H. SUBGROUP CALCULATION B, + 1 37HASED ON RESPONSE MATRICES. ISOTOPE='',A12,2H''.)') IG, + 2 TEXT12 + ENDIF +*---- +* SOLVE FOR THE MULTIBAND FLUX (VECTOR OF LENGTH NREG). +*---- + IPMACR=C_NULL_PTR + IPSOU=C_NULL_PTR + REBFLG=.FALSE. + CALL DOORFV(CDOOR,IPSYS,NPSYS,IPTRK,IFTRAK,IMPX,MI,NBMIX,IDIR, + 1 NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUN,FUN,IPMACR,IPSOU, + 2 REBFLG) + DEALLOCATE(NPSYS) +*---- +* INTEGRATE THE REGION-ORDERED FLUX OVER SUBGROUPS AND COMPUTE UNGAR, +* THE REGION-ORDERED FLUX. +*---- + UNGAR(:NREG,IRES,IG)=0.0 + DO 270 IM=1,MI + KPSYS=LCMGIL(IPSYS,IASM+IM) + IOF=(IM-1)*NUN + CALL LCMPUT(KPSYS,'FUNKNO$USS',NUN,2,FUN(IOF+1)) +* + DO 260 I=1,NREG + IOF=(IM-1)*NUN+KEYFLX(I) + UNGAR(I,IRES,IG)=UNGAR(I,IRES,IG)+REAL(WEIGHT_V%VECTOR(IM)* + 1 FUN(IOF),4) + 260 CONTINUE + 270 CONTINUE + DEALLOCATE(SUN,FUN) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(SIGG,SCAT_M,JPLIB2) + DEALLOCATE(MATRIX) + DEALLOCATE(AWPHI,PAV) + DEALLOCATE(NJJ) + RETURN + END |