File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/cooling_rate.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14 subroutine cooling_rate(s,mi,ni,n,m,Ti,T,chi,sigmai,alpha,rhoi,xvec)
15 Implicit NONE
16
17
18
19 integer, INTENT(IN) :: s
20 double precision :: mi(s)
21 double precision :: ni(s)
22 double precision :: n, m
23 double precision :: Ti(s)
24 double precision :: T
25 double precision :: chi(s,s)
26 double precision :: sigmai(s)
27 double precision :: alpha(s,s)
28 double precision :: rhoi(s)
29 double precision :: xvec(s)
30
31
32
33
34 integer :: L
35 integer :: ntrial
36 double precision :: tolx, tolf
37
38
39 = 100
40 tolx = 1d-14
41 tolf = 1d-14
42
43
44 DO L = 1, s
45 xvec(L) = mi(L) / m
46 ENDDO
47
48 CALL MNEWT(ntrial, xvec, s, tolx, tolf, &
49 mi,ni,n,m,chi,sigmai,alpha,rhoi)
50
51 return
52 end subroutine cooling_rate
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69 subroutine MNEWT(ntrial, x, s, tolx, tolf, &
70 mi,ni,n,m,chi,sigmai,alpha,rhoi)
71 Implicit NONE
72
73
74
75
76 INTEGER, intent(in) :: s
77 INTEGER, intent(in) :: ntrial
78 DOUBLE PRECISION, intent(in) :: tolx, tolf
79 DOUBLE PRECISION :: mi(s)
80 DOUBLE PRECISION :: ni(s)
81 DOUBLE PRECISION :: n, m
82 DOUBLE PRECISION :: chi(s,s)
83 DOUBLE PRECISION :: sigmai(s)
84 DOUBLE PRECISION :: alpha(s,s)
85 DOUBLE PRECISION :: rhoi(s)
86 DOUBLE PRECISION :: X(s)
87
88
89
90 INTEGER :: NP
91 PARAMETER (NP=15)
92
93
94
95 INTEGER :: I, K, INDX(s)
96 DOUBLE PRECISION :: D, ERRF, ERRX
97 DOUBLE PRECISION :: FJAC(s,s), FVEC(s)
98 DOUBLE PRECISION :: P(s)
99
100
101 DO K = 1, NTRIAL
102 CALL FUNC_JACOBI_EVAL(X, s, FVEC, FJAC, &
103 mi,ni,n,m,chi,sigmai,alpha,rhoi)
104 ERRF = 0d0
105 DO I = 1, s
106 ERRF = ERRF + DABS(FVEC(I))
107 ENDDO
108
109 IF(ERRF <= TOLF) RETURN
110 DO I = 1, s
111 P(I) = -FVEC(I)
112 ENDDO
113
114 CALL LUDCMP(fjac, s, NP, indx, d, 'MNewt')
115 CALL LUBKSB(fjac, s, NP, indx, p)
116
117 = 0d0
118 DO I = 1, s
119 ERRX = ERRX + DABS(P(I))
120 X(I) = X(I) + P(I)
121 ENDDO
122 IF(ERRX <= TOLX) RETURN
123 ENDDO
124
125 RETURN
126 END SUBROUTINE MNEWT
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142 SUBROUTINE FUNC_JACOBI_EVAL(X, s, FVEC, FJAC, &
143 mi,ni,n,m,chi,sigmai,alpha,rhoi)
144 Implicit NONE
145
146
147
148
149 INTEGER, intent(in) :: s
150 DOUBLE PRECISION :: X(s)
151
152 DOUBLE PRECISION :: FVEC(s), FJAC(s,s)
153
154 DOUBLE PRECISION :: mi(s)
155 DOUBLE PRECISION :: ni(s)
156 DOUBLE PRECISION :: n, m
157 DOUBLE PRECISION :: chi(s,s), sigmai(s), alpha(s,s), &
158 rhoi(s)
159
160
161
162 double precision :: pi
163 parameter (pi=3.14159265458979323846d0)
164 double precision :: one
165 parameter (one=1.d0)
166 double precision :: zero
167 parameter (zero=0.d0)
168
169
170
171
172 INTEGER :: I, J, K
173
174 DOUBLE PRECISION :: eta(s,s), lambda(s,s)
175 DOUBLE PRECISION :: gamma(s,s), sqrtXiXj(s,s)
176 DOUBLE PRECISION :: oneMGama(s,s)
177
178
179 DO i = 1, s
180 DO j = 1, s
181 ETA(i,j) = (one + alpha(i,j))/2.d0
182 lambda(i,j) = chi(i,j)* ni(j)*(mi(j)/(mi(j)+mi(i)))* &
183 ((sigmai(i)+sigmai(j))/2d0)**2 * (ETA(I,J)*2d0)
184
185 (i,j) = (mi(j)/(mi(j)+mi(i)))* (ETA(I,J))
186
187 (i,j) = dsqrt((x(i)+x(j))/(x(i)*x(j)))
188 oneMGama(i,j) = (one-gamma(i,j)*((x(i)+x(j))/x(j)))
189 ENDDO
190 ENDDO
191
192
193
194 DO i = 1, s
195 IF(i==1) THEN
196 FVEC(1) = ZERO
197 DO k = 1, s
198 FVEC(1) = FVEC(1) + (rhoi(k)/x(k))
199 ENDDO
200 FVEC(1) = FVEC(1) - n*m
201 ELSE
202 FVEC(i) = ZERO
203 DO k = 1, s
204 FVEC(i) = FVEC(i) &
205 + lambda(i,k) * sqrtXiXj(i,k) * oneMGama(i,k) &
206 - lambda(1,k) * sqrtXiXj(1,k) * oneMGama(1,k)
207 ENDDO
208 ENDIF
209 ENDDO
210
211
212
213
214
215
216 DO i = 1, s
217 DO j = 1, s
218 IF(i==1) THEN
219 FJAC(1,j) = -rhoi(j)/x(j)**2
220
221 ELSEIF(i==j) THEN
222 FJAC(j,j) = ZERO
223
224 DO k = 1, s
225 IF(k==i) THEN
226 FJAC(i,j) = FJAC(i,j) &
227 - lambda(1,k)*(one/(x(1)*x(j))-(x(j)+x(1))/ &
228 (x(1)*x(j)**2))*(one-gamma(1,k)*(x(j)+x(1))/ &
229 x(j))/dsqrt((x(j)+x(1))/(x(1)*x(j)))/2.d0 &
230 -lambda(1,k)*dsqrt((x(j)+x(1))/ &
231 (x(1)*x(j)))*(gamma(1,k)*(x(j)+x(1))/x(j)**2 &
232 -gamma(1,k)/ x(j))-dsqrt(2.d0)*(one-2.d0*gamma(i,k))* &
233 lambda(i,k)*x(j)**((-3.d0)/2.d0)/2.d0
234
235 ELSE
236 FJAC(i,j) = FJAC(i,j) &
237 + lambda(i,k)*(one-gamma(i,k)*(x(k)+x(j))/x(k))* &
238 (one/(x(j)*x(k))-(x(k)+x(j))/ &
239 (x(j)**2*x(k)))/ dsqrt((x(k)+x(j))/(x(j)*x(k)))/2.d0 &
240 - gamma(i,k)*lambda(i,k)* &
241 dsqrt((x(k)+x(j))/(x(j)*x(k)))/x(k)
242 ENDIF
243 ENDDO
244
245 ELSEIF(j==1) THEN
246 = 1
247 (i,j) = -lambda(i,k)*oneMGama(i,1)/ &
248 (x(1)**2*sqrtXiXj(i,1)*2.0d0) &
249 + lambda(i,k)*sqrtXiXj(i,1)*gamma(i,k)*x(i)/x(1)**2 &
250 + dsqrt(2d0)*(one-2d0*gamma(1,k))* &
251 lambda(1,k)*x(1)**((-3.0d0)/2.0d0)/2.0d0
252 DO k = 2, s
253 FJAC(i,j) = FJAC(i,j) &
254 + gamma(1,k)*lambda(1,k)*sqrtXiXj(k,1)/x(k) &
255 + lambda(1,k)*oneMGama(1,k)/ &
256 (x(1)**2*sqrtXiXj(k,1)*2.0d0)
257 ENDDO
258
259 ELSE
260
261 = j
262 FJAC(i,j) = - lambda(i,k)*oneMGama(i,k)/ &
263 (x(j)**2*sqrtXiXj(j,i)*2d0) &
264 + lambda(i,k)*sqrtXiXj(j,i)*gamma(i,k)*x(i)/x(j)**2 &
265 + lambda(1,k)*oneMGama(1,k)/(x(j)**2*sqrtXiXj(j,1)*2d0) &
266 - lambda(1,k)*sqrtXiXj(j,1)*gamma(1,k)*x(1)/x(j)**2
267
268 ENDIF
269 ENDDO
270 ENDDO
271
272
273 RETURN
274 END SUBROUTINE FUNC_JACOBI_EVAL
275
276
277
278
279
280
281
282
283
284
285
286
287
288 subroutine ludcmp(a,n,np,indx,d,calledFrom)
289 USE compar
290
291 implicit none
292
293
294
295 integer, intent(in) :: n
296 double precision, intent(inout) :: a(n,n)
297 integer :: np
298 integer, intent(out) :: indx(n)
299 double precision, intent(out) :: d
300 CHARACTER(len=*), intent(in) :: calledFrom
301
302
303
304 integer :: nmax
305 double precision :: TINY
306 parameter (NMAX=500, TINY=1.0D-20)
307
308
309
310 integer :: i, j, k, imax
311 double precision :: vv(NMAX)
312 double precision :: aamax, sum, dum
313
314
315 = 1.0d0
316 do i = 1,n
317 aamax=0.0d0
318 do j = 1,n
319 if (dabs(a(i,j)).gt.aamax) aamax = dabs(a(i,j))
320 enddo
321 if (aamax.eq.0.0d0) then
322 if(myPE==PE_IO) write(*,*) &
323 'Singular Matrix in ludcmp called from ', calledFrom
324 call mfix_exit(myPE)
325 endif
326 vv(i) = 1.d0/aamax
327 enddo
328 do j = 1,n
329 do i = 1,j-1
330 sum = a(i,j)
331 do k = 1,i-1
332 sum = sum-a(i,k)*a(k,j)
333 enddo
334 a(i,j) = sum
335 enddo
336 aamax = 0.0d0
337 do i = j,n
338 sum = a(i,j)
339 do k = 1,j-1
340 sum = sum-a(i,k)*a(k,j)
341 enddo
342 a(i,j) = sum
343 dum = vv(i)*dabs(sum)
344 if (dum.ge.aamax) then
345 imax = i
346 aamax = dum
347 endif
348 enddo
349 if (j.ne.imax) then
350 do k = 1,n
351 dum = a(imax,k)
352 a(imax,k) = a(j,k)
353 a(j,k) = dum
354 enddo
355 d = -d
356 vv(imax) = vv(j)
357 endif
358 indx(j) = imax
359 if (a(j,j).eq.0.0d0) a(j,j) = TINY
360 if (j.ne.n) then
361 dum = 1.0d0/a(j,j)
362 do i = j+1,n
363 a(i,j) = a(i,j)*dum
364 enddo
365 endif
366 enddo
367
368 return
369 end subroutine ludcmp
370
371
372
373
374
375
376
377
378
379
380
381
382
383 subroutine lubksb (a,n,np,indx,b)
384
385 implicit none
386
387
388
389 integer, intent(in) :: n
390 double precision, intent(in) :: a(n,n)
391 integer :: np
392 integer, intent(in) :: indx(n)
393 double precision, intent(inout) :: b(n)
394
395
396
397 integer :: i, ii, j, ll
398 double precision :: sum
399
400
401 =0
402 do i=1,n
403 ll=indx(i)
404 sum=b(ll)
405 b(ll)=b(i)
406 if (ii.ne.0) then
407 do j=ii,i-1
408 sum=sum-a(i,j)*b(j)
409 enddo
410 elseif (sum.ne.0.d0) then
411 ii=i
412 endif
413 b(i)=sum
414 enddo
415 do i=n,1,-1
416 sum=b(i)
417 if (i.lt.n) then
418 do j=i+1,n
419 sum=sum-a(i,j)*b(j)
420 enddo
421 endif
422 b(i)=sum/a(i,i)
423 enddo
424 return
425 end subroutine lubksb
426
427