Initial commit from Polytechnique Montreal

1 files changed, 454 insertions, 0 deletions

diff --git a/Trivac/src/KINSLB.f b/Trivac/src/KINSLB.f
new file mode 100755
index 0000000..6bb1b20
--- /dev/null
+++ b/Trivac/src/KINSLB.f
@@ -0,0 +1,454 @@
+*DECK KINSLB
+      SUBROUTINE KINSLB (IPTRK,IPSYS,IPKIN,LL4,ITY,NUN,NGR,IFL,IPR,IEXP,
+     1 NBM,NBFIS,NDG,ICL1,ICL2,IMPX,IMPH,TITR,EPS2,MAXINR,EPSINR,MAXX0,
+     2 PDC,TTF,TTP,DT,OVR,CHI,CHD,SGF,SGD,OMEGA,EVECT,SRC)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of the kinetics multigroup linear systems for the transient
+* neutron fluxes in Bivac. Use the inverse power method with a
+* two-parameter SVAT acceleration technique.
+*
+*Copyright:
+* Copyright (C) 2010 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
+* IPTRK   L_TRACK pointer to the tracking information.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IPKIN   L_KINET pointer to the KINET object.
+* LL4     order of the system matrices.
+* ITY     type of solution (1: classical Bivac/diffusion;
+*         11: Bivac/SPN).
+* NUN     number of unknowns in each energy group.
+* NGR     number of energy groups.
+* IFL     integration scheme for fluxes: =1 implicit;
+*         =2 Crank-Nicholson; =3 theta.
+* IPR     integration scheme for precursors: =1 implicit;
+*         =2 Crank-Nicholson; =3 theta; =4 exponential.
+* IEXP    exponential transformation flag (=1 to activate).
+* NBM     number of material mixtures.
+* NBFIS   number of fissile isotopes.
+* NDG     number of delayed-neutron groups.
+* ICL1    number of free iterations in one cycle of the inverse power
+*         method
+* ICL2    number of accelerated iterations in one cycle
+* IMPX    print parameter. =0: no print ; =1: minimum printing ;
+*         =2: iteration history is printed. =3: solution is printed
+* IMPH    =0: no action is taken
+*         =1: the flux is compared to a reference flux stored on lcm
+*         =2: the convergence histogram is printed
+*         =3: the convergence histogram is printed with axis and
+*            titles. The plotting file is completed
+*         =4: the convergence histogram is printed with axis, acce-
+*            leration factors and titles. The plotting file is
+*            completed
+* TITR    character*72 title
+* EPS2    convergence criteria for the flux
+* MAXINR  maximum number of thermal iterations.
+* EPSINR  thermal iteration epsilon.
+* MAXX0   maximum number of outer iterations
+* PDC     precursor decay constants.
+* TTF     value of theta-parameter for fluxes.
+* TTP     value of theta-parameter for precursors.
+* DT      current time increment.
+* OVR     reciprocal neutron velocities/DT.
+* CHI     steady-state fission spectrum.
+* CHD     delayed fission spectrum
+* SGF     nu*fission macroscopic x-sections/keff.
+* SGD     delayed nu*fission macroscopic x-sections/keff.
+* OMEGA   exponential transformation parameter.
+* SRC     fixed source
+*
+*Parameters: output
+* EVECT    converged solution
+*
+*References:
+* A. H\'ebert, 'Preconditioning the power method for reactor
+* calculations', Nucl. Sci. Eng., 94, 1 (1986).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      CHARACTER TITR*72
+      TYPE(C_PTR) IPTRK,IPSYS,IPKIN
+      INTEGER LL4,ITY,NUN,NGR,IFL,IPR,IEXP,NBM,NBFIS,NDG,ICL1,ICL2,IMPX,
+     1 IMPH,MAXINR,MAXX0
+      REAL EPS2,EPSINR,PDC(NDG),TTF,TTP,DT,OVR(NBM,NGR),
+     1 CHI(NBM,NBFIS,NGR),CHD(NBM,NBFIS,NGR,NDG),SGF(NBM,NBFIS,NGR),
+     2 SGD(NBM,NBFIS,NGR,NDG),OMEGA(NBM,NGR),EVECT(NUN,NGR)
+      DOUBLE PRECISION SRC(NUN,NGR)
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER*12 TEXT12
+      LOGICAL LOGTES,LMPH
+      DOUBLE PRECISION D2F(2,3),ALP,BET,DTF,DTP,DARG,DK
+      REAL ERR(250),ALPH(250),BETA(250),TKT,TKB
+      INTEGER  ITITR(18)
+      REAL, DIMENSION(:,:), ALLOCATABLE :: GRAD1,GRAD2
+      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: GAR1,GAR2,GAR3
+      REAL, DIMENSION(:), ALLOCATABLE :: WORK1,WORK2,WORK3,WORK4
+      DATA EPS1,MMAXX/1.0E-4,250/
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(GRAD1(NUN,NGR),GRAD2(NUN,NGR),GAR1(NUN,NGR),
+     1 GAR2(NUN,NGR),GAR3(NUN,NGR),WORK1(LL4),WORK2(LL4),WORK3(NBM))
+*
+      CALL MTOPEN(IMPX,IPTRK,LL4)
+      IF(LL4.GT.NUN) CALL XABORT('KINSLB: INVALID NUMBER OF UNKNOWNS.')
+*----
+*  INVERSE POWER METHOD.
+*----
+      DTF=9999.0D0
+      DTP=9999.0D0
+      TEST=0.0
+      IF(IFL.EQ.1)THEN
+        DTF=1.0D0
+      ELSEIF(IFL.EQ.2)THEN
+        DTF=0.5D0
+      ELSEIF(IFL.EQ.3)THEN
+        DTF=DBLE(TTF)
+      ENDIF
+      IF(IPR.EQ.2)THEN
+        DTP=0.5D0
+      ELSEIF(IPR.EQ.3)THEN
+        DTP=DBLE(TTP)
+      ENDIF
+      DCRIT=MINVAL(DT*PDC(:))
+*
+      ISTART=1
+      IF(IMPX.GE.1) WRITE (6,600)
+      IF(IMPX.GE.2) WRITE (6,610)
+      M=0
+   10 M=M+1
+*
+      DO 84 IGR=1,NGR
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,IGR),WORK1)
+      DO 15 IND=1,LL4
+      GAR1(IND,IGR)=DTF*WORK1(IND)
+   15 CONTINUE
+      IF(IEXP.EQ.0) THEN
+        DO 16 IBM=1,NBM
+        WORK3(IBM)=OVR(IBM,IGR)
+   16   CONTINUE
+      ELSE
+        DO 17 IBM=1,NBM
+        WORK3(IBM)=OVR(IBM,IGR)*(1.0+OMEGA(IBM,IGR)*DT)
+   17   CONTINUE
+      ENDIF
+      CALL KINBLM(IPTRK,NBM,LL4,WORK3,EVECT(1,IGR),WORK1)
+      DO 20 IND=1,LL4
+      GAR1(IND,IGR)=GAR1(IND,IGR)+WORK1(IND)
+   20 CONTINUE
+      DO 83 JGR=1,NGR
+      IF(JGR.EQ.IGR) GO TO 40
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 40
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK1)
+      DO 30 IND=1,LL4
+      GAR1(IND,IGR)=GAR1(IND,IGR)-DTF*WORK1(IND)
+   30 CONTINUE
+   40 DO 82 IFIS=1,NBFIS
+      DO 50 IBM=1,NBM
+      WORK3(IBM)=CHI(IBM,IFIS,IGR)*SGF(IBM,IFIS,JGR)
+   50 CONTINUE
+      CALL KINBLM(IPTRK,NBM,LL4,WORK3,EVECT(1,JGR),WORK1)
+      DO 60 IND=1,LL4
+      GAR1(IND,IGR)=GAR1(IND,IGR)-DTF*WORK1(IND)
+   60 CONTINUE
+      DO 81 IDG=1,NDG
+      DARG=PDC(IDG)*DT
+      IF(IPR.EQ.1)THEN
+        DK=1.0D0/(1.0D0+DARG)
+      ELSEIF(IPR.EQ.4)THEN
+        DK=(1.0D0-DEXP(-DARG))/DARG
+      ELSE
+        DK=1.0D0/(1.0D0+DTP*DARG)
+      ENDIF
+      DO 70 IBM=1,NBM
+      WORK3(IBM)=CHD(IBM,IFIS,IGR,IDG)*SGD(IBM,IFIS,JGR,IDG)
+   70 CONTINUE
+      CALL KINBLM(IPTRK,NBM,LL4,WORK3,EVECT(1,JGR),WORK1)
+      DO 80 IND=1,LL4
+      GAR1(IND,IGR)=GAR1(IND,IGR)+DTF*DK*WORK1(IND)
+   80 CONTINUE
+   81 CONTINUE
+   82 CONTINUE
+   83 CONTINUE
+   84 CONTINUE
+*----
+*  DIRECTION EVALUATION.
+*----
+      DO 120 IGR=1,NGR
+      DO 90 IND=1,LL4
+      GRAD1(IND,IGR)=REAL(SRC(IND,IGR)-GAR1(IND,IGR))
+   90 CONTINUE
+      DO 110 JGR=1,IGR-1
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 110
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK1)
+      DO 100 IND=1,LL4
+      GRAD1(IND,IGR)=GRAD1(IND,IGR)+REAL(DTF)*WORK1(IND)
+  100 CONTINUE
+  110 CONTINUE
+      CALL KDRCPU(TK2)
+      TKB=TKB+(TK2-TK1)
+*
+      CALL KDRCPU(TK1)
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL MTLDLS(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,IGR))
+      CALL KDRCPU(TK2)
+      DO 115 IND=1,LL4
+      GRAD1(IND,IGR)=GRAD1(IND,IGR)/REAL(DTF)
+  115 CONTINUE
+      TKT=TKT+(TK2-TK1)
+  120 CONTINUE
+*----
+*  PERFORM THERMAL (UP-SCATTERING) ITERATIONS
+*----
+      KTER=0
+      NADI=5 ! used with SPN approximations
+      IF(MAXINR.GT.1) THEN
+         CALL FLDBHR(IPTRK,IPSYS,.FALSE.,LL4,ITY,NUN,NGR,ICL1,ICL2,IMPX,
+     1   NADI,MAXINR,EPSINR,KTER,TKT,TKB,GRAD1)
+      ENDIF
+*----
+*  EVALUATION OF THE DISPLACEMENT AND OF THE TWO ACCELERATION PARAMETERS
+*  ALP AND BET.
+*----
+      DO 204 IGR=1,NGR
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,IGR),WORK1)
+      DO 130 IND=1,LL4
+      GAR2(IND,IGR)=DTF*WORK1(IND)
+  130 CONTINUE
+      IF(IEXP.EQ.0) THEN
+        DO 135 IBM=1,NBM
+        WORK3(IBM)=OVR(IBM,IGR)
+  135   CONTINUE
+      ELSE
+        DO 136 IBM=1,NBM
+        WORK3(IBM)=OVR(IBM,IGR)*(1.0+OMEGA(IBM,IGR)*DT)
+  136   CONTINUE
+      ENDIF
+      CALL KINBLM(IPTRK,NBM,LL4,WORK3,GRAD1(1,IGR),WORK1)
+      DO 140 IND=1,LL4
+      GAR2(IND,IGR)=GAR2(IND,IGR)+WORK1(IND)
+  140 CONTINUE
+      DO 203 JGR=1,NGR
+      IF(JGR.EQ.IGR) GO TO 160
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 160
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK1)
+      DO 150 IND=1,LL4
+      GAR2(IND,IGR)=GAR2(IND,IGR)-DTF*WORK1(IND)
+  150 CONTINUE
+  160 DO 202 IFIS=1,NBFIS
+      DO 170 IBM=1,NBM
+      WORK3(IBM)=CHI(IBM,IFIS,IGR)*SGF(IBM,IFIS,JGR)
+  170 CONTINUE
+      CALL KINBLM(IPTRK,NBM,LL4,WORK3,GRAD1(1,JGR),WORK1)
+      DO 180 IND=1,LL4
+      GAR2(IND,IGR)=GAR2(IND,IGR)-DTF*WORK1(IND)
+  180 CONTINUE
+      DO 201 IDG=1,NDG
+      DARG=PDC(IDG)*DT
+      IF(IPR.EQ.1)THEN
+        DK=1.0D0/(1.0D0+DARG)
+      ELSEIF(IPR.EQ.4)THEN
+        DK=(1.0D0-DEXP(-DARG))/DARG
+      ELSE
+        DK=1.0D0/(1.0D0+DTP*DARG)
+      ENDIF
+      DO 190 IBM=1,NBM
+      WORK3(IBM)=CHD(IBM,IFIS,IGR,IDG)*SGD(IBM,IFIS,JGR,IDG)
+  190 CONTINUE
+      CALL KINBLM(IPTRK,NBM,LL4,WORK3,GRAD1(1,JGR),WORK1)
+      DO 200 IND=1,LL4
+      GAR2(IND,IGR)=GAR2(IND,IGR)+DTF*DK*WORK1(IND)
+  200 CONTINUE
+  201 CONTINUE
+  202 CONTINUE
+  203 CONTINUE
+  204 CONTINUE
+*
+  270 ALP=1.0D0
+      BET=0.0D0
+      D2F(:2,:3)=0.0D0
+      IF(1+MOD(M-ISTART,ICL1+ICL2).GT.ICL1) THEN
+         IF(DCRIT.GT.1.0E-6) THEN
+*           TWO-PARAMETER ACCELERATION. SOLUTION OF A LINEAR SYSTEM.
+            DO 285 IGR=1,NGR
+            DO 280 I=1,LL4
+            D2F(1,1)=D2F(1,1)+GAR2(I,IGR)**2
+            D2F(1,2)=D2F(1,2)+GAR2(I,IGR)*GAR3(I,IGR)
+            D2F(2,2)=D2F(2,2)+GAR3(I,IGR)**2
+            D2F(1,3)=D2F(1,3)-(GAR1(I,IGR)-SRC(I,IGR))*GAR2(I,IGR)
+            D2F(2,3)=D2F(2,3)-(GAR1(I,IGR)-SRC(I,IGR))*GAR3(I,IGR)
+  280       CONTINUE
+  285       CONTINUE
+            D2F(2,1)=D2F(1,2)
+            CALL ALSBD(2,1,D2F,IER,2)
+            IF(IER.NE.0) THEN
+               DCRIT=1.0E-6
+               GO TO 270
+            ENDIF
+            ALP=D2F(1,3)
+            BET=D2F(2,3)/ALP
+            IF((ALP.LT.1.0D0).AND.(ALP.GT.0.0D0)) THEN
+               ALP=1.0D0
+               BET=0.0D0
+            ELSE IF(ALP.LE.0.0D0) THEN
+               ISTART=M+1
+               ALP=1.0D0
+               BET=0.0D0
+            ENDIF
+         ELSE
+*           ONE-PARAMETER ACCELERATION.
+            DO 295 IGR=1,NGR
+            DO 290 I=1,LL4
+            D2F(1,1)=D2F(1,1)+GAR2(I,IGR)**2
+            D2F(1,3)=D2F(1,3)-(GAR1(I,IGR)-SRC(I,IGR))*GAR2(I,IGR)
+  290       CONTINUE
+  295       CONTINUE
+            IF(D2F(1,1).NE.0.0D0) THEN
+               ALP=D2F(1,3)/D2F(1,1)
+            ELSE
+               ISTART=M+1
+            ENDIF
+         ENDIF
+         DO 305 IGR=1,NGR
+         DO 300 I=1,LL4
+         GRAD1(I,IGR)=REAL(ALP)*(GRAD1(I,IGR)+REAL(BET)*GRAD2(I,IGR))
+         GAR2(I,IGR)=ALP*(GAR2(I,IGR)+BET*GAR3(I,IGR))
+  300    CONTINUE
+  305    CONTINUE
+      ENDIF
+*
+      LOGTES=(M.LT.ICL1).OR.(MOD(M-ISTART,ICL1+ICL2).EQ.ICL1-1)
+      IF(LOGTES) THEN
+        ALLOCATE(WORK4(LL4))
+        DELT=0.0
+        DO 350 IGR=1,NGR
+        WORK1(:LL4)=0.0
+        WORK2(:LL4)=0.0
+        DO 320 JGR=1,NGR
+        WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+        CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+        IF(ILONG.EQ.0) GO TO 320
+        CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK4)
+        DO 310 I=1,LL4
+        WORK1(I)=WORK1(I)+WORK4(I)
+  310   CONTINUE
+        CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK4)
+        DO 315 I=1,LL4
+        WORK2(I)=WORK2(I)+WORK4(I)
+  315   CONTINUE
+  320   CONTINUE
+        DELN=0.0
+        DELD=0.0
+        DO 340 I=1,LL4
+        EVECT(I,IGR)=EVECT(I,IGR)+GRAD1(I,IGR)
+        GAR1(I,IGR)=GAR1(I,IGR)+GAR2(I,IGR)
+        GRAD2(I,IGR)=GRAD1(I,IGR)
+        GAR3(I,IGR)=GAR2(I,IGR)
+        DELN=MAX(DELN,ABS(WORK2(I)))
+        DELD=MAX(DELD,ABS(WORK1(I)))
+  340   CONTINUE
+        IF(DELD.NE.0.0) DELT=MAX(DELT,DELN/DELD)
+  350   CONTINUE
+        DEALLOCATE(WORK4)
+        IF(IMPX.GE.2) WRITE (6,620) M,ALP,BET,DELT
+*       COMPUTE THE CONVERGENCE HISTOGRAM.
+        IF((IMPH.GE.1).AND.(M.LE.250)) THEN
+          LMPH=IMPH.GE.1
+          CALL FLDXCO(IPKIN,LL4,NUN,EVECT(1,NGR),LMPH,ERR(M))
+          ALPH(M)=REAL(ALP)
+          BETA(M)=REAL(BET)
+        ENDIF
+        IF(DELT.LT.EPS2) GO TO 370
+      ELSE
+        DO 365 IGR=1,NGR
+        DO 360 I=1,LL4
+        EVECT(I,IGR)=EVECT(I,IGR)+GRAD1(I,IGR)
+        GAR1(I,IGR)=GAR1(I,IGR)+GAR2(I,IGR)
+        GRAD2(I,IGR)=GRAD1(I,IGR)
+        GAR3(I,IGR)=GAR2(I,IGR)
+  360   CONTINUE
+  365   CONTINUE
+        IF(IMPX.GE.2) WRITE (6,620) M,ALP,BET
+*       COMPUTE THE CONVERGENCE HISTOGRAM.
+        IF((IMPH.GE.1).AND.(M.LE.250)) THEN
+          LMPH=IMPH.GE.1
+          CALL FLDXCO(IPKIN,LL4,NUN,EVECT(1,NGR),LMPH,ERR(M))
+          ALPH(M)=REAL(ALP)
+          BETA(M)=REAL(BET)
+        ENDIF
+      ENDIF
+      IF(M.EQ.1) TEST=DELT
+      IF((M.GT.30).AND.(DELT.GT.TEST)) CALL XABORT('KINSLB: CONVERGENC'
+     1 //'E FAILURE.')
+      IF(M.GE.MIN(MAXX0,MMAXX)) THEN
+         WRITE (6,710)
+         GO TO 370
+      ENDIF
+      GO TO 10
+*----
+*  SOLUTION EDITION.
+*----
+  370 IF(IMPX.EQ.1) WRITE (6,640) M
+      IF(IMPX.GE.3) THEN
+         DO 380 IGR=1,NGR
+         WRITE (6,690) IGR,(EVECT(I,IGR),I=1,LL4)
+  380    CONTINUE
+      ENDIF
+      IF(IMPH.GE.2) THEN
+         IGRAPH=0
+  390    IGRAPH=IGRAPH+1
+         WRITE (TEXT12,'(5HHISTO,I3)') IGRAPH
+         CALL LCMLEN (IPKIN,TEXT12,ILENG,ITYLCM)
+         IF(ILENG.EQ.0) THEN
+            MDIM=MIN(250,M)
+            READ (TITR,'(18A4)') ITITR
+            CALL LCMSIX (IPKIN,TEXT12,1)
+            CALL LCMPUT (IPKIN,'HTITLE',18,3,ITITR)
+            CALL LCMPUT (IPKIN,'ALPHA',MDIM,2,ALPH)
+            CALL LCMPUT (IPKIN,'BETA',MDIM,2,BETA)
+            CALL LCMPUT (IPKIN,'ERROR',MDIM,2,ERR)
+            CALL LCMPUT (IPKIN,'IMPH',1,1,IMPH)
+            CALL LCMSIX (IPKIN,' ',2)
+         ELSE
+            GO TO 390
+         ENDIF
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(GRAD1,GRAD2,GAR1,GAR2,GAR3,WORK1,WORK2,WORK3)
+      RETURN
+*
+  600 FORMAT(1H1/50H KINSLB: ITERATIVE PROCEDURE BASED ON INVERSE POWE,
+     1 8HR METHOD/9X,30HSPACE-TIME KINETICS EQUATIONS.)
+  610 FORMAT(/11X,5HALPHA,3X,4HBETA,6X,8HACCURACY,12(1H.))
+  620 FORMAT(1X,I3,4X,2F8.3,1PE13.2)
+  640 FORMAT(/23H KINSLB: CONVERGENCE IN,I4,12H ITERATIONS.)
+  690 FORMAT(//52H KINSLB: SPACE-TIME KINETICS SOLUTION CORRESPONDING ,
+     1 12HTO THE GROUP,I4//(5X,1P,8E14.5))
+  710 FORMAT(/53H KINSLB: ***WARNING*** THE MAXIMUM NUMBER OF OUTER IT,
+     1 20HERATIONS IS REACHED.)
+      END