1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
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
|
