diff options
Diffstat (limited to 'Trivac/src/KINDRV.f')
| -rwxr-xr-x | Trivac/src/KINDRV.f | 295 |
1 files changed, 295 insertions, 0 deletions
diff --git a/Trivac/src/KINDRV.f b/Trivac/src/KINDRV.f new file mode 100755 index 0000000..8ad451d --- /dev/null +++ b/Trivac/src/KINDRV.f @@ -0,0 +1,295 @@ +*DECK KINDRV + SUBROUTINE KINDRV(NEN,KEN,CMOD,NGR,NBM,NBFIS,NDG,NLF,ITY,NEL, + 1 LL4,NUN,NUP,TTF,TTP,DT,IMPH,ICL1,ICL2,NADI,ADJ,MAXOUT,EPSOUT, + 2 MAXINR,EPSINR,IFL,IPR,IEXP,INORM,IMPX,POWTOT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Driver to perform the space-time kinetics calculations. +* +*Copyright: +* Copyright (C) 2008 Ecole Polytechnique de Montreal. +* +*Author(s): D. Sekki +* +*Parameters: input +* NEN number of LCM objects used in the module. +* KEN addresses of LCM objects: (1) L_KINET; (2) L_MACROLIB; +* (3) L_TRACK; (4) L_SYSTEM; (5) L_MACROLIB. +* CMOD name of the assembly door (BIVAC or TRIVAC). +* NGR number of energy groups. +* NBM number of material mixtures. +* NBFIS number of fissile isotopes. +* NDG number of delayed-neutron groups. +* NLF number of Legendre orders for fluxes. +* ITY type of finite elements and tracking. +* NEL total number of finite elements. +* LL4 number of flux unknowns per energy group. +* NUN total number of unknowns per energy group. +* NUP number of precursor unknowns per delayed group. +* TTF value of theta-parameter for fluxes. +* TTP value of theta-parameter for precursors. +* DT current time increment. +* IMPH management of convergence histogram. +* ICL1 number of free iterations in one cycle of the inverse power +* method +* ICL2 number of accelerated iterations in one cycle +* NADI number of inner adi iterations per outer iteration +* ADJ flag for adjoint space-time kinetics calculation +* MAXOUT maximum number of outer iterations +* EPSOUT convergence criteria for the flux +* MAXINR maximum number of thermal iterations. +* EPSINR thermal iteration epsilon. +* IFL temporal integration scheme for fluxes. +* IPR temporal integration scheme for precursors. +* IEXP exponential transformation flag (=1 to activate). +* INORM type of flux normalization (=0: no normalization; =1: imposed +* factor; =2: maximum flux; =3 initial power). +* IMPX printing parameter (=0 for no print). +* +*Parameter: output +* POWTOT power. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NEN,NGR,NBM,NBFIS,NDG,NLF,ITY,NEL,LL4,NUN,NUP,IMPH,ICL1, + 1 ICL2,NADI,MAXOUT,MAXINR,IFL,IPR,IEXP,INORM,IMPX + TYPE(C_PTR) KEN(NEN) + REAL TTF,TTP,DT,EPSOUT,EPSINR,POWTOT + CHARACTER CMOD*12 + LOGICAL ADJ +*---- +* LOCAL VARIABLES +*---- + PARAMETER(IOS=6) + INTEGER MAT(NEL),IDL(NEL),IDLPC(NEL) + REAL VOL(NEL),PDC(NDG),PMAX(NDG,NBFIS) + LOGICAL LNUD,LCHD + TYPE(C_PTR) IPMAC,IPSYS + REAL, DIMENSION(:), ALLOCATABLE :: DNF,AVG1,AVG2,WORK1,RM + REAL, DIMENSION(:,:), ALLOCATABLE :: EVECT,DNS,PHO,OVR,OMEGA + REAL, DIMENSION(:,:,:), ALLOCATABLE :: PC,CHI,SGF + REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: SGD,CHD,SGO + DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: SRC +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(EVECT(NUN,NGR),PC(NUP,NDG,NBFIS),SGD(NBM,NBFIS,NGR,NDG), + 1 OMEGA(NBM,NGR)) +*---- +* RECOVER INFORMATION +*---- + CALL KDRCPU(TK1) + TA1=TK1 + IF(IMPX.GT.0) WRITE(IOS,1001) + CALL LCMGET(KEN(1),'E-KEFF',EVL) + CALL LCMGET(KEN(1),'LAMBDA-D',PDC) + CALL LCMGET(KEN(1),'E-IDLPC',IDLPC) + CALL LCMLEN(KEN(1),'OMEGA',ILONG,ITYLCM) + IF((IEXP.EQ.0).OR.(ILONG.EQ.0)) THEN + OMEGA(:NBM,:NGR)=0.0 + ELSE + CALL LCMGET(KEN(1),'OMEGA',OMEGA) + ENDIF + CALL LCMGET(KEN(3),'VOLUME',VOL) + CALL LCMGET(KEN(3),'MATCOD',MAT) + CALL LCMGET(KEN(3),'KEYFLX',IDL) +*---- +* RECOVER CROSS SECTIONS (BEGINNING-OF-STEP) +*---- + ALLOCATE(DNF(NDG),DNS(NGR,NDG)) + CALL LCMLEN(KEN(1),'BETA-D',LEN,ITYL) + LNUD=(LEN.EQ.NDG) + IF(LNUD) CALL LCMGET(KEN(1),'BETA-D',DNF) + CALL LCMLEN(KEN(1),'CHI-D',LEN,ITYL) + LCHD=(LEN.EQ.NGR*NDG) + IF(LCHD) CALL LCMGET(KEN(1),'CHI-D',DNS) + ALLOCATE(OVR(NBM,NGR),CHI(NBM,NBFIS,NGR),CHD(NBM,NBFIS,NGR,NDG), + 1 SGF(NBM,NBFIS,NGR),SGO(NBM,NBFIS,NGR,NDG)) + IF(NEN.EQ.4) THEN + IPMAC=KEN(2) + IPSYS=KEN(4) + ELSE IF(NEN.EQ.6) THEN + IPMAC=KEN(5) + IPSYS=KEN(6) + ENDIF + CALL KINXSD(IPMAC,NGR,NBM,NBFIS,NDG,EVL,DT,DNF,DNS,LNUD,LCHD,OVR, + 1 CHI,CHD,SGF,SGO) +*---- +* COMPUTE THE SOURCE TERM +*---- + LL4B=LL4 + IF((ITY.EQ.11).OR.(ITY.EQ.13)) LL4B=LL4*NLF/2 + ALLOCATE(PHO(NUN,NGR),SRC(NUN,NGR)) + CALL LCMGET(KEN(1),'E-PREC',PC) + CALL LCMGET(KEN(1),'E-VECTOR',PHO) + CALL KINSRC(KEN(3),IPSYS,CMOD,IMPX,IFL,IPR,IEXP,NGR,NBM,NBFIS,NDG, + 1 ITY,LL4B,NUN,NUP,PDC,TTF,TTP,DT,ADJ,OVR,CHI,CHD,SGF,SGO,OMEGA, + 2 PHO,PC,SRC) +*---- +* RECOVER CROSS SECTIONS (END-OF-STEP) +*---- + CALL KINXSD(KEN(2),NGR,NBM,NBFIS,NDG,EVL,DT,DNF,DNS,LNUD,LCHD, + 1 OVR,CHI,CHD,SGF,SGD) + DEALLOCATE(DNS,DNF) +*---- +* RECOVER THE BEGINNING-OF-STEP FLUX +*---- + IF(IMPX.GT.0)THEN + CALL KDRCPU(TA2) + WRITE(IOS,1002) TA2-TA1 + WRITE(IOS,1003) + ENDIF + DO 15 IGR=1,NGR + DO 10 IND=1,NUN + EVECT(IND,IGR)=PHO(IND,IGR) + 10 CONTINUE + 15 CONTINUE +*---- +* COMPUTE THE FLUX SOLUTION +*---- + IF(CMOD.EQ.'BIVAC')THEN + IF(ADJ) CALL XABORT('KINDRV: ADJOINT CALCULATION NOT IMPLEMENT' + 1 //'ED WITH BIVAC.') + CALL KINSLB(KEN(3),KEN(4),KEN(1),LL4B,ITY,NUN,NGR,IFL,IPR, + 1 IEXP,NBM,NBFIS,NDG,ICL1,ICL2,IMPX,IMPH,TITR,EPSOUT,MAXINR, + 2 EPSINR,MAXOUT,PDC,TTF,TTP,DT,OVR,CHI,CHD,SGF,SGD,OMEGA,EVECT, + 3 SRC) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL KINSLT(KEN(3),KEN(4),KEN(1),LL4B,ITY,NUN,NGR,IFL,IPR, + 1 IEXP,NBM,NBFIS,NDG,ICL1,ICL2,IMPX,IMPH,TITR,EPSOUT,MAXINR, + 2 EPSINR,NADI,ADJ,MAXOUT,PDC,TTF,TTP,DT,OVR,CHI,CHD,SGF,SGD, + 3 OMEGA,EVECT,SRC) + ENDIF + DEALLOCATE(SRC) +*---- +* COMPUTE THE PRECURSOR SOLUTION +*---- + CALL KINPRC(KEN(3),KEN(4),CMOD,NGR,NBM,NBFIS,NDG,NEL,LL4,NUN,NUP, + 1 MAT,VOL,IDLPC,EVECT,PHO,CHD,CHO,SGD,SGO,PDC,DT,ADJ,TTP,PC,IPR, + 2 IEXP,OMEGA,IMPX) + CALL LCMPUT(KEN(1),'E-PREC',NDG*NUP*NBFIS,2,PC) +*---- +* COMPUTE THE EXPONENTIAL TRANSFORMATION FACTORS +*---- + IF(IEXP.EQ.1) THEN + ALLOCATE(WORK1(LL4),RM(LL4),AVG1(NBM),AVG2(NBM)) + DO 35 IGR=1,NGR + DO 20 IBM=1,NBM + AVG1(IBM)=EXP(OMEGA(IBM,IGR)*DT) + 20 CONTINUE + IF(CMOD.EQ.'BIVAC')THEN + CALL KINBLM(KEN(3),NBM,LL4,AVG1,EVECT(1,IGR),WORK1) + CALL MTLDLS('RM',KEN(3),KEN(4),LL4,1,WORK1) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL KINTLM(KEN(3),NBM,LL4,AVG1,EVECT(1,IGR),WORK1) + CALL LCMLEN(KEN(4),'RM',ILONG,ITYLCM) + CALL LCMGET(KEN(4),'RM',RM) + DO 25 IND=1,ILONG + FACT=RM(IND) + IF(FACT.EQ.0.0) CALL XABORT('KINDRV: SINGULAR RM.') + WORK1(IND)=WORK1(IND)/FACT + 25 CONTINUE + ENDIF + DO 30 IND=1,LL4 + EVECT(IND,IGR)=WORK1(IND) + 30 CONTINUE + 35 CONTINUE + CALL LCMPUT(KEN(1),'E-VECTOR',NGR*NUN,2,EVECT) +* + DO 60 IGR=1,NGR + AVG1(:NBM)=0.0 + AVG2(:NBM)=0.0 + DO 40 IEL=1,NEL + IBM=MAT(IEL) + IF(IBM.GT.0) THEN + AVG1(IBM)=AVG1(IBM)+VOL(IEL)*PHO(IDL(IEL),IGR) + AVG2(IBM)=AVG2(IBM)+VOL(IEL)*EVECT(IDL(IEL),IGR) + ENDIF + 40 CONTINUE + DO 50 IBM=1,NBM + RATIO=MIN(10.0,ABS(AVG2(IBM)/AVG1(IBM))) + OMEGA(IBM,IGR)=LOG(RATIO)/DT + 50 CONTINUE + IF(IMPX.GT.1) THEN + WRITE(IOS,1006) (OMEGA(IBM,IGR),IBM=1,NBM) + ENDIF + 60 CONTINUE + CALL LCMPUT(KEN(1),'OMEGA',NBM*NGR,2,OMEGA) + DEALLOCATE(AVG2,AVG1,RM,WORK1) + ENDIF +*---- +* COMPUTE AVERAGED FLUX VALUES. +*---- + DO 70 IGR=1,NGR + IF(CMOD.EQ.'BIVAC')THEN + CALL FLDBIV(KEN(3),NEL,NUP,EVECT(1,IGR),MAT,VOL,IDL) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL FLDTRI(KEN(3),NEL,NUP,EVECT(1,IGR),MAT,VOL,IDL) + ENDIF + 70 CONTINUE + CALL LCMPUT(KEN(1),'E-VECTOR',NGR*NUN,2,EVECT) +*---- +* FIND THE MAXIMUM FLUX VALUE +*---- + FMAX=0.0 + IDMX=0 + DO 85 IGR=1,NGR + DO 80 IEL=1,NEL + IND=IDL(IEL) + IF(IND.EQ.0) GO TO 80 + IF(ABS(EVECT(IND,IGR)).GT.FMAX) THEN + FMAX=EVECT(IND,IGR) + IDMX=IEL + IGMX=IGR + ENDIF + 80 CONTINUE + 85 CONTINUE + IF(IDMX.EQ.0) CALL XABORT('KINDRV: UNABLE TO SET FMAX.') + IND=IDLPC(IDMX) + IF(IND.EQ.0) CALL XABORT('KINDRV: UNABLE TO SET PMAX.') + DO 95 IFIS=1,NBFIS + DO 90 IDG=1,NDG + PMAX(IDG,IFIS)=PC(IND,IDG,IFIS) + 90 CONTINUE + 95 CONTINUE + IF(IMPX.GT.0) THEN + WRITE(IOS,1004) FMAX,IDMX,IGMX + CALL KDRCPU(TK2) + WRITE(IOS,1005)TK2-TK1 + ENDIF + CALL LCMPUT(KEN(1),'CTRL-FLUX',1,2,FMAX) + CALL LCMPUT(KEN(1),'CTRL-PREC',NDG*NBFIS,2,PMAX) + CALL LCMPUT(KEN(1),'CTRL-IDL',1,1,IDMX) + CALL LCMPUT(KEN(1),'CTRL-IGR',1,1,IGMX) +*---- +* COMPUTE REACTOR POWER +*---- + IF(INORM.EQ.3) THEN + CALL KINPOW(KEN(2),NGR,NBM,NUN,NEL,MAT,VOL,IDL,EVECT,POWTOT) + CALL LCMPUT(KEN(1),'E-POW',1,2,POWTOT) + IF(IMPX.GT.0) WRITE(6,*) 'REACTOR POWER (MW) =',POWTOT + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(PHO,SGO,SGF,CHD,CHI,OVR) + DEALLOCATE(OMEGA,SGD,PC,EVECT) + RETURN +* + 1001 FORMAT(/1X,'=> ASSEMBLY OF THE SYSTEM MATRICES'/) + 1002 FORMAT(/1X,'TOTAL CPU TIME USED FOR THE ASSEMBLING', + 1 1X,'OF ALL SYSTEM MATRICES =',F6.3/) + 1003 FORMAT(/1X,'=> COMPUTING THE KINETICS SOLUTION'/) + 1004 FORMAT(/1X,'CONTROLLING PARAMETERS:',2X,'MAX-VA', + 1 'L',1X,1PE12.5,3X,'IDL #',I5.5,3X,'IGR #',I2.2/) + 1005 FORMAT(/1X,'TOTAL CPU TIME USED FOR KINETICS CALC', + 1 'ULATIONS =',F10.3//1X,'=> SPACE-TIME',1X, + 2 'KINETICS CALCULATION IS DONE.') + 1006 FORMAT(39H KINDRV: MIXTURE-ORDERED OMEGA FACTORS:/(1P,10E14.6)) + END |