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/BREERA.f | 443 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 443 insertions(+) create mode 100644 Dragon/src/BREERA.f (limited to 'Dragon/src/BREERA.f') diff --git a/Dragon/src/BREERA.f b/Dragon/src/BREERA.f new file mode 100644 index 0000000..336032f --- /dev/null +++ b/Dragon/src/BREERA.f @@ -0,0 +1,443 @@ +*DECK BREERA + SUBROUTINE BREERA(IPMAC1,NC,NG,NL,LX1,NMIX1,IMIX,ICODE,ISPH,ZKEFF, + 1 B2,ENER,XXX1,VOL1,FLX1,DC1,TOT1,CHI1,SIGF1,SCAT1,JXM,JXP,FHETXM, + 2 FHETXP,ADF1,NGET,ADFREF,IPRINT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Implement the 1D ERM-ANM reflector model. +* +*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 +* IPMAC1 nodal macrolib. +* NC number of sn macrolibs. +* NG number of energy groups. +* NL Legendre order of TOT1 and SCAT1 arrays (=1 for isotropic +* scattering in LAB). +* LX1 number of nodes in the reflector model. +* NMIX1 number of mixtures in the nodal calculation. +* IMIX mix index of each node. +* ICODE physical albedo index on each side of the domain. +* ISPH SPH flag (=0: use discontinuity factors; =1: use SPH factors). +* ZKEFF effective multiplication factor. +* B2 buckling. +* ENER energy limits. +* XXX1 spatial mesh. +* VOL1 volumes. +* FLX1 averaged fluxes +* DC1 diffusion coefficients. +* TOT1 total cross sections. +* CHI1 fission spectra. +* SIGF1 nu*fission cross sections. +* SCAT1 scattering P0 cross sections. +* JXM left boundary currents. +* JXP right boundary currents. +* FHETXM left boundary fluxes. +* FHETXP right boundary fluxes. +* ADF1 assembly discontinuity factors from macrolib. +* NGET type of NGET normalization if discontinuity factors +* (=0: simple; =1: imposed ADF on fuel assembly; =2: recover +* fuel assembly ADF from input macrolib). +* ADFREF imposed ADF values on fuel assembly side. +* IPRINT edition flag. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPMAC1 + INTEGER NC,NG,NL,LX1,NMIX1,IMIX(LX1),ICODE(2),ISPH,NGET,IPRINT + REAL ZKEFF(NC),B2(NC),ENER(NG+1),XXX1(LX1+1),VOL1(NMIX1,NC), + 1 FLX1(NMIX1,NG,NC),DC1(NMIX1,NG,NC),TOT1(NMIX1,NG,NL,NC), + 2 CHI1(NMIX1,NG,NC),SIGF1(NMIX1,NG,NC),SCAT1(NMIX1,NG,NG,NL,NC), + 3 JXM(NMIX1,NG,NC),JXP(NMIX1,NG,NC),FHETXM(NMIX1,NG,NL,NC), + 4 FHETXP(NMIX1,NG,NL,NC),ADF1(NMIX1,NG,NC),ADFREF(NG) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (NSTATE=40) + INTEGER ISTATE(NSTATE) + CHARACTER(LEN=8) HADF(2) + TYPE(C_PTR) JPMAC1,KPMAC1 +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS + REAL, ALLOCATABLE, DIMENSION(:) :: WORK1D,WORK1,WORK2,WORK4,WORK5, + 1 VOLTOT + REAL, ALLOCATABLE, DIMENSION(:,:) :: FLX,DC,TOT,CHI,SIGF, + 1 ADF,AFACTOR,BETA,WORK3 + REAL, ALLOCATABLE, DIMENSION(:,:,:) ::FDXM,FDXP,SCAT + REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: TAU,B,X + REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: WORK2D + REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: FHOMM,FHOMP,L,R +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(WORK1(NG),WORK2(NG),WORK4(NG),WORK5(NG),VOLTOT(NMIX1), + 1 FLX(NMIX1,NG),DC(NMIX1,NG),TOT(NMIX1,NG),CHI(NMIX1,NG), + 1 SIGF(NMIX1,NG),ADF(NMIX1,NG),AFACTOR(NG,NG),BETA(NG,NG)) + ALLOCATE(FDXM(NMIX1,NG,NG),FDXP(NMIX1,NG,NG),SCAT(NMIX1,NG,NG), + 1 WORK3(NG,NG)) + ALLOCATE(FHOMM(NC,NG,NMIX1),FHOMP(NC,NG,NMIX1),L(NG,2*NG,LX1), + 1 R(NG,2*NG,LX1)) +*---- +* AVERAGE THE OUTPUT NODAL MACROLIB +*---- + VOLTOT(:)=0.0 + FLX(:,:)=0.0 + DC(:,:)=0.0 + TOT(:,:)=0.0 + CHI(:,:)=0.0 + SIGF(:,:)=0.0 + SCAT(:,:,:)=0.0 + ADF(:,:)=0.0 + FHOMM(:NC,:NG,:NMIX1)=0.0D0 + FHOMP(:NC,:NG,:NMIX1)=0.0D0 + DO IC=1,NC + DO IBM=1,NMIX1 + VOLTOT(IBM)=VOLTOT(IBM)+VOL1(IBM,IC) + DO IGR=1,NG + FLX(IBM,IGR)=FLX(IBM,IGR)+FLX1(IBM,IGR,IC) + DC(IBM,IGR)=DC(IBM,IGR)+DC1(IBM,IGR,IC) + TOT(IBM,IGR)=TOT(IBM,IGR)+TOT1(IBM,IGR,1,IC) + CHI(IBM,IGR)=CHI(IBM,IGR)+CHI1(IBM,IGR,IC) + SIGF(IBM,IGR)=SIGF(IBM,IGR)+SIGF1(IBM,IGR,IC) + DO JGR=1,NG + SCAT(IBM,IGR,JGR)=SCAT(IBM,IGR,JGR)+SCAT1(IBM,IGR,JGR,1,IC) + ENDDO + ADF(IBM,IGR)=ADF(IBM,IGR)+ADF1(IBM,IGR,IC) + ENDDO + ENDDO + ENDDO + VOLTOT(:)=VOLTOT(:)/REAL(NC) + FLX(:,:)=FLX(:,:)/REAL(NC) + DC(:,:)=DC(:,:)/REAL(NC) + TOT(:,:)=TOT(:,:)/REAL(NC) + CHI(:,:)=CHI(:,:)/REAL(NC) + SIGF(:,:)=SIGF(:,:)/REAL(NC) + SCAT(:,:,:)=SCAT(:,:,:)/REAL(NC) + ADF(:,:)=ADF(:,:)/REAL(NC) +*---- +* LOOP OVER CASES +*---- + IF(ISPH.EQ.1) CALL XABORT('BREERA: SPH OPTION NOT IMPLEMENTED.') + J_FUEL=0 + DO IC=1,NC +*---- +* SET AND SOLVE ANALYTIC NODAL SYSTEM +*---- + DO I=1,LX1 + IBM=IMIX(I) + IF(IBM.EQ.0) CYCLE + WORK1(:NG)=DC1(IBM,:NG,IC) + WORK3(:NG,:NG)=SCAT1(IBM,:NG,:NG,1,IC) + WORK4(:NG)=CHI1(IBM,:NG,IC) + WORK5(:NG)=SIGF1(IBM,:NG,IC) + DO IGR=1,NG + IF(SIGF1(IBM,IGR,IC).GT.0.0) J_FUEL=I + WORK2(IGR)=TOT1(IBM,IGR,1,IC)+B2(IC)*DC1(IBM,IGR,IC)- + 1 SCAT1(IBM,IGR,IGR,1,IC) + ENDDO + VOL=XXX1(I+1)-XXX1(I) + CALL NSSLR1(ZKEFF(IC),NG,VOL,WORK1,WORK2,WORK3,WORK4,WORK5, + 1 L(1,1,I),R(1,1,I)) + ! + FHOMM(IC,:NG,IBM)=FHOMM(IC,:NG,IBM)+REAL(MATMUL(L(:NG,:NG,I), + 1 FLX1(IBM,:NG,IC))+ + 2 MATMUL(L(:NG,NG+1:2*NG,I),JXM(IBM,:NG,IC)),4)*VOL + FHOMP(IC,:NG,IBM)=FHOMP(IC,:NG,IBM)+REAL(MATMUL(R(:NG,:NG,I), + 1 FLX1(IBM,:NG,IC))+ + 2 MATMUL(R(:NG,NG+1:2*NG,I),JXP(IBM,:NG,IC)),4)*VOL + ENDDO + DO IBM=1,NMIX1 + FHOMM(IC,:NG,IBM)=FHOMM(IC,:NG,IBM)/VOLTOT(IBM) + FHOMP(IC,:NG,IBM)=FHOMP(IC,:NG,IBM)/VOLTOT(IBM) + ENDDO + IF(IPRINT.GT.0) THEN + WRITE(6,'(/39H BREERA: NODAL SURFACE FLUXES FOR CASE=,I5)') IC + DO IBM=1,NMIX1 + WRITE(6,'(/9H MIXTURE=,I5)') IBM + WRITE(6,20) 'FHOMM',FHOMM(IC,:NG,IBM) + WRITE(6,20) 'FHOMP',FHOMP(IC,:NG,IBM) + ENDDO + ENDIF +*---- +* END OF LOOP OVER CASES +*---- + ENDDO +*---- +* COMPUTE DISCONTINUITY AND ALBEDO FACTORS +*---- + AFACTOR(:,:)=0.0 + DO IBM=1,NMIX1 + IF(NC.EQ.1) THEN + ! DF-NEM approach + FDXM(IBM,:,:)=0.0 + FDXP(IBM,:,:)=0.0 + DO IGR=1,NG + FDXM(IBM,IGR,IGR)=FHETXM(IBM,IGR,1,1)/REAL(FHOMM(1,IGR,IBM)) + FDXP(IBM,IGR,IGR)=FHETXP(IBM,IGR,1,1)/REAL(FHOMP(1,IGR,IBM)) + ENDDO + IF(IBM.EQ.NMIX1) THEN + DO IGR=1,NG + AFACTOR(IGR,IGR)=JXP(IBM,IGR,1)/REAL(FHOMP(1,IGR,IBM)) + ENDDO + ENDIF + ELSE IF(NC.LT.NG) THEN + CALL XABORT('BREERA: DEGENERATE SYSTEM') + ELSE IF(NC.EQ.NG) THEN + ! ERM-ANM approach: linear system resolution + ALLOCATE(WORK2D(NC,2*NG)) + DO IGR=1,NG + DO IC=1,NC + WORK2D(IC,IGR)=FHOMM(IC,IGR,IBM) + WORK2D(IC,NG+IGR)=FHETXM(IBM,IGR,1,IC) + ENDDO + ENDDO + CALL ALSBD(NC,NG,WORK2D,IER,NC) + IF(IER.NE.0) CALL XABORT('BREERA: SINGULAR MATRIX(1).') + DO IGR=1,NG + DO IC=1,NC + FDXM(IBM,IGR,IC)=REAL(WORK2D(IC,NG+IGR)) + ENDDO + ENDDO + DO IGR=1,NG + DO IC=1,NC + WORK2D(IC,IGR)=FHOMP(IC,IGR,IBM) + WORK2D(IC,NG+IGR)=FHETXP(IBM,IGR,1,IC) + ENDDO + ENDDO + CALL ALSBD(NC,NG,WORK2D,IER,NC) + IF(IER.NE.0) CALL XABORT('BREERA: SINGULAR MATRIX(2).') + DO IGR=1,NG + DO IC=1,NC + FDXP(IBM,IGR,IC)=REAL(WORK2D(IC,NG+IGR)) + ENDDO + ENDDO + IF(IBM.EQ.NMIX1) THEN + DO IGR=1,NG + DO IC=1,NC + WORK2D(IC,IGR)=FHOMP(IC,IGR,IBM) + WORK2D(IC,NG+IGR)=JXP(IBM,IGR,IC) + ENDDO + ENDDO + CALL ALSBD(NC,NG,WORK2D,IER,NC) + IF(IER.NE.0) CALL XABORT('BREERA: SINGULAR MATRIX(3).') + DO IGR=1,NG + DO JGR=1,NG + AFACTOR(IGR,JGR)=REAL(WORK2D(JGR,NG+IGR)) + ENDDO + ENDDO + ENDIF + DEALLOCATE(WORK2D) + ELSE IF(NC.GE.NG) THEN + ! ERM-ANM approach: pseudo inversion + ALLOCATE(TAU(NG),B(NC),X(NG)) + CALL ALST2F(NC,NC,NG,FHOMM(1,1,IBM),TAU) + DO IGR=1,NG + B(:)=FHETXM(IBM,IGR,1,:) + CALL ALST2S(NC,NC,NG,FHOMM(1,1,IBM),TAU,B,X) + FDXM(IBM,IGR,:)=REAL(X(:)) + ENDDO + CALL ALST2F(NC,NC,NG,FHOMP(1,1,IBM),TAU) + DO IGR=1,NG + B(:)=FHETXP(IBM,IGR,1,:) + CALL ALST2S(NC,NC,NG,FHOMP(1,1,IBM),TAU,B,X) + FDXP(IBM,IGR,:)=REAL(X(:)) + ENDDO + IF(IBM.EQ.NMIX1) THEN + DO IGR=1,NG + B(:)=JXP(IBM,IGR,:) + CALL ALST2S(NC,NC,NG,FHOMP(1,1,IBM),TAU,B,X) + AFACTOR(IGR,:)=REAL(X(:)) + ENDDO + ENDIF + DEALLOCATE(X,B,TAU) + ENDIF + ENDDO + IF(IPRINT.GT.0) THEN + WRITE(6,'(/48H BREERA: DISCONTINUITY FACTORS BEFORE NORMALIZAT, + 1 3HION)') + DO IBM=1,NMIX1 + WRITE(6,'(/9H MIXTURE=,I5)') IBM + WRITE(6,20) 'FDXM',FDXM(IBM,:NG,:NG) + WRITE(6,20) 'FDXP',FDXP(IBM,:NG,:NG) + ENDDO + ENDIF +*---- +* COMPUTE ALBEDOS +*---- + IF(ICODE(2).NE.0) THEN + BETA(:,:)=0.0 + DO IGR=1,NG + DO JGR=1,NG + BETA(IGR,JGR)=(1.0-2.0*AFACTOR(IGR,JGR))/(1.0+2.0* + 1 AFACTOR(IGR,JGR)) + ENDDO + ENDDO + IF(IPRINT.GT.0) THEN + WRITE(6,'(/16H BREERA: ALBEDOS)') + WRITE(6,20) 'BETA',BETA(:NG,:NG) + ENDIF + ENDIF +*---- +* NGET NORMALIZATION OF THE DISCONTINUITY FACTORS +*---- + ALLOCATE(WORK2D(NG,2*NG)) + DO J=1,LX1-1 + IBM=IMIX(J) + IBMP=IMIX(J+1) + IF((IBM.EQ.0).OR.(IBMP.EQ.0)) CYCLE + DO IGR=1,NG + DO JGR=1,NG + WORK2D(IGR,JGR)=FDXP(IBM,IGR,JGR) + WORK2D(IGR,NG+JGR)=FDXM(IBMP,IGR,JGR) + ENDDO + ENDDO + CALL ALSBD(NG,NG,WORK2D,IER,NG) + IF(IER.NE.0) CALL XABORT('BREERA: SINGULAR MATRIX(3).') + DO IGR=1,NG + ! impose the adf on the fuel assembly side + IF((J.EQ.J_FUEL).AND.(NGET.EQ.1)) THEN + FNORM=ADFREF(IGR) + ELSE IF((J.EQ.J_FUEL).AND.(NGET.EQ.2)) THEN + FNORM=ADF(IBM,IGR) + ELSE + FNORM=FDXP(IBM,IGR,IGR) + ENDIF + FDXP(IBM,IGR,:)=0.0 + FDXP(IBM,IGR,IGR)=FNORM + DO JGR=1,NG + FDXM(IBMP,IGR,JGR)=REAL(WORK2D(IGR,NG+JGR))*FNORM + ENDDO + ENDDO + ENDDO + DEALLOCATE(WORK2D) + IF(J_FUEL.GT.0) THEN + DO J=J_FUEL,1,-1 + IBM=IMIX(J) + IF(IBM.EQ.0) CYCLE + DO IGR=1,NG + FNORM=FDXP(IBM,IGR,IGR)/FDXM(IBM,IGR,IGR) + DO JGR=1,NG + IF(J>1) THEN + IBMM=IMIX(J-1) + IF(IBMM.GT.0) FDXP(IBMM,IGR,JGR)=FDXP(IBMM,IGR,JGR)*FNORM + ENDIF + FDXM(IBM,IGR,JGR)=FDXM(IBM,IGR,JGR)*FNORM + ENDDO + ENDDO + ENDDO + ENDIF + DO J=J_FUEL+1,LX1 + IBM=IMIX(J) + IF(IBM.EQ.0) CYCLE + DO IGR=1,NG + FNORM=FDXM(IBM,IGR,IGR)/FDXP(IBM,IGR,IGR) + DO JGR=1,NG + IF(J