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/B1DIF.f | 381 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 381 insertions(+) create mode 100644 Dragon/src/B1DIF.f (limited to 'Dragon/src/B1DIF.f') diff --git a/Dragon/src/B1DIF.f b/Dragon/src/B1DIF.f new file mode 100644 index 0000000..8739ccc --- /dev/null +++ b/Dragon/src/B1DIF.f @@ -0,0 +1,381 @@ +*DECK B1DIF + SUBROUTINE B1DIF(OPTION,TYPE,NGRO,ST,SFNU,XHI,IJJ0,IJJ1,NJJ0,NJJ1, + 1 SCAT0,SCAT1,REFKEF,LFISSI,IMPX,DHOM,GAMMA,B2,ALAM1,CAET,A2,PHI) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solve the B-n equations and perform a buckling search if required. +* +*Copyright: +* Copyright (C) 2002 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 +* OPTION type of leakage coefficient. Can be 'LKRD', 'RHS', 'B0', 'P0', +* 'B1', 'P1', 'B0TR' or 'P0TR'. 'LKRD' and 'RHS' are used to +* impose a leakage coefficient. +* TYPE type of buckling search. Can be 'DIFF', 'K', 'B' or 'L'. +* NGRO number of energy groups. +* ST macroscopic total cross sections. +* SFNU nu*macroscopic fission cross sections. +* XHI fission spectrum normalized to one. +* IJJ0 most thermal group in band for P0 scattering. +* NJJ0 number of groups in band for P0 scattering. +* IJJ1 most thermal group in band for P1 scattering. +* NJJ1 number of groups in band for P1 scattering. +* SCAT0 packed diffusion P0 macroscopic cross sections. +* SCAT1 packed diffusion P1 macroscopic cross sections. +* REFKEF target K-effective for type B or type L calculations. +* LFISSI fissile isotope flag (=.TRUE. if present). +* IMPX print flag. +* +*Parameters: input/output +* PHI homogeneous flux from heterogeneous calculation on input and +* fundamental flux at output. +* +*Parameters: output +* DHOM homogeneous leakage coefficients. +* GAMMA gamma factors. +* B2 buckling. +* ALAM1 effective multiplication factor. +* CAET infinite multiplication factor. +* A2 migration area. +* +*----------------------------------------------------------------------- +* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + CHARACTER OPTION*4,TYPE*4 + INTEGER NGRO,IJJ0(NGRO),IJJ1(NGRO),NJJ0(NGRO),NJJ1(NGRO),IMPX + REAL ST(NGRO),SFNU(NGRO),XHI(NGRO),SCAT0(*),SCAT1(*),DHOM(NGRO), + > GAMMA(NGRO) + DOUBLE PRECISION B2,ALAM1,CAET,A2,PHI(NGRO),REFKEF + LOGICAL LFISSI +*---- +* LOCAL VARIABLES +*---- + PARAMETER (EPS=1.0D-6,MAXIT=50) + DOUBLE PRECISION FFITX(MAXIT),B2ITX(MAXIT) + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: CSTOC,SA + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: ASTOC,BSTOC +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(ASTOC(NGRO,NGRO+1),BSTOC(NGRO,NGRO),CSTOC(NGRO),SA(NGRO)) +* + IF((IMPX.GT.0).AND.(TYPE.EQ.'DIFF')) THEN + WRITE (6,400) + ELSE IF(IMPX.GT.0) THEN + WRITE (6,410) OPTION,TYPE + ENDIF + IAPROX=2 + IF((OPTION.EQ.'P0').OR.(OPTION.EQ.'P1').OR.(OPTION.EQ.'P0TR')) + > IAPROX=1 + IF((OPTION.EQ.'LKRD').OR.(OPTION.EQ.'RHS')) IAPROX=0 + BIL1=0.0D0 + DO 5 I=1,NGRO + BIL1=BIL1+XHI(I) + SA(I)=ST(I) + 5 CONTINUE + IF((BIL1.GT.0.9D0).AND.(ABS(BIL1-1.0D0).GT.EPS)) THEN + IF(IMPX.GT.0) + > WRITE(6,'(46H B1DIF: WARNING INCONSISTENT FISSION SPECTRUM.)') + ENDIF + IGAR=0 + DO 11 I=1,NGRO + DO 10 J=IJJ0(I),IJJ0(I)-NJJ0(I)+1,-1 + IGAR=IGAR+1 + SA(J)=SA(J)-SCAT0(IGAR) + 10 CONTINUE + 11 CONTINUE + IF(TYPE.EQ.'DIFF') THEN + DO 15 I=1,NGRO + ST2=DBLE(ST(I)) + GAMMA(I)=REAL(B1GAMA(IAPROX,B2,ST2)) + IF(IAPROX.NE.0) DHOM(I)=REAL(1.0D0/(3.0D0*GAMMA(I)*ST2)) + 15 CONTINUE + RETURN + ELSE IF((TYPE.EQ.'B').OR.(TYPE.EQ.'L')) THEN +* COMPUTE THE INITIAL BUCKLING. + BIL1=0.0D0 + BIL2=0.0D0 + IF(LFISSI) THEN + DO 20 I=1,NGRO + BIL1=BIL1+(SFNU(I)/REFKEF)*PHI(I) + 20 CONTINUE + ENDIF + DO 21 I=1,NGRO + BIL1=BIL1-SA(I)*PHI(I) + BIL2=BIL2+DHOM(I)*PHI(I) + 21 CONTINUE + B2=BIL1/BIL2 + DO 25 I=1,NGRO + ST2=-0.7D0*(ST(I)**2) + IF(B2.LT.ST2) THEN + IF(IMPX.GT.0) WRITE (6,415) B2,ST2 + B2=ST2 + ENDIF + PHI(I)=1.0D0 + 25 CONTINUE + ENDIF + IF(IMPX.GT.1) WRITE(6,420) B2 + IF(TYPE.EQ.'L') GO TO 160 +*---- +* COMPUTE THE FUNDAMENTAL FLUX WITH TYPE K OR TYPE B +*---- + ITEX=0 + 30 ITEX=ITEX+1 + IF(ITEX.GT.MAXIT) CALL XABORT('B1DIF: UNABLE TO CONVERGE(1).') + IGAR=0 + DO 55 I=1,NGRO + ST2=DBLE(ST(I)) + DD=DBLE(DHOM(I)) + BETA=B1BETA(IAPROX,B2,ST2,DD) + DO 40 J=1,NGRO + ASTOC(I,J)=0.0D0 + BSTOC(I,J)=0.0D0 + 40 CONTINUE + ASTOC(I,I)=ST2 + ASTOC(I,NGRO+1)=(1.0D0-B2*BETA)*XHI(I) + DO 50 J=IJJ0(I),IJJ0(I)-NJJ0(I)+1,-1 + IGAR=IGAR+1 + ASTOC(I,J)=ASTOC(I,J)-(1.0D0-B2*BETA)*SCAT0(IGAR) + BSTOC(I,J)=SCAT0(IGAR) + 50 CONTINUE + 55 CONTINUE + IF((OPTION.EQ.'P1').OR.(OPTION.EQ.'B1')) THEN + DO 72 J=1,NGRO + DO 60 K=1,NGRO + CSTOC(K)=BSTOC(K,J) + BSTOC(K,J)=0.0D0 + 60 CONTINUE + IGAR=0 + DO 71 I=1,NGRO + DO 70 K=IJJ1(I),IJJ1(I)-NJJ1(I)+1,-1 + IGAR=IGAR+1 + BSTOC(I,J)=BSTOC(I,J)+3.0D0*SCAT1(IGAR)*CSTOC(K) + 70 CONTINUE + 71 CONTINUE + 72 CONTINUE + IGAR=0 + DO 95 I=1,NGRO + ST2=DBLE(ST(I)) + DD=DBLE(DHOM(I)) + BETA=B1BETA(IAPROX,B2,ST2,DD)*ST2 + DO 80 J=1,NGRO + ASTOC(I,J)=ASTOC(I,J)+BETA*BSTOC(I,J) + 80 CONTINUE + DO 90 J=IJJ1(I),IJJ1(I)-NJJ1(I)+1,-1 + IGAR=IGAR+1 + ASTOC(I,NGRO+1)=ASTOC(I,NGRO+1)-3.0D0*BETA*SCAT1(IGAR)*XHI(J) + ASTOC(I,J)=ASTOC(I,J)-3.0D0*BETA*SCAT1(IGAR)*ST(J) + 90 CONTINUE + 95 CONTINUE + ENDIF + CALL B1SOL(NGRO,ASTOC,IER) + IF(IER.NE.0) CALL XABORT('B1DIF: SINGULAR MATRIX(1).') + ALAM1=0.0D0 + CAET=0.0D0 + DO 130 I=1,NGRO + ALAM1=ALAM1+SFNU(I)*ASTOC(I,NGRO+1) + CAET=CAET+SA(I)*ASTOC(I,NGRO+1) + 130 CONTINUE + IF(IMPX.GT.1) WRITE (6,430) ITEX,ALAM1,B2 + B2ITX(ITEX)=B2 + FFITX(ITEX)=REFKEF-ALAM1 + IF(TYPE.EQ.'K') THEN + DO 140 I=1,NGRO + PHI(I)=ASTOC(I,NGRO+1)/ALAM1 + 140 CONTINUE + ELSE IF(TYPE.EQ.'B') THEN +* COMPUTE THE EXTRAPOLATED BUCKLING. + IF(ITEX.LE.5) THEN +* USE A BALANCE RELATION. + B2=B2*(ALAM1/REFKEF-CAET)/(1.0D0-CAET) + ELSE + IF(ITEX.EQ.6) THEN +* SORT THE ROOT CONVERGENCE HISTORY. + DO I=2,ITEX-1 + WORKF=FFITX(ITEX-I) + WORKB=B2ITX(ITEX-I) + J=I + DO WHILE((J.GT.0).AND. + > (ABS(FFITX(ITEX-J+1)).GT.ABS(WORKF))) + FFITX(ITEX-J)=FFITX(ITEX-J+1) + B2ITX(ITEX-J)=B2ITX(ITEX-J+1) + J=J-1 + ENDDO + FFITX(ITEX-J)=WORKF + B2ITX(ITEX-J)=WORKB + ENDDO + ENDIF + J=0 + DO I=ITEX-1,2,-1 + IF(FFITX(I)*FFITX(ITEX).LT.0.0) THEN + J=I + EXIT + ENDIF + ENDDO + IF(J.NE.0) THEN +* USE A BISSECTION METHOD. + B2=0.5D0*(B2ITX(J)+B2ITX(ITEX)) + ELSE +* USE THE SECANT METHOD. + AA=FFITX(ITEX)-FFITX(ITEX-1) + B2=(B2ITX(ITEX-1)*FFITX(ITEX)-B2ITX(ITEX)*FFITX(ITEX-1))/AA + ENDIF + ENDIF +* CHECK THE CONVERGENCE. + BIL1=0.0D0 + BIL2=0.0D0 + DO 150 I=1,NGRO + ST2=ST(I)**2 + BIL1=MAX(BIL1,ABS(ASTOC(I,NGRO+1)/ALAM1)) + BIL2=MAX(BIL2,ABS(PHI(I)-ASTOC(I,NGRO+1)/ALAM1)) + PHI(I)=ASTOC(I,NGRO+1)/ALAM1 + 150 CONTINUE + ERR3=ABS(REFKEF-ALAM1) + IF((BIL2.GE.10*EPS*BIL1).OR.(ERR3.GE.EPS)) GO TO 30 + ENDIF + GO TO 300 +*---- +* COMPUTE THE FUNDAMENTAL FLUX WITH TYPE L +*---- + 160 ITEX=0 + 170 ITEX=ITEX+1 + IF(ITEX.GT.MAXIT) CALL XABORT('B1DIF: UNABLE TO CONVERGE(2).') + IF((OPTION.EQ.'P1').OR.(OPTION.EQ.'B1')) THEN + IGAR=0 + DO 200 I=1,NGRO + ST2=DBLE(ST(I)) + DD=DBLE(DHOM(I)) + BETA=B1BETA(IAPROX,B2,ST2,DD) + BETA=BETA*ST2/(1.0D0-B2*BETA) + DO 180 J=1,NGRO + ASTOC(I,J)=0.0D0 + 180 CONTINUE + ASTOC(I,I)=1.0D0 + ASTOC(I,NGRO+1)=BETA + DO 190 J=IJJ1(I),IJJ1(I)-NJJ1(I)+1,-1 + IGAR=IGAR+1 + ASTOC(I,J)=ASTOC(I,J)-3.0D0*BETA*SCAT1(IGAR)*PHI(J)/PHI(I) + 190 CONTINUE + 200 CONTINUE + CALL B1SOL(NGRO,ASTOC,IER) + IF(IER.NE.0) CALL XABORT('B1DIF: SINGULAR MATRIX(2).') + DO 210 I=1,NGRO + DHOM(I)=REAL(ASTOC(I,NGRO+1)) + 210 CONTINUE + ELSE IF((OPTION.NE.'LKRD').AND.(OPTION.NE.'RHS')) THEN + DO 220 I=1,NGRO + ST2=ST(I) + GAMMA(I)=REAL(B1GAMA(IAPROX,B2,ST2)) + IF(IAPROX.NE.0) DHOM(I)=REAL(1.0D0/(3.0D0*GAMMA(I)*ST2)) + 220 CONTINUE + ENDIF + ASTOC(:NGRO,:NGRO)=0.0D0 + BSTOC(:NGRO,:NGRO)=0.0D0 + ASTOC(:NGRO,NGRO+1)=PHI(:NGRO) + IGAR=0 + DO 250 I=1,NGRO + ASTOC(I,I)=ASTOC(I,I)+ST(I) + BSTOC(I,I)=BSTOC(I,I)-DHOM(I) + DO 230 J=IJJ0(I),IJJ0(I)-NJJ0(I)+1,-1 + IGAR=IGAR+1 + ASTOC(I,J)=ASTOC(I,J)-SCAT0(IGAR) + 230 CONTINUE + IF(LFISSI) THEN + DO 240 J=1,NGRO + ASTOC(I,J)=ASTOC(I,J)-XHI(I)*SFNU(J)/REFKEF + 240 CONTINUE + ENDIF + 250 CONTINUE + B2OLD=B2 + CALL ALEIGD(ASTOC,BSTOC,NGRO,B2,ASTOC(1,NGRO+1),EPS,IT) + IF(IMPX.GT.1) WRITE (6,440) ITEX,IT,B2 + BIL1=SQRT(DOT_PRODUCT(ASTOC(:NGRO,NGRO+1),ASTOC(:NGRO,NGRO+1))) + BIL2=0.0D0 + DO 260 I=1,NGRO + BIL2=MAX(BIL2,ABS(PHI(I)-ASTOC(I,NGRO+1)/BIL1)) + PHI(I)=0.5*(PHI(I)+ASTOC(I,NGRO+1)/BIL1) + 260 CONTINUE + ERR3=ABS(B2-B2OLD) + IF((BIL2.GE.10.0*EPS*BIL1).OR.(ERR3.GE.EPS)) GO TO 170 +*---- +* COMPUTE THE LEAKAGE COEFFICIENTS +*---- + 300 IF((OPTION.EQ.'P1').OR.(OPTION.EQ.'B1')) THEN + IGAR=0 + DO 325 I=1,NGRO + ST2=DBLE(ST(I)) + DD=DBLE(DHOM(I)) + BETA=B1BETA(IAPROX,B2,ST2,DD) + BETA=BETA*ST2/(1.0D0-B2*BETA) + DO 310 J=1,NGRO + ASTOC(I,J)=0.0D0 + 310 CONTINUE + ASTOC(I,I)=1.0D0 + ASTOC(I,NGRO+1)=BETA + DO 320 J=IJJ1(I),IJJ1(I)-NJJ1(I)+1,-1 + IGAR=IGAR+1 + ASTOC(I,J)=ASTOC(I,J)-3.0D0*BETA*SCAT1(IGAR)*PHI(J)/PHI(I) + 320 CONTINUE + 325 CONTINUE + CALL B1SOL(NGRO,ASTOC,IER) + IF(IER.NE.0) CALL XABORT('B1DIF: SINGULAR MATRIX(2).') + DO 330 I=1,NGRO + DHOM(I)=REAL(ASTOC(I,NGRO+1)) + 330 CONTINUE + ELSE IF((OPTION.NE.'LKRD').AND.(OPTION.NE.'RHS')) THEN + DO 340 I=1,NGRO + ST2=ST(I) + GAMMA(I)=REAL(B1GAMA(IAPROX,B2,ST2)) + IF(IAPROX.NE.0) DHOM(I)=REAL(1.0D0/(3.0D0*GAMMA(I)*ST2)) + 340 CONTINUE + ENDIF + A2=0.0D0 + CAET=0.0D0 + ZXC=0.0D0 + DO 350 I=1,NGRO + A2=A2+DHOM(I)*PHI(I) + CAET=CAET+SA(I)*PHI(I) + IF(LFISSI) ZXC=ZXC+SFNU(I)*PHI(I)/REFKEF + ST2=DBLE(ST(I)) + GAMMA(I)=REAL(B1GAMA(IAPROX,B2,ST2)) + 350 CONTINUE + A2=A2/CAET + CAET=REFKEF*ZXC/CAET + IF(CAET.EQ.0.0) THEN + ALAM2=0.0 + ELSE + ALAM2=CAET/(1.0D0+A2*B2) + ENDIF + IF(TYPE.EQ.'L') ALAM1=ALAM2 + IF(IMPX.GT.0) WRITE (6,450) ITEX,B2,ALAM1,ALAM2,CAET,A2 +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(SA,CSTOC,BSTOC,ASTOC) + RETURN +* + 400 FORMAT(/42H B1DIF: DIFFUSION COEFFICIENT CALCULATION.) + 410 FORMAT(/20H B1DIF: SOLUTION OF ,A4,21H EQUATIONS WITH TYPE ,A4) + 415 FORMAT(47H B1DIF: THE INITIAL BUCKLING WAS INCREASED FROM,1P, + 1 E13.5,3H TO,E13.5) + 420 FORMAT(26H B1DIF: INITIAL BUCKLING =,1P,E13.5) + 430 FORMAT(33H B1DIF: K-EFFECTIVE ITERATION NO.,I3,13H. K-EFFECTIVE, + 1 2H =,F10.6,11H BUCKLING =,1P,E13.5) + 440 FORMAT(30H B1DIF: BUCKLING ITERATION NO.,I3,13H CONVERGED IN,I5, + 1 29H INNER ITERATIONS. BUCKLING =,1P,E13.5) + 450 FORMAT(8X,22HNUMBER OF ITERATIONS =,I3/8X,10HBUCKLING =,1P,E13.5, + 1 0P/8X,13HK-EFFECTIVE =,F10.6,3H (,F10.6,2H )/8X,12HK-INFINITE =, + 2 F10.6/8X,16HMIGRATION AREA =,1P,E13.5/) + END -- cgit v1.2.3