diff options
Diffstat (limited to 'Dragon/src/FLUDRV.f')
| -rw-r--r-- | Dragon/src/FLUDRV.f | 345 |
1 files changed, 345 insertions, 0 deletions
diff --git a/Dragon/src/FLUDRV.f b/Dragon/src/FLUDRV.f new file mode 100644 index 0000000..428113d --- /dev/null +++ b/Dragon/src/FLUDRV.f @@ -0,0 +1,345 @@ +*DECK FLUDRV + SUBROUTINE FLUDRV(IPRT,IPFLUX,IPTRK,IPMACR,IPSOU,IFTRAK,IPSYS, + 1 IPHASE,ITPIJ,CXDOOR,ITRANC,TITLE,B2,INITFL,LFORW,LEAKSW,IREBAL, + 2 NGRP,NMAT,NIFIS,NANIS,NLF,NLIN,NFUNL,OPTION,NUN,MAXINR,EPSINR, + 3 MAXOUT,EPSUNK,EPSOUT,IFRITR,IACITR,ITYPEC,ILEAK,NREG,NSOUT, + 4 MATCOD,KEYFLX,VOL,REFKEF,NMERG,IMERG) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Driver for Boltzmann equation solvers. +* +*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): R. Roy +* +*Parameters: input +* IPRT print flag. +* IPFLUX pointer to the flux LCM object. +* IPTRK pointer to the tracking LCM object. +* IPMACR pointer to the macrolib LCM object. +* IPSOU pointer to the fixed source LCM object. +* IFTRAK tracking file unit number. +* IPSYS pointer to the system LCM object (=0 for the method of +* characteristics). +* IPHASE 1 for asm 2 for pij. +* ITPIJ type of collision probability information available: +* =1 scattering modified pij (wij); +* =2 standard pij; +* =3 scattering modified pij+pijk (wij,wijk); +* =4 standard pij+pijk. +* CXDOOR name of the flux solution door. +* ITRANC type of transport correction (>0 to perform a correction). +* TITLE title. +* B2 initial or imposed directional bucklings. +* INITFL flux initialization flag (=0/1/2: uniform flux/LCM/DSA). +* LFORW flag set to .false. to solve an adjoint problem. +* LEAKSW leakage flag (=.true. if leakage is present on the outer +* surface). +* IREBAL flux rebalancing flag (=1: perform rebalancing). +* NGRP number of energy groups. +* NMAT number of mixtures. +* NIFIS number of fissile isotopes. +* NANIS maximum cross section Legendre order. +* NLF number of Legendre orders for the flux. +* NLIN number of polynomial components in flux spatial expansion. +* NFUNL number of spherical harmonics components. +* OPTION type of leakage coefficients: +* 'LKRD' (recover leakage coefficients in Macrolib); +* 'RHS' (recover leakage coefficients in RHS flux object); +* 'B0' (B-0), 'P0' (P-0), 'B1' (B-1), +* 'P1' (P-1), 'B0TR' (B-0 with transport correction) or 'P0TR' +* (P-0 with transport correction). +* NUN number of unknowns per energy group including spherical +* harmonic terms, interface currents and fundamental +* currents. +* MAXINR maximum number of thermal iterations. +* EPSINR thermal iterations epsilon. +* MAXOUT maximum number of outer iterations. +* EPSUNK outer iterations eigenvector epsilon. +* EPSOUT outer iterations eigenvalue epsilon. +* IFRITR number of free iterations in an acceleration cycle. +* IACITR number of accelerated iterations in an acceleration cycle. +* ITYPEC type of flux evaluation: +* =-1 skip the flux calculation; +* =0 fixed sources; +* =1 fixed source eigenvalue problem (GPT type); +* =2 fission sources/K-eff convergence; +* =3 fission sources/K-eff convergence/db2 buckling evaluation; +* =4 fission sources/db2 buckling convergence; +* =5 b2 sources/db2 buckling convergence. +* ILEAK method used to include DB2 effect: +* =1 the scattering modified cp matrix is multiplied by PNLR; +* =2 the reduced cp matrix is multiplied by PNL; +* =3 sigs0-db2 approximation; +* =4 albedo approximation; +* =5 Todorova-type isotropic streaming model; +* =6 Ecco-type isotropic streaming model; +* >6 Tibere type anisotropic streaming model. +* NREG number of regions. +* NSOUT number of outer surfaces. +* MATCOD mixture indices. +* KEYFLX index of L-th order flux components in unknown vector. +* VOL volumes. +* REFKEF target effective multiplication factor (K-eff). +* NMERG number of leakage zones. +* IMERG leakage zone index in each material mixture zone. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + CHARACTER CXDOOR*12,TITLE*72,OPTION*4 + TYPE(C_PTR) IPFLUX,IPTRK,IPMACR,IPSOU,IPSYS + INTEGER IPRT,IFTRAK,IPHASE,ITPIJ,ITRANC,INITFL,IREBAL,NGRP, + > NMAT,NIFIS,NANIS,NLF,NLIN,NFUNL,NUN,MAXINR,MAXOUT, + > IFRITR,IACITR,ITYPEC,ILEAK,NREG,NSOUT,MATCOD(NREG), + > KEYFLX(NREG,NLIN,NFUNL),NMERG,IMERG(NMAT) + REAL EPSUNK,EPSINR,B2(4),VOL(NREG) + LOGICAL LFORW,LEAKSW + DOUBLE PRECISION REFKEF +*---- +* LOCAL VARIABLES +*---- + PARAMETER (IUNOUT=6,NSTATE=40,NDIMO=2,IGPT=0) + TYPE(C_PTR) JPMACR,KPMACR,JPFLUX,IPFLUP,JPSYS,KPSYS + LOGICAL LREBAL + INTEGER ISTATE(NSTATE) + REAL EPSCON(5) + CHARACTER CAN(0:19)*2 +*---- +* ALLOCATABLE ARRAYS +*---- + REAL, ALLOCATABLE, DIMENSION(:) :: SIGT,SIGS0 + REAL, ALLOCATABLE, DIMENSION(:,:) :: FLUXO,XSTRC,XSTK + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XSDIA,XSCHI,XSNUF +*---- +* DATA STATEMENTS +*---- + SAVE CAN + DATA CAN /'00','01','02','03','04','05','06','07','08','09', + > '10','11','12','13','14','15','16','17','18','19'/ +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(FLUXO(NUN,NGRP),XSTRC(0:NMAT,NGRP), + > XSDIA(0:NMAT,0:NANIS,NGRP),XSCHI(0:NMAT,NIFIS,NGRP), + > XSNUF(0:NMAT,NIFIS,NGRP),XSTK(NMAT,NIFIS)) +* + IF(IPRT.GE.3) THEN + WRITE(IUNOUT,6000) + WRITE(IUNOUT,6001) (IREGIO,VOL(IREGIO),MATCOD(IREGIO), + > IREGIO=1,NREG) + ENDIF + IPFLUP=C_NULL_PTR +*---- +* RECOVER CROSS SECTIONS. +*---- + JPMACR=LCMGID(IPMACR,'GROUP') + DO 110 IGR=1,NGRP + KPMACR=LCMGIL(JPMACR,IGR) + DO 10 IANI=0,NANIS + CALL LCMLEN(KPMACR,'NJJS'//CAN(IANI),ILNLCM,ITYLCM) + IF(ILNLCM.NE.NMAT) THEN + CALL LCMLIB(KPMACR) + CALL XABORT('FLUDRV: FLUX CALCULATION ERROR, SCATTERING'// + > ' MATRIX OF ORDER ANIS ='//CAN(IANI)//' NOT ON LCM') + ENDIF + 10 CONTINUE +*---- +* RECOVER CORRECTED TOTAL AND WITHIN-GROUP SCATTERING CROSS SECTIONS +*---- + JPSYS=LCMGID(IPSYS,'GROUP') + KPSYS=LCMGIL(JPSYS,IGR) + CALL LCMLEN(KPSYS,'DRAGON-TXSC',ILCTX,ITYLCM) + CALL LCMLEN(KPSYS,'DRAGON-S0XSC',ILCS0X,ITYLCM) + IF(ILCTX.NE.NMAT+1) THEN + CALL XABORT('FLUDRV: INVALID LENGTH FOR DRAGON-TXSC.') + ELSE IF(MOD(ILCS0X,NMAT+1).NE.0) THEN + CALL XABORT('FLUDRV: INVALID LENGTH FOR DRAGON-S0XSC.') + ENDIF + NANI_ASM=ILCS0X/(NMAT+1)-1 + ALLOCATE(SIGT(ILCTX),SIGS0(ILCS0X)) + SIGS0(:ILCS0X)=0.0 + CALL LCMGET(KPSYS,'DRAGON-TXSC',SIGT) + CALL LCMGET(KPSYS,'DRAGON-S0XSC',SIGS0) + XSTRC(0:NMAT,IGR)=SIGT(:NMAT+1) + XSDIA(0:NMAT,0:NANIS,IGR)=0.0 + DO IANI=0,MIN(NANIS,NANI_ASM) + DO IMAT=1,NMAT + XSDIA(IMAT,IANI,IGR)=SIGS0(IANI*(NMAT+1)+IMAT+1) + ENDDO + ENDDO + DEALLOCATE(SIGS0,SIGT) + IF(IPRT.GE.3) THEN + WRITE(IUNOUT,6002) IGR + WRITE(IUNOUT,6003) (IMAT,XSTRC(IMAT,IGR), + > XSDIA(IMAT,0,IGR),IMAT=1,NMAT) + ENDIF +*---- +* RECOVER FISSION CROSS SECTIONS +*---- + CALL LCMLEN(KPMACR,'CHI',ILONG,ITYLCM) + IF( ILONG.EQ.0 )THEN + IF(( ITYPEC.NE.0 ).AND.( ITYPEC.NE.5 ) + 1 .AND.( ITYPEC.NE.-2 ))THEN + CALL XABORT('FLUDRV: TYPE S, F OR L REQUESTED') + ENDIF + XSCHI(0:NMAT,:NIFIS,:NGRP)=0.0 + XSNUF(0:NMAT,:NIFIS,:NGRP)=0.0 + ELSE + CALL LCMGET(KPMACR,'CHI',XSTK) + DO 60 IFIS= 1, NIFIS + XSCHI(0,IFIS,IGR)= 0.0 + DO 50 IMAT= 1, NMAT + XSCHI(IMAT,IFIS,IGR)= XSTK(IMAT,IFIS) + 50 CONTINUE + 60 CONTINUE + CALL LCMGET(KPMACR,'NUSIGF',XSTK) + DO 80 IFIS= 1, NIFIS + XSNUF(0,IFIS,IGR)= 0.0 + DO 70 IMAT= 1, NMAT + XSNUF(IMAT,IFIS,IGR)= XSTK(IMAT,IFIS) + 70 CONTINUE + 80 CONTINUE + IF( IPRT.GT.3 )THEN + WRITE(IUNOUT,6004) (IMAT,XSNUF(IMAT,1,IGR), + > XSCHI(IMAT,1,IGR),IMAT=1,NMAT) + ENDIF + ENDIF + DO 100 IANI=0,NANIS + CALL LCMLEN(KPMACR,'NJJS'//CAN(IANI),ILNLCM,ITYLCM) + IF(ILNLCM.NE.NMAT) THEN + CALL LCMLIB(KPMACR) + CALL XABORT('FLUDRV: FLUX CALCULATION ERROR, SCATTERING'// + > ' MATRIX OF ORDER ANIS ='//CAN(IANI)//' NOT ON LCM') + ENDIF + 100 CONTINUE + 110 CONTINUE +*---- +* FLUX INITIALIZATION +*---- + IF(LFORW) THEN + CALL LCMLEN(IPFLUX,'FLUX',ILINIT,ITYLCM) + ELSE + CALL LCMLEN(IPFLUX,'AFLUX',ILINIT,ITYLCM) + ENDIF + IF((ILINIT.EQ.0).OR.(INITFL.EQ.0)) THEN + IF(LFORW) THEN + JPFLUX=LCMLID(IPFLUX,'FLUX',NGRP) + ELSE + JPFLUX=LCMLID(IPFLUX,'AFLUX',NGRP) + ENDIF + DO 130 IGR=1,NGRP + FLUXO(:NUN,IGR)=0.0 + IF(ITYPEC.GT.0) THEN + IF((CXDOOR.EQ.'BIVAC').OR.(CXDOOR.EQ.'TRIVAC')) THEN + FLUXO(:NUN,IGR)=1.0 + ELSE + DO 120 IREGIO=1,NREG + IND=KEYFLX(IREGIO,1,1) + IF(IND.GT.0) FLUXO(IND,IGR)=1.0 + 120 CONTINUE + ENDIF + ENDIF + IF(LFORW) THEN + CALL LCMPDL(JPFLUX,IGR,NUN,2,FLUXO(1,IGR)) + ELSE + CALL LCMPDL(JPFLUX,NGRP-IGR+1,NUN,2,FLUXO(1,IGR)) + ENDIF + 130 CONTINUE + ENDIF + IF(ITYPEC.GE.2) THEN + EIGENK=REAL(REFKEF) + CALL LCMPUT(IPFLUX,'K-EFFECTIVE',1,2,EIGENK) + ENDIF + IF(ITYPEC.GE.3) THEN + CALL LCMPUT(IPFLUX,'B2 B1HOM',1,2,B2(4)) + IF(ILEAK.GE.7) CALL LCMPUT(IPFLUX,'B2 HETE',3,2,B2) + ENDIF + IF(ITYPEC.EQ.-1) GO TO 1001 +* + IF(ILEAK.GE.7.AND.ITYPEC.GE.3) THEN + IF(ITRANC.EQ.0) THEN + IF(OPTION.EQ.'B0TR'.OR.OPTION.EQ.'P0TR'.OR.OPTION.EQ.'LKRD' + > .OR.OPTION.EQ.'RHS') + > CALL XABORT('FLUDRV: ILLEGAL OPTION = '//OPTION// + > ' FOR HETEROGENEOUS LEAKAGE CALCULATION'// + > ' WITHOUT TRANSPORT CORRECTED CROSS SECTIONS') + ENDIF + ENDIF +* + IF (CXDOOR.EQ.'MCCG') THEN + CALL LCMLEN(IPTRK,'KEYCUR$MCCG',ICREB,ITYLCM) + CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) + IF ((ICREB.GT.0).AND.(ISTATE(24).EQ.0)) THEN + LREBAL=(IREBAL.EQ.1) + ELSE + LREBAL=(IREBAL.EQ.1).AND.(.NOT.LEAKSW) + ENDIF + ELSE + LREBAL=(IREBAL.EQ.1).AND.(.NOT.LEAKSW) + ENDIF + CALL FLU2DR(IPRT,IPMACR,IPFLUX,IPSYS,IPTRK,IPFLUP,IPSOU,IGPT, + 1 IFTRAK,CXDOOR,TITLE,NUN,NREG,NSOUT,NANIS,NLF,NLIN,NFUNL,NGRP, + 2 NMAT,NIFIS,LFORW,LEAKSW,MAXINR,EPSINR,MAXOUT,EPSUNK,EPSOUT, + 3 IFRITR,IACITR,ITYPEC,IPHASE,ITPIJ,ILEAK,OPTION,REFKEF,MATCOD, + 4 KEYFLX,VOL,XSTRC,XSDIA,XSNUF,XSCHI,LREBAL,INITFL,NMERG,IMERG) +* + 1001 CALL LCMLEN(IPFLUX,'FLUX',ILON1,ITYLCM) + CALL LCMLEN(IPFLUX,'AFLUX',ILON2,ITYLCM) + ISTATE(:NSTATE)=0 + ISTATE(1)=NGRP + ISTATE(2)=NUN + IF((ILON1.GT.0).AND.(ILON2.GT.0)) THEN + ISTATE(3)=11 + ELSE IF(ILON1.GT.0) THEN + ISTATE(3)=1 + ELSE IF(ILON2.GT.0) THEN + ISTATE(3)=10 + ENDIF + ISTATE(4)=0 + ISTATE(5)=0 + ISTATE(6)=ITYPEC + ISTATE(7)=ILEAK + ISTATE(8)=IFRITR + ISTATE(9)=IACITR + ISTATE(10)=IREBAL + ISTATE(11)=MAXINR + ISTATE(12)=MAXOUT + ISTATE(17)=NMAT + ISTATE(18)=NMERG + CALL LCMPUT(IPFLUX,'STATE-VECTOR',NSTATE,1,ISTATE) + EPSCON(1)=EPSINR + EPSCON(2)=EPSUNK + EPSCON(3)=EPSOUT + EPSCON(4:5)=0.0 + CALL LCMPUT(IPFLUX,'EPS-CONVERGE',5,2,EPSCON) + CALL LCMPUT(IPFLUX,'KEYFLX',NREG,1,KEYFLX) + CALL LCMPTC(IPFLUX,'OPTION',4,OPTION) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(XSTK,XSNUF,XSCHI,XSDIA,XSTRC,FLUXO) + RETURN +*---- +* FORMATS +*---- + 6000 FORMAT(//30X,' EDITION REGION/VOL/MIXTURE '// + >3(5X,'REGION',5X,'VOL ',5X,'MIXTURE')/) + 6001 FORMAT(1P,3(1X,I8,4X,E12.5,I8,4X)) + 6002 FORMAT(//30X,' G R O U P : ',I5// + >30X,' TOTAL MACROSCOPIC CROSS SECTIONS PER MIXTURE '/) + 6003 FORMAT(3(1X,'MIXTURE',4X,'NTOT0',11X,'SIGW',3X)/ + >1P,3(1X,I4,3X,E12.5,3X,E12.5)) + 6004 FORMAT(3(1X,'MIXTURE',4X,'NUSIGF',11X,'CHI ',3X)/ + >1P,3(1X,I4,3X,E12.5,3X,E12.5)) + END |