MFIX  2016-1
cooling_rate.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2 !
3 ! subroutine name: cooling_rate(s,mi,ni,n,m,T,chi,sigmai,alpha,rhoi,xvec)
4 !
5 ! author: C. Hrenya, Jan 2009
6 !
7 ! Purpose: find theta, Ti, and zeta0 according to GHD polydisperse KT
8 !
9 ! Literature/References:
10 ! C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
11 !
12 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
13 
14  subroutine cooling_rate(s,mi,ni,n,m,Ti,T,chi,sigmai,alpha,rhoi,xvec)
15  Implicit NONE
16 !-----------------------------------------------
17 ! Dummy arguments
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 ! Local variables
33 !-----------------------------------------------
34  integer :: L
35  integer :: ntrial
36  double precision :: tolx, tolf
37 !-----------------------------------------------
38 
39  ntrial = 100
40  tolx = 1d-14
41  tolf = 1d-14
42 
43 ! Initial guess for theta
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 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
56 !
57 ! subroutine name: MNEWT(ntrial, x, mmax, tolx, tolf, mi,ni,n,m,chi,sigmai,alpha,rhoi)
58 ! Purpose: use Newton-Raphson method to solve set of non-linear
59 ! equations to be used for GHD polydisperse KT
60 !
61 ! Literature/Document References:
62 ! Numerical Recipies in Fortran 77, page 374
63 !
64 ! Note: this subroutine has been modified by CMH from Sofiane original
65 ! to accommodate the inputs/outputs of the stand-along program ghd
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 ! Dummy arguments
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 ! Local parameters
89 !-----------------------------------------------
90  INTEGER :: NP
91  parameter(np=15) ! solves up to NP variable/equations;
92 !-----------------------------------------------
93 ! Local variables
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) ! RHS of linear equations.
112  ENDDO
113 
114  CALL ludcmp(fjac, s, np, indx, d, 'MNewt') ! solve system of s linear equations using
115  CALL lubksb(fjac, s, np, indx, p) ! LU decomposition
116 
117  errx = 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 ! for K trials
124 
125  RETURN
126  END SUBROUTINE mnewt
127 
128 
129 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
130 !
131 ! Module name: FUNC_JACOBI_EVAL(X, s, FVEC, FJAC, mi,ni,n,m,chi,sigmai,alpha,rhoi)
132 ! Purpose: computes values of functions and Jacobians of homogenous
133 ! cooling rates to get values of theta_i (or Ti)
134 ! Literature/Document References:
135 ! C. Hrenya handnotes and S. Benyahia derivation, Dec-2008
136 !
137 ! Note: this subroutine has been modified by CMH from Sofiane original
138 ! to accommodate the inputs/outputs of the stand-along program ghd
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 ! Dummy arguments
148 !-----------------------------------------------
149  INTEGER, intent(in) :: s
150  DOUBLE PRECISION :: X(s)
151 ! vector function and matrix jacobian
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 ! Local parameters
161 !-----------------------------------------------
162  double precision :: one
163  parameter(one=1.d0)
164  double precision :: zero
165  parameter(zero=0.d0)
166 !-----------------------------------------------
167 ! Local variables
168 !-----------------------------------------------
169 ! indices
170  INTEGER :: I, J, K
171 ! various quantities
172  DOUBLE PRECISION :: eta(s,s), lambda(s,s)
173  DOUBLE PRECISION :: gamma(s,s), sqrtXiXj(s,s)
174  DOUBLE PRECISION :: oneMGama(s,s)
175 !-----------------------------------------------
176 
177  DO i = 1, s
178  DO j = 1, s
179  eta(i,j) = (one + alpha(i,j))/2.d0
180  lambda(i,j) = chi(i,j)* ni(j)*(mi(j)/(mi(j)+mi(i)))* &
181  ((sigmai(i)+sigmai(j))/2d0)**2 * (eta(i,j)*2d0)
182 ! eta = (1+e)/2
183  gamma(i,j) = (mi(j)/(mi(j)+mi(i)))* (eta(i,j))
184 ! The "2" in mu_ji/2 and 2 * eta will cancel
185  sqrtxixj(i,j) = dsqrt((x(i)+x(j))/(x(i)*x(j)))
186  onemgama(i,j) = (one-gamma(i,j)*((x(i)+x(j))/x(j)))
187  ENDDO
188  ENDDO
189 
190 
191 ! Start computing values of the function FVEC(i)
192  DO i = 1, s
193  IF(i==1) THEN
194  fvec(1) = zero
195  DO k = 1, s
196  fvec(1) = fvec(1) + (rhoi(k)/x(k))
197  ENDDO
198  fvec(1) = fvec(1) - n*m
199  ELSE
200  fvec(i) = zero
201  DO k = 1, s
202  fvec(i) = fvec(i) &
203  + lambda(i,k) * sqrtxixj(i,k) * onemgama(i,k) &
204  - lambda(1,k) * sqrtxixj(1,k) * onemgama(1,k)
205  ENDDO
206  ENDIF
207  ENDDO
208 ! End of computing values of the function FVEC(i)
209 
210 ! Evaluation of functions FVEC(i) is done above, and their Jacobian is computed below.
211 
212 
213 ! Start computing jacobian (J_ij) of the functions FVEC(i).
214  DO i = 1, s
215  DO j = 1, s
216  IF(i==1) THEN
217  fjac(1,j) = -rhoi(j)/x(j)**2
218 
219  ELSEIF(i==j) THEN
220  fjac(j,j) = zero
221 
222  DO k = 1, s
223  IF(k==i) THEN
224  fjac(i,j) = fjac(i,j) &
225  - lambda(1,k)*(one/(x(1)*x(j))-(x(j)+x(1))/ &
226  (x(1)*x(j)**2))*(one-gamma(1,k)*(x(j)+x(1))/ &
227  x(j))/dsqrt((x(j)+x(1))/(x(1)*x(j)))/2.d0 &
228  -lambda(1,k)*dsqrt((x(j)+x(1))/ &
229  (x(1)*x(j)))*(gamma(1,k)*(x(j)+x(1))/x(j)**2 &
230  -gamma(1,k)/ x(j))-dsqrt(2.d0)*(one-2.d0*gamma(i,k))* &
231  lambda(i,k)*x(j)**((-3.d0)/2.d0)/2.d0
232 
233  ELSE
234  fjac(i,j) = fjac(i,j) &
235  + lambda(i,k)*(one-gamma(i,k)*(x(k)+x(j))/x(k))* &
236  (one/(x(j)*x(k))-(x(k)+x(j))/ &
237  (x(j)**2*x(k)))/ dsqrt((x(k)+x(j))/(x(j)*x(k)))/2.d0 &
238  - gamma(i,k)*lambda(i,k)* &
239  dsqrt((x(k)+x(j))/(x(j)*x(k)))/x(k)
240  ENDIF
241  ENDDO
242 
243  ELSEIF(j==1) THEN ! for i /= 1 and i /= j and j == 1
244  k = 1 ! k =1 is a special case for computing Jacobioan
245  fjac(i,j) = -lambda(i,k)*onemgama(i,1)/ &
246  (x(1)**2*sqrtxixj(i,1)*2.0d0) &
247  + lambda(i,k)*sqrtxixj(i,1)*gamma(i,k)*x(i)/x(1)**2 &
248  + dsqrt(2d0)*(one-2d0*gamma(1,k))* &
249  lambda(1,k)*x(1)**((-3.0d0)/2.0d0)/2.0d0
250  DO k = 2, s
251  fjac(i,j) = fjac(i,j) &
252  + gamma(1,k)*lambda(1,k)*sqrtxixj(k,1)/x(k) &
253  + lambda(1,k)*onemgama(1,k)/ &
254  (x(1)**2*sqrtxixj(k,1)*2.0d0)
255  ENDDO
256 
257  ELSE ! for i /= 1 and i /= j and j /= 1
258 ! no do loop for k here, just one term k = j, the rest will be zero.
259  k = j
260  fjac(i,j) = - lambda(i,k)*onemgama(i,k)/ &
261  (x(j)**2*sqrtxixj(j,i)*2d0) &
262  + lambda(i,k)*sqrtxixj(j,i)*gamma(i,k)*x(i)/x(j)**2 &
263  + lambda(1,k)*onemgama(1,k)/(x(j)**2*sqrtxixj(j,1)*2d0) &
264  - lambda(1,k)*sqrtxixj(j,1)*gamma(1,k)*x(1)/x(j)**2
265 
266  ENDIF ! for checks on i
267  ENDDO ! for j
268  ENDDO ! for i
269 ! End of computing jacobian (J_ij).
270 
271  RETURN
272  END SUBROUTINE func_jacobi_eval
273 
274 
275 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
276 !
277 ! subroutine name: ludcmp(a,n,np,indx,d, calledFrom)
278 ! Purpose: Replaces matrix a (n,n) by the LU decomposition of a rowwise
279 ! permutation of itself. Used in combination with lubksb.
280 !
281 ! Literature/Document References:
282 ! Numerical Recipies in Fortran 77, page 38-39
283 !
284 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
285 !
286  subroutine ludcmp(a,n,np,indx,d,calledFrom)
287  USE compar
288  USE exit, only: mfix_exit
289 
290  implicit none
291 !-----------------------------------------------
292 ! Dummy arguments
293 !-----------------------------------------------
294  integer, intent(in) :: n
295  double precision, intent(inout) :: a(n,n)
296  integer :: np
297  integer, intent(out) :: indx(n)
298  double precision, intent(out) :: d
299  CHARACTER(len=*), intent(in) :: calledFrom
300 !-----------------------------------------------
301 ! Local parameters
302 !-----------------------------------------------
303  integer :: nmax
304  double precision :: TINY
305  parameter(nmax=500, tiny=1.0d-20)
306 !-----------------------------------------------
307 ! Local variables
308 !-----------------------------------------------
309  integer :: i, j, k, imax
310  double precision :: vv(nmax)
311  double precision :: aamax, sum, dum
312 !-----------------------------------------------
313 
314  d = 1.0d0
315  do i = 1,n
316  aamax=0.0d0
317  do j = 1,n
318  if (dabs(a(i,j)).gt.aamax) aamax = dabs(a(i,j))
319  enddo
320  if (aamax.eq.0.0d0) then
321  if(mype==pe_io) write(*,*) &
322  'Singular Matrix in ludcmp called from ', calledfrom
323  call mfix_exit(mype)
324  endif
325  vv(i) = 1.d0/aamax
326  enddo
327  do j = 1,n
328  do i = 1,j-1
329  sum = a(i,j)
330  do k = 1,i-1
331  sum = sum-a(i,k)*a(k,j)
332  enddo
333  a(i,j) = sum
334  enddo
335  aamax = 0.0d0
336  do i = j,n
337  sum = a(i,j)
338  do k = 1,j-1
339  sum = sum-a(i,k)*a(k,j)
340  enddo
341  a(i,j) = sum
342  dum = vv(i)*dabs(sum)
343  if (dum.ge.aamax) then
344  imax = i
345  aamax = dum
346  endif
347  enddo
348  if (j.ne.imax) then
349  do k = 1,n
350  dum = a(imax,k)
351  a(imax,k) = a(j,k)
352  a(j,k) = dum
353  enddo
354  d = -d
355  vv(imax) = vv(j)
356  endif
357  indx(j) = imax
358  if (a(j,j).eq.0.0d0) a(j,j) = tiny
359  if (j.ne.n) then
360  dum = 1.0d0/a(j,j)
361  do i = j+1,n
362  a(i,j) = a(i,j)*dum
363  enddo
364  endif
365  enddo
366 
367  return
368  end subroutine ludcmp
369 
370 
371 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
372 !
373 ! subroutine name: lubksb (a,n,np,indx,b)
374 ! Purpose: solves the set of n linear equations A(n,n).X(n) = B(n).
375 !
376 !
377 ! Literature/Document References:
378 ! Numerical Recipies in Fortran 77, page 39.
379 !
380 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
381 
382  subroutine lubksb (a,n,np,indx,b)
383 
384  implicit none
385 !-----------------------------------------------
386 ! Dummy arguments
387 !-----------------------------------------------
388  integer, intent(in) :: n
389  double precision, intent(in) :: a(n,n)
390  integer :: np
391  integer, intent(in) :: indx(n)
392  double precision, intent(inout) :: b(n)
393 !-----------------------------------------------
394 ! Local variables
395 !-----------------------------------------------
396  integer :: i, ii, j, ll
397  double precision :: sum
398 !-----------------------------------------------
399 
400  ii=0
401  do i=1,n
402  ll=indx(i)
403  sum=b(ll)
404  b(ll)=b(i)
405  if (ii.ne.0) then
406  do j=ii,i-1
407  sum=sum-a(i,j)*b(j)
408  enddo
409  elseif (sum.ne.0.d0) then
410  ii=i
411  endif
412  b(i)=sum
413  enddo
414  do i=n,1,-1
415  sum=b(i)
416  if (i.lt.n) then
417  do j=i+1,n
418  sum=sum-a(i,j)*b(j)
419  enddo
420  endif
421  b(i)=sum/a(i,i)
422  enddo
423  return
424  end subroutine lubksb
425 
mnewt
subroutine mnewt(ntrial, x, s, tolx, tolf, mi, ni, n, m, chi, sigmai, alpha, rhoi)
Definition: cooling_rate.f:71
func_jacobi_eval
subroutine func_jacobi_eval(X, s, FVEC, FJAC, mi, ni, n, m, chi, sigmai, alpha, rhoi)
Definition: cooling_rate.f:144
cooling_rate
subroutine cooling_rate(s, mi, ni, n, m, Ti, T, chi, sigmai, alpha, rhoi, xve
Definition: cooling_rate.f:15
compar::pe_io
integer pe_io
Definition: compar_mod.f:30
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
exit
Definition: exit.f:2
ludcmp
subroutine ludcmp(a, n, np, indx, d, calledFrom)
Definition: cooling_rate.f:287
compar::mype
integer mype
Definition: compar_mod.f:24
lubksb
subroutine lubksb(a, n, np, indx, b)
Definition: cooling_rate.f:383