File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/cooling_rate.f

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 :: 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     ! Local variables
170     !-----------------------------------------------
171     ! indices
172           INTEGER :: I, J, K
173     ! various quantities
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     ! eta = (1+e)/2
185                  gamma(i,j) = (mi(j)/(mi(j)+mi(i)))* (ETA(I,J))
186     ! The "2" in mu_ji/2 and 2 * eta will cancel
187                  sqrtXiXj(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     ! Start computing values of the function FVEC(i)
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     ! End of computing values of the function FVEC(i)
211     
212     ! Evaluation of functions FVEC(i) is done above, and their Jacobian is computed below.
213     
214     
215     ! Start computing jacobian (J_ij) of the functions FVEC(i).
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 ! for i /= 1 and i /= j and j == 1
246                 k = 1  ! k =1 is a special case for computing Jacobioan
247                 FJAC(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 ! for i /= 1 and i /= j and j /= 1
260     ! no do loop for k here, just one term k = j, the rest will be zero.
261                 k = 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 ! for checks on i
269             ENDDO ! for j
270           ENDDO ! for i
271     ! End of computing jacobian (J_ij).
272     
273           RETURN
274           END SUBROUTINE FUNC_JACOBI_EVAL
275     
276     
277     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
278     !
279     !  subroutine name: ludcmp(a,n,np,indx,d, calledFrom)
280     !  Purpose: Replaces matrix a (n,n) by the LU decomposition of a rowwise
281     !           permutation of itself. Used in combination with lubksb.
282     !
283     !  Literature/Document References:
284     !     Numerical Recipies in Fortran 77, page 38-39
285     !
286     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
287     !
288           subroutine ludcmp(a,n,np,indx,d,calledFrom)
289           USE compar
290     
291           implicit none
292     !-----------------------------------------------
293     ! Dummy arguments
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     ! Local parameters
303     !-----------------------------------------------
304           integer :: nmax
305           double precision :: TINY
306           parameter (NMAX=500, TINY=1.0D-20)
307     !-----------------------------------------------
308     ! Local variables
309     !-----------------------------------------------
310           integer :: i, j, k, imax
311           double precision :: vv(NMAX)
312           double precision :: aamax, sum, dum
313     !-----------------------------------------------
314     
315           d = 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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
373     !
374     !  subroutine name: lubksb (a,n,np,indx,b)
375     !  Purpose: solves the set of n linear equations A(n,n).X(n) = B(n).
376     !
377     !
378     !  Literature/Document References:
379     !     Numerical Recipies in Fortran 77, page 39.
380     !
381     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
382     
383           subroutine lubksb (a,n,np,indx,b)
384     
385           implicit none
386     !-----------------------------------------------
387     ! Dummy arguments
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     ! Local variables
396     !-----------------------------------------------
397           integer :: i, ii, j, ll
398           double precision :: sum
399     !-----------------------------------------------
400     
401           ii=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