summaryrefslogtreecommitdiff
path: root/Dragon/src/BRERT.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/BRERT.f')
-rw-r--r--Dragon/src/BRERT.f414
1 files changed, 414 insertions, 0 deletions
diff --git a/Dragon/src/BRERT.f b/Dragon/src/BRERT.f
new file mode 100644
index 0000000..3dd87fb
--- /dev/null
+++ b/Dragon/src/BRERT.f
@@ -0,0 +1,414 @@
+*DECK BRERT
+ SUBROUTINE BRERT(IPMAC1,IELEM,ICOL,NG,NL,LX1,NMIX1,IMIX,ICODE,
+ 1 ISPH,IDIFF,ZKEFF,B2,ENER,XXX1,VOL1,FLX1,DC1,TOT1,CHI1,SIGF1,
+ 2 SCAT1,JXM,JXP,FHETXM,FHETXP,ADF1,NGET,ADFREF,IMPX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Implement the 1D DF-RT (Raviart-Thomas) reflector model.
+*
+*Copyright:
+* Copyright (C) 2025 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.
+* IELEM Raviart-Thomas polynomial order.
+* ICOL Raviart-Thomas polynomial integration type.
+* NG number of energy groups.
+* NL Legendre order of TOT1 and SCAT1 arrays (=1 for isotropic
+* scattering in LAB). (NL-1) is the SPN order (if IDIFF>1,
+* NL is an even integer).
+* 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).
+* IDIFF PN calculation option (=0: diffusion theory; =1: SPN theory
+* with 'NTOT1'; =2: SPN theory with 1/(3*D)).
+* 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.
+* IMPX edition flag.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPMAC1
+ INTEGER IELEM,ICOL,NG,NL,LX1,NMIX1,IMIX(LX1),ICODE(2),ISPH,IDIFF,
+ 1 NGET,IMPX
+ REAL ZKEFF,B2,ENER(NG+1),XXX1(LX1+1),VOL1(NMIX1),FLX1(NMIX1,NG),
+ 1 DC1(NMIX1,NG),TOT1(NMIX1,NG,NL),CHI1(NMIX1,NG),SIGF1(NMIX1,NG),
+ 2 SCAT1(NMIX1,NG,NG,NL),JXM(NMIX1,NG),JXP(NMIX1,NG),
+ 3 FHETXM(NMIX1,NG,NL),FHETXP(NMIX1,NG,NL),ADF1(NMIX1,NG),ADFREF(NG)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (NSTATE=40)
+ INTEGER ISTATE(NSTATE)
+ CHARACTER CM*2,HADF*8,TEXT12*12
+ TYPE(C_PTR) JPMAC1,KPMAC1
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS
+ REAL, ALLOCATABLE, DIMENSION(:) :: WORK,AFACTOR,BETA,WORK1,WORK2,
+ 1 VOLTOT
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: FDXM,FDXP,WORK3,WORK4,WORK5,
+ 1 WORK6,WORK7
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: FHOMM,FHOMP
+ REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:) :: RCAT
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ J_FUEL=0
+ NLMAX=1
+ IF(IDIFF.GT.0) THEN
+ IF(NL.LT.2) CALL XABORT('BRERT: EVEN NL>=2 EXPECTED WITH SPN.')
+ NLMAX=NL/2
+ ENDIF
+ ALLOCATE(FHOMM(NMIX1,NG,NLMAX),FHOMP(NMIX1,NG,NLMAX),
+ 1 FDXM(NMIX1,NG),FDXP(NMIX1,NG),AFACTOR(NG),BETA(NG),VOLTOT(NMIX1),
+ 2 WORK1(NG),WORK2(NG),WORK3(NG,NG),WORK6(NG,NLMAX),WORK7(NG,NLMAX))
+*----
+* COMPUTE BOUNDARY FLUXES
+*----
+ FDXM(:NMIX1,:NG)=0.0
+ FDXP(:NMIX1,:NG)=0.0
+ FHOMM(:NMIX1,:NG,:NLMAX)=0.0
+ FHOMP(:NMIX1,:NG,:NLMAX)=0.0
+ VOLTOT(:NMIX1)=0.0
+ DO I=1,LX1
+ IBM=IMIX(I)
+ IF(IBM.EQ.0) CYCLE
+ WORK6(:NG,:NLMAX)=0.0
+ WORK7(:NG,:NLMAX)=0.0
+ DELX=XXX1(I+1)-XXX1(I)
+ IF(IMPX.GT.0) WRITE(6,'(/15H BRERT: REGION=,I5)') I
+ IF(IDIFF.EQ.0) THEN
+* diffusion theory
+ ALLOCATE(WORK4(NG,1),WORK5(NG,1))
+ WORK1(:NG)=DC1(IBM,:NG)
+ WORK4(:NG,1)=JXM(IBM,:NG)
+ WORK5(:NG,1)=JXP(IBM,:NG)
+ DO IG=1,NG
+ IF(SIGF1(IBM,IG).GT.0.0) J_FUEL=I
+ WORK2(IG)=TOT1(IBM,IG,1)+B2*DC1(IBM,IG)-SCAT1(IBM,IG,IG,1)
+ DO JG=1,NG
+ WORK3(IG,JG)=CHI1(IBM,IG)*SIGF1(IBM,JG)/ZKEFF
+ IF(JG.NE.IG) WORK3(IG,JG)=WORK3(IG,JG)+SCAT1(IBM,IG,JG,1)
+ ENDDO
+ ENDDO
+ CALL BRERTD(IELEM,ICOL,NG,DELX,WORK1,WORK2,WORK3,WORK4,WORK5,
+ 1 IMPX,WORK6,WORK7)
+ DEALLOCATE(WORK5,WORK4)
+ ELSE
+* SPN theory
+ ALLOCATE(WORK4(NG,NL/2),WORK5(NG,NL/2),RCAT(NG,NG,NL))
+ DO IL=1,NL/2
+ WORK4(:NG,IL)=FHETXM(IBM,:NG,2*IL)
+ WORK5(:NG,IL)=FHETXP(IBM,:NG,2*IL)
+ ENDDO
+ RCAT(:NG,:NG,:NL)=0.0
+ DO IG=1,NG
+ IF(SIGF1(IBM,IG).GT.0.0) J_FUEL=I
+ DO JG=1,NG
+ RCAT(IG,JG,1)=-CHI1(IBM,IG)*SIGF1(IBM,JG)/ZKEFF
+ ENDDO
+ RCAT(IG,IG,1)=RCAT(IG,IG,1)+B2*DC1(IBM,IG)
+ DO IL=1,NL,2
+ RCAT(IG,IG,IL)=RCAT(IG,IG,IL)+TOT1(IBM,IG,IL)
+ DO JG=1,NG
+ RCAT(IG,JG,IL)=RCAT(IG,JG,IL)-SCAT1(IBM,IG,JG,IL)
+ ENDDO
+ ENDDO
+ DO IL=2,NL,2
+ IF(IDIFF.EQ.1) THEN
+ DO JG=1,NG
+ RCAT(IG,JG,IL)=RCAT(IG,JG,IL)-SCAT1(IBM,IG,JG,IL)
+ ENDDO
+ ELSE
+ TOT1(IBM,IG,IL)=1.0/(3.0*DC1(IBM,IG))
+ SCAT1(IBM,IG,:NG,IL)=0.0
+ ENDIF
+ RCAT(IG,IG,IL)=RCAT(IG,IG,IL)+TOT1(IBM,IG,IL)
+ ENDDO
+ ENDDO
+ DO IL=1,NL
+ RCAT(:NG,:NG,IL)=RCAT(:NG,:NG,IL)*REAL(2*IL-1)
+ ENDDO
+ CALL BRERTS(IELEM,ICOL,NG,NL,DELX,RCAT,WORK4,WORK5,IMPX,
+ 1 WORK6,WORK7)
+ DEALLOCATE(RCAT,WORK5,WORK4)
+ ENDIF
+ FHOMM(IBM,:NG,:NLMAX)=FHOMM(IBM,:NG,:NLMAX)+WORK6(:NG,:NLMAX)*
+ 1 DELX
+ FHOMP(IBM,:NG,:NLMAX)=FHOMP(IBM,:NG,:NLMAX)+WORK7(:NG,:NLMAX)*
+ 1 DELX
+ VOLTOT(IBM)=VOLTOT(IBM)+DELX
+ ENDDO
+ DEALLOCATE(WORK7,WORK6,WORK3,WORK2,WORK1)
+ DO IBM=1,NMIX1
+ DO IGR=1,NG
+ IF(NL.LE.2) THEN
+ FDXM(IBM,IGR)=VOLTOT(IBM)*FHETXM(IBM,IGR,1)/FHOMM(IBM,IGR,1)
+ FDXP(IBM,IGR)=VOLTOT(IBM)*FHETXP(IBM,IGR,1)/FHOMP(IBM,IGR,1)
+ ELSE
+ ! Yamamoto formula
+ FDXM(IBM,IGR)=VOLTOT(IBM)*(FHETXM(IBM,IGR,1)+2.0*
+ 1 FHETXM(IBM,IGR,2))/(FHOMM(IBM,IGR,1)+2.0*FHOMM(IBM,IGR,2))
+ FDXP(IBM,IGR)=VOLTOT(IBM)*(FHETXP(IBM,IGR,1)+2.0*
+ 1 FHETXP(IBM,IGR,2))/(FHOMP(IBM,IGR,1)+2.0*FHOMP(IBM,IGR,2))
+ ENDIF
+ ENDDO
+ ENDDO
+ IF(IMPX.GT.0) THEN
+ WRITE(6,'(/48H BRERT: DISCONTINUITY FACTORS BEFORE NORMALIZATI,
+ 1 2HON)')
+ DO IBM=1,NMIX1
+ WRITE(6,'(/9H MIXTURE=,I5)') IBM
+ WRITE(6,20) 'FDXM',FDXM(IBM,:NG)
+ WRITE(6,20) 'FDXP',FDXP(IBM,:NG)
+ ENDDO
+ ENDIF
+*----
+* COMPUTE ALBEDOS
+*----
+ IF(ICODE(2).NE.0) THEN
+ BETA(:)=0.0
+ IBM=IMIX(LX1)
+ DO IGR=1,NG
+ IF(IBM.EQ.0) CYCLE
+ AFACTOR(IGR)=FDXP(IBM,IGR)*JXP(IBM,IGR)/FHETXP(IBM,IGR,1)
+ BETA(IGR)=(1.0-2.0*AFACTOR(IGR))/(1.0+2.0*AFACTOR(IGR))
+ ENDDO
+ IF(IMPX.GT.0) THEN
+ WRITE(6,'(/15H BRERT: ALBEDOS)')
+ WRITE(6,20) 'BETA',BETA(:NG)
+ ENDIF
+ ENDIF
+*----
+* NGET NORMALIZATION OF THE DISCONTINUITY FACTORS
+*----
+ IF(J_FUEL.GT.0) THEN
+ IF(NGET.GT.0) THEN
+ IBM=IMIX(J_FUEL)
+ DO IGR=1,NG
+ ! impose the adf on the fuel assembly side
+ IF(IBM.EQ.0) CYCLE
+ IF(NGET.EQ.1) THEN
+ FNORM=ADFREF(IGR)/FDXP(IBM,IGR)
+ ELSE
+ FNORM=ADF1(IBM,IGR)/FDXP(IBM,IGR)
+ ENDIF
+ FDXP(IBM,IGR)=FDXP(IBM,IGR)*FNORM
+ IF(J_FUEL<LX1) THEN
+ IBMP=IMIX(J_FUEL+1)
+ IF(IBMP.GT.0) FDXM(IBMP,IGR)=FDXM(IBMP,IGR)*FNORM
+ ENDIF
+ ENDDO
+ ENDIF
+ DO J=J_FUEL,1,-1
+ IBM=IMIX(J)
+ IF(IBM.EQ.0) CYCLE
+ DO IGR=1,NG
+ IF(J>1) THEN
+ IBMM=IMIX(J-1)
+ IF(IBMM.GT.0) FDXP(IBMM,IGR)=FDXP(IBMM,IGR)*FDXP(IBM,IGR)/
+ 1 FDXM(IBM,IGR)
+ ENDIF
+ FDXM(IBM,IGR)=FDXP(IBM,IGR)
+ ENDDO
+ ENDDO
+ ENDIF
+ DO J=J_FUEL+1,LX1
+ IBM=IMIX(J)
+ IF(IBM.EQ.0) CYCLE
+ DO IGR=1,NG
+ IF(J<LX1) THEN
+ IBMP=IMIX(J+1)
+ IF(IBMP.GT.0) FDXM(IBMP,IGR)=FDXM(IBMP,IGR)*FDXM(IBM,IGR)/
+ 1 FDXP(IBM,IGR)
+ ENDIF
+ FDXP(IBM,IGR)=FDXM(IBM,IGR)
+ ENDDO
+ ENDDO
+ IF(IMPX.GT.0) THEN
+ WRITE(6,'(/48H BRERT: DISCONTINUITY FACTORS AFTER NGET NORMALI,
+ 1 6HZATION)')
+ DO IBM=1,NMIX1
+ WRITE(6,'(/9H MIXTURE=,I5)') IBM
+ WRITE(6,20) 'FDX',FDXM(IBM,:NG)
+ ENDDO
+ ENDIF
+*----
+* APPLY SPH FACTORS
+*----
+ IF(ISPH.EQ.1) THEN
+ DO IGR=1,NG
+ DO IBM=1,NMIX1
+ IF(FDXM(IBM,IGR).LE.0) CALL XABORT('BRERT: NEGATIVE SPH F'
+ 1 //'ACTOR.')
+ ENDDO
+ DO J=1,LX1
+ IBM=IMIX(J)
+ IF(IBM.EQ.0) CYCLE
+ DC1(IBM,IGR)=DC1(IBM,IGR)/FDXM(IBM,IGR)
+ SIGF1(IBM,IGR)=SIGF1(IBM,IGR)/FDXM(IBM,IGR)
+ DO IL=1,NL,2
+ TOT1(IBM,IGR,IL)=TOT1(IBM,IGR,IL)/FDXM(IBM,IGR)
+ DO JGR=1,NG
+ SCAT1(IBM,IGR,JGR,IL)=SCAT1(IBM,IGR,JGR,IL)/FDXM(IBM,JGR)
+ ENDDO
+ ENDDO
+ DO IL=2,NL,2
+ TOT1(IBM,IGR,IL)=TOT1(IBM,IGR,IL)*FDXM(IBM,IGR)
+ DO JGR=1,NG
+ SCAT1(IBM,IGR,JGR,IL)=SCAT1(IBM,IGR,JGR,IL)*FDXM(IBM,IGR)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ IF(ICODE(2).NE.0) THEN
+ BETA(:)=0.0
+ IF(ICODE(2).NE.0) THEN
+ IBM=IMIX(LX1)
+ DO IGR=1,NG
+ IF(IBM.EQ.0) CYCLE
+ AFACTOR(IGR)=AFACTOR(IGR)/FDXM(IBM,IGR)
+ BETA(IGR)=(1.0-2.0*AFACTOR(IGR))/(1.0+2.0*AFACTOR(IGR))
+ ENDDO
+ ENDIF
+ IF(IMPX.GT.0) THEN
+ WRITE(6,'(/29H BRERT: SPH CORRECTED ALBEDOS)')
+ WRITE(6,20) 'BETA',BETA(:NG)
+ ENDIF
+ ENDIF
+ ENDIF
+ IF(IMPX.GT.0) THEN
+ WRITE(6,'(/30H BRERT: DIFFUSION COEFFICIENTS)')
+ DO IBM=1,NMIX1
+ WRITE(6,'(/9H MIXTURE=,I5)') IBM
+ WRITE(6,20) 'DIFF',DC1(IBM,:NG)
+ ENDDO
+ ENDIF
+*----
+* SAVE THE OUTPUT NODAL MACROLIB
+*----
+ ALLOCATE(IJJ(NMIX1),NJJ(NMIX1),IPOS(NMIX1),WORK(NMIX1*NG))
+ ISTATE(:)=0
+ ISTATE(1)=NG
+ ISTATE(2)=NMIX1
+ ISTATE(3)=NL
+ IF(J_FUEL.GT.0) ISTATE(4)=1
+ IF(ICODE(2).NE.0) ISTATE(8)=1 ! physical albedo information
+ ISTATE(9)=1 ! diffusion coefficient information
+ IF(ISPH.EQ.0) ISTATE(12)=3 ! discontinuity factor information
+ IF(ISPH.EQ.1) ISTATE(14)=1 ! SPH factor information
+ CALL LCMPUT(IPMAC1,'STATE-VECTOR',NSTATE,1,ISTATE)
+ CALL LCMPUT(IPMAC1,'ENERGY',NG+1,2,ENER)
+ CALL LCMPUT(IPMAC1,'VOLUME',NMIX1,2,VOL1)
+ CALL LCMPUT(IPMAC1,'B2 B1HOM',1,2,B2)
+ IF(ICODE(2).NE.0) CALL LCMPUT(IPMAC1,'ALBEDO',NG,2,BETA)
+ IF(ISPH.EQ.0) THEN
+ CALL LCMSIX(IPMAC1,'ADF',1)
+ NTYPE=1
+ HADF='FD_B'
+ CALL LCMPUT(IPMAC1,'NTYPE',1,1,NTYPE)
+ CALL LCMPTC(IPMAC1,'HADF',8,HADF)
+ CALL LCMPUT(IPMAC1,HADF,NMIX1*NG,2,FDXM)
+ CALL LCMSIX(IPMAC1,' ',2)
+ ELSE IF(ISPH.EQ.1) THEN
+ CALL LCMSIX(IPMAC1,'SPH',1)
+ ISTATE(:)=0
+ ISTATE(1)=4
+ ISTATE(2)=1
+ ISTATE(6)=1
+ ISTATE(7)=1
+ ISTATE(8)=NG
+ CALL LCMPUT(IPMAC1,'STATE-VECTOR',NSTATE,1,ISTATE)
+ CALL LCMSIX(IPMAC1,' ',2)
+ ENDIF
+ JPMAC1=LCMLID(IPMAC1,'GROUP',NG)
+ DO IGR=1,NG
+ KPMAC1=LCMDIL(JPMAC1,IGR)
+ DO IBM=1,NMIX1
+ WORK(IBM)=VOL1(IBM)*FLX1(IBM,IGR)
+ ENDDO
+ CALL LCMPUT(KPMAC1,'FLUX-INTG',NMIX1,2,WORK)
+ DO IL=1,NL
+ WRITE(TEXT12,'(4HNTOT,I1)') IL-1
+ CALL LCMPUT(KPMAC1,TEXT12,NMIX1,2,TOT1(:,IGR,IL))
+ ENDDO
+ CALL LCMPUT(KPMAC1,'DIFF',NMIX1,2,DC1(:,IGR))
+ CALL LCMPUT(KPMAC1,'CHI',NMIX1,2,CHI1(:,IGR))
+ CALL LCMPUT(KPMAC1,'NUSIGF',NMIX1,2,SIGF1(:,IGR))
+ IF(ISPH.EQ.1) THEN
+ DO IBM=1,NMIX1
+ WORK(IBM)=1.0/FDXM(IBM,IGR)
+ ENDDO
+ CALL LCMPUT(KPMAC1,'NSPH',NMIX1,2,WORK)
+ ENDIF
+ DO IL=1,NL
+ WRITE(CM,'(I2.2)') IL-1
+ WORK(:NMIX1)=SCAT1(:NMIX1,IGR,IGR,IL)
+ CALL LCMPUT(KPMAC1,'SIGW'//CM,NMIX1,2,WORK)
+ IPOSDE=0
+ DO IBM=1,NMIX1
+ J2=IGR
+ J1=IGR
+ DO JGR=1,NG
+ IF(SCAT1(IBM,IGR,JGR,IL).NE.0.0) THEN
+ J2=MAX(J2,JGR)
+ J1=MIN(J1,JGR)
+ ENDIF
+ ENDDO
+ NJJ(IBM)=J2-J1+1
+ IJJ(IBM)=J2
+ IPOS(IBM)=IPOSDE+1
+ DO JGR=J2,J1,-1
+ IPOSDE=IPOSDE+1
+ IF(IPOSDE.GT.NG*NMIX1) CALL XABORT('BRERT: SCAT OVERFLOW.')
+ WORK(IPOSDE)=SCAT1(IBM,IGR,JGR,IL)
+ ENDDO
+ ENDDO
+ CALL LCMPUT(KPMAC1,'SCAT'//CM,IPOSDE,2,WORK)
+ CALL LCMPUT(KPMAC1,'NJJS'//CM,NMIX1,1,NJJ)
+ CALL LCMPUT(KPMAC1,'IJJS'//CM,NMIX1,1,IJJ)
+ CALL LCMPUT(KPMAC1,'IPOS'//CM,NMIX1,1,IPOS)
+ ENDDO
+ ENDDO
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(WORK,IPOS,NJJ,IJJ)
+ DEALLOCATE(VOLTOT,BETA,AFACTOR,FDXP,FDXM,FHOMP,FHOMM)
+ RETURN
+ 20 FORMAT(1X,A9,1P,10E12.4,/(10X,10E12.4))
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:24:30 +0000