MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
kij-fitting.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! MODULE KIJ_FITTING
3 !
4 ! vectors of experimental data needed for fitting kij-parameters. A vextor
5 ! giving deviations between experiment and calculations is also included
6 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
7 
8 module kij_fitting
9 
10  implicit none
11  save
12 
13  PRIVATE
14  INTEGER :: n_exp
15  REAL, DIMENSION(400) :: x_exp, y_exp, t_exp, p_exp
16  REAL, DIMENSION(400,2) :: deviation
17 
18 
19  PUBLIC :: kij_fit_lij, kij_fit
20 
21 contains
22 
23 
24 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
25 ! SUBROUTINE KIJ_FIT_LIJ
26 !
27 ! routine for fitting kij (and lij, if needed) to binary systems
28 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
29 
30 SUBROUTINE kij_fit_lij
31 
32  USE basic_variables
33  use cg_minimization
34  use optimizer_2d
35  use utilities, only: file_open, paus
36 
37  !-----------------------------------------------------------------------------
38  INTEGER :: i, nadjust, prin
39  REAL, ALLOCATABLE :: optpars(:)
40  REAL :: fmin, t0, h0, machep
41  REAL :: dummy2, dummy3, dummy4, dummy5
42  REAL :: rms1, rms2, aad_m1, aad_m2
43  CHARACTER :: fitin*50, fitin2*71, fitout*80, dummy1*30
44  LOGICAL :: go_up
45 
46  ! INTEGER :: STATUS, iter, nfunc, ngrad
47  ! REAL :: gnorm
48  ! REAL, ALLOCATABLE :: d(:), g(:), xtemp(:), gtemp(:)
49  REAL, ALLOCATABLE :: g(:)
50  !-----------------------------------------------------------------------------
51 
52  ncomp = 2
53  nphas = 2
54  u_in_p = 1.0
55  u_out_p = 1.e6
56 
57 
58  write (*,*) ' CHOOSE THE INPUT FILE (./input_file/kij-dat/...) :'
59  write (*,*) ' SPECIFY FULL NAME:'
60  READ (*,*) fitin
61  ! fitin = 'co2-methanol.dat'
62  ! fitin = 'co2-heptane.dat'
63  ! fitin = 'methane-tetracontane.dat'
64  ! fitin = 'n2-acetone.dat'
65  ! fitin = 'meoh-octane-lle.dat'
66  ! fitin = 'ethygly-hexane-lle.dat'
67  ! fitin = 'co2-benzene_all.dat'
68  ! fitin = '2-propanol-benzene_318K.dat'
69  fitin2 = './input_file/kij-dat/'//fitin
70 
71  CALL file_open(fitin2,87)
72 
73  i = 0
74  READ (87,*) compna(1)
75  READ (87,*) compna(2)
76  go_up = .true.
77 
78  DO WHILE ( go_up )
79  READ (87,*) dummy1,dummy2,dummy3,dummy4,dummy5
80  go_up = .false.
81  IF (dummy1 /= 'end') THEN
82  if ( .NOT.( dummy2 == 1.0 .OR. dummy3 == 1.0 ) .AND. &
83  .NOT.(dummy2 == 0.0 .AND. dummy3 == 0.0 ) ) then
84  i = i + 1
85  x_exp(i) = dummy2
86  y_exp(i) = dummy3
87  p_exp(i) = dummy4 * 1.e5
88  t_exp(i) = dummy5 + 273.15
89  end if
90  go_up = .true.
91  END IF
92  END DO
93  CLOSE (87)
94  n_exp = i
95 
96  write (*,*) ' CHOOSE THE EOS :'
97  write (*,*) ' 0: SAFT EOS '
98  write (*,*) ' 1: PC-SAFT EOS '
99  write (*,*) ' 2: SRK EOS '
100  write (*,*) ' 3: PR EOS '
101  read (*,*) eos
102 
103  IF (eos == 1) THEN
104  write (*,*) ' FURTHER SPECIFY THE EOS :'
105  write (*,*) ' 0: PC-SAFT EOS '
106  write (*,*) ' 1: PCP-SAFT EOS '
107  write (*,*) ' 2: PCIP-SAFT EOS '
108  read (*,*) pol
109  END IF
110 
111  write (*,*) ' CHOOSE THE BINARY PARAMETERS TO BE ADJUSTED :'
112  write (*,*) ' 1: kij '
113  write (*,*) ' 2: kij AND lij'
114  read (*,*) nadjust
115  IF ( nadjust == 2 ) THEN
116  num = 1
117  CALL set_default_eos_numerical
118  write (*,*) 'calculation conducted with numerical derivatives (press enter)'
119  call paus
120  END IF
121  ALLOCATE ( optpars(nadjust) )
122 
123  fitout = fitin
124  IF ( eos == 1 ) fitout = './output_file/kij-pc-saft/'//fitout
125  IF ( eos == 0 ) fitout = './output_file/kij-saft/'//fitout
126  IF ( eos == 2 ) fitout = './output_file/kij-srk/'//fitout
127  IF ( eos == 3 ) fitout = './output_file/kij-pr/'//fitout
128  OPEN ( 88, file=fitout )
129 
130 
131  CALL para_input
132 
133  write (*,'(a,F12.5)') ' current value of kij=',kij(1,2)
134  write (*,'(a)') ' do you want to give another starting value? ( 1 : yes, 0 : no)'
135  read(*,*) i
136  if ( i == 1 ) read (*,*) kij(1,2)
137 
138  !-----------------------------------------------------------------------------
139  t0 = 1.e-4
140  h0 = 0.01
141  prin = 0
142  machep = 1.e-15
143 
144  optpars(1) = kij(1,2) + 1.0
145  IF (nadjust == 2) optpars(2) = lij(1,2) + 1.0
146  !optpars(1) = kij(1,2)
147  !IF (nadjust == 2) optpars(2) = lij(1,2)
148 
149 
150  IF (nadjust > 1) THEN
151  !call PRAXIS2( t0, MACHEP, h0, nadjust, PRIN, optpars, BINDEV_KIJ_p, fmin )
152  ALLOCATE ( g(nadjust) )
153  CALL newton_opt_2d ( bindev_kij_p, optpars, nadjust, 1.e-8, 1.e-8, g, fmin)
154  DEALLOCATE ( g )
155  !ALLOCATE( d(nadjust) )
156  !ALLOCATE( g(nadjust) )
157  !ALLOCATE( xtemp(nadjust) )
158  !ALLOCATE( gtemp(nadjust) )
159  !CALL cg_descent (1.d-5, optpars, nadjust, BINDEV_KIJ, DEV_KIJ_grad, STATUS, &
160  ! gnorm, fmin, iter, nfunc, ngrad, d, g, xtemp, gtemp)
161  ELSE
162  call dim1min( t0, h0, nadjust, optpars, bindev_kij_p, fmin )
163  END IF
164 
165 
166  write (*,*)' Ergebnis: kij=',kij(1,2),' AAD=',fmin / real( n_exp )
167  write(88,*)' Ergebnis: kij=',kij(1,2),' AAD=',fmin / real( n_exp )
168  write (*,*)' Ergebnis: lij=',lij(1,2)
169  write(88,*)' Ergebnis: lij=',lij(1,2)
170 
171 
172  rms1 = 0.0
173  rms2 = 0.0
174  aad_m1 = 0.0
175  aad_m2 = 0.0
176  write (88,*) 'experimental data'
177  DO i = 1, n_exp
178  write (88,'(i4,4(2x,G13.6))') i,x_exp(i),y_exp(i), p_exp(i)/1.e5,t_exp(i)-273.15
179  END DO
180  write (88,*) ' '
181  write (88,*) 'calculated values'
182  DO i = 1, n_exp
183  write (88,'(i4,4(2x,G13.6))') i,x_exp(i)+deviation(i,1), &
184  y_exp(i)+deviation(i,2), &
185  p_exp(i)/1.e5,t_exp(i)-273.15
186  END DO
187  write (88,*) ' '
188  DO i = 1, n_exp
189  write (*,*) ' DEV% ph_1=',100.0*deviation(i,1),' DEV% ph_2=',100.0*deviation(i,2)
190  write (88,*) ' DEVIATION% phase_1=',100.0*deviation(i,1), &
191  ' DEVIATION% phase_2=',100.0*deviation(i,2)
192  rms1 = rms1 + deviation(i,1)**2
193  rms2 = rms2 + deviation(i,2)**2
194  aad_m1 = aad_m1 + abs( deviation(i,1) )
195  aad_m2 = aad_m2 + abs( deviation(i,2) )
196  END DO
197  rms1 = 100.0 * ( rms1 / real( n_exp ) )**0.5
198  rms2 = 100.0 * ( rms2 / real( n_exp ) )**0.5
199  aad_m1 = 100.0 * aad_m1 / real( n_exp )
200  aad_m2 = 100.0 * aad_m2 / real( n_exp )
201 
202  ! write(*,*)' '
203  ! write(*,*)' MAX.DEV% ph_1=',maxerr(1,1),' AT',(maxerr(1,i),i=2,3)
204  ! write(*,*)' MAX.DEV% ph_2=',maxerr(2,1),' AT',(maxerr(2,i),i=2,3)
205  write(*,*)' '
206  write(*,*)' AAD% phase_1=',aad_m1
207  write(*,*)' AAD% phase_2=',aad_m2
208  write(*,*)' '
209  write(*,*)' RMS% phase_1=',rms1
210  write(*,*)' RMS% phase_2=',rms2
211 
212  ! write(88,*)' '
213  ! write(88,*)' MAX. DEVIATION% phase_1=',maxerr(1,1),' AT',(maxerr(1,i),i=2,3)
214  ! write(88,*)' MAX. DEVIATION% phase_2=',maxerr(2,1),' AT',(maxerr(2,i),i=2,3)
215  write(88,*)' '
216  write(88,*)' AAD% phase_1=',aad_m1
217  write(88,*)' AAD% phase_2=',aad_m2
218  write(88,*)' '
219  write(88,*)' RMS% phase_1=',rms1
220  write(88,*)' RMS% phase_2=',rms2
221 
222  write (*,*) ' '
223  write(*,*)' Result: kij=',kij(1,2)
224  IF ( nadjust == 2 ) write(*,*)' lij=',lij(1,2)
225  write(88,*)' Result: kij=',kij(1,2)
226  IF ( nadjust == 2 ) write(88,*)' lij=',lij(1,2)
227  write (*,*) ' '
228  write (88,*) ' '
229 
230  CLOSE (88)
231  DEALLOCATE ( optpars )
232  ! DEALLOCATE( d, g, xtemp, gtemp )
233 
234 END SUBROUTINE kij_fit_lij
235 
236 
237 
238 
239 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
240 !
241 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
242 
243 SUBROUTINE bindev_kij_p ( fmin, optpars, nadjust )
244 
245  USE basic_variables
246  USE starting_values
247  USE utilities
248 
249  !-----------------------------------------------------------------------------
250  INTEGER, INTENT(IN) :: nadjust
251  REAL, INTENT(IN) :: optpars(:)
252  REAL, INTENT(IN OUT) :: fmin
253 
254  !-----------------------------------------------------------------------------
255  INTEGER :: i, converg, iterate_t
256  REAL :: error
257  REAL :: x_sta, y_sta
258  REAL :: dev_x_i, dev_y_i, weigh_x_p, weigh_x_y
259  REAL, DIMENSION(400) :: dev
260  REAL :: xi_save
261  real, dimension(nc) :: rhoi1, rhoi2
262  CHARACTER (LEN=1) :: report_success_step
263  CHARACTER (LEN=30) :: report_character
264  !-----------------------------------------------------------------------------
265 
266 
267  kij(1,2) = optpars(1) - 1.0
268  kij(2,1) = kij(1,2)
269 
270  IF (nadjust == 2) lij(1,2) = optpars(2)-1.0
271  IF (nadjust == 2) lij(2,1) = - lij(1,2)
272 
273  weigh_x_p = 0.25 ! weight of errors in x to errors in p
274  weigh_x_y = 0.8 ! weight of errors in x to errors in y
275  ! weigh_x_p = 0.0 ! weight of errors in x to errors in p
276  ! weigh_x_y = 1.0 ! weight of errors in x to errors in y
277 
278  !-----------------------------------------------------------------------------
279  ! loop for experimental mixture data
280  !-----------------------------------------------------------------------------
281 
282  DO i = 1, n_exp
283 
284  report_character = ' '
285  report_success_step = ' '
286 
287  t = t_exp(i)
288  p = p_exp(i)
289 
290  n_unkw = ncomp ! number of quantities to be iterated
291  it(1) = 'x11' ! iteration of mol fraction of comp.1 phase 1
292  it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
293  sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
294  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
295 
296  x_sta = x_exp(i)
297  y_sta = y_exp(i)
298  if ( x_sta == 0.0 ) x_sta = 0.001 ! these are only starting values
299  if ( x_sta == 1.0 ) x_sta = 0.999 ! these are only starting values
300  if ( y_sta == 0.0 ) y_sta = 0.001 ! these are only starting values
301  if ( y_sta == 1.0 ) y_sta = 0.999 ! these are only starting values
302 
303  xi(1,2) = x_sta
304  xi(2,2) = y_sta
305  xi(1,1) = 1.0 - x_sta
306  xi(2,1) = 1.0 - y_sta
307 
308  !--------------------------------------------------------------------------
309  ! step A for solving the pT-flash problem: Rachford-Rice
310  !--------------------------------------------------------------------------
311 
312  xif(1) = 0.5 * ( xi(1,1) + xi(2,1) )
313  xif(2) = 1.0 - xif(1)
314 
315  dense(1) = 0.4
316  dense(2) = 0.001
317 
318  !outp = 2
319  call rachford_rice ( converg, rhoi1, rhoi2 )
320  if ( converg == 1 ) then
321  CALL occupy_val_init ( rhoi1, rhoi2 )
322  val_conv = val_init
323  report_success_step = 'A'
324  end if
325 
326 
327  !--------------------------------------------------------------------------
328  ! step B for solving the pT-flash problem: VLE-scan
329  !--------------------------------------------------------------------------
330 
331  IF ( converg /= 1 ) THEN
332 
333  call vle_min
334  xif( 1:ncomp ) = xi( 1, 1:ncomp )
335  call start_var_fixed_composition ( converg )
336  if ( converg == 1 ) then
337  val_conv = val_init
338  call restore_converged
339  report_success_step = 'B'
340  end if
341 
342  END IF
343 
344  !--------------------------------------------------------------------------
345  ! step C for solving the pT-flash problem: x-scan and flash
346  !--------------------------------------------------------------------------
347 
348  IF ( converg /= 1 ) THEN
349 
350  outp = 0 ! output to terminal
351  call scan_compositions ( converg )
352  if ( converg == 1 ) then
353  val_conv = val_init
354  call restore_converged
355  report_success_step = 'C'
356  end if
357 
358  END IF
359 
360 
361  !==========================================================================
362  ! calculate entry to the objective function
363  !==========================================================================
364 
365  IF (converg == 1 .AND. abs(val_conv(1)/val_conv(2)-1.0) > 0.25 ) THEN
366 
367  !-----------------------------------------------------------------------
368  ! converged to VLE, first calculate the x-deviations and y-deviations
369  !-----------------------------------------------------------------------
370 
371  if ( val_init(1)/val_init(2) < 0.6 ) then
372  write (*,*) 'swap A',i
373  xi_save = xi(1,2)
374  xi(1,2) = xi(2,2)
375  xi(2,2) = xi_save
376  xi_save = dense(1)
377  dense(1) = dense(2)
378  dense(2) = xi_save
379  end if
380 
381  dev(i) = 0.0
382  dev_x_i = 0.0
383  dev_y_i = 0.0
384  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) dev_x_i = (xi(1,2)-x_exp(i))**2
385  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) dev_y_i = (xi(2,2)-y_exp(i))**2
386 
387  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) deviation(i,1) = xi(1,2)-x_exp(i)
388  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) deviation(i,2) = xi(2,2)-y_exp(i)
389 
390  !-----------------------------------------------------------------------
391  ! secondly, calculate the deviation in bubble point pressure
392  !-----------------------------------------------------------------------
393 
394  xi(1,2) = x_exp(i)
395  xi(1,1) = 1.0 - xi(1,2)
396  iterate_t = 0
397 
398  call bubble_point_rachford_rice ( iterate_t, converg, rhoi1, rhoi2 )
399  IF ( converg == 1 ) CALL occupy_val_init ( rhoi1, rhoi2 )
400  IF ( converg == 1 ) val_conv = val_init
401 
402  !-----------------------------------------------------------------------
403  ! now calculate the contribution of point i to the objective function
404  !-----------------------------------------------------------------------
405 
406  if ( converg == 1 ) then
407  dev(i) = ( weigh_x_p * 1.0/dev_x_i &
408  + ( 1.0-weigh_x_p )*( val_conv(4)/p_exp(i)-1.0 )**(-2) )**(-1)
409  dev(i) = weigh_x_y * dev(i) + ( 1.0 - weigh_x_y ) * dev_y_i
410  else
411  dev(i) = weigh_x_y * dev_x_i + ( 1.0 - weigh_x_y ) * dev_y_i
412  end if
413  if ( converg /= 1 ) write (*,*) 'second step not converged, i=',i
414 
415  ELSE IF (converg == 1 .AND. abs(val_conv(1)/val_conv(2)-1.0) <= 0.25 ) THEN
416 
417  !-----------------------------------------------------------------------
418  ! converged to LLE, first check, whether the phase index has to be exchanged
419  !-----------------------------------------------------------------------
420 
421  if ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 .AND. abs(xi(1,2)-x_exp(i)) > abs(xi(2,2)-x_exp(i)) ) then
422  write (*,*) 'swap1',i
423  xi_save = xi(1,2)
424  xi(1,2) = xi(2,2)
425  xi(2,2) = xi_save
426  end if
427  if ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 .AND. abs(xi(2,2)-y_exp(i)) > abs(xi(1,2)-y_exp(i)) ) then
428  write (*,*) 'swap2',i
429  xi_save = xi(1,2)
430  xi(1,2) = xi(2,2)
431  xi(2,2) = xi_save
432  end if
433 
434  !-----------------------------------------------------------------------
435  ! secondly, calculate the x- and y-deviations as the objectiv function
436  !-----------------------------------------------------------------------
437 
438  dev(i) = 0.0
439  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) dev(i) = (xi(1,2)-x_exp(i))**2
440  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) dev(i) = dev(i) + (xi(2,2)-y_exp(i))**2
441 
442  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) deviation(i,1) = xi(1,2)-x_exp(i)
443  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) deviation(i,2) = xi(2,2)-y_exp(i)
444 
445  ELSE
446 
447  !-----------------------------------------------------------------------
448  ! no solution found for (p,T) of point i. Now calculate delta(p)
449  !-----------------------------------------------------------------------
450 
451  dense(1) = 0.4
452  dense(2) = 1.e-5
453  xi(1,2) = x_exp(i)
454  xi(1,1) = 1.0 - xi(1,2)
455  xi(2,1:2) = xi(1,1:2)
456  t = t_exp(i)
457  p = p_exp(i)
458  iterate_t = 0
459  call bubble_point_rachford_rice ( iterate_t, converg, rhoi1, rhoi2 )
460  if ( converg == 1 ) then
461  CALL occupy_val_init ( rhoi1, rhoi2 )
462  val_conv = val_init
463  end if
464 
465  !-----------------------------------------------------------------------
466  ! calculate the contribution of point i to the objective function
467  !-----------------------------------------------------------------------
468 
469  if ( converg == 1 ) dev(i) = 5.0*( val_conv(4)/p_exp(i)-1.0 )**2
470  if ( converg == 1 ) report_character = 'only delta(p) converged'
471  if ( converg /= 1 ) write (*,'(a,i4,a,2G18.11)') ' point number ',i, &
472  ' did not converge, x,y=',x_exp(i),y_exp(i)
473  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) deviation(i,1) = 0.0 ! temporary
474  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) deviation(i,2) = 0.0 ! temporary
475 
476 
477  END IF
478  !write (*,'(a,i3,G20.12,2x,a,x,a)') 'deviation',i,dev(i), report_success_step, report_character
479  !read (*,*)
480 
481  END DO
482 
483  error = sum( dev(1:n_exp) )
484  fmin = error
485 
486  write(*,'(a,f14.8,a,3(f14.8))') 'deviation=',error,' kij=',optpars(1:nadjust) - 1.0
487 
488 END SUBROUTINE bindev_kij_p
489 
490 
491 
492 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
493 ! SUBROUTINE KIJ_FIT
494 !
495 ! routine for fitting kij to binary systems
496 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
497 
498 SUBROUTINE kij_fit
499 
500  USE basic_variables
501  use utilities, only: file_open, paus
502 
503  !-----------------------------------------------------------------------------
504  INTEGER :: i, j, k, opt
505  REAL :: dummy2, dummy3, dummy4
506  REAL :: kijorg, devges(5)
507  REAL :: dist, kijsav(5), devopt, dummy5, rms1, rms2
508  REAL :: aad1(400,5), aad2(400,5), aad_m1, aad_m2, maxerr(2,3)
509  REAL :: max_t, min_t, max_p, min_p, aver_y
510  REAL :: fmin, optpars(1)
511  CHARACTER :: fitin*50, fitin2*71, fitout*80, dummy1*30
512  LOGICAL :: go_up
513  !-----------------------------------------------------------------------------
514 
515  ncomp = 2
516  nphas = 2
517  n_unkw = ncomp
518  it(1) = 'x11' ! iteration of mol fraction of comp.1 phase 1
519  it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
520  sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
521  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
522 
523  dist = 0.01
524 
525  ! write (*,*) ' CHOOSE YOUR OPTIMIZATION STRATEGY :'
526  ! write (*,*) ' OPTIMIZE THE x AND y VALUES ...... CHOOSE: 1'
527  ! write (*,*) ' OPTIMIZE THE K1 AND K2 VALUES ...... CHOOSE: 2'
528  ! READ (*,*) c_option
529 
530  write (*,*) ' CHOOSE THE INPUT FILE (./input_file/kij-dat/...) :'
531  write (*,*) ' SPECIFY FULL NAME:'
532  READ (*,*) fitin
533  ! fitin='h2s-n-decane.dat'
534  ! fitin='methane-tetracontane.dat'
535  fitin2 = './input_file/kij-dat/'//fitin
536 
537  CALL file_open(fitin2,87)
538 
539  n_exp = -1
540  i = 0
541  READ (87,*) compna(1)
542  READ (87,*) compna(2)
543 
544  go_up = .true.
545  DO WHILE ( go_up )
546  READ (87,*) dummy1,dummy2,dummy3,dummy4,dummy5
547  n_exp = n_exp + 1
548  IF (dummy1 == 'end') THEN
549  go_up = .false.
550  ELSE
551  i = i + 1
552  x_exp(i) = dummy2
553  y_exp(i) = dummy3
554  p_exp(i) = dummy4 * 1.e5
555  t_exp(i) = dummy5 + 273.15
556  ENDIF
557  END DO
558  CLOSE (87)
559 
560  write (*,*) ' CHOOSE THE EOS :'
561  write (*,*) ' 0: SAFT EOS '
562  write (*,*) ' 1: PC-SAFT EOS '
563  write (*,*) ' 2: SRK EOS '
564  write (*,*) ' 3: PR EOS '
565  read (*,*) eos
566 
567  IF (eos == 1) THEN
568  write (*,*) ' FURTHER SPECIFY THE EOS :'
569  write (*,*) ' 0: PC-SAFT EOS '
570  write (*,*) ' 1: PCP-SAFT EOS '
571  write (*,*) ' 2: PCIP-SAFT EOS '
572  read (*,*) pol
573  IF(pol == 2 .AND. num == 0) write(*,*)'switch to num. derivatives'
574  IF(pol == 2 .AND. num == 0) num=1
575  ENDIF
576 
577  fitout = fitin
578  IF (eos == 1) fitout = './output_file/kij-pc-saft/'//fitout
579  IF (eos == 0) fitout = './output_file/kij-saft/'//fitout
580  IF (eos == 2) fitout = './output_file/kij-srk/'//fitout
581  IF (eos == 3) fitout = './output_file/kij-pr/'//fitout
582  OPEN (88,file=fitout)
583 
584  CALL para_input
585 
586 
587  kij_adjust: DO
588 
589  kijorg = kij(1,2)
590 
591  DO j = 1,5
592 
593  IF (j == 1) THEN
594  kij(1,2) = kijorg - 1.0 * dist
595  ELSEIF (j == 2) THEN
596  kij(1,2) = kijorg - 0.5 * dist
597  ELSEIF (j == 3) THEN
598  kij(1,2) = kijorg
599  ELSEIF (j == 4) THEN
600  kij(1,2) = kijorg + 0.5 * dist
601  ELSE
602  kij(1,2) = kijorg + 1.0 * dist
603  ENDIF
604  kijsav(j) = kij(1,2)
605  kij(2,1) = kij(1,2)
606  optpars(1) = kij(1,2) + 1.0
607  ! optpars(1) = kij(1,2)
608 
609  call bindev_kij ( fmin, optpars, 1 )
610  devges(j) = fmin
611  write (*,'(a,i3,2(f14.8))') ' deviation',j,devges(j),kij(1,2)
612 
613  DO i = 1, n_exp
614  aad1(i,j) = deviation(i,1)
615  aad2(i,j) = deviation(i,2)
616  ENDDO
617 
618  END DO
619 
620  devopt = devges(1)
621  kij(1,2) = kijsav(1)
622  kij(2,1) = kijsav(1)
623  opt = 1
624  DO k = 2,5
625  IF (devges(k) < devopt) THEN
626  kij(1,2) = kijsav(k)
627  kij(2,1) = kijsav(k)
628  devopt = devges(k)
629  opt = k
630  ENDIF
631  END DO
632 
633  IF ( opt /= 1 .AND. opt /= 5 ) dist = 0.25*dist
634 
635  write(*,'(a,f14.8,a,2(f14.8))')' new kij',kij(1,2),' +/-',dist,devopt
636 
637 
638 
639  IF (dist <= 0.00005) EXIT kij_adjust
640 
641  END DO kij_adjust
642 
643  rms1 = 0.0
644  rms2 = 0.0
645  aad_m1 = 0.0
646  aad_m2 = 0.0
647  maxerr(1,1) = 100.0 * aad1(1,opt)
648  maxerr(1,2) = t_exp(1)
649  maxerr(1,3) = p_exp(1)
650  maxerr(2,1) = 100.0 * aad2(1,opt)
651  maxerr(2,2) = t_exp(1)
652  maxerr(2,3) = p_exp(1)
653  DO i = 1, n_exp
654  write (88,'(i4,4(2x,G13.6))') i,x_exp(i),y_exp(i),p_exp(i)/1.e5,t_exp(i)-273.15
655  ENDDO
656  write (88,*) ' '
657  DO i = 1, n_exp
658  IF( (100.0 * abs(aad1(i,opt))) > maxerr(1,1) ) THEN
659  maxerr(1,1) = 100.0 * abs(aad1(i,opt))
660  maxerr(1,2) = t_exp(i)
661  maxerr(1,3) = p_exp(i)
662  ENDIF
663  IF( (100.0 * abs(aad2(i,opt))) > maxerr(2,1) ) THEN
664  maxerr(2,1) = 100.0 * abs(aad2(i,opt))
665  maxerr(2,2) = t_exp(i)
666  maxerr(2,3) = p_exp(i)
667  ENDIF
668  write (*,*) ' DEV% ph_1=', 100.0 * aad1(i,opt),' DEV% ph_2=',100.0 * aad2(i,opt)
669  write (88,*) ' DEVIATION% phase_1=',100.0 * aad1(i,opt), &
670  ' DEVIATION% phase_2=',100.0 * aad2(i,opt)
671  rms1 = rms1 + aad1(i,opt)**2
672  rms2 = rms2 + aad2(i,opt)**2
673  aad_m1 = aad_m1 + abs(aad1(i,opt))
674  aad_m2 = aad_m2 + abs(aad2(i,opt))
675  END DO
676  rms1 = 100.0 * ( rms1 / real( n_exp ) )**0.5
677  rms2 = 100.0 * ( rms2 / real( n_exp ) )**0.5
678  aad_m1 = 100.0 * aad_m1 / real( n_exp )
679  aad_m2 = 100.0 * aad_m2 / real( n_exp )
680 
681  write(*,*)' '
682  write(*,*)' MAX.DEV% ph_1=',maxerr(1,1),' AT',(maxerr(1,i),i=2,3)
683  write(*,*)' MAX.DEV% ph_2=',maxerr(2,1),' AT',(maxerr(2,i),i=2,3)
684  write(*,*)' '
685  write(*,*)' AAD% phase_1=',aad_m1
686  write(*,*)' AAD% phase_2=',aad_m2
687  write(*,*)' '
688  write(*,*)' RMS% phase_1=',rms1
689  write(*,*)' RMS% phase_2=',rms2
690 
691  write(88,*)' '
692  write(88,*)' MAX. DEVIATION% phase_1=',maxerr(1,1),' AT',(maxerr(1,i),i=2,3)
693  write(88,*)' MAX. DEVIATION% phase_2=',maxerr(2,1),' AT',(maxerr(2,i),i=2,3)
694  write(88,*)' '
695  write(88,*)' AAD% phase_1=',aad_m1
696  write(88,*)' AAD% phase_2=',aad_m2
697  write(88,*)' '
698  write(88,*)' RMS% phase_1=',rms1
699  write(88,*)' RMS% phase_2=',rms2
700 
701  write (*,*) ' '
702  write(*,*)' Result: kij=',kij(1,2)
703  write(88,*)' Result: kij=',kij(1,2)
704  write (*,*) ' '
705  write (88,*) ' '
706 
707  max_t=t_exp(1)
708  min_t=t_exp(1)
709  max_p=p_exp(1)
710  min_p=p_exp(1)
711  aver_y = y_exp(1) / real( n_exp )
712  DO i = 2, n_exp
713  IF (t_exp(i) > max_t) max_t = t_exp(i)
714  IF (t_exp(i) < min_t) min_t = t_exp(i)
715  IF (p_exp(i) > max_p) max_p = p_exp(i)
716  IF (p_exp(i) < min_p) min_p = p_exp(i)
717  aver_y = aver_y + y_exp(i) / real( n_exp )
718  END DO
719  write (88,'(a,1x,12(1x,e12.5),1x,a)') fitin,kij(1,2),aad_m1 &
720  ,aad_m2,rms1,rms2,maxerr(1,1),maxerr(2,1),min_t-273.15 &
721  ,max_t-273.15,min_p/1.e5,max_p/1.e5,aver_y
722  write (*,*) ' '
723 
724  CLOSE (88)
725 
726 END SUBROUTINE kij_fit
727 
728 
729 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
730 !
731 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
732 
733 SUBROUTINE bindev_kij ( fmin, optpars, nadjust )
734 
735  USE basic_variables
736  USE starting_values
737 
738  !-----------------------------------------------------------------------------
739  INTEGER, INTENT(IN) :: nadjust
740  REAL, INTENT(IN) :: optpars(:)
741  REAL, INTENT(IN OUT) :: fmin
742 
743  !-----------------------------------------------------------------------------
744  INTEGER :: i, k, ph, ph_split, converg
745  REAL :: error
746  REAL :: x_sta, y_sta
747  REAL, DIMENSION(400) :: dev
748  REAL :: end_x, steps
749  !-----------------------------------------------------------------------------
750 
751 
752  kij(1,2) = optpars(1) - 1.0
753  ! kij(1,2) = optpars(1)
754  kij(2,1) = kij(1,2)
755  IF (nadjust == 2) lij(1,2) = optpars(2)-1.0
756  ! IF (nadjust == 2) lij(1,2) = optpars(2)
757  IF (nadjust == 2) lij(2,1) = - lij(1,2)
758 
759  DO i = 1, n_exp
760 
761  t = t_exp(i)
762  p = p_exp(i)
763 
764  n_unkw = ncomp ! number of quantities to be iterated
765  it(1) = 'x11' ! iteration of mol fraction of comp.1 phase 1
766  it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
767  sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
768  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
769 
770  x_sta = x_exp(i)
771  y_sta = y_exp(i)
772  if ( x_sta == 0.0 ) x_sta = 0.01 ! these are only starting values
773  if ( x_sta == 1.0 ) x_sta = 0.99 ! these are only starting values
774  if ( y_sta == 0.0 ) y_sta = 0.01 ! these are only starting values
775  if ( y_sta == 1.0 ) y_sta = 0.99 ! these are only starting values
776 
777  lnx(1,2) = log( x_sta )
778  lnx(2,2) = log( y_sta )
779  lnx(1,1) = log( 1.0 - x_sta )
780  lnx(2,1) = log( 1.0 - y_sta )
781  val_init(1) = 0.4
782  val_init(2) = 0.001
783  val_init(3) = t
784  val_init(4) = p
785  DO ph = 1,nphas
786  DO k = 1,ncomp
787  val_init(4+k+(ph-1)*ncomp) = lnx(ph,k)
788  ENDDO
789  ENDDO
790 
791  CALL objective_ctrl (converg)
792 
793  IF (converg /= 1) THEN
794 
795  outp = 0 ! output to terminal
796  call vle_min
797  xif( 1:ncomp ) = xi( 1, 1:ncomp )
798  call start_var_fixed_composition ( converg )
799 
800  IF (converg /= 1) THEN
801  write (*,*) ' NO SOLUTION FOUND FOR REG. PRESSURE',i
802  flashcase = .true.
803  xif(2) = ( x_exp(i) + y_exp(i) ) / 2.0
804  xif(1) = 1.0 - xif(2)
805  p = p_exp(i) - 0.1*p_exp(i)
806  val_init(4) = p
807  ph_split = 0
808  call start_var_fixed_composition ( converg )
809  IF (converg == 1) THEN
810  steps = 10.0
811  end_x = p_exp(i)
812  running = 'p'
813  call phase_equilib (end_x,steps,converg)
814  END IF
815  IF (converg == 1) write (*,*) ' NOW A SOLUTION WAS FOUND ?, CHECK !',i,p_exp(i)
816  END IF
817  ! IF (converg /= 1) write (*,*) ' NO SOLUTION FOUND',i,p_exp(i)
818  ENDIF
819 
820 
821  dev(i) = 0.0
822  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) dev(i) = (xi(1,2)-x_exp(i))**2
823  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) dev(i) = dev(i) + (xi(2,2)-y_exp(i))**2
824 
825  IF ( x_exp(i) /= 0.0 .AND. x_exp(i) /= 1.0 ) deviation(i,1) = xi(1,2)-x_exp(i)
826  IF ( y_exp(i) /= 0.0 .AND. y_exp(i) /= 1.0 ) deviation(i,2) = xi(2,2)-y_exp(i)
827 
828  ! write (*,*) 'converg A',converg, xi(1,2)-x_exp(i), deviation(i,1)
829 
830  END DO
831  !call paus
832 
833  error = sum( dev(1:n_exp) )
834  fmin = error
835 
836  write(*,'(a,f14.8,a,3(f14.8))') 'deviation=',error,' kij=',optpars(1:nadjust) - 1.0
837 
838 END SUBROUTINE bindev_kij
839 
840 
841 end module kij_fitting
842 
starting_values
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module STARTING_VALUES This m...
Definition: modules.f90:251
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29