diff options
Diffstat (limited to 'Dragon/src/FLUGPI.f')
| -rw-r--r-- | Dragon/src/FLUGPI.f | 439 |
1 files changed, 439 insertions, 0 deletions
diff --git a/Dragon/src/FLUGPI.f b/Dragon/src/FLUGPI.f new file mode 100644 index 0000000..91c5ef2 --- /dev/null +++ b/Dragon/src/FLUGPI.f @@ -0,0 +1,439 @@ +*DECK FLUGPI + SUBROUTINE FLUGPI(IPFLUX,IPMACR,ITYPEC,MAXOUT,MAXINR,EPSOUT, + > EPSUNK,EPSINR,IREBAL,IFRITR,IACITR,COPTIO, + > ILEAK,B2,NGROUP,NREGIO,NMAT,NIFISS,LEAKSW, + > REFKEF,ITPIJ,IPRINT,REC,INITFL,NMERG,IMERG) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Read data for flux solution operator. +* +*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): G. Marleau +* +*Parameters: input +* IPFLUX pointer to the flux LCM object. +* IPMACR pointer to the macrolib LCM object. +* NGROUP number of energy groups. +* NREGIO number of regions. +* NMAT number of material mixtures. +* NIFISS number of fissile isotopes. +* LEAKSW leakage flag (=.true. if leakage is present on the outer +* surface). +* 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. +* REC flux recovery flag: +* =.true. recover the existing solution as initial estimate; +* =.false. use a new initial estimate. +* +*Parameters: output +* ITYPEC type of flux evaluation: +* =-2 Fourier analysis; +* =-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; +* MAXOUT maximum number of outer iterations. +* MAXINR maximum number of thermal iterations. +* EPSOUT outer iterations eigenvalue epsilon. +* EPSUNK outer iterations eigenvector epsilon. +* EPSINR thermal iterations epsilon. +* IREBAL flux rebalancing flag (=1: perform rebalancing). +* IFRITR number of free iterations in an acceleration cycle. +* IACITR number of accelerated iterations in an acceleration cycle. +* COPTIO 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). +* 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; +* =17,27,37,47,57,67 for heterogeneous method with pijk and +* fixed b (17) or search b_x(27), b_y (37), b_z (47), +* b_r (57) or b_x=b_y=b_z (67). +* B2 initial or imposed directional bucklings. +* REFKEF target effective multiplication factor (K-eff). +* IPRINT print selection for FLU: module (= 0/1/2/3 no print/short +* print/long print). +* INITFL flux initialisation flags: +* = 0 flux initialisation (=0.0 or 1.0); +* = 1 flux read on LCM; +* = 2 initialization from DSA flux read on LCM. +* NMERG number of leakage zones. +* IMERG leakage zone index in each material mixture zone. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPFLUX,IPMACR + INTEGER ITYPEC,MAXOUT,MAXINR,IREBAL,IFRITR,IACITR,ILEAK, + > NGROUP,NREGIO,NMAT,NIFISS,ITPIJ,IPRINT,INITFL,NMERG, + > IMERG(NMAT) + REAL EPSOUT,EPSUNK,EPSINR,B2(4) + CHARACTER COPTIO*4 + LOGICAL LEAKSW,REC + DOUBLE PRECISION REFKEF +*---- +* LOCAL VARIABLES +*---- + PARAMETER (NBUCKN=7,NLEAK=7,NSDIR=5,ROUMIN=1.0E-5,RINMIN=5.0E-5, + > NSTATE=40) + CHARACTER CARLIR*4,CTYPEC*1,CBUCKN(0:NBUCKN)*4,CLEAK(NLEAK)*6, + > CSDIR(NSDIR)*1 + SAVE CBUCKN,CLEAK,CSDIR + INTEGER ITYPLU,INTLIR,IRSDIR(NSDIR),ISTATE(NSTATE) + REAL REALIR,BSDIR(NSDIR),EPSCON(5) + DOUBLE PRECISION DBLINP + DATA (CBUCKN(JJ),JJ=0,NBUCKN) + > /'LKRD','RHS','B0','P0','B1','P1','B0TR','P0TR'/ + DATA (CLEAK(JJ),JJ=1,NLEAK) + > /'PNLR','PNL','SIGS','ALBS','HETE','ECCO','TIBERE'/ + DATA (CSDIR(III),III=1,NSDIR) + > /'X','Y','Z','R','G'/ +*---- +* INITIALIZE TO DEFAULT VALUE +*---- + ISDIR=0 + COPTIO='B0' + IF((ITPIJ.EQ.1).OR.(ITPIJ.GE.3)) ISDIR=5 + DO 10 III=1,NSDIR + BSDIR(III)=0.0 + IRSDIR(III)=0 + 10 CONTINUE + REFKEF=1.0D0 + IPRINT=1 + IF(REC) THEN + CALL LCMGET(IPFLUX,'STATE-VECTOR',ISTATE) + ITYPEC=ISTATE(6) + ILEAK=ISTATE(7) + IFRITR=ISTATE(8) + IACITR=ISTATE(9) + IREBAL=ISTATE(10) + MAXOUT=ISTATE(12) + NMERG=ISTATE(18) + CALL LCMGET(IPFLUX,'EPS-CONVERGE',EPSCON) + EPSINR=EPSCON(1) + EPSUNK=EPSCON(2) + EPSOUT=EPSCON(3) + INITFL=1 + CALL LCMLEN(IPFLUX,'B2 HETE',ILCML1,ITYLCM) + CALL LCMLEN(IPFLUX,'B2 B1HOM',ILCML2,ITYLCM) + IF(ILCML1.EQ.3) THEN + CALL LCMGET(IPFLUX,'B2 HETE',BSDIR) + IRSDIR(1)=1 + IRSDIR(2)=1 + IRSDIR(3)=1 + ELSE IF(ILCML2.EQ.1) THEN + CALL LCMGET(IPFLUX,'B2 B1HOM',BSDIR(5)) + IRSDIR(5)=1 + ENDIF + IF(NMERG.GT.0) CALL LCMGET(IPFLUX,'IMERGE-LEAK',IMERG) + ELSE + ITYPEC=-99 + ILEAK=0 + IFRITR=3 + IACITR=3 + IREBAL=1 + MAXOUT=0 + EPSINR=RINMIN + EPSUNK=0.0 + EPSOUT=ROUMIN + INITFL=0 + NMERG=1 + IMERG(:NMAT)=1 + CALL LCMPUT(IPFLUX,'IMERGE-LEAK',NMAT,1,IMERG) + ENDIF + IF(NGROUP.EQ.1) MAXINR=1 + IF(MOD(ITPIJ,2).EQ.0) THEN + MAXINR=4*NGROUP + ELSE + MAXINR=2*NGROUP + ENDIF +*---- +* READ OPTION NAME +*---- + 20 CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.EQ.10) GO TO 140 + 30 IF(ITYPLU.NE.3) + > CALL XABORT('FLUGPI: READ ERROR - CHARACTER VARIABLE EXPECTED') + IF(CARLIR.EQ.';') THEN + GO TO 140 + ELSE IF(CARLIR.EQ.'EDIT') THEN + CALL REDGET(ITYPLU,IPRINT,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.1) CALL XABORT('FLUGPI: READ ERROR - INTEGER VA' + > //'RIABLE EXPECTED') + ELSE IF(CARLIR.EQ.'TYPE') THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CTYPEC,DBLINP) + IF(ITYPLU.NE.3) CALL XABORT('FLUGPI: READ ERROR - CHARACTER ' + > //'VARIABLE EXPECTED') + IF(CTYPEC.EQ.'N') THEN + ITYPEC=-1 + ELSE IF(CTYPEC.EQ.'S') THEN + ITYPEC=0 + ELSE IF(CTYPEC.EQ.'P') THEN + ITYPEC=1 + ELSE IF(CTYPEC.EQ.'K') THEN + ITYPEC=2 + ELSE IF(CTYPEC.EQ.'B') THEN + ITYPEC=4 + ILEAK=3 + ELSE IF(CTYPEC.EQ.'L') THEN + ITYPEC=5 + ILEAK=3 + ELSE IF(CTYPEC.EQ.'F') THEN + ITYPEC=-2 + ILEAK=0 + MAXOUT=1 + MAXINR=1 + ELSE + CALL XABORT('FLUGPI: READ ERROR - INVALID TYPE KEYWORD= ' + > //CTYPEC//' -- ONLY VALUES ALLOWED ARE: N,S,K,B,L OR F') + ENDIF + IF(ITYPEC.GE.2) THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.3) CALL XABORT('FLUGPI: READ ERROR - CHARACTE' + > //'R VARIABLE EXPECTED') + DO 40 JBUC=0,NBUCKN + IF(CARLIR.EQ.CBUCKN(JBUC)) THEN + COPTIO=CARLIR + GO TO 50 + ENDIF + 40 CONTINUE + GO TO 30 + 50 IF(ITYPEC.EQ.2) ITYPEC=3 + ILEAK=3 + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.3) CALL XABORT('FLUGPI: READ ERROR - CHARACTE' + > //'R VARIABLE EXPECTED') + DO 70 JLEA=1,NLEAK + IF(CARLIR.EQ.CLEAK(JLEA)(:4)) THEN + ILEAK=JLEA + IF(LEAKSW.AND.(ILEAK.NE.3).AND.(ILEAK.NE.5)) THEN + CALL XABORT('FLUGPI: FUNDAMENTAL MODE EXPECTED WITH A ' + > //'LEAKAGE MODEL OTHER THAN SIGS OR HETE.') + ENDIF + IF(ILEAK.EQ.5) THEN + DO IBM=1,NMAT + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF((IBM.EQ.1).AND.(ITYPLU.EQ.3)) GO TO 30 + IF(ITYPLU.NE.1) CALL XABORT('FLUGPI: READ ERROR - IN' + > //'TEGER VARIABLE EXPECTED') + IMERG(IBM)=INTLIR + NMERG=MAX(NMERG,IMERG(IBM)) + ENDDO + CALL LCMPUT(IPFLUX,'IMERGE-LEAK',NMAT,1,IMERG) + ELSE IF(ILEAK.EQ.7) THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.3) CALL XABORT('FLUGPI: READ ERROR - ' + > //'CHARACTER VARIABLE EXPECTED') + DO 60 III=1,NSDIR + IF(CARLIR.EQ.CSDIR(III)) THEN + ISDIR=III + GO TO 20 + ENDIF + 60 CONTINUE + GO TO 30 + ENDIF + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + GO TO 30 + ENDIF + 70 CONTINUE + GO TO 30 + ENDIF + ELSE IF(CARLIR.EQ.'REBA') THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.3) CALL XABORT('FLUGPI: READ ERROR - CHARACTE' + > //'R VARIABLE EXPECTED') + IF(CARLIR.EQ.'OFF ') THEN + IREBAL=0 + ELSE + IREBAL=1 + GO TO 30 + ENDIF + ELSE IF(CARLIR.EQ.'INIT') THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.3) CALL XABORT('FLUGPI: READ ERROR - CHARACTE' + > //'R VARIABLE EXPECTED') + IF(CARLIR.EQ.'OFF') THEN +* initial flat distribution + INITFL=0 + ELSE IF(CARLIR.EQ.'ON') THEN +* use LCM flux + INITFL=1 + ELSE IF(CARLIR.EQ.'DSA') THEN +* use DSA flux + INITFL=2 + ELSE + CALL XABORT('FLUGPI: OFF/ON/DSA KEYWORD EXPECTED') + ENDIF + ELSE IF(CARLIR.EQ.'EXTE') THEN + DO 90 II=1,3 + CALL REDGET(ITYPLU,NITMA,EPSOUT,CARLIR,DBLINP) + IF(ITYPLU.EQ.1) MAXOUT=NITMA + IF(ITYPLU.EQ.3) GO TO 30 + 90 CONTINUE + ELSE IF(CARLIR.EQ.'UNKT') THEN + CALL REDGET(ITYPLU,INTLIR,EPSUNK,CARLIR,DBLINP) + IF(ITYPLU.NE.2) + > CALL XABORT('FLUGPI: REAL VALUE OF EPSUNK MUST FOLLOW' + > //' UNKT') + GO TO 20 + ELSE IF(CARLIR.EQ.'THER') THEN + DO 100 II=1,3 + CALL REDGET(ITYPLU,NITMA,EPSINR,CARLIR,DBLINP) + IF(ITYPLU.EQ.1) MAXINR=NITMA + IF(ITYPLU.EQ.3) GO TO 30 + 100 CONTINUE + ELSE IF(CARLIR.EQ.'ACCE') THEN + CALL REDGET(ITYPLU,IFRITR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.1) CALL XABORT('FLUGPI: READ ERROR - INTEGER VA' + > //'RIABLE EXPECTED') + CALL REDGET(ITYPLU,IACITR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.1) CALL XABORT('FLUGPI: READ ERROR - INTEGER VA' + > //'RIABLE EXPECTED') + ELSE IF(CARLIR.EQ.'KEFF') THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.2) CALL XABORT('FLUGPI: READ ERROR - REAL VA' + > //'RIABLE EXPECTED FOLLOWING KEFF KEYWORD') + REFKEF=REALIR + ELSE IF(CARLIR.EQ.'BUCK') THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IRSDIR(:NSDIR)=0 + IF(ITYPLU.EQ.2) THEN + BSDIR(5)=REALIR + IRSDIR(5)=1 + GO TO 20 + ELSE IF(ITYPLU.EQ.1) THEN + CALL XABORT('FLUGPI: READ ERROR - INTEGER ' + > //'VARIABLE FOUND FOLLOWING BUCK KEYWORD') + ENDIF + 110 CONTINUE + DO 120 III=1,NSDIR + IF(CARLIR.EQ.CSDIR(III)) THEN + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.2) + > CALL XABORT('FLUGPI: READ ERROR - REAL VARIABLE '// + > 'EXPECTED FOLLOWING BUCKLING DIRECTION KEYWORD') + BSDIR(III)=REALIR + IRSDIR(III)=1 + GO TO 130 + ENDIF + 120 CONTINUE + GO TO 30 + 130 CONTINUE + CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLINP) + IF(ITYPLU.NE.3) + > CALL XABORT('FLUGPI: READ ERROR - CHARACTER ' + > //'VARIABLE EXPECTED') + GO TO 110 + ELSE IF(CARLIR.EQ.'IDEM') THEN + IRSDIR(:NSDIR)=0 + CALL LCMLEN(IPMACR,'B2 HETE',ILCMLN,ITYLCM) + IF(ILCMLN.EQ.3) THEN + CALL LCMGET(IPMACR,'B2 HETE',BSDIR) + IRSDIR(1)=1 + IRSDIR(2)=1 + IRSDIR(3)=1 + ENDIF + CALL LCMLEN(IPMACR,'B2 B1HOM',ILCMLN,ITYLCM) + IF(ILCMLN.EQ.1) THEN + CALL LCMGET(IPMACR,'B2 B1HOM',BSDIR(5)) + IRSDIR(5)=1 + ENDIF + ELSE + CALL XABORT('FLUGPI: READ ERROR - ILLEGAL KEYWORD '//CARLIR) + ENDIF + GO TO 20 +*---- + 140 CONTINUE + IF(ITYPEC.EQ.3) THEN + ISDIR=0 + IF(ILEAK.EQ.7) THEN + DO 150 III=1,NSDIR + IF(IRSDIR(III).EQ.1) GO TO 160 + 150 CONTINUE + ELSE + GO TO 160 + ENDIF + CALL XABORT('FLUGPI: NO BUCKLING READ FOR TYPE K '// + > 'CALCULATION WITH IMPOSED BUCKLING') + 160 CONTINUE + ENDIF + IF(ILEAK.EQ.7) THEN + IF(IRSDIR(5).EQ.1) THEN + DO 210 III=1,NSDIR-1 + IF(IRSDIR(III).EQ.1) THEN + CALL XABORT('FLUGPI: GLOBAL INITIAL BUCKLING '// + > 'INCONSISTENT WITH X, Y, Z, R BUCKLING') + ENDIF + 210 CONTINUE + B2(1)=BSDIR(5)/3.0 + B2(2)=B2(1) + B2(3)=B2(1) + B2(4)=BSDIR(5) + ELSE IF(IRSDIR(4).EQ.1) THEN + DO 220 III=1,2 + IF(IRSDIR(III).EQ.1) THEN + CALL XABORT('FLUGPI: RADIAL INITIAL BUCKLING '// + > 'INCONSISTENT WITH X, Y BUCKLING') + ENDIF + 220 CONTINUE + B2(1)=BSDIR(4)/2.0 + B2(2)=B2(1) + B2(3)=BSDIR(3) + B2(4)=BSDIR(3)+BSDIR(4) + ELSE + B2(1)=BSDIR(1) + B2(2)=BSDIR(2) + B2(3)=BSDIR(3) + B2(4)=BSDIR(1)+BSDIR(2)+BSDIR(3) + ENDIF + ILEAK=(ISDIR+1)*10+ILEAK + ELSE + IF(IRSDIR(4).NE.0.OR.IRSDIR(3).NE.0.OR. + > IRSDIR(2).NE.0.OR.IRSDIR(1).NE.0) + > CALL XABORT('FLUGPI: FOR HOMOGENEOUS LEAKAGE METHOD'// + > 'DIRECTIONS X ,Y, Z AND R BUCKLING ARE ILLEGAL') + B2(4)=BSDIR(5) + ENDIF + IF(EPSOUT.LT.1.0E-10) THEN + CALL XABORT('FLUGPI: ERROR -- EPSOUT MUST BE GREATER '// + > 'THAN 1.0E-10') + ENDIF + IF(EPSINR.LT.1.0E-10) THEN + CALL XABORT('FLUGPI: ERROR -- EPSINR MUST BE GREATER '// + > 'THAN 1.0E-10') + ENDIF + IF(EPSUNK.LE.1.0E-10) EPSUNK=EPSINR + IF(ITYPEC.EQ.-99) CALL XABORT('FLUGPI: TYPE NOT DEFINED.') + IF(MAXOUT.EQ.0)THEN + IF(ITYPEC.LE.2) THEN + MAXOUT=MAX(2*NREGIO-1,2*NIFISS-1) + ELSE + MAXOUT=MAX(10*NREGIO+1,10*NIFISS+1) + ENDIF + ENDIF + RETURN + END |