Initial commit from Polytechnique Montreal

1 files changed, 182 insertions, 0 deletions

diff --git a/Donjon/src/PKIDRV.f b/Donjon/src/PKIDRV.f
new file mode 100644
index 0000000..fc3dcad
--- /dev/null
+++ b/Donjon/src/PKIDRV.f
@@ -0,0 +1,182 @@
+*DECK PKIDRV
+      SUBROUTINE PKIDRV(IPMAP,NALPHA,NGROUP,LAMBDA,EPSILON,BETAI,
+     1 LAMBDAI,DT,PARAMI,PARAMB,T,Y)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of the point kinetic equations using the Runge-Kutta method.
+*
+*Copyright:
+* Copyright (C) 2017 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
+* IPMAP   pointer to the point kinetic directory
+* NALPHA  number of feedback parameters
+* NGROUP  number of delayed precursor groups
+* LAMBDA  prompt neutron generation time
+* EPSILON Runge-Kutta epsilon
+* BETAI   delayed neutron fraction vector
+* LAMBDAI delayed neutron time constant vector
+* DT      stage duration (double precision value)
+* PARAMI  initial values of the global parameters corresponding to
+*         RHO=0
+* PARAMB  values of global parameters at beginning of stage
+* T       time at beggining of stage (double precision value)
+* Y       solution of the point kinetic equations at beginning of stage
+*
+*Parameters: ouput
+* PARAMB  values of global parameters at end of stage
+* T       time at end of stage
+* Y       solution of the point kinetic equations at end of stage
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  Subroutine arguments
+*----
+      TYPE(C_PTR) IPMAP
+      INTEGER NALPHA,NGROUP
+      REAL LAMBDA,EPSILON,BETAI(NGROUP),LAMBDAI(NGROUP),PARAMI(NALPHA),
+     1 PARAMB(NALPHA)
+      DOUBLE PRECISION DT,T,Y(NGROUP+1)
+*----
+*  Local variables
+*----
+      PARAMETER(NRKMIN=100,NRKMAX=100000)
+      DOUBLE PRECISION DH,DPP,P0,T1,BETA,MAXI,RHO(3)
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: YSAV,YSUM,Y1,Y2,
+     1 Y3,Y4
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: A
+*----
+*  Scratch storage allocation
+*----
+      ALLOCATE(YSAV(NGROUP+1),YSUM(NGROUP+1),Y1(NGROUP+1),Y2(NGROUP+1),
+     1 Y3(NGROUP+1),Y4(NGROUP+1),A(NGROUP+1,NGROUP+1))
+*----
+*   Runge-Kutta and calcul parameters initialisation
+*----
+      DPP=1.D0
+      NRK=NRKMIN
+      P0=-1.D0
+*----
+*  Set the Runge-Kutta evolution matrix
+*----
+      BETA=0.D0
+      DO I=1,NGROUP
+         BETA=BETA+BETAI(I)
+      ENDDO
+      A(:NGROUP+1,:NGROUP+1)=0.0D0
+      DO I=2,NGROUP+1
+         A(I,1)=BETAI(I-1)/LAMBDA
+      ENDDO
+      DO I=2,NGROUP+1
+         A(1,I)=LAMBDAI(I-1)
+      ENDDO
+      DO I=2,NGROUP+1
+         A(I,I)=-LAMBDAI(I-1)
+      ENDDO
+*----
+*  Runge-Kutta convergence loop
+*----
+      RHO(:3)=0.0D0
+      DO WHILE ((DPP.GE.EPSILON).AND.(NRK.LE.NRKMAX))
+*       time and time-step initialisation
+        DH=DT/REAL(NRK)
+        T1=T
+*
+*       save of the working vector
+        DO I=1,NGROUP+1
+          YSAV(I)=Y(I)
+        ENDDO       
+*
+*       Runge-Kutta iteration loop
+        DO I=1,NRK
+*         total reactivity calculation with feedback
+          IF(NALPHA.GT.0) CALL PKIRHO(IPMAP,NALPHA,T,DH,PARAMI,PARAMB,
+     1    RHO)
+*
+*         Runge-Kutta procedure
+          A(1,1)=(RHO(1)-BETA)/LAMBDA
+          DO J=1,NGROUP+1
+             Y1(J)=0.0D0
+             DO K=1,NGROUP+1
+                Y1(J)=Y1(J)+A(J,K)*Y(K)
+             ENDDO
+          ENDDO
+          DO J=1,NGROUP+1
+             YSUM(J)=Y(J)+(DH/6.D0)*Y1(J)
+             Y1(J)=Y(J)+DH/2.D0*Y1(J)
+          ENDDO
+          A(1,1)=(RHO(2)-BETA)/LAMBDA
+          DO J=1,NGROUP+1
+             Y2(J)=0.0D0
+             DO K=1,NGROUP+1
+                Y2(J)=Y2(J)+A(J,K)*Y1(K)
+             ENDDO
+          ENDDO
+          DO J=1,NGROUP+1
+             YSUM(J)=YSUM(J)+(DH/3.D0)*Y2(J)
+             Y2(J)=Y(J)+DH/2.D0*Y2(J)
+          ENDDO
+          A(1,1)=(RHO(2)-BETA)/LAMBDA
+          DO J=1,NGROUP+1
+             Y3(J)=0.0D0
+             DO K=1,NGROUP+1
+                Y3(J)=Y3(J)+A(J,K)*Y2(K)
+             ENDDO
+          ENDDO
+          DO J=1,NGROUP+1
+             YSUM(J)=YSUM(J)+(DH/3.D0)*Y3(J)
+             Y3(J)=Y(J)+DH*Y3(J)
+          ENDDO
+          A(1,1)=(RHO(3)-BETA)/LAMBDA
+          DO J=1,NGROUP+1
+             Y4(J)=0.0D0
+             DO K=1,NGROUP+1
+                Y4(J)=Y4(J)+A(J,K)*Y3(K)
+             ENDDO
+          ENDDO
+          DO J=1,NGROUP+1
+             YSUM(J)=YSUM(J)+(DH/6.D0)*Y4(J)
+             Y(J)=YSUM(J)
+          ENDDO
+          T=T+DH
+*
+*         convergence test initialisation
+          MAXI=0.D0
+          DO J=1,NGROUP+1
+             MAXI=MAX(ABS(Y(J)),MAXI)
+          ENDDO
+          IF(MAXI.GT.1.0D30) GOTO 100
+        ENDDO
+* 
+*       convergence test
+ 100    IF(P0.NE.-1.D0) DPP=ABS(Y(1)-P0)/ABS(P0)
+        P0=Y(1)
+*
+*       reinitialisation of the number of Runge-Kutta time-steps
+        NRK=2*NRK
+*
+*       reinitialisation of the working vector if not converged
+        IF((DPP.GE.EPSILON).AND.(NRK.LE.NRKMAX)) THEN
+          DO I=1,NGROUP+1
+             Y(I)=YSAV(I)
+          ENDDO
+          T=T1
+        ENDIF
+      ENDDO
+*----
+*  Scratch storage deallocation
+*----
+      DEALLOCATE(A,Y4,Y3,Y2,Y1,YSUM,YSAV)
+      RETURN
+      END