File: RELATIVE:/../../../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 :: 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     
289           implicit none
290     !-----------------------------------------------
291     ! Dummy arguments
292     !-----------------------------------------------
293           integer, intent(in) :: n
294           double precision, intent(inout) :: a(n,n)
295           integer :: np
296           integer, intent(out) :: indx(n)
297           double precision, intent(out) :: d
298           CHARACTER(len=*), intent(in) :: calledFrom
299     !-----------------------------------------------
300     ! Local parameters
301     !-----------------------------------------------
302           integer :: nmax
303           double precision :: TINY
304           parameter (NMAX=500, TINY=1.0D-20)
305     !-----------------------------------------------
306     ! Local variables
307     !-----------------------------------------------
308           integer :: i, j, k, imax
309           double precision :: vv(NMAX)
310           double precision :: aamax, sum, dum
311     !-----------------------------------------------
312     
313           d = 1.0d0
314           do i = 1,n
315              aamax=0.0d0
316              do j = 1,n
317                 if (dabs(a(i,j)).gt.aamax) aamax = dabs(a(i,j))
318              enddo
319              if (aamax.eq.0.0d0) then
320                if(myPE==PE_IO) write(*,*) &
321                   'Singular Matrix in ludcmp called from ', calledFrom
322                call mfix_exit(myPE)
323              endif
324              vv(i) = 1.d0/aamax
325           enddo
326           do j = 1,n
327              do i = 1,j-1
328                 sum = a(i,j)
329                 do k = 1,i-1
330                    sum = sum-a(i,k)*a(k,j)
331                 enddo
332                 a(i,j) = sum
333              enddo
334              aamax = 0.0d0
335              do i = j,n
336                 sum = a(i,j)
337                 do k = 1,j-1
338                    sum = sum-a(i,k)*a(k,j)
339                 enddo
340                 a(i,j) = sum
341                 dum = vv(i)*dabs(sum)
342                 if (dum.ge.aamax) then
343                    imax = i
344                    aamax = dum
345                 endif
346              enddo
347              if (j.ne.imax) then
348                 do k = 1,n
349                    dum = a(imax,k)
350                    a(imax,k) = a(j,k)
351                    a(j,k) = dum
352                 enddo
353                 d = -d
354                 vv(imax) = vv(j)
355              endif
356              indx(j) = imax
357              if (a(j,j).eq.0.0d0) a(j,j) = TINY
358              if (j.ne.n) then
359                 dum = 1.0d0/a(j,j)
360                 do i = j+1,n
361                    a(i,j) = a(i,j)*dum
362                 enddo
363              endif
364           enddo
365     
366           return
367           end subroutine ludcmp
368     
369     
370     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
371     !
372     !  subroutine name: lubksb (a,n,np,indx,b)
373     !  Purpose: solves the set of n linear equations A(n,n).X(n) = B(n).
374     !
375     !
376     !  Literature/Document References:
377     !     Numerical Recipies in Fortran 77, page 39.
378     !
379     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
380     
381           subroutine lubksb (a,n,np,indx,b)
382     
383           implicit none
384     !-----------------------------------------------
385     ! Dummy arguments
386     !-----------------------------------------------
387           integer, intent(in) :: n
388           double precision, intent(in) :: a(n,n)
389           integer :: np
390           integer, intent(in) :: indx(n)
391           double precision, intent(inout) :: b(n)
392     !-----------------------------------------------
393     ! Local variables
394     !-----------------------------------------------
395           integer :: i, ii, j, ll
396           double precision :: sum
397     !-----------------------------------------------
398     
399           ii=0
400           do i=1,n
401             ll=indx(i)
402             sum=b(ll)
403             b(ll)=b(i)
404             if (ii.ne.0) then
405               do j=ii,i-1
406                 sum=sum-a(i,j)*b(j)
407               enddo
408              elseif (sum.ne.0.d0) then
409                ii=i
410              endif
411              b(i)=sum
412           enddo
413            do i=n,1,-1
414              sum=b(i)
415              if (i.lt.n) then
416                do j=i+1,n
417                  sum=sum-a(i,j)*b(j)
418                enddo
419              endif
420              b(i)=sum/a(i,i)
421            enddo
422            return
423            end subroutine lubksb
424     
425