diff options
Diffstat (limited to 'Dragon/src/FLUGPT.f')
| -rw-r--r-- | Dragon/src/FLUGPT.f | 361 |
1 files changed, 361 insertions, 0 deletions
diff --git a/Dragon/src/FLUGPT.f b/Dragon/src/FLUGPT.f new file mode 100644 index 0000000..de7d160 --- /dev/null +++ b/Dragon/src/FLUGPT.f @@ -0,0 +1,361 @@ +*DECK FLUGPT + SUBROUTINE FLUGPT(IPRT,IPFLUX,IPTRK,IPMACR,IPFLUP,IPSOU,IFTRAK, + 1 IPSYS,IPHASE,ITPIJ,CXDOOR,TITLE,INITFL,LFORW,LEAKSW,IREBAL, + 2 NGRP,NMAT,NIFIS,NANIS,NLF,NLIN,NFUNL,OPTION,NUN,MAXINR,EPSINR, + 3 MAXOUT,EPSUNK,EPSOUT,IFRITR,IACITR,ILEAK,NREG,NSOUT,MATCOD, + 4 KEYFLX,VOL,REFKEF,NMERG,IMERG) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Driver for Boltzmann equation solvers. Solution of a fixed source +* eigenvalue problem. +* +*Copyright: +* Copyright (C) 2008 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 +* IPRT print flag. +* IPFLUX pointer to the flux LCM object. +* IPTRK pointer to the tracking LCM object. +* IPMACR pointer to the macrolib LCM object. +* IPFLUP pointer to the unperturbed flux LCM object. +* IPSOU pointer to the GPT 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. +* TITLE title. +* 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. +* 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,IPFLUP,IPSOU,IPSYS + INTEGER IPRT,IFTRAK,IPHASE,ITPIJ,INITFL,IREBAL,NGRP,NMAT, + > NIFIS,NANIS,NLF,NLIN,NFUNL,NUN,MAXINR,MAXOUT, + > IFRITR,IACITR,ILEAK,NREG,NSOUT,MATCOD(NREG), + > KEYFLX(NREG,NLIN,NFUNL),NMERG,IMERG(NMAT) + REAL EPSUNK,EPSINR,VOL(NREG) + LOGICAL LFORW,LEAKSW + DOUBLE PRECISION REFKEF +*---- +* LOCAL VARIABLES +*---- + PARAMETER (IUNOUT=6,NSTATE=40,NDIMO=2,ITYPEC=1) + TYPE(C_PTR) JPMACR,KPMACR,JPFLUX,JPFLUP1,JPFLUP2,JPGPT,KPFLUX, + 1 KPGPT,JPSYS,KPSYS + LOGICAL LREBAL + INTEGER ISTATE(NSTATE) + DOUBLE PRECISION AIL,BIL,GAZ + REAL EPSCON(5) + CHARACTER CAN(0:19)*2 +*---- +* ALLOCATABLE ARRAYS +*---- + REAL, ALLOCATABLE, DIMENSION(:) :: SIGT,SIGS0 + REAL, ALLOCATABLE, DIMENSION(:,:) :: FLUXO,SOURO,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),SOURO(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 +*---- +* RECOVER CROSS SECTIONS. +*---- + JPMACR=LCMGID(IPMACR,'GROUP') + DO 100 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('FLUGPT: 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('FLUGPT: INVALID LENGTH FOR DRAGON-TXSC.') + ELSE IF(MOD(ILCS0X,NMAT+1).NE.0) THEN + CALL XABORT('FLUGPT: 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 + CALL XABORT('FLUGPT: NO FISSION SPECTRA FOUND ON MACROLIB.') + 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 + ENDIF + IF( IPRT.GT.3 )THEN + WRITE(IUNOUT,6004) (IMAT,XSNUF(IMAT,1,IGR), + > XSCHI(IMAT,1,IGR),IMAT=1,NMAT) + ENDIF + DO 90 IANI=0,NANIS + CALL LCMLEN(KPMACR,'NJJS'//CAN(IANI),ILNLCM,ITYLCM) + IF(ILNLCM.NE.NMAT) THEN + CALL LCMLIB(KPMACR) + CALL XABORT('FLUGPT: FLUX CALCULATION ERROR, '// + > 'SCATTERING MATRIX OF ORDER ANIS ='//CAN(IANI)//' NOT ON LCM') + ENDIF + 90 CONTINUE + 100 CONTINUE +*---- +* GPT FLUX INITIALIZATION +*---- + CALL LCMGET(IPSOU,'STATE-VECTOR',ISTATE) + IF(LFORW) THEN + CALL LCMLEN(IPFLUX,'DLUX',ILINIT,ITYLCM) + ELSE + CALL LCMLEN(IPFLUX,'ADFLUX',ILINIT,ITYLCM) + ENDIF + IF((ILINIT.EQ.0).OR.(INITFL.EQ.0)) THEN + IF(LFORW) THEN + MAXGPT=ISTATE(3) + JPFLUX=LCMLID(IPFLUX,'DFLUX',MAXGPT) + ELSE + MAXGPT=ISTATE(4) + JPFLUX=LCMLID(IPFLUX,'ADFLUX',MAXGPT) + ENDIF + ENDIF +*---- +* MAIN LOOP OVER EIGENVALUE FIXED SOURCE SOLUTIONS. +*---- + CALL LCMGET(IPFLUP,'STATE-VECTOR',ISTATE) + IF(ISTATE(3).NE.11) CALL XABORT('FLUGPT: MISSING UNPERTURBED DIR' + 1 //'ECT AND/OR ADJOINT FLUXES.') + JPFLUP1=LCMGID(IPFLUP,'FLUX') + JPFLUP2=LCMGID(IPFLUP,'AFLUX') + IF(LFORW) THEN + JPGPT=LCMGID(IPSOU,'DSOUR') + ELSE + JPGPT=LCMGID(IPSOU,'ASOUR') + ENDIF + DO 1000 IGPT=1,MAXGPT + CALL LCMLEL(JPGPT,IGPT,ILONG,ITYLCM) + IF(ILONG.EQ.0) GO TO 1000 + IF(IPRT.GT.0) THEN + WRITE(IUNOUT,'(1X,29(1H-)/25H FLUGPT: GPT EQUATION NB.,I5/ + 1 1X,29(1H-))') IGPT + ENDIF + IF((ILINIT.EQ.0).OR.(INITFL.EQ.0)) THEN + KPFLUX=LCMLIL(JPFLUX,IGPT,NGRP) + DO 120 IGR=1,NGRP + FLUXO(:NUN,IGR)=0.0 + DO 110 IREGIO=1,NREG + IND=KEYFLX(IREGIO,1,1) + IF(IND.GT.0) FLUXO(IND,IGR)=1.0 + 110 CONTINUE + IF(LFORW) THEN + CALL LCMPDL(KPFLUX,IGR,NUN,2,FLUXO(1,IGR)) + ELSE + CALL LCMPDL(KPFLUX,NGRP-IGR+1,NUN,2,FLUXO(1,IGR)) + ENDIF + 120 CONTINUE + ENDIF +*---- +* RECOVER UNPERTURBED FLUXES AND VALIDATION OF THE FIXED SOURCE TERM. +*---- + IF(LFORW) THEN + KPGPT=LCMGIL(JPGPT,IGPT) + DO 130 IGR=1,NGRP + CALL LCMGDL(JPFLUP2,IGR,FLUXO(1,IGR)) + CALL LCMGDL(KPGPT,IGR,SOURO(1,IGR)) + 130 CONTINUE + ELSE + KPGPT=LCMGIL(JPGPT,IGPT) + DO 140 IGR=1,NGRP + CALL LCMGDL(JPFLUP1,IGR,FLUXO(1,IGR)) + CALL LCMGDL(KPGPT,IGR,SOURO(1,IGR)) + 140 CONTINUE + ENDIF + AIL=0.0D0 + BIL=0.0D0 + DO 155 IGR=1,NGRP + DO 150 IUN=1,NUN + GAZ=FLUXO(IUN,IGR)*SOURO(IUN,IGR) + AIL=AIL+GAZ + BIL=BIL+GAZ**2 + 150 CONTINUE + 155 CONTINUE + IF(REAL(NUN)*SQRT(BIL).LT.EPSINR) GO TO 1000 + GAZ=ABS(AIL)/SQRT(REAL(NUN)*BIL) + IF(GAZ.GT.EPSINR) CALL XABORT('FLUGPT: THE SOURCE TERM IS NOT OR' + 1 //'THOGONAL TO THE ADJOINT REFERENCE FLUX.') +* + 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,NFUNL,NGRP,NMAT, + 2 NIFIS,LFORW,LEAKSW,MAXINR,EPSINR,MAXOUT,EPSUNK,EPSOUT,IFRITR, + 3 IACITR,ITYPEC,IPHASE,ITPIJ,ILEAK,OPTION,REFKEF,MATCOD,KEYFLX, + 4 VOL,XSTRC,XSDIA,XSNUF,XSCHI,LREBAL,INITFL,NMERG,IMERG) + 1000 CONTINUE +*---- +* END OF MAIN LOOP OVER EIGENVALUE FIXED SOURCE SOLUTIONS. +*---- + ISTATE(:NSTATE)=0 + ISTATE(1)=NGRP + ISTATE(2)=NUN + IF(LFORW) THEN + ISTATE(3)=100 + ELSE + ISTATE(3)=1000 + ENDIF + ISTATE(4)=0 + ISTATE(5)=MAXGPT + 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 LCMPTC(IPFLUX,'OPTION',4,OPTION) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(XSTK,XSNUF,XSCHI,XSDIA,XSTRC,SOURO,FLUXO) + RETURN +*---- +* FORMATS +*---- + 6000 FORMAT(//30X,' EDITION REGION/VOL/MIXTURE '// + >3(5X,'REGION',5X,'VOL ',5X,'MIXTURE')/) + 6001 FORMAT(1P,3(5X,I4,4X,E12.5,4X,I4,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 |