MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
crit_point_mixtures.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! Module Module_Heidemann_Khalil
3 !
4 ! This module ....
5 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
6 
7 Module module_heidemann_khalil
8 
9  implicit none
10  save
11 
12  real :: error_condition2
13 
14 End Module module_heidemann_khalil
15 
16 
17 
18 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
19 ! SUBROUTINE Heidemann_Khalil
20 !
21 ! This subroutine ....
22 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
23 
24 SUBROUTINE heidemann_khalil
25 
26  USE basic_variables
27  USE module_heidemann_khalil
28  USE solve_nonlin
29  IMPLICIT NONE
30 
31  !-----------------------------------------------------------------------------
32  INTERFACE
33  SUBROUTINE heidemann_khalil_obj ( iter_no, y, residu, dummy )
34  INTEGER, INTENT(IN) :: iter_no
35  REAL, INTENT(IN) :: y(iter_no)
36  REAL, INTENT(OUT) :: residu(iter_no)
37  INTEGER, INTENT(IN OUT) :: dummy
38  END SUBROUTINE heidemann_khalil_obj
39  END INTERFACE
40 
41  INTEGER :: info
42  REAL, ALLOCATABLE :: y(:),diag(:),residu(:)
43 
44  INTEGER :: i, count, nphas_save
45  REAL :: error0, error1, dense0, dfdx, delta_rho
46  CHARACTER (LEN=2) :: ensemble_save
47  !-----------------------------------------------------------------------------
48 
49  ensemble_save = ensemble_flag
50  ensemble_flag = 'tv'
51 
52  nphas_save = nphas
53  nphas = 1
54 
55  ! xiF(2) = 0.5 * ( 0.71928411 + 0.72025642 )
56  ! xiF(1) = 1.0 - xiF(2)
57  ! dense(1) = 0.5 * ( 0.159315 + 0.158817 ) + 0.02
58  ! t = 500.0
59 
60  dense0 = dense(1)
61 
62  info = 1
63  n_unkw = ncomp + 1
64  acc_a = 5.e-8
65  step_a = 1.e-8
66 
67 
68  ALLOCATE( y(n_unkw), diag(n_unkw), residu(n_unkw) )
69  DO i = 1,ncomp
70  y(i) = 1.0 / sqrt( REAL( ncomp ) )
71  END DO
72  y(ncomp+1) = t
73  count = 0
74  error0 = 1.0
75 
76  DO WHILE ( abs( error0) > 0.001 .AND. count < 20 )
77  count = count + 1
78 
79  dense(1) = dense0 + 0.0001
80  CALL hbrd (heidemann_khalil_obj, n_unkw, y, residu, step_a, acc_a, info, diag)
81  error1 = error_condition2
82  IF (sum( abs( residu(1:n_unkw) ) ) > 1.e-5) write (*,*) 'caution: error 1st inner loop', &
83  sum( abs( residu(1:n_unkw) ) )
84 
85  dense(1) = dense0
86  CALL hbrd (heidemann_khalil_obj, n_unkw, y, residu, step_a, acc_a, info, diag)
87  IF (sum( abs( residu(1:n_unkw) ) ) > 1.e-5) write (*,*) 'caution: error 2nd inner loop', &
88  sum( abs( residu(1:n_unkw) ) )
89  error0 = error_condition2
90 
91  ! write (*,'(a,4G18.10)') ' t, p, eta error', t, p, dense(1), error0
92  ! read (*,*)
93  dfdx = ( error1 - error0 ) / 0.0001
94  delta_rho = min( error0 / dfdx, 0.02)
95  delta_rho = max( delta_rho, -0.02)
96  dense0 = dense0 - delta_rho
97  dense(1) = dense0
98 
99  END DO
100 
101  !tc = t
102  !pc = p
103 
104  DEALLOCATE( y, diag, residu )
105 
106  ensemble_flag = ensemble_save
107  nphas = nphas_save
108  IF ( abs( error_condition2 ) > 1.e-1 ) write (*,*) 'caution: error outer loop', abs( error_condition2 )
109 
110 END SUBROUTINE heidemann_khalil
111 
112 
113 
114 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
115 ! SUBROUTINE Heidemann_Khalil
116 !
117 ! This subroutine ....
118 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
119 
120 SUBROUTINE heidemann_khalil_obj ( iter_no, y, residu, dummy )
121 
122  USE parameters, ONLY: pi
123  USE basic_variables
124  USE module_heidemann_khalil
125  IMPLICIT NONE
126 
127  !-----------------------------------------------------------------------------
128  INTEGER, INTENT(IN) :: iter_no
129  REAL, INTENT(IN) :: y(iter_no)
130  REAL, INTENT(OUT) :: residu(iter_no)
131  INTEGER, INTENT(IN OUT) :: dummy
132 
133  !-----------------------------------------------------------------------------
134  INTEGER :: i, j, k
135  REAL, DIMENSION(nc) :: dn
136  REAL, DIMENSION(nc) :: dhs, rhoi00, rhoi0
137  REAL :: rho, d_rho
138  ! REAL :: lnphi(np,nc)
139  REAL :: qij(nc,nc), qij0(nc,nc), qijk(nc,nc,nc)
140  ! CHARACTER (LEN=3) :: char_len
141  !-----------------------------------------------------------------------------
142 
143  dn(1:ncomp) = y(1:ncomp)
144  t = y(ncomp+1)
145  !dense(1) = y(ncomp+2)
146 
147  dhs(1:ncomp) = parame(1:ncomp,2) * ( 1.0 - 0.12 *exp( -3.0*parame(1:ncomp,3)/t ) )
148  rho = dense(1) / sum( pi/6.0*xif(1:ncomp)*parame(1:ncomp,1)* dhs(1:ncomp)**3 )
149  rhoi00(1:ncomp) = xif(1:ncomp)*rho
150 
151  d_rho = 0.000001
152 
153  DO k = 1, ncomp
154 
155  rhoi0(1:ncomp) = rhoi00(1:ncomp)
156  rhoi0(k) = rhoi00(k) + d_rho
157  CALL qij_matrix ( rhoi0, dhs, d_rho, qij )
158  qij0(:,:) = qij(:,:)
159 
160  rhoi0(1:ncomp) = rhoi00(1:ncomp)
161  CALL qij_matrix ( rhoi0, dhs, d_rho, qij )
162 
163  DO i = 1, ncomp
164  DO j = 1, ncomp
165  qijk(i,j,k) = ( qij0(i,j) - qij(i,j) ) / d_rho
166  ! write(*,*) i,j,k,qijk(i,j,k)
167  END DO
168  END DO
169 
170  END DO
171 
172  ! write (*,'(a,4G18.10)') 'det',qij(2,2)*qij(1,1) - qij(2,1)*qij(1,2)
173  ! write (*,'(a,3G18.10)') ' t,p,eta ', t, p, dense(1)
174  DO j = 1, ncomp
175  residu(j) = sum( qij(1:ncomp,j)*dn(1:ncomp) )
176  END DO
177  residu(ncomp+1) = 1.0 - sum( dn(1:ncomp)*dn(1:ncomp) )
178 
179  error_condition2 = 0.0
180  DO k = 1, ncomp
181  DO i = 1, ncomp
182  DO j = 1, ncomp
183  error_condition2 = error_condition2 + qijk(i,j,k) * dn(i) * dn(j) * dn(k)
184  END DO
185  END DO
186  END DO
187  error_condition2 = error_condition2 * 1.e-4 ! the values are scaled down to prevent numerical dominance of this error
188  !residu(ncomp+2) = error_condition2
189 
190  ! write (char_len,'(I3)') ncomp+2
191  ! write (*,'(a,'//char_len//'G18.10)') ' error',residu(1:ncomp+1)
192  ! write (*,*) ' '
193  ! read (*,*)
194 
195 END SUBROUTINE heidemann_khalil_obj
196 
197 
198 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
199 ! SUBROUTINE Heidemann_Khalil
200 !
201 ! This subroutine calculates the spinodal for a binary mixture to given
202 ! T, p and for a given starting value of the density vector rhoi
203 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
204 
205 SUBROUTINE binary_tp_spinodal ( rhoi )
206 
207  USE basic_variables
208  USE module_heidemann_khalil
209  USE solve_nonlin
210  IMPLICIT NONE
211 
212  !-----------------------------------------------------------------------------
213  INTERFACE
214  SUBROUTINE binary_spinodal_obj ( iter_no, y, residu, dummy )
215  INTEGER, INTENT(IN) :: iter_no
216  REAL, INTENT(IN) :: y(iter_no)
217  REAL, INTENT(OUT) :: residu(iter_no)
218  INTEGER, INTENT(IN OUT) :: dummy
219  END SUBROUTINE binary_spinodal_obj
220  END INTERFACE
221 
222  REAL, dimension(nc) :: rhoi
223  !-----------------------------------------------------------------------------
224 
225  INTEGER :: info
226  REAL, ALLOCATABLE :: y(:), diag(:), residu(:)
227 
228  INTEGER :: i, nphas_save
229  CHARACTER (LEN=2) :: ensemble_save
230  !-----------------------------------------------------------------------------
231 
232  ensemble_save = ensemble_flag
233  ensemble_flag = 'tv'
234 
235  nphas_save = nphas
236  nphas = 1
237 
238  info = 1
239  n_unkw = ncomp + 2
240  acc_a = 5.e-8
241  step_a = 1.e-8
242 
243 
244  ALLOCATE( y(n_unkw), diag(n_unkw), residu(n_unkw) )
245  DO i = 1, ncomp
246  y(i) = 1.0 / sqrt( REAL( ncomp ) )
247  END DO
248  DO i = 1, ncomp
249  y( ncomp + i ) = rhoi( i )
250  END DO
251 
252  CALL hbrd (binary_spinodal_obj, n_unkw, y, residu, step_a, acc_a, info, diag)
253 
254  DO i = 1, ncomp
255  rhoi( i ) = y( ncomp + i )
256  END DO
257  write (*,*) 'info',info
258  write (*,*) 'rhoi',rhoi(1:ncomp)
259  !write (*,*) 'eta',PI / 6.0 * sum( rhoi(1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
260  write (*,*) 'x',rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
261 
262  DEALLOCATE( y, diag, residu )
263 
264  ensemble_flag = ensemble_save
265  nphas = nphas_save
266 
267 END SUBROUTINE binary_tp_spinodal
268 
269 
270 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
271 ! SUBROUTINE Heidemann_Khalil
272 !
273 ! This subroutine ....
274 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
275 
276 SUBROUTINE binary_spinodal_obj ( iter_no, y, residu, dummy )
277 
278  USE parameters, ONLY: pi
279  USE basic_variables
280  USE module_heidemann_khalil
281  USE utilities
282  IMPLICIT NONE
283 
284  !-----------------------------------------------------------------------------
285  INTEGER, INTENT(IN) :: iter_no
286  REAL, INTENT(IN) :: y(iter_no)
287  REAL, INTENT(OUT) :: residu(iter_no)
288  INTEGER, INTENT(IN OUT) :: dummy
289 
290  !-----------------------------------------------------------------------------
291  INTEGER :: j, k
292  REAL :: p_calculated, zges, p_sp
293  REAL :: d_rho
294  REAL, DIMENSION(nc) :: dn
295  REAL, DIMENSION(nc) :: dhs, rhoi
296  REAL, DIMENSION(nc,nc) :: qij
297  !-----------------------------------------------------------------------------
298 
299  dn(1:ncomp) = y(1:ncomp)
300  rhoi(1:ncomp) = y( (ncomp+1) : (ncomp+ncomp) )
301  dhs(1:ncomp) = parame(1:ncomp,2) * ( 1.0 - 0.12 *exp( -3.0*parame(1:ncomp,3)/t ) ) ! is this calc. needed?
302 
303  call p_calc ( p_calculated, zges )
304 
305  d_rho = 0.000001
306 
307  DO k = 1, ncomp
308 
309  CALL qij_matrix ( rhoi, dhs, d_rho, qij )
310 
311  END DO
312 
313  !write (*,'(a,4G18.10)') 'det',qij(2,2)*qij(1,1) - qij(2,1)*qij(1,2)
314  write (*,'(a,3G18.10)') ' t,p,eta ', t, p, dense(1)
315  DO j = 1, ncomp
316  residu(j) = sum( qij(1:ncomp,j)*dn(1:ncomp) )
317  END DO
318  residu(ncomp+1) = 1.0 - sum( dn(1:ncomp)*dn(1:ncomp) )
319  write (*,*) 'p_sp has to be handed over to the obj fct properly!'
320  write (*,*) 'Can I simply set p = p_sp? In other words is p altered during the calculation?'
321  stop
322  residu(ncomp+2) = p_sp - p_calculated
323 
324  write (*,'(a,4G18.10)') ' error',residu(1:4)
325  call paus (' ')
326 
327 END SUBROUTINE binary_spinodal_obj
328 
329 
330 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
331 ! SUBROUTINE Heidemann_Khalil
332 !
333 ! This subroutine ....
334 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
335 
336 SUBROUTINE qij_matrix ( rhoi0, dhs, d_rho, qij )
337 
338  USE parameters, ONLY: pi, kbol
339  USE basic_variables
340  USE eos_variables, ONLY: pges
341  IMPLICIT NONE
342 
343  !-----------------------------------------------------------------------------
344  REAL, DIMENSION(nc), INTENT(IN) :: rhoi0
345  REAL, DIMENSION(nc), INTENT(IN) :: dhs
346  REAL, INTENT(IN) :: d_rho
347  REAL, DIMENSION(nc,nc), INTENT(OUT) :: qij
348 
349  !-----------------------------------------------------------------------------
350  INTEGER :: i
351  REAL, DIMENSION(nc) :: rhoi, lnf0, lnf1, lnf2
352  REAL :: lnphi(np,nc), zges
353  !-----------------------------------------------------------------------------
354 
355  DO i = 1, ncomp
356  rhoi(1:ncomp) = rhoi0(1:ncomp)
357  rhoi(i) = rhoi0(i) + d_rho
358  dense(1) = sum( pi/6.0*rhoi(1:ncomp)*parame(1:ncomp,1)* dhs(1:ncomp)**3 )
359  densta(1) = dense(1)
360  xi(1,1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
361  CALL fugacity ( lnphi )
362  p = pges
363  zges = (pges * 1.e-30) / ( kbol*t*sum(rhoi(1:ncomp)) )
364  lnf1(1:ncomp) = lnphi(1,1:ncomp) + log( rhoi(1:ncomp) )
365 
366  IF (rhoi0(i) - d_rho > 0.0 ) THEN
367  rhoi(1:ncomp) = rhoi0(1:ncomp)
368  rhoi(i) = rhoi0(i) - d_rho
369  dense(1) = sum( pi/6.0*rhoi(1:ncomp)*parame(1:ncomp,1)* dhs(1:ncomp)**3 )
370  densta(1) = dense(1)
371  xi(1,1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
372  CALL fugacity ( lnphi )
373  p = pges
374  zges = (pges * 1.e-30) / ( kbol*t*sum(rhoi(1:ncomp)) )
375  lnf2(1:ncomp) = lnphi(1,1:ncomp) + log( rhoi(1:ncomp) )
376  END IF
377 
378  rhoi(1:ncomp) = rhoi0(1:ncomp)
379  dense(1) = sum( pi/6.0*rhoi(1:ncomp)*parame(1:ncomp,1)* dhs(1:ncomp)**3 )
380  densta(1) = dense(1)
381  xi(1,1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
382  CALL fugacity ( lnphi )
383  p = pges
384  zges = (pges * 1.e-30) / ( kbol*t*sum(rhoi(1:ncomp)) )
385  lnf0(1:ncomp) = lnphi(1,1:ncomp) + log( rhoi(1:ncomp) )
386 
387  IF (rhoi0(i) - d_rho > 0.0 ) THEN
388  qij(i,1:ncomp) = ( lnf1(1:ncomp) - lnf2(1:ncomp) ) / (2.0*d_rho) ! qij = d(F/V) / (d_rho_i*d_rho_j)
389  ELSE
390  qij(i,1:ncomp) = ( lnf1(1:ncomp) - lnf0(1:ncomp) ) / d_rho ! qij = d(F/V) / (d_rho_i*d_rho_j)
391  END IF
392  ! write (*,*) i,1,qij(i,1)
393  ! write (*,*) i,2,qij(i,2)
394  END DO
395 
396 END SUBROUTINE qij_matrix
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29