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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: ordinary_diff(s,mi,sigmai,alpha,phii,T,phi,ni,n,rhoi,rho,
4     !             m,mu,sigma,chi,zeta0,theta,Ti,DT,nu,Dij,I_ilj,dTl_dnj,
5     !             dzeta0_dnj,dchi0il_dnj)
6     !
7     !  author:  C. Hrenya, Feb 2009
8     !
9     !  Purpose: find ordinary diffusion coefficient according to GHD polydisperse KT
10     !
11     !  Literature/References:
12     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
13     !
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
15     
16           subroutine ordinary_diff(s,mi,sigmai,alpha,phii,T,phi,ni,n,rhoi, &
17               rho,m,mu,sigma,chi,zeta0,theta,Ti,DT,nu,Dij,I_ilj,dTl_dnj, &
18               dzeta0_dnj,dchi0il_dnj)
19           Implicit NONE
20     
21           integer s, indx(s)
22     
23           double precision mi(s),sigmai(s),alpha(s,s),phii(s),T,phi, &
24                           ni(s),n,rhoi(s),rho,m,mu(s,s),sigma(s,s), &
25                           chi(s,s),zeta0, &
26                           theta(s),Ti(s),DT(s),nu(s,s),Dij(s,s), &
27                           I_ilj(s,s,s),dTl_dnj(s,s),dzeta0_dnj(s), &
28                           dchi0il_dnj(s,s,s)
29     
30           integer i,j,k,l,kk
31           double precision kronecker(s,s),M1,M2,M3,dphi_dnl(s), &
32                           dM1_dnl(s),dM2_dnl(s),dM3_dnl(s), &
33                           dni_dnl(s,s), &
34                           dmuioverT_dnl(s,s),perturbation,delta_ni(s), &
35                           delta_phii(s),delta_rhoi(s),chi_ij, &
36                           theta_pos(s),Ti_pos(s), &
37                           zeta0_pos,group1(s,s),group2(s,s), &
38                           p_pos,niTi_pos(s),theta_neg(s), &
39                           Ti_neg(s),zeta0_neg,p_neg,niTi_neg(s), &
40                           dp_dnj(s),dniTi_dnj(s,s),sum1(s,s), &
41                           Amat(s,s),bmat(s,s), &
42                           Amat0(s,s),bmat0(s)
43           double precision d, twoDeltaNi
44     
45           integer NP
46           parameter (NP=15)     !max no. of linear equations to solve
47     
48           double precision pi
49           parameter (pi=3.14159265458979323846d0)
50     
51           do i=1,s
52              do j=1,s
53                 if (i.eq.j) then
54                    kronecker(i,j) = 1.d0
55                 else
56                    kronecker(i,j) = 0.d0
57                 endif
58              enddo
59           enddo
60     
61     !calculate Mn's - p 6.1 CMH notes
62           M1=0.d0
63           M2=0.d0
64           M3=0.d0
65           do i=1,s
66              M1 = M1 + ni(i)*sigmai(i)
67              M2 = M2 + ni(i)*sigmai(i)**2
68              M3 = M3 + ni(i)*sigmai(i)**3
69           enddo
70           M1 = M1/n
71           M2 = M2/n
72           M3 = M3/n
73     
74     !calculate partial derivative of phi and Mn's wrt n_l - p 6.2 CMH notes
75           do l=1,s
76              dphi_dnl(l) = pi/6.d0*sigmai(l)**3
77              dM1_dnl(l)  = (sigmai(l)-M1)/n
78              dM2_dnl(l)  = (sigmai(l)**2-M2)/n
79              dM3_dnl(l)  = (sigmai(l)**3-M3)/n
80              do i=1,s
81                 dni_dnl(i,l)  = kronecker(i,l)
82              enddo
83           enddo
84     
85     !calculate partial derivative of mu_i/T wrt n_l - p 6.4 CMH notes;
86     !this Matlab code was generated by Mathematica, where the symbolic derivative was evaluated
87     ! note that 1/ni term will cancel with -1/ni defined in I_ilj
88           do i=1,s
89              do l=1,s
90                 dmuioverT_dnl(i,l) = &
91                  (3.d0*sigmai(i)*phi*dM2_dnl(l))/(M3*(1.d0 - phi)) -  &
92            (3.d0*sigmai(i)*M2*phi*dM3_dnl(l))/(M3**2*(1.d0 - phi)) +  &
93            dphi_dnl(l)/(1.d0 - phi) + (3.d0*sigmai(i)*M2* &
94            dphi_dnl(l))/(M3*(1.d0 - phi)) +  &
95            (3.d0*sigmai(i)*M2*phi*dphi_dnl(l))/(M3*(1.d0 - phi)**2)+  &
96            3.d0*sigmai(i)**2*((phi*dM1_dnl(l))/(M3*(1.d0 - phi)) +  &
97               (2.d0*DLog(1.d0 - phi)*M2*dM2_dnl(l))/M3**2 +  &
98               (2.d0*M2*phi*dM2_dnl(l))/(M3**2*(1.d0 - phi)**2) -  &
99               (2.d0*DLog(1.d0 - phi)*M2**2*dM3_dnl(l))/M3**3 -  &
100               (2.d0*M2**2*phi*dM3_dnl(l))/(M3**3*(1.d0 - phi)**2) -  &
101               (M1*phi*dM3_dnl(l))/(M3**2*(1.d0 - phi)) +  &
102               (M2**2*dphi_dnl(l))/(M3**2*(1.d0 - phi)**2) -  &
103               (M2**2*dphi_dnl(l))/(M3**2*(1.d0 - phi)) +  &
104               (M1*dphi_dnl(l))/(M3*(1.d0 - phi)) +  &
105               (2.d0*M2**2*phi*dphi_dnl(l))/(M3**2*(1.d0 - phi)**3) +  &
106               (M1*phi*dphi_dnl(l))/(M3*(1.d0 - phi)**2)) +  &
107            sigmai(i)**3*((3.d0*M2*phi**2*dM1_dnl(l))/ &
108               (M3**2*(1.d0 - phi)**2) -  &
109               (6.d0*DLog(1.d0 - phi)*M2**2*dM2_dnl(l))/M3**3 +  &
110               (3.d0*M1*phi**2*dM2_dnl(l))/(M3**2*(1.d0 - phi)**2) -  &
111               (3.d0*M2**2*phi*(2.d0 - 5.d0*phi + phi**2)* &
112               dM2_dnl(l))/(M3**3*(1.d0 - phi)**3) +  &
113               (6.d0*DLog(1.d0 - phi)*M2**3*dM3_dnl(l))/M3**4 -  &
114               (phi*dM3_dnl(l))/(M3**2*(1.d0 - phi)) -  &
115               (6.d0*M1*M2*phi**2*dM3_dnl(l))/ &
116               (M3**3*(1.d0 - phi)**2) +  &
117               (3.d0*M2**3*phi*(2.d0 - 5.d0*phi + phi**2)* &
118               dM3_dnl(l))/(M3**4*(1.d0 - phi)**3) +  &
119               (2.d0*M2**3*dphi_dnl(l))/(M3**3*(1.d0 - phi)) +  &
120                dphi_dnl(l)/(M3*(1.d0 - phi)) +  &
121               (6.d0*M1*M2*phi*dphi_dnl(l))/(M3**2*(1.d0 - phi)**2) +  &
122               (phi*dphi_dnl(l))/(M3*(1.d0 - phi)**2) +  &
123               (6.d0*M1*M2*phi**2*dphi_dnl(l))/ &
124               (M3**2*(1.d0 - phi)**3) -  &
125               (M2**3*(2.d0 - 5.d0*phi + phi**2)*dphi_dnl(l))/ &
126               (M3**3*(1.d0 - phi)**3) -  &
127               (3.d0*M2**3*phi*(2.d0 - 5.d0*phi +  &
128               phi**2)*dphi_dnl(l))/(M3**3*(1.d0 - phi)**4) -  &
129               (M2**3*phi*(-5.d0*dphi_dnl(l) +  &
130               2.d0*phi*dphi_dnl(l)))/(M3**3*(1.d0 - phi)**3))
131              enddo
132           enddo
133     
134     !calculate partial derivative of pair correlation function wrt n_j - p 6.5 CMH notes;
135     !this Matlab code was generated by Mathematica, where the symbolic derivative was evaluated
136           do i=1,s
137              do l=1,s
138                 do j=1,s
139                   dchi0il_dnj(i,l,j) = &
140                  (1.5d0*sigmai(i)*sigmai(l)*phi*dM2_dnl(j))/ &
141            (sigma(i,l)*M3*(1.d0 - phi)**2) +  &
142            (1.d0*sigmai(i)**2*sigmai(l)**2*M2*phi**2*dM2_dnl(j))/ &
143             (sigma(i,l)**2*M3**2*(1.d0 - phi)**3) -  &
144            (1.5d0*sigmai(i)*sigmai(l)*M2*phi*dM3_dnl(j))/ &
145            (sigma(i,l)*M3**2*(1.d0 - phi)**2) -  &
146            (1.d0*sigmai(i)**2*sigmai(l)**2*M2**2*phi**2* &
147            dM3_dnl(j))/(sigma(i,l)**2*M3**3*(1.d0 - phi)**3) +  &
148            dphi_dnl(j)/(1.d0 - phi)**2 + (1.5d0*sigmai(i)* &
149            sigmai(l)*M2*dphi_dnl(j))/ &
150             (sigma(i,l)*M3*(1.d0 - phi)**2) + (1.d0*sigmai(i)**2* &
151            sigmai(l)**2*M2**2*phi*dphi_dnl(j))/ &
152             (sigma(i,l)**2*M3**2*(1.d0 - phi)**3) +  &
153            (3.d0*sigmai(i)*sigmai(l)*M2*phi*dphi_dnl(j))/ &
154            (sigma(i,l)*M3*(1.d0 - phi)**3) + (1.5d0*sigmai(i)**2* &
155            sigmai(l)**2*M2**2*phi**2*dphi_dnl(j))/ &
156             (sigma(i,l)**2*M3**2*(1.d0 - phi)**4)
157                 enddo
158              enddo
159           enddo
160     
161     !determination of I_ilj - p 6.1 of CMH notes
162     !      do i=1,s
163     !         do l=1,s
164     !            do j=1,s
165     !               I_ilj(i,l,j) = (1.5d0/pi*(kronecker(j,l)*ni(l)* &
166     !                 dmuioverT_dnl(i,l)  &
167     !                 -ni(l)/ni(i)*kronecker(j,l)*kronecker(i,j)  &
168     !                 -4.d0*pi/3.d0*kronecker(j,l)*ni(l)*chi(i,j)* &
169     !                 sigma(i,j)**3)-ni(l)**2*sigma(i,l)**3* &
170     !                 dchi0il_dnj(i,l,j))
171     !            enddo
172     !         enddo
173     !      enddo
174     !
175     !CMH correction (10/9/09) to previous error in notes - this I_ilj for binary mix only
176           I_ilj(1,2,1) = ( 3.d0/2.d0/pi*(dmuioverT_dnl(1,1)-4.d0*pi/3.d0*chi(1,1)*sigma(1,1)**3) &
177               -ni(1)*sigma(1,1)**3*dchi0il_dnj(1,1,1) &
178               -ni(2)*sigma(1,2)**3*dchi0il_dnj(1,2,1) ) &
179               /(chi(1,2)*sigma(1,2)**3)
180     
181           I_ilj(1,2,2) = ( 3.d0/2.d0/pi*(dmuioverT_dnl(1,2)- 4.d0*pi/3.d0*chi(1,2)*sigma(1,2)**3) &
182               -ni(1)*chi(1,1)*sigma(1,1)**3/chi(1,1)*dchi0il_dnj(1,1,2) &
183               -ni(2)*chi(1,2)*sigma(1,2)**3/chi(1,2)*dchi0il_dnj(1,2,2) ) &
184               /(chi(1,2)*sigma(1,2)**3)
185     
186           I_ilj(2,1,1) = ( 3.d0/2.d0/pi*(dmuioverT_dnl(2,1)- 4.d0*pi/3.d0*chi(2,1)*sigma(2,1)**3) &
187              -ni(2)*chi(2,2)*sigma(2,2)**3/chi(2,2)*dchi0il_dnj(2,2,1) &
188              -ni(1)*chi(2,1)*sigma(2,1)**3/chi(2,1)*dchi0il_dnj(2,1,1) ) &
189               /(chi(2,1)*sigma(2,1)**3)
190     
191           I_ilj(2,1,2) = ( 3.d0/2.d0/pi*(dmuioverT_dnl(2,2)-4.d0*pi/3.d0*chi(2,2)*sigma(2,2)**3) &
192               -ni(2)*chi(2,2)*sigma(2,2)**3/chi(2,2)*dchi0il_dnj(2,2,2) &
193               -ni(1)*chi(2,1)*sigma(2,1)**3/chi(2,1)*dchi0il_dnj(2,1,2) ) &
194               /(chi(2,1)*sigma(2,1)**3)
195     
196           I_ilj(1,1,1) = 0.d0
197           I_ilj(1,1,2) = 0.d0
198           I_ilj(2,2,1) = 0.d0
199           I_ilj(2,2,2) = 0.d0
200     !
201     ! for smax > 2, use GHD theory in conjunction with standard Enskog theory (i.e. I_ilj's = 0)
202     ! This fix is only temporary until these terms are derived by C.M. Hrenya and co-workers.
203     !
204           if(s > 2) then
205             do i=1,s
206                do l=1,s
207                   do j=1,s
208                      I_ilj(i,l,j) = 0d0
209                   enddo
210                enddo
211             enddo
212           endif
213     !
214     
215     !Numerical evaulation of several partial derivatives - p 6.6 of CMH notes.
216     !Each of the partial derivatives is with respect to nj, where T and other
217     !ni are constant, so need to
218     !   1) perturb nj a small positive amount (and correspondingly n,phii,phi)
219     !   2) evaulate dependent quantity of interest
220     !   3) perturb nj a small negative amount (and correspondingly n,phii,phi)
221     !   4) evaluate dependent quantity of interest
222     !   5) calculate partial derivative
223     !   6) return nj & n to original values
224     
225           perturbation = 0.000001d0    ! small perturb of current value in each direction
226           do j=1,s
227              !calculate perturbation amounts
228                 if(ni(j)/=0d0) then
229                   delta_ni(j)  = perturbation*ni(j)
230                   delta_phii(j)= perturbation*ni(j)*pi/6.d0*sigmai(j)**3
231                   delta_rhoi(j)   = perturbation*ni(j)*mi(j) !;
232                   twoDeltaNi        = 2d0*delta_ni(j)
233                 else
234                   delta_ni(j)   = perturbation
235                   delta_phii(j) = perturbation*pi/6.d0*sigmai(j)**3
236                   delta_rhoi(j) = perturbation*mi(j)
237                   twoDeltaNi      = delta_ni(j)
238                 endif
239              !perturb current nj values in + direction
240                 ni(j)   = ni(j) + delta_ni(j)
241                 n       = n     + delta_ni(j)
242                 phii(j) = phii(j)+ delta_phii(j)
243                 phi     = phi    + delta_phii(j)
244                 rhoi(j) = rhoi(j) + delta_rhoi(j) !;
245                 do i=1,s
246                    do k=1,s
247                       call chi_ij_GHD(s,i,k,sigmai,phi,ni,chi_ij)
248                       chi(i,k) = chi_ij
249                    enddo
250                 enddo
251              !evaluate dependent quantites at increased value of nj et al.
252                 call cooling_rate(s,mi,ni,n,m,Ti,T,chi,sigmai,alpha,rhoi, &
253                                  theta_pos)
254                 do i=1,s
255                    Ti_pos(i) = mi(i)*T/(m*theta_pos(i))
256                 enddo
257                 do l=1,s
258                    group1(1,l)  = chi(1,l)*ni(l)*mu(l,1)* &
259                                  sigma(1,l)**2*(1d0+alpha(1,l))
260                    group2(1,l)  = mu(l,1)/2.d0*(1d0+alpha(1,l))
261                 enddo
262                 zeta0_pos = 0d0 !;
263                 do k=1,s
264                    zeta0_pos = zeta0_pos + group1(1,k)* &
265                            dsqrt((theta_pos(1)+theta_pos(k)) &
266                            /(theta_pos(1)*theta_pos(k))) &
267                            * ((1.d0-group2(1,k)*(theta_pos(1)+ &
268                            theta_pos(k))/theta_pos(k)))
269                 enddo
270                 zeta0_pos = 8.d0/3.d0*dsqrt(2.d0*pi*T/m)*zeta0_pos
271                 call pressure (s,alpha,ni,n,mu,sigma,chi,T,Ti_pos, &
272                               p_pos)
273                 do i=1,s
274                    niTi_pos(i) = ni(i)*Ti_pos(i)
275                 enddo
276              !perturb current nj values in - direction
277                 if(ni(j) /= delta_ni(j)) then
278                   ni(j)   = ni(j) - 2.d0*delta_ni(j)
279                   n       = n     - 2.d0*delta_ni(j)
280                   phii(j) = phii(j) - 2.d0*delta_phii(j)
281                   phi     = phi     - 2.d0*delta_phii(j)
282                   rhoi(j) = rhoi(j) - 2.d0*delta_rhoi(j) !;
283                 else ! bring ni back to zero and do a one-sided (1st order) derivative
284                   ni(j)   = ni(j) - delta_ni(j)
285                   n       = n     - delta_ni(j)
286                   phii(j) = phii(j) - delta_phii(j)
287                   phi     = phi     - delta_phii(j)
288                   rhoi(j) = rhoi(j) - delta_rhoi(j)
289                 endif
290                 do i=1,s
291                    do k=1,s
292                       call chi_ij_GHD(s,i,k,sigmai,phi,ni,chi_ij)
293                       chi(i,k) = chi_ij
294                    enddo
295                 enddo
296              !evaluate dependent quantites at decreased value of nj et al.
297                 call cooling_rate(s,mi,ni,n,m,Ti,T,chi,sigmai,alpha,rhoi, &
298                                  theta_neg)
299                 do i=1,s
300                    Ti_neg(i) = mi(i)*T/(m*theta_neg(i))
301                 enddo
302                 do l=1,s
303                    group1(1,l)  = chi(1,l)*ni(l)*mu(l,1)* &
304                                  sigma(1,l)**2*(1d0+alpha(1,l))
305                    group2(1,l)  = mu(l,1)/2.d0*(1d0+alpha(1,l))
306                 enddo
307                 zeta0_neg = 0d0 !;
308                 do k=1,s
309                    zeta0_neg = zeta0_neg + group1(1,k)* &
310                            dsqrt((theta_neg(1)+theta_neg(k)) &
311                            /(theta_neg(1)*theta_neg(k))) &
312                            * ((1.d0-group2(1,k)*(theta_neg(1)+ &
313                            theta_neg(k))/theta_neg(k)))
314                 enddo
315                 zeta0_neg = 8.d0/3.d0*dsqrt(2.d0*pi*T/m)*zeta0_neg
316                 call pressure (s,alpha,ni,n,mu,sigma,chi,T,Ti_neg, &
317                               p_neg)
318                 do i=1,s
319                    niTi_neg(i) = ni(i)*Ti_neg(i)
320                 enddo
321              !find partial derivatives
322                 dzeta0_dnj(j) = (zeta0_pos - zeta0_neg)/ &
323                                ( twoDeltaNi)
324                 dp_dnj(j) = (p_pos - p_neg)/( twoDeltaNi)
325                 do l=1,s
326                    dTl_dnj(l,j) = (Ti_pos(l) - Ti_neg(l))/ &
327                                  ( twoDeltaNi)
328                 enddo
329                 do i=1,s
330                    dniTi_dnj(i,j) = (niTi_pos(i) - niTi_neg(i))/ &
331                                    ( twoDeltaNi)
332                 enddo
333              !reset nj et al. to original values
334                 if(ni(j) /= 0d0) then
335                   ni(j)   = ni(j) + delta_ni(j)
336                   n       = n     + delta_ni(j)
337                   phii(j) = phii(j) + delta_phii(j)
338                   phi     = phi     + delta_phii(j)
339                   rhoi(j) = rhoi(j) + delta_rhoi(j) !;
340                 endif
341                 do i=1,s
342                    do k=1,s
343                       call chi_ij_GHD(s,i,k,sigmai,phi,ni,chi_ij)
344                       chi(i,k) = chi_ij
345                    enddo
346                 enddo
347           enddo
348     
349           do i=1,s
350              do j=1,s
351                 sum1(i,j) = 0.d0    !calculate summation used in b vector - p 6 CMH notes
352              enddo
353           enddo
354           do i=1,s
355              do j=1,s
356                 do l=1,s
357                    sum1(i,j) = sum1(i,j) + chi(i,l)*sigma(i,l)**3* &
358                      mu(l,i)*(1.d0+alpha(i,l))  &
359                      *((Ti(i)/mi(i)+Ti(l)/mi(l))*(kronecker(j,l)  &
360                      +0.5d0*(ni(l)/chi(i,l)* &
361                      dchi0il_dnj(i,l,j)+I_ilj(i,l,j)))  &
362                      +ni(l)/mi(l)*dTl_dnj(l,j))
363                 enddo
364              enddo
365           enddo
366     
367           do i=1,s
368              do j=1,s
369                 Amat(i,j) = (nu(i,j)-0.5d0*zeta0*kronecker(i,j))*  &        !A matrix for solution of Dij (p 7 CMH notes)
370                             mi(j)/mi(i)
371                 bmat(i,j) = rho**2/mi(i)/mi(j)*dzeta0_dnj(j)*DT(i)  &       !b matrix for solution of Dij (p 7 CMH notes)
372                   + rho/mi(i)/mi(j)*dniTi_dnj(i,j)  &
373                   - ni(i)/mi(j)*dp_dnj(j)  &
374                   + 2.d0*pi/3.d0*rho*ni(i)/mi(j)*sum1(i,j)
375              enddo
376           enddo
377     
378     ! this extra kk loop and addition of Amat0 and bmat0 is necessary
379     ! since x & b in Ax=b are s by s matrices rather than vectors of length s,
380     ! whereas LUBSKB is specific to x & b vectors of length s
381           do kk=1,s
382              do i=1,s
383                 do j=1,s
384                     Amat0(i,j) = Amat(i,j)
385                 enddo
386              enddo
387              do i=1,s
388                 bmat0(i) = bmat(i,kk)
389              enddo
390     
391              CALL LUDCMP(Amat0, s, NP, indx, d, 'ordinary_diff') ! solve system of s linear equations using
392              CALL LUBKSB(Amat0, s, NP, indx, bmat0) ! LU decomposition
393     
394              do i=1,s
395                 Dij(i,kk) = bmat0(i)
396              enddo
397           enddo
398     
399           if(s==2) then ! for binary only
400             Dij(1,1) = -mi(2)/mi(1)*Dij(2,1)
401             Dij(1,2) = -mi(2)/mi(1)*Dij(2,2)
402           endif
403     
404           return
405           end subroutine ordinary_diff
406     
407