MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
starting_value.f90
Go to the documentation of this file.
1 
3 
4 
5 
6 
7 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
8 ! module STARTING_VALUES
9 !
10 ! This module contains parameters and variables for a phase stability
11 ! analyis as part of a flash calculation.
12 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
13 
14 module starting_values
15 
16  use parameters, only: nc
17  implicit none
18  save
19 
20  integer :: scan_index
21  real, dimension(nc), public :: rhoi_trial
22  real :: fdenf
23  logical :: first_call = .true.
24  logical, public :: generate_starting_val
25  logical, public :: flashcase
26 
27  real, dimension(nc) :: rhoi_best
28  real :: fmin_best
29 
30  PRIVATE
31  PUBLIC :: start_var, start_var_fixed_composition, phase_stability, vle_min, &
32  tangent_plane_line_search, rachford_rice, occupy_val_init, &
33  scan_compositions, bubble_point_rachford_rice, stability_hessian
34 
35 CONTAINS
36 
37 
38 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
62 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
63 
64 subroutine start_var (converg)
65 
66  use basic_variables
67  implicit none
68 
69  !-----------------------------------------------------------------------------
70  integer, intent(in out) :: converg
71 
72  !-----------------------------------------------------------------------------
73  logical :: renormalize
74  !-----------------------------------------------------------------------------
75 
76  converg = 0
77  nphas = 2
78  n_unkw = ncomp
79  outp = 0 ! output to terminal
80  generate_starting_val = .true.
81 
82  ensemble_flag = 'tp'
83 
84  if ( first_call .AND. xif(1) == 0.0 ) call read_mode_staring_value ( scan_index, outp )
85 
86  renormalize = .false. ! for renormalization group theory (RGT)
87  if (num == 2) renormalize = .true.
88  if (num == 2) num = 0 ! if RGT: initial phase equilibr. is for non-renormalized model
89 
90  if ( xif(1) /= 0.0 ) then
91 
92  flashcase = .true.
93 
94  call start_var_fixed_composition ( converg )
95 
96  else
97 
98  flashcase = .true.
99 
100  if ( ncomp > 2 ) then
101  write (*,*) 'Solubility Code: SR starting_value requires a defined feed composition for all but binary mixtures'
102  stop 5
103  end if
104 
105  if ( scan_index == 2 ) then
106  call vle_min
107  xif( 1:ncomp ) = xi( 1, 1:ncomp )
108  call start_var_fixed_composition ( converg )
109  else
110  call scan_compositions ( converg )
111  end if
112 
113  xif(:) = 0.0
114 
115  end if
116 
117  IF (renormalize) num = 2
118  generate_starting_val = .false.
119 
120 end subroutine start_var
121 
122 
123 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
124 ! subroutine scan_compositions
125 !
126 ! this SR initiates initiates flash calculations along a grid of compositions.
127 ! Only for binary mixtures and for given t, p.
128 !
129 ! input: t, p
130 ! output: converg
131 ! val_init ( !!! ) if a converged solution is found
132 !
133 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
134 
135 subroutine scan_compositions ( converg )
136 
137  use basic_variables
138  use utilities
139  implicit none
140 
141  !-----------------------------------------------------------------------------
142  integer, intent(in out) :: converg
143 
144  !-----------------------------------------------------------------------------
145  integer :: i, i_conv, steps
146  integer :: converg_overall
147  real :: x1_right
148  real :: diff_para
149  real :: xiF_in_between_x
150  real, dimension( 4, 0:nc*np+6 ) :: val_save
151  logical :: different_equilibrium_as_previous
152  real :: density(np),w(np,nc)
153  !-----------------------------------------------------------------------------
154 
155  steps = 100
156  x1_right = 0.0
157  converg_overall = 0
158  if ( ncomp > 2 ) write (*,*) 'SR scan_compositions only suited for binary systems'
159  if ( ncomp > 2 ) stop 5
160 
161  do i = 0, steps
162 
163  xif(1) = REAL(i) / REAL(steps)
164  if ( xif(1) <= 1.e-30 ) xif(1) = 1.e-30
165  if ( xif(1) >= (1.0 - 1.e-12) ) xif(1) = 1.0 - 1.e-12
166  xif(2) = 1.0 - xif(1)
167 
168  if ( xif(1) > x1_right ) then
169 
170  if ( outp >= 1 ) write (*,'(a,G19.11)') ' scan x, with x(1)=',xif(1)
171 
172  call start_var_fixed_composition ( converg )
173 
174  xif_in_between_x = ( exp(val_init(5)) - xif(1) ) / (xif(1) - exp(val_init(7)) )
175 
176  if ( converg == 1 .AND. xif_in_between_x > 0.0 ) then
177 
178  different_equilibrium_as_previous = .true.
179  if ( converg_overall > 0 ) then
180  diff_para = abs( exp(val_save(converg_overall,5))-exp(val_init(5)) ) &
181  + abs( exp(val_save(converg_overall,7))-exp(val_init(7)) )
182  if ( diff_para < 1.e-5 ) different_equilibrium_as_previous = .false.
183  end if
184 
185  if ( different_equilibrium_as_previous ) then
186  converg_overall = converg_overall + 1
187  val_save( converg_overall, : ) = val_init(:)
188  x1_right = max( xi(1,1), xi(2,1) )
189  if ( outp >= 1 ) write (*,'(a,2i4,5G19.11)') ' found equil.',i,converg_overall, &
190  exp(val_init(5) ),xif(1),exp(val_init(7) ),val_init(1:2)
191  if ( outp >= 1 ) write (*,*) ' '
192  ! call paus ('end of scan_compositions cycle')
193 
194  end if
195  end if
196 
197  end if
198 
199  end do
200 
201  if ( converg_overall > 0) converg = 1
202 
203  if ( outp >= 1 .OR. converg_overall > 1 ) then
204  if ( converg_overall >= 1) write (*,*) ' '
205  if ( converg_overall >= 1) write (*,*) '================================='
206  if ( converg_overall == 1) write (*,*) 'phase equilibrium:'
207  if ( converg_overall == 2) write (*,*) 'two phase equilibria found:'
208  if ( converg_overall >= 1) write (*,*) '================================='
209  do i_conv = 1, converg_overall
210  dense(1:2) = val_save( i_conv, 1:2 )
211  lnx(1,1:2) = val_save( i_conv, 5:6 )
212  lnx(2,1:2) = val_save( i_conv, 7:8 )
213  xi(1,1:2) = exp( lnx(1,1:2) )
214  xi(2,1:2) = exp( lnx(2,1:2) )
215  call si_dens ( density, w )
216  write (*,*) ' '
217  if ( converg_overall > 1) write (*,*) 'EQUILIBRIUM', i_conv
218  if ( converg_overall > 1) write (*,*) '---------------------------------------------'
219  write(*,'(t20,a,f7.2,a,f10.5,a)') 'T =',t-u_out_t,' P =',p/u_out_p
220  write (*,'(t26,a)') 'PHASE I PHASE II '
221  write (*,'(a,t20,F13.4,t39,F13.4)') ' Density ',density(1),density(2)
222  write (*,'(x,a,t20,a3,E14.6,t42,E14.6)') compna(1),'x =',exp( lnx(1:2,1) )
223  write (*,'(x,a,t20,a3,E14.6,t42,E14.6)') compna(2),'x =',exp( lnx(1:2,2) )
224  write (*,'(x,a,t20,a3,E14.6,t42,E14.6)') compna(1),'w =',w(1:2,1)
225  write (*,'(x,a,t20,a3,E14.6,t42,E14.6)') compna(2),'w =',w(1:2,2)
226  end do
227  end if
228 
229  val_init( : ) = val_save( 1, : )
230  if ( converg_overall > 1) then
231  write (*,*) ' '
232  write (*,*) 'choose a phase equilibrium'
233  write (*,*) '---------------------------------------------'
234  read (*,*) i_conv
235  val_init( : ) = val_save( i_conv, : )
236  end if
237 
238 end subroutine scan_compositions
239 
240 
241 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
242 ! subroutine calculate_equilibrium_feed
243 !
244 ! input: t, p
245 ! rhoi1, rhoi2 ( intent(in), starting values )
246 !
247 ! output: converg, val_init (meaningful, if converg=1)
248 !
249 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
250 
251 subroutine calculate_equilibrium_feed ( rhoi1, rhoi2, converg )
252 
253  use basic_variables
254  use eos_variables, only: pi, mseg, dhs
255  use utilities
256  implicit none
257 
258  !-----------------------------------------------------------------------------
259  integer, intent(in out) :: converg
260  real, dimension(nc), intent(in out) :: rhoi1
261  real, dimension(nc), intent(in out) :: rhoi2
262 
263  !-----------------------------------------------------------------------------
264  integer :: promising_min
265  real :: phi_2_start
266  real, dimension(np,nc) :: rho_enter
267  !-----------------------------------------------------------------------------
268 
269  rho_enter(1,1:ncomp) = rhoi1(1:ncomp)
270  rho_enter(2,1:ncomp) = rhoi2(1:ncomp)
271 
272  xi(1,1:ncomp) = rhoi1(1:ncomp) / sum(rhoi1(1:ncomp))
273  xi(2,1:ncomp) = rhoi2(1:ncomp) / sum(rhoi2(1:ncomp))
274  dense(1) = pi/6.0 * sum( rhoi1(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
275  dense(2) = pi/6.0 * sum( rhoi2(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
276 
277  CALL rachford_rice ( converg, rhoi1, rhoi2 )
278 
279  if ( converg == 0 ) then
280 
281  rhoi1(1:ncomp) = rho_enter(1,1:ncomp)
282  rhoi2(1:ncomp) = rho_enter(2,1:ncomp)
283  call tangent_plane_line_search ( rhoi1, rhoi2, phi_2_start, promising_min )
284 
285  if ( promising_min == 1 .AND. maxval( abs( xi(1,1:ncomp) - xi(2,1:ncomp) ) ) > 0.0001 ) then
286 
287  call tangent_plane_2 ( phi_2_start, rhoi1, rhoi2 )
288  xi( 1, 1:ncomp ) = rhoi1( 1:ncomp ) / sum( rhoi1( 1:ncomp ) )
289  xi( 2, 1:ncomp ) = rhoi2( 1:ncomp ) / sum( rhoi2( 1:ncomp ) )
290  if ( maxval( abs( xi( 1, 1:ncomp ) - xi( 2, 1:ncomp ) ) ) > 0.0001 ) then
291 
292  CALL occupy_val_init ( rhoi1, rhoi2 )
293  if ( .NOT. flashcase ) then
294  CALL select_sum_rel (1,0,1)
295  CALL select_sum_rel (2,0,2)
296  else
297  CALL determine_flash_it2
298  end if
299 
300  CALL objective_ctrl ( converg )
301 
302  if ( converg == 1 ) val_init = val_conv
303 
304  end if
305 
306  end if
307 
308  else
309 
310  CALL occupy_val_init ( rhoi1, rhoi2 )
311  val_conv = val_init
312 
313  end if
314 
315  if ( outp >= 2 ) write (*,*) ' '
316  if ( outp >= 2 ) write (*,*) 'leaving calculate_equilibrium_feed: convergence=',converg
317 
318 end subroutine calculate_equilibrium_feed
319 
320 
321 
322 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
323 ! subroutine start_var_fixed_composition
324 !
325 ! This subroutine performs a flash calculation for a mixtue with defined compo-
326 ! sition. The density is not a priori defined.
327 !
328 ! input: t, p, xiF(:)
329 ! output: converg
330 ! val_init ( !!! ) if a converged solution is found
331 !
332 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
333 
334 subroutine start_var_fixed_composition ( converg )
335 
336  use basic_variables
337  use eos_variables, only: pi, kbol
338  use utilities
339  implicit none
340 
341  !-----------------------------------------------------------------------------
342  integer, intent(in out) :: converg
343 
344  !-----------------------------------------------------------------------------
345  integer :: ph_split
346  real, dimension(np,nc) :: lnphi
347  real, dimension(nc) :: rhoi_feed
348  real, dimension(nc) :: rhoi1, rhoi2
349  real, dimension(nc) :: rhoi_trial
350  real :: eta_trial
351  real :: phi_1, fmin
352  !-----------------------------------------------------------------------------
353 
354  if ( sum( xif(1:ncomp) ) /= 1.0 ) write (*,*) 'Solubility Code: feed composition /= 1.0', sum( xif(1:ncomp) )
355  if ( sum( xif(1:ncomp) ) /= 1.0 ) stop 5
356 
357  xi(1,1:ncomp) = xif(1:ncomp)
358  xi(2,1:ncomp) = xif(1:ncomp)
359 
360  !-----------------------------------------------------------------------------
361  ! start with a Rachford-Rice iteration that searches for VLE. The attempt is
362  ! naive, because standard starting values for dense(1), dense(2), and
363  ! xi(2,:)=xiF(:) are used.
364  !-----------------------------------------------------------------------------
365 
366  dense(1) = 0.4
367  dense(2) = 1.e-6
368  CALL rachford_rice ( converg, rhoi1, rhoi2 )
369 
370  !-----------------------------------------------------------------------------
371  ! distinguish: case (1), RR did not converge. Continue with a conventional
372  ! stability analysis
373  ! case (2), RR converged. Then we have to verify that the liquid
374  ! phase is stable.
375  !-----------------------------------------------------------------------------
376 
377  if ( converg /= 1 ) then
378 
379  !--------------------------------------------------------------------------
380  ! No solution was found with the Rachford-Rice procedure. Now do a
381  ! traditional stability test first. If two densities are possible at
382  ! xiF(1:ncomp), then do the phase stability for both "feed" conditions.
383  !--------------------------------------------------------------------------
384  if ( outp >= 2 ) write (*,*) ' '
385  if ( outp >= 2 ) write (*,*) '--------------------------------------------------------'
386  if ( outp >= 1 ) write (*,*) 'do phase stability analysis'
387  if ( outp >= 2 ) write (*,*) '--------------------------------------------------------'
388  if ( outp >= 2 ) write (*,*) ' '
389 
390  densta(1) = 0.4
391  densta(2) = 1.e-6
392  xi( 1, 1:ncomp ) = xif( 1:ncomp )
393  xi( 2, 1:ncomp ) = xif( 1:ncomp )
394  CALL fugacity (lnphi)
395 
396  !--------------------------------------------------------------------------
397  ! phase stability based on liquid "feed"
398  !--------------------------------------------------------------------------
399  my_f( 1:ncomp ) = my_cal(1,1:ncomp)
400  rhoi_feed( 1:ncomp ) = rhoi_cal( 1, 1:ncomp )
401 
402  eta_trial = 0.4
403  ! to do: check wheather eta_trial = dense(1), or (2) is better choice
404  CALL phase_stability ( rhoi_feed, eta_trial, ph_split, rhoi_trial )
405  rhoi1( 1:ncomp ) = rhoi_feed( 1:ncomp )
406  rhoi2( 1:ncomp ) = rhoi_trial( 1:ncomp )
407  if ( outp > 0 ) write (*,*) 'phase split =', ph_split
408  if ( ph_split == 0 .AND. outp > 0 ) write (*,*) 'no phase split detected, stable phase is likely'
409  if ( ph_split == 2 .AND. outp > 0 ) write (*,*) 'condition could be close to critical point'
410 
411  if ( ph_split >= 1 ) then
412 
413  !-----------------------------------------------------------------------
414  ! If phase_stability gave ph_split = 1 or 2: do tangent plane iteration
415  !-----------------------------------------------------------------------
416 
417  call calculate_equilibrium_feed ( rhoi1, rhoi2, converg )
418 
419  if ( converg == 1 .AND. outp > 0 ) write (*,*) 'end of starting value - solution found'
420 
421  end if
422 
423  if ( converg == 0 ) then
424 
425  !-----------------------------------------------------------------------
426  ! test phase stability based on vapor "feed"
427  !-----------------------------------------------------------------------
428 !!$ densta(1) = 0.4
429 !!$ densta(2) = 1.E-6
430 !!$ xi( 1, 1:ncomp ) = xiF( 1:ncomp )
431 !!$ xi( 2, 1:ncomp ) = xiF( 1:ncomp )
432 !!$ CALL fugacity (lnphi)
433 !!$
434 !!$ my_f( 1:ncomp ) = my_cal( 2, 1:ncomp )
435 !!$ eta_trial = dense(2)
436 !!$ rhoi_feed( 1:ncomp ) = rhoi_cal( 2, 1:ncomp )
437 !!$
438 !!$ CALL phase_stability ( rhoi_feed, eta_trial, ph_split, rhoi_trial )
439 !!$ rhoi1( 1:ncomp ) = rhoi_feed( 1:ncomp )
440 !!$ rhoi2( 1:ncomp ) = rhoi_trial( 1:ncomp )
441 !!$ if ( outp > 0 ) write (*,*) 'phase split =', ph_split
442 !!$ if ( ph_split == 0 .AND. outp > 0 ) write (*,*) 'no phase split', &
443 !!$ ' detected, most probably a stable phase'
444 !!$ if ( ph_split == 2 .AND. outp > 0 ) write (*,*) 'condition could be close to critical point'
445 !!$
446 !!$ if ( ph_split >= 1 ) then
447 !!$
448 !!$ call calculate_equilibrium_feed ( rhoi1, rhoi2, converg )
449 !!$
450 !!$ end if
451 !!$
452 !!$ if ( converg == 0 .AND. ph_split == 1 .AND. outp > 0 ) &
453 !!$ write (*,*) 'vapor phase is unstable, but phase split did not converge'
454  end if
455 
456  else
457 
458  !--------------------------------------------------------------------------
459  ! RR was successful. VLE was calculated. Now, test stability of liquid phase
460  !--------------------------------------------------------------------------
461  if ( outp >= 2 ) write (*,*) ' now test phase stability'
462  if ( outp >= 2 ) write (*,*) ' '
463 
464  my_f( 1:ncomp ) = my_cal(1,1:ncomp)
465  eta_trial = 0.45
466  rhoi_feed( 1:ncomp ) = rhoi1( 1:ncomp )
467  if ( dense(2) > dense(1) ) rhoi_feed( 1:ncomp ) = rhoi2( 1:ncomp )
468 
469  CALL phase_stability ( rhoi_feed, eta_trial, ph_split, rhoi_trial )
470  if ( outp > 0 ) write (*,*) 'result of stability-test for L-phase, ph_split =', ph_split
471 
472  if ( ph_split /= 1 ) then
473 
474  !-----------------------------------------------------------------------
475  ! for the calculated VLE: liquid phase is stable
476  !-----------------------------------------------------------------------
477  if ( outp >= 2 ) write (*,*) 'the calculated VLE is stable'
478  CALL occupy_val_init ( rhoi1, rhoi2 )
479 
480  else
481 
482  !-----------------------------------------------------------------------
483  ! for the calculated VLE: liquid phase is not stable.
484  ! Determine stable equilibrium. Use the density vector rhoi_trial (de-
485  ! termined from phase_stability) as a starting value. The starting value
486  ! of second phase is determined from the component balance (using the
487  ! L-phase and afterwards the V-phase of the previously determined VLE)
488  !-----------------------------------------------------------------------
489  if ( outp >= 2 ) write (*,*) '========================================================'
490  if ( outp > 0 ) write (*,*) ' converged VLE (or LLE) not stable: searching (other) LLE'
491  if ( outp >= 2 ) write (*,*) '========================================================'
492  if ( outp >= 2 ) write (*,*) ' '
493 
494  !-----------------------------------------------------------------------
495  ! temporarily save the converged VLE solution into val_init
496  !-----------------------------------------------------------------------
497  call occupy_val_init ( rhoi1, rhoi2 )
498 
499  rhoi1( 1:ncomp ) = rhoi_feed( 1:ncomp )
500  rhoi2( 1:ncomp ) = rhoi_trial( 1:ncomp )
501  ! I am not sure the following line (and the condition 0<phi_1<1) is helpful. Why did I introduce this?
502  phi_1 = ( xif(1) - rhoi2(1)/sum(rhoi2(1:ncomp)) ) &
503  / ( rhoi1(1)/sum(rhoi1(1:ncomp)) - rhoi2(1)/sum(rhoi2(1:ncomp)) )
504 
505  if ( outp > 0 ) write (*,*) 'initial phase fraction', phi_1, kij(1,2)
506 
507  !if ( phi_1 > 0.0 .AND. phi_1 < 1.0 ) then
508 
509  call calculate_equilibrium_feed ( rhoi1, rhoi2, converg )
510 
511  if ( converg == 1 ) then
512  phi_1 = ( xif(1) - xi(2,1) ) / ( xi(1,1) - xi(2,1) )
513  fmin = phi_1 * gibbs(1) + ( 1.0 - phi_1 ) * gibbs(2)
514  if ( outp >= 2 ) write (*,*) ' '
515  if ( outp >= 2 ) write (*,*) '--------------------------------------------------------'
516  if ( outp > 0 ) write (*,*) 'first phase equilibrium trial converged ',fmin
517  if ( outp >= 2 ) write (*,*) '--------------------------------------------------------'
518  if ( outp >= 2 ) write (*,*) ' '
519  else
520  if ( outp >= 2 ) write (*,*) 'search for additional phase split gave no solution'
521  end if
522  !else
523  ! if ( outp >= 2 ) write (*,*) 'first phase equilibrium trial did not qualify'
524  !end if
525 
526 !!$ !rhoi_test1( 1:ncomp ) = rhoi_trial( 1:ncomp )
527 !!$ !rhoi_test2( 1:ncomp ) = rhoi_save( 1:ncomp )
528 !!$ !phi_1 = ( xif(1) - rhoi_test2(1)/sum(rhoi_test2(1:ncomp)) ) &
529 !!$ ! / ( rhoi_test1(1)/sum(rhoi_test1(1:ncomp)) - rhoi_test2(1)/sum(rhoi_test2(1:ncomp)) )
530 !!$ if( outp > 0 ) write (*,*) 'initial phase fraction', phi_1
531 !!$
532 !!$ if ( phi_1 > 0.0 .AND. phi_1 < 1.0 ) then
533 !!$
534 !!$ !rhoi1( 1:ncomp ) = rhoi_test1( 1:ncomp )
535 !!$ !rhoi2( 1:ncomp ) = rhoi_test2( 1:ncomp )
536 !!$ call calculate_equilibrium_feed ( rhoi1, rhoi2, converg ) ! rhoi_test1, rhoi_test2
537 !!$
538 !!$ converg_all = converg_all + converg
539 !!$
540 !!$ if ( converg == 1 ) then
541 !!$ phi_1 = ( xif(1) - xi(2,1) ) / ( xi(1,1) - xi(2,1) )
542 !!$ fmin = phi_1 * gibbs(1) + ( 1.0 - phi_1 ) * gibbs(2)
543 !!$ if ( outp >= 2 ) write (*,*) ' '
544 !!$ if ( outp >= 2 ) write (*,*) '--------------------------------------------------------'
545 !!$ if ( outp > 0 ) write (*,*) 'second phase equilibrium trial converged ',fmin
546 !!$ if ( outp >= 2 ) write (*,*) '--------------------------------------------------------'
547 !!$ if ( outp >= 2 ) write (*,*) ' '
548 !!$ else
549 !!$ if ( outp >= 2 ) write (*,*) 'second phase equilibrium trial gave no solution'
550 !!$ end if
551 !!$ else
552 !!$ if ( outp >= 2 ) write (*,*) 'second phase equilibrium trial did not qualify'
553 !!$ end if
554 
555  if ( converg == 0 ) then ! this is done to recover the VLE solution of Rachford-Rice
556  !--------------------------------------------------------------------
557  ! the search for additional phase split did not deliver stable phases.
558  ! Restore the initially converged results -- no particular action
559  ! needed because the results are already stored in "val_init".
560  !--------------------------------------------------------------------
561  converg = 1
562  end if
563 
564  end if
565 
566  end if
567  if ( outp > 0 ) write (*,*) ' '
568 
569 end subroutine start_var_fixed_composition
570 
571 
572 
573 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
574 !
575 ! input: rhoi1, rhoi2
576 !
577 ! output: rhoi1, rhoi2, xi(1,:), xi(2,:)
578 !
579 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
580 
581 subroutine tangent_plane_line_search ( rhoi1, rhoi2, phi_2_opt, promising_min )
582 
583  use basic_variables
584  use eos_variables, only: pi, mseg, dhs
585  implicit none
586 
587  !-----------------------------------------------------------------------------
588  real, dimension(nc), intent(in out) :: rhoi1
589  real, dimension(nc), intent(in out) :: rhoi2
590  real, intent(out) :: phi_2_opt
591  integer, intent(out) :: promising_min
592 
593  !-----------------------------------------------------------------------------
594  integer :: n, i, steps
595  real :: fmin
596  real :: phi_2, phi_2_min, phi_2_max, fmin_scan
597  real, allocatable :: optpara(:)
598  !-----------------------------------------------------------------------------
599 
600  phi_2_opt = 0.5
601  promising_min = 0
602 
603  n = ncomp
604  ALLOCATE( optpara(n) )
605 
606  xi( 1, 1:ncomp ) = rhoi1( 1:ncomp ) / sum( rhoi1( 1:ncomp ) )
607  xi( 2, 1:ncomp ) = rhoi2( 1:ncomp ) / sum( rhoi2( 1:ncomp ) )
608  lnx( 1, 1:ncomp ) = log( xi( 1, 1:ncomp ) )
609  lnx( 2, 1:ncomp ) = log( xi( 2, 1:ncomp ) )
610 
611  densta(1) = pi / 6.0 * sum( rhoi1( 1:ncomp ) * mseg( 1:ncomp ) * dhs( 1:ncomp )**3 )
612  densta(2) = pi / 6.0 * sum( rhoi2( 1:ncomp ) * mseg( 1:ncomp ) * dhs( 1:ncomp )**3 )
613  !pause
614 
615  phi_2_min = 0.001
616  phi_2_max = 1.0
617  do i = 1, ncomp
618  phi_2_max = min( phi_2_max, xif(i) / xi(2,i) )
619  end do
620 
621  steps = 50
622  fmin_scan = 1.e10
623 
624  do i = 0, steps
625 
626  phi_2 = phi_2_min + ( phi_2_max - phi_2_min ) * real(i) / real(steps)
627  if ( phi_2 < 1.e-6 ) phi_2 = 1.e-6
628  if ( phi_2 > (1.0 - 1.e-6) ) phi_2 = 1.0 - 1.e-6
629 
630  optpara( 1:ncomp ) = log( xi( 2, 1:ncomp ) * phi_2 )
631 
632  call tangent_value ( fmin, optpara, n )
633 
634  if ( fmin < fmin_scan ) then
635  fmin_scan = fmin
636  phi_2_opt = phi_2
637  rhoi1( 1:ncomp ) = rhoi_cal( 1, 1:ncomp )
638  rhoi2( 1:ncomp ) = rhoi_cal( 2, 1:ncomp )
639  if ( i > 0 ) promising_min = 1
640  end if
641 
642  end do
643 
644  DEALLOCATE( optpara )
645 
646  xi( 1, 1:ncomp ) = rhoi1( 1:ncomp ) / sum( rhoi1( 1:ncomp ) )
647  xi( 2, 1:ncomp ) = rhoi2( 1:ncomp ) / sum( rhoi2( 1:ncomp ) )
648 
649  if ( outp >= 1 ) write (*,*) 'finished line search'
650  if ( outp >= 2 ) write (*,*) ' '
651 
652 end subroutine tangent_plane_line_search
653 
654 
655 
656 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
657 !
658 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
659 
660 subroutine tangent_value ( fmin, optpara, n )
661 
662  use basic_variables, only: ncomp, np, outp, gibbs, xif, xi, lnx, my_f
663  implicit none
664 
665  !-----------------------------------------------------------------------------
666  integer, intent(IN) :: n
667  real, intent(IN) :: optpara(:)
668  real, intent(IN OUT) :: fmin
669 
670  !-----------------------------------------------------------------------------
671  integer :: i
672  real :: lnphi(np,nc)
673  real :: ph_1_frac, punish
674  real, dimension(nc) :: ni_1, ni_2
675  !-----------------------------------------------------------------------------
676 
677  !-----------------------------------------------------------------------------
678  ! saveguard that the entering parameters are reasonable
679  !-----------------------------------------------------------------------------
680  if ( maxval( optpara(1:n) ) > 5.0 ) then
681  fmin = 1.e8 * maxval( optpara(1:n) )
682  return
683  end if
684  if ( maxval( optpara(1:n) ) < -100.0 ) then
685  fmin = 1.e8
686  return
687  end if
688  do i = 1, n
689  if ( optpara(i) /= optpara(i) ) then
690  fmin = 1.e8
691  return
692  end if
693  end do
694 
695  !-----------------------------------------------------------------------------
696  ! setting of mole numbers of all species in second phase
697  !-----------------------------------------------------------------------------
698 
699  do i = 1, ncomp
700  if ( optpara(i) < -100.0 ) then
701  ni_2(i) = 0.0
702  else
703  ni_2(i) = exp( optpara(i) )
704  end if
705  end do
706 
707  punish = 0.0
708  do i = 1, ncomp
709  ni_1(i) = xif(i) - ni_2(i)
710  if ( ni_2(i) > xif(i) ) then
711  punish = punish - 1.e8*ni_1(i) + 0.1
712  !ni_2(i) = xif(i)
713  ni_1(i) = xif(i) * 1.e-50
714  end if
715  if ( ni_1(i) <= 0.0 ) ni_1(i) = xif(i) * 1.e-50
716  end do
717 
718 
719  !-----------------------------------------------------------------------------
720  ! calculate mole fractions of second and first phase (from component balance)
721  !-----------------------------------------------------------------------------
722 
723  xi(2,1:ncomp) = ni_2(1:ncomp) / sum( ni_2(1:ncomp) )
724  lnx(2,1:ncomp) = optpara(1:ncomp) - log( sum( ni_2(1:ncomp) ) )
725 
726  ph_1_frac = sum( ni_1(1:ncomp) )
727  xi( 1, 1:ncomp ) = ni_1( 1:ncomp ) / ph_1_frac
728  lnx( 1, 1:ncomp ) = log( ni_1( 1:ncomp ) ) - log( ph_1_frac )
729 
730  !-----------------------------------------------------------------------------
731  ! calculate total gibbs energy (for given T.p, and x(2,:) of trial phase)
732  !-----------------------------------------------------------------------------
733 
734  ! write (*,*) 'f_TP x', xi( 1, 1:ncomp )
735  ! write (*,*) 'f_TP x', xi( 2, 1:ncomp )
736  CALL fugacity (lnphi)
737 
738  fmin = gibbs(1) * ph_1_frac + gibbs(2) * ( 1.0 - ph_1_frac ) + punish
739  ! write (*,*) 'f_TP f=',gibbs(1) , gibbs(2)
740  fmin = fmin - sum( xif(1:ncomp)*my_f( 1:ncomp ))
741  if ( outp >= 2 ) write (*,'(a,4G18.8)') 'TP',fmin, ni_1(1:ncomp) !,xi(1,1:ncomp)
742  ! write (*,'(a,4G18.8)') 'al', ph_1_frac, (xi(2,i), i=1,ncomp)
743  ! write (*,*) ' '
744  ! pause
745 
746 end subroutine tangent_value
747 
748 
749 
750 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
751 !
752 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
753 
754 subroutine tangent_grad ( g, optpara, n )
755 
756  use basic_variables
757  use eos_variables, only: kbol
758  implicit none
759 
760  !-----------------------------------------------------------------------------
761  integer, intent(IN) :: n
762  real, intent(IN) :: optpara(:)
763  real, intent(in out) :: g(:)
764 
765  !-----------------------------------------------------------------------------
766  integer :: i
767  real :: lnphi(np,nc), ph_1_frac
768  real, dimension(nc) :: ni_1, ni_2
769  !-----------------------------------------------------------------------------
770 
771 
772  !-----------------------------------------------------------------------------
773  ! saveguard that the entering parameters are reasonable
774  !-----------------------------------------------------------------------------
775  if ( maxval( optpara(1:n) ) < -100.0 ) then
776  g(:) = 0.0
777  return
778  end if
779  do i = 1, n
780  if ( optpara(i) /= optpara(i) .OR. optpara(i) > 5.0 ) then
781  g(:) = 0.0
782  return
783  end if
784  end do
785 
786  !-----------------------------------------------------------------------------
787  ! setting of mole numbers of all species in second phase
788  !-----------------------------------------------------------------------------
789 
790  do i = 1, n
791  if ( optpara(i) < -300.0 ) then
792  ni_2(i) = 0.0
793  else
794  ni_2(i) = exp( optpara(i) )
795  end if
796  end do
797 
798  do i = 1, ncomp
799  ni_1(i) = xif(i) - ni_2(i)
800  if ( ni_2(i) > xif(i) ) then
801  !write (*,*) 'now active bound ================================'
802  !ni_2(i) = xif(i)
803  ni_1(i) = xif(i) * 1.e-100
804  end if
805  end do
806 
807  !-----------------------------------------------------------------------------
808  ! calculate mole fractions of second and first phase (from component balance)
809  !-----------------------------------------------------------------------------
810 
811  xi(2,1:ncomp) = ni_2(1:ncomp) / sum( ni_2(1:ncomp) )
812  lnx(2,1:ncomp) = optpara(1:ncomp) - log( sum( ni_2(1:ncomp) ) )
813 
814  ph_1_frac = sum( ni_1(1:ncomp) )
815  xi( 1, 1:ncomp ) = ni_1( 1:ncomp ) / ph_1_frac
816  lnx( 1, 1:ncomp ) = log( ni_1( 1:ncomp ) ) - log( ph_1_frac )
817 
818  !-----------------------------------------------------------------------------
819  ! calculate chemical potentials (for given T.p, and x(2,:) of trial phase)
820  !-----------------------------------------------------------------------------
821 
822  CALL fugacity (lnphi)
823 
824  !g( 1:ncomp ) = - ( lnphi(1,1:ncomp) + log( rhoi_cal(1,1:ncomp) ) &
825  ! + LOG(p/(kbol*t*1.E30*sum(rhoi_cal(1,1:ncomp)))) ) &
826  ! + ( lnphi(2,1:ncomp) + log( rhoi_cal(2,1:ncomp) ) &
827  ! + LOG(p/(kbol*t*1.E30*sum(rhoi_cal(2,1:ncomp)))) )
828  !g( 1:ncomp ) = g( 1:ncomp ) * EXP( optpara(1:ncomp) )
829 
830  g( 1:ncomp ) = - ( lnphi(1,1:ncomp) + log( xi(1,1:ncomp) ) ) &
831  + ( lnphi(2,1:ncomp) + log( xi(2,1:ncomp) ) )
832  g( 1:ncomp ) = g( 1:ncomp ) * exp( optpara(1:ncomp) )
833 
834  ! write (*,'(a,4G18.8)') 'g', g(1:ncomp)
835 
836 end subroutine tangent_grad
837 
838 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
839 !
840 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
841 
842 subroutine tangent_hessian (hessian, gtrans, optpara, n)
843 
844  use basic_variables
845  implicit none
846 
847  integer, intent(IN) :: n
848  real, intent(IN) :: optpara(:)
849  real, intent(IN OUT) :: gtrans(:)
850  real, intent(IN OUT) :: hessian(:,:)
851 
852  !-----------------------------------------------------------------------------
853  integer :: i, j
854  real :: delta
855  real :: optpara_mod(n), g(n), gi(n,n), gi_left(n,n)
856  !-----------------------------------------------------------------------------
857 
858 
859  delta = 1.0e-4
860 
861  optpara_mod = optpara
862 
863  do i = 1, n
864 
865  optpara_mod = optpara
866  optpara_mod(i) = optpara(i)*(1.0+delta)
867 
868  call tangent_grad (g, optpara_mod, n)
869  gi(i,1:n) = g(1:n)
870 
871  optpara_mod(i) = optpara(i)*(1.0-delta)
872 
873  call tangent_grad (g, optpara_mod, n)
874  gi_left(i,1:n) = g(1:n)
875 
876  end do
877 
878  call tangent_grad (g, optpara, n)
879 
880 
881  do i = 1, n
882  do j = 1, n
883 
884  hessian(i,j) = ( gi(i,j) - gi_left(i,j) ) / ( 2.0*optpara(i)*delta )
885  ! hessian(j,i) = hessian(i,j)
886  ! write (*,*) i,j,hessian(i,j)
887 
888  end do
889  end do
890 
891  gtrans = g
892  ! call tangent_value (fmin, optpara, n)
893  ! write (*,'(4G20.12)') g(1:n),xi(1,1),xi(2,1)
894 
895 end subroutine tangent_hessian
896 
897 
898 
899 
900 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
901 !
902 ! input: xi(1,:), xi(2,:), densta(1:2)
903 !
904 ! output: rhoi1, rhoi2, xi(1,:), xi(2,:), dense(:)
905 !
906 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
907 
908 subroutine tangent_plane_2 ( phi_2_start, rhoi1, rhoi2 )
909 
910  use basic_variables
911  use cg_minimization
912  implicit none
913 
914  !-----------------------------------------------------------------------------
915  real, intent(in) :: phi_2_start
916  real, dimension(nc), intent(out) :: rhoi1
917  real, dimension(nc), intent(out) :: rhoi2
918 
919  !-----------------------------------------------------------------------------
920  integer :: n
921  integer :: small_i, min_ph, other_ph
922  real :: lnphi(np,nc)
923  real, allocatable :: optpara(:)
924  integer :: STATUS, iter, nfunc, ngrad
925  real :: gnorm, fmin
926  real, allocatable :: d(:), g(:), xtemp(:), gtemp(:)
927  !-----------------------------------------------------------------------------
928 
929  n = ncomp
930 
931  ALLOCATE( optpara(n) )
932  ALLOCATE( d(n) )
933  ALLOCATE( g(n) )
934  ALLOCATE( xtemp(n) )
935  ALLOCATE( gtemp(n) )
936 
937  lnx( 1, 1:ncomp ) = log( xi( 1, 1:ncomp ) )
938  lnx( 2, 1:ncomp ) = log( xi( 2, 1:ncomp ) )
939 
940  optpara( 1:ncomp ) = log( xi( 2, 1:ncomp ) * phi_2_start )
941 
942  !if ( n == 2 ) then
943  ! CALL Newton_Opt_2D ( tangent_hessian, tangent_value, optpara, n, 1.E-8, 1.E-8, gtemp, fmin )
944  !else
945  CALL cg_descent ( 1.e-6, optpara, n, tangent_value, tangent_grad, status, &
946  gnorm, fmin, iter, nfunc, ngrad, d, g, xtemp, gtemp )
947  !end if
948 
949  !-----------------------------------------------------------------------------
950  ! If one component is a polymer (indicated by a low component-density)
951  ! then get an estimate of the polymer-lean composition, by solving for
952  ! xi_p1 = ( xi_p2 * phii_p2) / phii_p1 (phase equilibrium condition,
953  ! with p1 for phase 1)
954  !-----------------------------------------------------------------------------
955  if ( minval( lnx( 1, 1:ncomp ) ) < minval( lnx( 2, 1:ncomp ) ) ) then
956  min_ph = 1
957  other_ph = 2
958  else
959  min_ph = 2
960  other_ph = 1
961  end if
962  small_i = minloc( lnx(min_ph,1:ncomp), 1 )
963  !---- if one component is a polymer ------------------------------------------
964  if ( minval( lnx(min_ph,1:ncomp) ) < -20.0 ) then
965  CALL fugacity ( lnphi )
966  lnx(min_ph,small_i) = lnx(other_ph,small_i)+lnphi(other_ph,small_i) - lnphi(min_ph,small_i)
967  xi(min_ph,1:ncomp) = exp( lnx(min_ph,1:ncomp) ) / sum( exp( lnx(min_ph,1:ncomp) ) )
968  call fugacity ( lnphi )
969  end if
970 
971  rhoi1( 1:ncomp ) = rhoi_cal( 1, 1:ncomp )
972  rhoi2( 1:ncomp ) = rhoi_cal( 2, 1:ncomp )
973 
974  DEALLOCATE( optpara, d, g, xtemp, gtemp )
975 
976 
977 end subroutine tangent_plane_2
978 
979 
980 
981 
982 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
983 !
984 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
985 
986 subroutine tangent_plane ( phi_2_start, rhoi1, rhoi2 )
987 
988  use basic_variables
989  implicit none
990 
991  !-----------------------------------------------------------------------------
992  real, intent(in) :: phi_2_start
993  real, dimension(nc), intent(out) :: rhoi1
994  real, dimension(nc), intent(out) :: rhoi2
995 
996  !-----------------------------------------------------------------------------
997  integer :: n
998  integer :: small_i, min_ph, other_ph
999  integer :: PRIN
1000  real :: fmin , t0, h0, MACHEP
1001  real :: lnphi(np,nc)
1002  real, allocatable :: optpara(:)
1003  !-----------------------------------------------------------------------------
1004 
1005  n = ncomp
1006  t0 = 1.e-4
1007  h0 = 0.1
1008  prin = 0
1009  machep = 1.e-15
1010 
1011  ALLOCATE( optpara(n) )
1012 
1013  lnx( 1, 1:ncomp ) = log( xi( 1, 1:ncomp ) )
1014  lnx( 2, 1:ncomp ) = log( xi( 2, 1:ncomp ) )
1015 
1016  optpara( 1:ncomp ) = log( xi( 2, 1:ncomp ) * phi_2_start )
1017 
1018 
1019  call praxis( t0, machep, h0, n, prin, optpara, tangent_value, fmin )
1020 
1021  ! The optimal optpara-vector is not necessarily the one that was last evaluated.
1022  ! tangent_value is reexecuted with the optimal vector optpara, in order to update the ln(x) values
1023  call tangent_value( fmin, optpara, n )
1024 
1025 
1026  !-----------------------------------------------------------------------------
1027  ! If one component is a polymer (indicated by a low component-density)
1028  ! then get an estimate of the polymer-lean composition, by solving for
1029  ! xi_p1 = ( xi_p2 * phii_p2) / phii_p1 (phase equilibrium condition,
1030  ! with p1 for phase 1)
1031  !-----------------------------------------------------------------------------
1032  if ( minval( lnx( 1, 1:ncomp ) ) < minval( lnx( 2, 1:ncomp ) ) ) then
1033  min_ph = 1
1034  other_ph = 2
1035  else
1036  min_ph = 2
1037  other_ph = 1
1038  end if
1039  small_i = minloc( lnx(min_ph,1:ncomp), 1 )
1040  !---- if one component is a polymer ------------------------------------------
1041  if ( minval( lnx(min_ph,1:ncomp) ) < -20.0 ) then
1042  CALL fugacity ( lnphi )
1043  lnx(min_ph,small_i) = lnx(other_ph,small_i)+lnphi(other_ph,small_i) - lnphi(min_ph,small_i)
1044  optpara(small_i) = lnx(2,small_i) + log( sum( exp( optpara(1:ncomp) ) ) )
1045  end if
1046 
1047  !-----------------------------------------------------------------------------
1048  ! caution: these initial values are for a flashcase overwritten in
1049  ! SUBROUTINE determine_flash_it2, because in that case, the lnx-values
1050  ! treated as ln(mole_number).
1051  !-----------------------------------------------------------------------------
1052  rhoi1( 1:ncomp ) = rhoi_cal( 1, 1:ncomp )
1053  rhoi2( 1:ncomp ) = rhoi_cal( 2, 1:ncomp )
1054 
1055  DEALLOCATE( optpara )
1056 
1057 end subroutine tangent_plane
1058 
1059 
1060 
1061 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1062 !
1063 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1064 
1065 subroutine helmholtz_flash_grads ( n, x_in, f_tpd, grad, hessian, diagonal )
1066 
1067  use parameters, only: kbol, pi
1068  use basic_variables, only: np, nc, ncomp, t, p, xi, xif, densta, my_cal, parame
1069  use eos_variables, only: dhs
1070  use utilities, only: fden_calc
1071  implicit none
1072 
1073  !-----------------------------------------------------------------------------
1074  integer, intent(in) :: n
1075  real, dimension(n), intent(in) :: x_in
1076  real, intent(out) :: f_tpd
1077  real, dimension(n), intent(out) :: grad
1078  real, dimension(n,n), intent(out) :: hessian
1079  real, dimension(n,n), intent(out) :: diagonal
1080 
1081  !-----------------------------------------------------------------------------
1082  integer :: i, j
1083  real :: phi
1084  real :: V
1085  real, dimension(ncomp) :: w
1086  real, dimension(ncomp) :: rhoi1, rhoi2
1087  real :: f1, f2
1088  real, dimension(ncomp) :: my1, my2
1089  real :: lnphi(np,nc)
1090  real :: df_phi, df_V
1091  real, dimension(ncomp) :: df_wi
1092  real, allocatable, dimension(:,:) :: A1_rr, Aig1_rr, A2_rr, Aig2_rr
1093  real, dimension(ncomp,ncomp) :: df_wi_wj
1094  real, dimension(ncomp) :: df_wi_V, df_wi_phi
1095  real :: df_phi_phi, df_phi_V, df_V_V
1096  real, dimension(2) :: sum_f_rhoi_rhoj_w_w, sum_f_rhoi_rhoj_w_r, sum_f_rhoi_rhoj_r_r
1097  !-----------------------------------------------------------------------------
1098 
1099  w(1:ncomp) = x_in(1:ncomp)
1100  phi = x_in(ncomp+1)
1101  v = x_in(ncomp+2)
1102 
1103  ! V1 = V * ( 0.5 + phi )
1104  ! V2 = V * ( 0.5 - phi )
1105 
1106  do i = 1, ncomp
1107  rhoi1(i) = xif(i) / v + ( phi - 0.5 ) * w(i)
1108  rhoi2(i) = xif(i) / v + ( phi + 0.5 ) * w(i)
1109  if ( rhoi1(i) < 0.0 ) rhoi1(i) = 0.0
1110  if ( rhoi2(i) < 0.0 ) rhoi2(i) = 0.0
1111  end do
1112 
1113  call fden_calc (f1, rhoi1)
1114  call fden_calc (f2, rhoi2)
1115 
1116  densta(1) = pi / 6.0 * sum( rhoi1(1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
1117  densta(2) = pi / 6.0 * sum( rhoi2(1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
1118  xi(1,1:ncomp) = rhoi1(1:ncomp) / sum( rhoi1(1:ncomp) )
1119  xi(2,1:ncomp) = rhoi2(1:ncomp) / sum( rhoi2(1:ncomp) )
1120  call fugacity ( lnphi )
1121 
1122  my1( 1:ncomp ) = my_cal(1,1:ncomp) ! lnphi(1,1:ncomp) + log( xi(1,1:ncomp) )
1123  my2( 1:ncomp ) = my_cal(2,1:ncomp) ! lnphi(2,1:ncomp) + log( xi(2,1:ncomp) )
1124 
1125  ! f_tpd = V1 * f1 + V2 * f2 + p / ( KBOL*1.E30 * t ) * V
1126  f_tpd = v * ( ( 0.5 + phi ) * f1 + ( 0.5 - phi ) * f2 + p / ( kbol*1.e30 * t ) )
1127  ! write (*,*) 'f_tpd',f_tpd
1128 
1129  df_phi = v * ( f1 -f2 + ( 0.5 + phi )*sum( my1(1:ncomp)*w(1:ncomp) ) &
1130  + ( 0.5 - phi )*sum( my2(1:ncomp)*w(1:ncomp) ) )
1131  df_v = ( f_tpd - ( 0.5 + phi )*sum( my1(1:ncomp)*xif(1:ncomp) ) &
1132  - ( 0.5 - phi )*sum( my2(1:ncomp)*xif(1:ncomp) ) ) / v
1133  do i = 1, ncomp
1134  df_wi(i) = v * ( phi*phi - 0.25 ) * ( my1(i) - my2(i) )
1135  end do
1136 
1137  grad(1:ncomp) = df_wi(1:ncomp)
1138  grad(ncomp+1) = df_phi
1139  grad(ncomp+2) = df_v
1140 
1141  allocate( a1_rr( ncomp, ncomp ), aig1_rr( ncomp, ncomp ) )
1142  allocate( a2_rr( ncomp, ncomp ), aig2_rr( ncomp, ncomp ) )
1143 
1144  call dda_drhoi_drhoj_eos ( ncomp, rhoi1, a1_rr, aig1_rr )
1145  a1_rr(:,:) = a1_rr(:,:) + aig1_rr(:,:)
1146 
1147  sum_f_rhoi_rhoj_w_w(1) = 0.0
1148  sum_f_rhoi_rhoj_w_r(1) = 0.0
1149  sum_f_rhoi_rhoj_r_r(1) = 0.0
1150  do i = 1, ncomp
1151  do j = 1, ncomp
1152  sum_f_rhoi_rhoj_w_w(1) = sum_f_rhoi_rhoj_w_w(1) + a1_rr(i,j)*w(i)*w(j)
1153  sum_f_rhoi_rhoj_w_r(1) = sum_f_rhoi_rhoj_w_r(1) + a1_rr(i,j)*w(i)*xif(j)
1154  sum_f_rhoi_rhoj_r_r(1) = sum_f_rhoi_rhoj_r_r(1) + a1_rr(i,j)*xif(i)*xif(j)
1155  end do
1156  end do
1157 
1158  call dda_drhoi_drhoj_eos ( ncomp, rhoi2, a2_rr, aig2_rr )
1159  a2_rr(:,:) = a2_rr(:,:) + aig2_rr(:,:)
1160 
1161 
1162  sum_f_rhoi_rhoj_w_w(2) = 0.0
1163  sum_f_rhoi_rhoj_w_r(2) = 0.0
1164  sum_f_rhoi_rhoj_r_r(2) = 0.0
1165  do i = 1, ncomp
1166  do j = 1, ncomp
1167  sum_f_rhoi_rhoj_w_w(2) = sum_f_rhoi_rhoj_w_w(2) + a2_rr(i,j)*w(i)*w(j)
1168  sum_f_rhoi_rhoj_w_r(2) = sum_f_rhoi_rhoj_w_r(2) + a2_rr(i,j)*w(i)*xif(j)
1169  sum_f_rhoi_rhoj_r_r(2) = sum_f_rhoi_rhoj_r_r(2) + a2_rr(i,j)*xif(i)*xif(j)
1170  end do
1171  end do
1172 
1173 
1174  do i = 1, ncomp
1175  do j = 1, ncomp
1176  df_wi_wj(i,j) = v * (phi*phi - 0.25 ) * ( (phi-0.5) * a1_rr(i,j) - (phi+0.5) * a2_rr(i,j) )
1177  end do
1178  end do
1179 
1180  do i = 1, ncomp
1181  df_wi_phi(i) = v * ( 2.0*phi*( my1(i)-my2(i) ) &
1182  + (phi*phi-0.25)*sum( w(1:ncomp)*(a1_rr(i,1:ncomp)-a2_rr(i,1:ncomp)) ) )
1183  df_wi_v(i) = ( df_wi(i) - (phi*phi-0.25)*sum( xif(1:ncomp)*(a1_rr(i,1:ncomp)-a2_rr(i,1:ncomp)) ) ) / v
1184  end do
1185 
1186  df_phi_phi = v * ( 2.0*sum( ( my1(1:ncomp)-my2(1:ncomp))*w(1:ncomp) ) &
1187  + ( 0.5 + phi ) * sum_f_rhoi_rhoj_w_w(1) &
1188  + ( 0.5 - phi ) * sum_f_rhoi_rhoj_w_w(2) )
1189  df_phi_v = ( df_phi - sum( ( my1(1:ncomp)-my2(1:ncomp))*xif(1:ncomp) ) &
1190  - ( 0.5 + phi ) * sum_f_rhoi_rhoj_w_r(1) &
1191  - ( 0.5 - phi ) * sum_f_rhoi_rhoj_w_r(2) ) / v
1192  df_v_v = ( (0.5+phi) * sum_f_rhoi_rhoj_r_r(1) + (0.5-phi) * sum_f_rhoi_rhoj_r_r(2) ) / v**3
1193 
1194  do i = 1, ncomp
1195  do j = 1, ncomp
1196  hessian(i,j) = df_wi_wj(i,j)
1197  end do
1198  hessian(i,ncomp+1) = df_wi_phi(i)
1199  hessian(i,ncomp+2) = df_wi_v(i)
1200  hessian(ncomp+1,i) = df_wi_phi(i)
1201  hessian(ncomp+2,i) = df_wi_v(i)
1202  end do
1203  hessian(ncomp+1,ncomp+1) = df_phi_phi
1204  hessian(ncomp+1,ncomp+2) = df_phi_v
1205  hessian(ncomp+2,ncomp+1) = df_phi_v
1206  hessian(ncomp+2,ncomp+2) = df_v_v
1207 
1208  do i = 1, ncomp
1209  diagonal(i,i) = v * (phi*phi - 0.25 ) * ( (phi-0.5)/rhoi1(i) - (phi+0.5)/rhoi2(i) )
1210  end do
1211  diagonal(ncomp+1,ncomp+1) = v * ( 2.0*sum( ( log( rhoi1(1:ncomp) )- log(rhoi2(1:ncomp)) )*w(1:ncomp) ) &
1212  + (0.5+phi) * sum(w(1:ncomp)*w(1:ncomp)/rhoi1(1:ncomp)) &
1213  + (0.5-phi) * sum(w(1:ncomp)*w(1:ncomp)/rhoi2(1:ncomp)) )
1214  diagonal(ncomp+2,ncomp+2) = ( (0.5+phi) * sum(xif(1:ncomp)*xif(1:ncomp)/rhoi1(1:ncomp)) &
1215  + (0.5-phi) * sum(xif(1:ncomp)*xif(1:ncomp)/rhoi2(1:ncomp)) ) / v**3
1216 
1217  !write (*,*) ' '
1218  !write (*,*) 'G1_i:',df_wi(1:ncomp)
1219  !write (*,*) 'G2 :',df_phi
1220  !write (*,*) 'G3 :',df_V
1221 
1222  deallocate( a1_rr, aig1_rr, a2_rr, aig2_rr )
1223 
1224 end subroutine helmholtz_flash_grads
1225 
1226 
1227 
1228 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1229 ! subroutine phase_stability
1230 !
1231 ! input: my_f (not in argument list) as the chem. potential of the feed-phase.
1232 ! rhoi_feed (argument list) density vector of feed phase
1233 ! eta_trial (argument list) as the starting density for trial phase.
1234 !
1235 ! output is: ph_split=0 for stable phase or ph_split=1 for a phase split. The
1236 ! array rhoi_trial( 1:ncomp ) is then an output, defining the minimum of the
1237 ! tangent plane distance.
1238 ! The value ph_split=2 is assigned for cases, where the TPD is not negative, but
1239 ! very close to zero and where the compositions are different to the feed-phase.
1240 ! This is encountered close to the critical point.
1241 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1242 
1243 subroutine phase_stability ( rhoi_feed, eta_trial, ph_split, rhoi_trial )
1244 
1245  use basic_variables
1246  use eos_variables, only: dhs, pi, mseg
1247  use eos_numerical_derivatives
1248  use cg_minimization
1249  use optimizer_2d
1250  implicit none
1251 
1252  !-----------------------------------------------------------------------------
1253  real, dimension(nc), intent(IN) :: rhoi_feed
1254  real, intent(IN) :: eta_trial
1255  integer, intent(OUT) :: ph_split
1256  real, dimension(nc), intent(OUT) :: rhoi_trial
1257  !-----------------------------------------------------------------------------
1258 
1259  integer :: n
1260  real :: fmin
1261  real, allocatable :: optpara(:)
1262 
1263  integer :: i
1264  integer :: i_trial_composition
1265  real :: rhoi(nc)
1266  real :: delta_rhoi(nc)
1267  real :: w(np,nc), mean_mass
1268  integer :: ph_split_i( 2*ncomp+1 )
1269  real :: fmin_vec( 2*ncomp+1 )
1270  real :: rhoi2_vec(2*ncomp+1, 1:ncomp)
1271  character (LEN=2) :: ensemble_flag_save
1272 
1273  integer :: STATUS, iter, nfunc, ngrad
1274  real :: gnorm
1275  real, allocatable :: d(:), g(:), xtemp(:), gtemp(:)
1276  !-----------------------------------------------------------------------------
1277 
1278 
1279  ph_split = 0
1280  ensemble_flag_save = ensemble_flag
1281  ensemble_flag = 'tv'
1282 
1283  n = ncomp
1284  ALLOCATE( optpara(n) )
1285  ALLOCATE( d(n) )
1286  ALLOCATE( g(n) )
1287  ALLOCATE( xtemp(n) )
1288  ALLOCATE( gtemp(n) )
1289 
1290 
1291  if ( outp >= 2 ) then
1292  write (*,'(a)') ' reference x_i, tested for stability:'
1293  write (*,'(a,4G21.12)') ' x_i feed:', rhoi_feed(1:ncomp) / sum( rhoi_feed(1:ncomp) )
1294  write (*,*) ' '
1295  end if
1296 
1297 
1298  do i_trial_composition = 0, (ncomp + ncomp)
1299 
1300  !--------------------------------------------------------------------------
1301  ! setting trial-phase mole-fractions
1302  !--------------------------------------------------------------------------
1303  if ( i_trial_composition == 0 ) w(2,1:ncomp) = 1.0 / REAL(ncomp)
1304  if ( i_trial_composition >= 1 .AND. i_trial_composition <= ncomp ) then
1305  w( 2, 1:ncomp ) = 1.0 / REAL( ncomp-1 ) * 0.05
1306  w( 2, i_trial_composition ) = 0.95
1307  else if ( i_trial_composition >= ncomp+1 ) then
1308  w( 2, 1:ncomp ) = 1.0 / REAL( ncomp-1 ) * 0.00001
1309  w( 2, i_trial_composition - ncomp ) = 0.99999
1310  end if
1311 
1312  mean_mass = 1.0 / sum( w(2,1:ncomp)/mm(1:ncomp) )
1313  xi(2,1:ncomp) = w(2,1:ncomp)/mm(1:ncomp) * mean_mass
1314  if ( outp >= 2 ) write (*,'(a,4G20.12)') ' trial phase x ---- ', xi(2,1:ncomp)
1315 
1316  !--------------------------------------------------------------------------
1317  ! starting values for iteration (optpara)
1318  !--------------------------------------------------------------------------
1319  do i = 1,ncomp
1320  rhoi(i) = xi(2,i)*eta_trial/sum( pi/6.0*xi(2,1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
1321  optpara(i) = log( rhoi(i) )
1322  !optpara(i) = rhoi(i)
1323  end do
1324 
1325  !--------------------------------------------------------------------------
1326  ! minimizing the objective fct. Phase split for values of fmin < 0.0
1327  !--------------------------------------------------------------------------
1328 
1329  if ( n == 2 ) then
1330  CALL newton_opt_2d ( f_stability, optpara, n, 1.e-8, 1.e-8, gtemp, fmin )
1331  fmin_best = fmin
1332  do i = 1, ncomp
1333  if ( optpara(i) > -50.0 ) then
1334  rhoi_best(i) = exp( optpara(i) )
1335  else
1336  rhoi_best(i) = 1.e-50
1337  end if
1338  end do
1339  else
1340  fmin_best = 1.e30
1341  CALL cg_descent ( 1.e-5, optpara, n, f_stability, stability_grad, status, &
1342  gnorm, fmin, iter, nfunc, ngrad, d, g, xtemp, gtemp)
1343  fmin = fmin_best
1344  end if
1345 
1346  !--------------------------------------------------------------------------
1347  ! save result to vector and determine stability of resulting phase
1348  !--------------------------------------------------------------------------
1349 
1350  delta_rhoi( 1:n ) = abs( 1.0 - rhoi_best(1:n) / rhoi_feed(1:n) )
1351 
1352  rhoi2_vec( i_trial_composition + 1, 1:ncomp ) = rhoi_best(1:ncomp)
1353  fmin_vec( i_trial_composition + 1 ) = 1.e8
1354 
1355  ph_split_i( i_trial_composition + 1 ) = 0
1356  if ( fmin < 1.e-3 .AND. maxval( delta_rhoi(1:n) ) > 0.1 ) then
1357  ph_split_i( i_trial_composition + 1 ) = 2
1358  fmin_vec( i_trial_composition + 1 ) = fmin
1359  end if
1360 
1361  if (fmin < -1.e-8 .AND. maxval( delta_rhoi(1:n) ) > 0.001 ) then
1362  ph_split_i( i_trial_composition + 1 ) = 1
1363  fmin_vec( i_trial_composition + 1 ) = fmin
1364  end if
1365 
1366 
1367  !--------------------------------------------------------------------------
1368  ! output to terminal
1369  !--------------------------------------------------------------------------
1370  if ( outp >= 2 ) then
1371  write (*,'(a, 5G20.12)') ' eta, xi_trialphase ',pi/6.0 &
1372  * sum( rhoi_best(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 ), &
1373  rhoi_best(1:ncomp) / sum( rhoi_best(1:ncomp) )
1374  write (*,'(a,G20.12,i4,G20.12,i3)') ' trial-result: fmin, N, d(x), ph-split:', &
1375  fmin, i_trial_composition, maxval( delta_rhoi( 1:n ) ), &
1376  ph_split_i( i_trial_composition+1 )
1377  write (*,*) ' '
1378  end if
1379 
1380 
1381  end do
1382 
1383  !-----------------------------------------------------------------------------
1384  ! write vector (rhoi_trial). Vector rhoi_trial gives the (overall) minimum of the TPD
1385  !-----------------------------------------------------------------------------
1386  i_trial_composition = minloc( fmin_vec( 1: (ncomp+ncomp)+1 ), 1 )
1387  rhoi_trial( 1:ncomp ) = rhoi2_vec( i_trial_composition, 1:ncomp )
1388  ph_split = ph_split_i( i_trial_composition )
1389 
1390  if ( outp >= 2 ) write (*,*) 'selected trial-result',i_trial_composition,fmin_vec( i_trial_composition )
1391 
1392  DEALLOCATE( optpara )
1393  DEALLOCATE( d, g, xtemp, gtemp )
1394  ensemble_flag = ensemble_flag_save
1395 
1396 end subroutine phase_stability
1397 
1398 
1399 
1400 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1401 !
1402 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1403 !
1404 subroutine f_stability ( f_tpd, optpara, n )
1405 
1406  use parameters, only: pi, kbol
1407  use basic_variables
1408  use eos_variables, only: dhs
1409  use utilities, only: fden_calc
1410  implicit none
1411 
1412  !-----------------------------------------------------------------------------
1413  integer, intent(in) :: n
1414  real, intent(in) :: optpara(:)
1415  real, intent(in out) :: f_tpd
1416 
1417  !-----------------------------------------------------------------------------
1418  integer :: i
1419  real :: rhoi(nc),gradterm
1420  real :: fden,punish
1421  real :: dens
1422  !-----------------------------------------------------------------------------
1423 
1424  !-----------------------------------------------------------------------------
1425  ! saveguard that the entering parameters are reasonable
1426  !-----------------------------------------------------------------------------
1427  if ( maxval( optpara(1:n) ) > 2.0 ) then
1428  f_tpd = 100.0 * maxval( optpara(1:n) )
1429  return
1430  end if
1431  do i = 1, n
1432  if ( optpara(i) /= optpara(i) ) then
1433  f_tpd = 1.e8
1434  return
1435  end if
1436  end do
1437  if ( maxval( optpara(1:n) ) < -50.0 ) then
1438  f_tpd = - 100.0 * sum( optpara(1:n) )
1439  return
1440  end if
1441 
1442  !-----------------------------------------------------------------------------
1443  ! assign density of each component
1444  !-----------------------------------------------------------------------------
1445 
1446  punish = 0.0
1447 
1448  do i = 1, n
1449  if ( optpara(i) > -50.0 ) then
1450  rhoi(i) = exp( optpara(i) )
1451  else
1452  rhoi(i) = 1.e-50
1453  end if
1454  !if ( optpara(i) < 0.0) then
1455  ! rhoi(i) = 1.E-50 / (1.0-optpara(i)*1.E10)
1456  ! punish = punish - optpara(i)*1.E6
1457  !end if
1458  !if ( optpara(i) >= 0.6) rhoi(i) = 0.6
1459  end do
1460 
1461  dens = pi/6.0 * sum( rhoi(1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
1462 
1463  !-----------------------------------------------------------------------------
1464  ! correcting too high densities
1465  !-----------------------------------------------------------------------------
1466 
1467  if (dens > 0.65) then
1468  punish = punish + (dens-0.65)*1000.0
1469  rhoi(1:n) = rhoi(1:n)*0.65/dens
1470  end if
1471 
1472  !-----------------------------------------------------------------------------
1473  ! calculate the Helmholtz energ and the tangent plane distance (TPD)
1474  !-----------------------------------------------------------------------------
1475 
1476  call fden_calc (fden, rhoi)
1477 
1478  gradterm = sum( my_f(1:n) * rhoi(1:n) )
1479 
1480  f_tpd = fden + (p * 1.e-30) / (kbol*t) - gradterm + punish
1481  f_tpd = f_tpd * 1000.0
1482  ! write (*,'(a,5G16.8)') 'f_tpd', f_tpd, optpara(1:n), punish
1483 
1484  if ( f_tpd < fmin_best ) then
1485  fmin_best = f_tpd
1486  rhoi_best(1:n) = rhoi(1:n)
1487  end if
1488 
1489 end subroutine f_stability
1490 
1491 
1492 
1493 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1494 !
1495 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1496 
1497 subroutine stability_grad (g, optpara, n)
1498 
1499  use basic_variables
1500  use eos_variables, only: dhs, pi
1501  implicit none
1502 
1503  integer, intent(in) :: n
1504  real, intent(in) :: optpara(:)
1505  real, intent(in out) :: g(:)
1506 
1507  !-----------------------------------------------------------------------------
1508  integer :: i
1509  integer :: nphas_save
1510  real :: rhoi(nc)
1511  real :: lnphi(np,nc)
1512  !-----------------------------------------------------------------------------
1513 
1514  nphas_save = nphas
1515  nphas = 1
1516 
1517  !-----------------------------------------------------------------------------
1518  ! saveguard that the entering parameters are reasonable
1519  !-----------------------------------------------------------------------------
1520 
1521  do i = 1, n
1522  if ( optpara(i) /= optpara(i) .OR. optpara(i) > 5.0 ) then
1523  g(:) = 0.0
1524  nphas = nphas_save
1525  return
1526  end if
1527  end do
1528  if ( maxval( optpara(1:n) ) < -50.0 ) then
1529  g(:) = 0.0
1530  nphas = nphas_save
1531  return
1532  end if
1533 
1534  !-----------------------------------------------------------------------------
1535  ! assign density of all species
1536  !-----------------------------------------------------------------------------
1537 
1538  do i = 1, n
1539  if ( optpara(i) > -50.0 ) then
1540  rhoi(i) = exp( optpara(i) )
1541  else
1542  rhoi(i) = 1.e-50
1543  end if
1544  !if ( optpara(i) < 0.0) rhoi(i) = 1.E-50 / (1.0-optpara(i)*1.E10)
1545  !if ( optpara(i) >= 0.6) rhoi(i) = 0.6
1546  end do
1547 
1548  densta(1) = pi/6.0 * sum( rhoi(1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
1549 
1550  !-----------------------------------------------------------------------------
1551  ! correcting too high densities
1552  !-----------------------------------------------------------------------------
1553 
1554  if (densta(1) > 0.65) then
1555  rhoi(1:n) = rhoi(1:n) * 0.65 / densta(1)
1556  densta(1) = pi/6.0 * sum( rhoi(1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
1557  end if
1558 
1559  xi(1,1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
1560 
1561  !-----------------------------------------------------------------------------
1562  ! calculate the chemical potential of all species
1563  !-----------------------------------------------------------------------------
1564 
1565  call fugacity (lnphi)
1566 
1567  g( 1:n ) = lnphi( 1, 1:n ) + log( rhoi( 1:n ) ) - my_f( 1:n )
1568  g( 1:n ) = g( 1:n ) * rhoi( 1:n ) * 1000.0
1569  ! write (*,'(a,5G16.8)') 'f_tpd_grad',g(1:n)
1570 
1571  nphas = nphas_save
1572 
1573 end subroutine stability_grad
1574 
1575 
1576 
1577 
1578 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1579 !
1580 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1581 
1582 subroutine stability_hessian (hessian, gtrans, fmin, optpara, n)
1583 
1584  use basic_variables
1585  implicit none
1586 
1587  integer, intent(IN) :: n
1588  real, intent(IN) :: optpara(:)
1589  real, intent(IN OUT) :: fmin
1590  real, intent(IN OUT) :: gtrans(:)
1591  real, intent(IN OUT) :: hessian(:,:)
1592 
1593  !-----------------------------------------------------------------------------
1594  integer :: i, j
1595  real :: delta
1596  real :: optpara_mod(n), g(n), gi_right(n,n), gi_left(n,n)
1597  real, allocatable, dimension(:,:) :: A_rr, Aig_rr
1598  !-----------------------------------------------------------------------------
1599 
1600 
1601  delta = 1.0e-4
1602 
1603  optpara_mod = optpara
1604 
1605  do i = 1, n
1606 
1607  optpara_mod = optpara
1608  optpara_mod(i) = optpara(i)*(1.0+delta)
1609 
1610  call stability_grad (g, optpara_mod, n)
1611  gi_right(i,1:n) = g(1:n)
1612 
1613  optpara_mod(i) = optpara(i)*(1.0-delta)
1614 
1615  call stability_grad (g, optpara_mod, n)
1616  gi_left(i,1:n) = g(1:n)
1617 
1618  end do
1619 
1620  call stability_grad (g, optpara, n)
1621 
1622  do i = 1, n
1623  do j = 1, n
1624 
1625  hessian(i,j) = ( gi_right(i,j) - gi_left(i,j) ) / ( 2.0*optpara(i)*delta )
1626  ! hessian(j,i) = hessian(i,j)
1627  write (*,*) i,j,hessian(i,j)
1628 
1629  end do
1630  end do
1631  write (*,*) ' '
1632 
1633  gtrans = g
1634  call f_stability (fmin, optpara, n)
1635 
1636  allocate( a_rr(ncomp,ncomp), aig_rr(ncomp,ncomp) )
1637  call dda_drhoi_drhoj_eos ( ncomp, exp( optpara(:) ), a_rr, aig_rr )
1638  do i = 1, n
1639  do j = 1, n
1640  hessian(i,j) = a_rr(i,j)
1641  if (i==j) hessian(i,j) = hessian(i,j) + 1.0 / exp( optpara(i) )
1642  !--- convert to second derivative ddA / dln(rho_k) dln(rho_l) ( from ... d(rho_l) )
1643  if (i/=j) hessian(i,j) = hessian(i,j)*exp( optpara(i) )*exp( optpara(j) )*1000.0
1644  if (i==j) hessian(i,j) = hessian(i,j)*exp( optpara(i) )*exp( optpara(j) )*1000.0 + g(i)
1645  if (i/=j) write (*,*) i,j,hessian(i,j)
1646  if (i==j) write (*,*) i,j,hessian(i,j)
1647  end do
1648  end do
1649  deallocate( a_rr, aig_rr )
1650  read (*,*)
1651 
1652 end subroutine stability_hessian
1653 
1654 
1655 
1656 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1657 !
1658 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1659 
1660 subroutine vle_min
1661 
1662  use parameters, only: rgas
1663  use basic_variables
1664  use utilities
1665  implicit none
1666 
1667  !-----------------------------------------------------------------------------
1668  integer, parameter :: steps = 40
1669  logical :: lle_check ! currently not any more used
1670  integer :: i, j, k, phasen(0:steps)
1671  real :: lnphi(np,nc)
1672  real, dimension(0:steps) :: vlemin, llemin, xval, start_xv, start_xl
1673  real :: x_sav,dg_dx2
1674  !-----------------------------------------------------------------------------
1675 
1676  j = 0
1677  k = 0
1678  nphas = 2
1679 
1680  start_xv(:) = 0.0
1681  start_xl(:) = 0.0
1682 
1683  x_sav = xi(1,1)
1684  sum_rel(1) = 'x12' ! summation relation
1685  sum_rel(2) = 'x22' ! summation relation
1686 
1687  do i = 0, steps
1688 
1689  densta(1) = 0.45
1690  densta(2) = 1.e-6
1691  xi(1,1) = 1.0 - REAL(i) / REAL(steps)
1692  if ( xi(1,1) <= 1.e-50 ) xi(1,1) = 1.e-50
1693  xi(2,1) = xi(1,1)
1694  lnx(1,1) = log(xi(1,1))
1695  lnx(2,1) = log(xi(2,1))
1696 
1697  CALL x_summation
1698  CALL fugacity ( lnphi )
1699 
1700  xval(i) = xi(1,1)
1701  llemin(i)= gibbs(1) * rgas * t
1702 
1703  if ( abs(1.0-dense(1)/dense(2)) > 0.0001 ) then
1704  vlemin(i) = ( gibbs(1) - gibbs(2) ) * rgas * t
1705  phasen(i) = 2
1706  else
1707  phasen(i) = 1
1708  end if
1709 
1710  if (i > 0 .AND. phasen(i) == 2) then
1711 
1712  if ( phasen(i-1) == 2 .AND. abs( vlemin(i) + vlemin(i-1) ) < &
1713  abs( vlemin(i) ) + abs( vlemin(i-1) ) ) then
1714  j = j + 1
1715  start_xv(j) = xval(i-1) + (xval(i)-xval(i-1)) &
1716  * abs( vlemin(i-1) ) / abs( vlemin(i) - vlemin(i-1) )
1717  end if
1718 
1719  end if
1720 
1721  end do
1722 
1723  do i = 2, steps - 2
1724 
1725  dg_dx2 = (-llemin(i-2) + 16.0*llemin(i-1) - 30.0*llemin(i) &
1726  + 16.0*llemin(i+1) - llemin(i+2)) / ( 12.0 * ((xval(i)-xval(i-1))**2) )
1727  if ( dg_dx2 < 0.0 ) then
1728  k = k + 1
1729  start_xl(k) = xval(i)
1730  end if
1731 
1732  end do
1733 
1734  if ( start_xl(1) == 0.0 .AND. start_xv(1) /= 0.0 ) then
1735  xi(1,1) = start_xv(1)
1736  xi(1,2) = 1.0 - xi(1,1)
1737  lle_check = .false.
1738  if ( outp >= 1) write (*,*) 'VLE is likely', xi(1,1),xi(1,2)
1739  else if ( start_xl(1) /= 0.0 .AND. start_xv(1) == 0.0 ) then
1740  xi(1,1) = start_xl(1)
1741  xi(1,2) = 1.0 - xi(1,1)
1742  if ( outp >= 1) write (*,*) 'LLE is likely', xi(1,1),xi(1,2)
1743  lle_check = .true.
1744  else if ( start_xl(1) /= 0.0 .AND. start_xv(1) /= 0.0 ) then
1745  xi(1,1) = start_xv(1)
1746  xi(1,2) = 1.0 - xi(1,1)
1747  if ( outp >= 1) write(*,*) 'starting with VLE and check for LLE'
1748  lle_check = .true.
1749  else
1750  xi(1,1) = x_sav
1751  xi(1,2) = 1.0 - xi(1,1)
1752  end if
1753 
1754  CALL x_summation
1755 
1756 end subroutine vle_min
1757 
1758 
1759 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1760 ! subroutine rachford_rice
1761 !
1762 ! This subroutine performs a p-T flash for a defined composition (xiF) in a
1763 ! feed-trial phase by iterating the Rachford-Rice equation.
1764 !
1765 ! input: xiF(:), xi(1,:), xi(2,:), dense(1), dense(2), t, p
1766 !
1767 ! output: rhoi1(:), rhoi2(:)
1768 ! xi(1,:), xi(2,:)
1769 ! values of: dense(:), lnphi(:,:) are also converged values, if converg=1
1770 ! ph_frac (phase fraction) could be defined as output, but is not yet.
1771 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1772 
1773 subroutine rachford_rice ( converg, rhoi1, rhoi2 )
1774 
1775  use basic_variables
1776  use optimizer_2d, only: newton_opt_analytic
1777  implicit none
1778 
1779  !-----------------------------------------------------------------------------
1780  integer, intent(in out) :: converg
1781  real, dimension(nc), intent(out) :: rhoi1
1782  real, dimension(nc), intent(out) :: rhoi2
1783 
1784  !-----------------------------------------------------------------------------
1785  integer :: k1, k2
1786  integer :: nr_inside_steps, nr_outside_steps
1787  real, parameter :: tol_out = 1.e-11
1788  real, parameter :: tol_in = 1.e-13
1789  real :: ph_frac, ph_frac0
1790  real :: lnphi(np,nc), Ki(nc)
1791  real :: f1, d_ph, dfdph
1792  real :: error1, error2 = 0.0
1793  real :: xi_compare(2,nc)
1794  logical :: comp_balance_violation = .false.
1795 
1796  integer :: i
1797  real, allocatable, dimension(:) :: w
1798  real, allocatable, dimension(:) :: x_in, grad
1799  real, allocatable, dimension(:,:) :: hessian
1800  real :: phi, V1, V2, V, phi_00, f_tpd
1801  !-----------------------------------------------------------------------------
1802  outp = 3
1803  converg = 0
1804  nphas = 2 ! number of phases
1805 
1806  nr_outside_steps = 80
1807  nr_inside_steps = 10
1808 
1809  xi_compare(1:2,:) = xi(1:2,:)
1810 
1811  ensemble_flag = 'tp'
1812  densta(1) = dense(1) ! Index 1 is for liquid density (here: packing fraction eta)
1813  densta(2) = dense(2) ! Index 2 is for vapour density (here: packing fraction eta)
1814 
1815  ph_frac = 0.5
1816  d_ph = 0.00001
1817  if ( outp >= 2 ) write (*,'(a,4G20.12)') 'RR, N_out, N_inner, xi(1,1), xiF(1), xi(2,1), ph_frac, error1'
1818 
1819  !-----------------------------------------------------------------------------
1820  ! start iteration
1821  !-----------------------------------------------------------------------------
1822 
1823  k1 = 0
1824  error1 = tol_out + 1.0
1825  do while ( error1 > tol_out .AND. k1 < nr_outside_steps )
1826 
1827  !---------------------------------------------------------------------------
1828  ! outer loop (converging the mole fractions)
1829  !---------------------------------------------------------------------------
1830 
1831  k1 = k1 + 1
1832 
1833  CALL fugacity ( lnphi )
1834 
1835  ki(1:ncomp) = exp( lnphi(1,1:ncomp) - lnphi(2,1:ncomp) )
1836  write (*,*) 'dense',dense(1),dense(2)
1837 
1838  if ( abs(1.0-dense(1)/dense(2)) < 1.e-6 .AND. sum(abs(ki(1:ncomp)-1.0)) < 1.e-8 ) exit
1839 
1840  k2 = 0
1841  error2 = 1.0
1842  do while ( error2 > tol_in .AND. k2 < nr_inside_steps )
1843 
1844  !------------------------------------------------------------------------
1845  ! inner loop (converging the fraction of one phase)
1846  !------------------------------------------------------------------------
1847 
1848  k2 = k2 + 1
1849 
1850  ph_frac0 = ph_frac
1851  f1 = sum( xif(1:ncomp)*(ki(1:ncomp)-1.0) / ( 1.0 + (ki(1:ncomp)-1.0)*ph_frac ) )
1852 
1853  dfdph = - sum( xif(1:ncomp) / ( 1.0/(ki(1:ncomp)-1.0) + ph_frac )**2 )
1854  if ( dfdph == 0.0 ) dfdph = 0.0000001
1855 
1856  ph_frac = ph_frac0 - f1 / dfdph
1857 
1858  error2 = abs( ph_frac - ph_frac0 )
1859  if ( outp >= 3 ) write (*,'(a,3G20.8)') 'phase_fraction, error', ph_frac, ph_frac0, error2
1860 
1861  end do
1862 
1863  !---------------------------------------------------------------------------
1864  ! constrain the phase fraction to ( 0 < ph_frac < 1 )
1865  !---------------------------------------------------------------------------
1866  comp_balance_violation = .false.
1867  if ( ph_frac > 1.0 ) comp_balance_violation = .true.
1868  if ( ph_frac < 0.0 ) comp_balance_violation = .true.
1869  if ( ph_frac > 1.0 ) ph_frac = 0.9999999
1870  if ( ph_frac < 0.0 ) ph_frac = 0.0000001
1871  !if ( k2 >= nr_inside_steps ) write (*,*) 'Rachford-Rice: Newton-loop not converged'
1872 
1873  !---------------------------------------------------------------------------
1874  ! determine x and error of outer loop
1875  !---------------------------------------------------------------------------
1876  xi(1,1:ncomp) = xif(1:ncomp) / ( 1.0 + ph_frac *( ki(1:ncomp)-1.0 ) )
1877  xi(2,1:ncomp) = ki(1:ncomp)* xif(1:ncomp) / ( 1.0 + ph_frac *( ki(1:ncomp)-1.0 ) )
1878 
1879  if ( sum( xi(2,1:ncomp) ) < 0.2 ) xi(2,1:ncomp) = xi(2,1:ncomp) / sum( xi(2,1:ncomp) )
1880  if ( sum( xi(1,1:ncomp) ) < 0.2 ) xi(1,1:ncomp) = xi(1,1:ncomp) / sum( xi(1,1:ncomp) )
1881 
1882  error1 = sum( abs( xi_compare(1,1:ncomp) - xi(1,1:ncomp) ) ) &
1883  + sum( abs( xi_compare(2,1:ncomp) - xi(2,1:ncomp) ) )
1884  xi_compare(1:2,:) = xi(1:2,:)
1885 
1886  if ( outp >= 2 ) write (*,'(a,2i5,5G20.12)') 'RR',k1,k2,xi(1,1), xif(1),xi(2,1), ph_frac, error1
1887 
1888  if ( k1 == 100 .AND. outp > 0 ) write (*,*) 'Rachford-Rice: outside loop slow convergence'
1889 
1890  end do
1891 
1892 
1893  !-----------------------------------------------------------------------------
1894  ! if convergence: accept solution
1895  !-----------------------------------------------------------------------------
1896 
1897  if ( error1 < tol_out .AND. error2 < tol_in .AND. .NOT. comp_balance_violation ) then
1898 
1899  rhoi1( 1:ncomp ) = rhoi_cal( 1, 1:ncomp )
1900  rhoi2( 1:ncomp ) = rhoi_cal( 2, 1:ncomp )
1901  converg = 1
1902  if ( outp >= 2 ) write (*,*) ' '
1903  if ( outp >= 2 ) write (*,*) '========================================================'
1904  if ( outp >= 1 ) write (*,*) 'Rachford-Rice converged'
1905  if ( outp >= 2 ) write (*,*) '========================================================'
1906  if ( outp >= 2 ) write (*,*) ' x p1 = ',rhoi1( 1:ncomp ) / sum( rhoi1( 1:ncomp ) )
1907  if ( outp >= 2 ) write (*,*) ' x p2 = ',rhoi2( 1:ncomp ) / sum( rhoi2( 1:ncomp ) )
1908  if ( outp >= 2 ) write (*,*) 'eta 1,2 = ',dense(1), dense(2)
1909  if ( outp >= 2 ) write (*,*) ' '
1910  alpha = ph_frac
1911 
1912  end if
1913 
1914  !-----------------------------------------------------------------------------
1915  ! for slow convergence, take second order approach.
1916  ! algorithm proposed by Nagarajan, Cullick, Griewank (FPE, 1991, 191-210).
1917  ! So far, a rough criterion for the approach to convergence is used.
1918  !-----------------------------------------------------------------------------
1919 
1920  if ( converg /= 1 .AND. error1 < 1.e-4 ) then
1921 
1922  if( outp >= 1) write (*,*) 'Rachf.-R.: not converged'
1923 
1924  phi_00 = 0.5
1925  ensemble_flag = 'tv'
1926  allocate ( w(ncomp), x_in(ncomp+2), grad(ncomp+2), hessian(ncomp+2,ncomp+2) )
1927  rhoi1(1:ncomp) = rhoi_cal(1,1:ncomp)
1928  rhoi2(1:ncomp) = rhoi_cal(2,1:ncomp)
1929  w(1:ncomp) = rhoi2(1:ncomp) - rhoi1(1:ncomp)
1930  v1 = (1.0-phi_00) * sum( xif(1:ncomp) ) / sum( rhoi1(1:ncomp) )
1931  v2 = phi_00 * sum( xif(1:ncomp) ) / sum( rhoi2(1:ncomp) )
1932  write (*,*) 'phi_00, xiF(1:ncomp)', phi_00, xif(1:ncomp)
1933  write (*,*) 'V1, V2', v1, v2
1934  write (*,*) 'check:',1, (1.0-ph_frac)*xi(1,1) + ph_frac*xi(2,1), xif(1)
1935  write (*,*) 'check:',2, (1.0-ph_frac)*xi(1,2) + ph_frac*xi(2,2), xif(2)
1936  v = v1 + v2
1937  phi = 0.5 * ( v1 - v2 ) / v
1938  x_in(1:ncomp) = w(1:ncomp)
1939  x_in(ncomp+1) = phi
1940  x_in(ncomp+2) = v
1941  call newton_opt_analytic ( helmholtz_flash_grads, (ncomp+2), x_in, 1.e-7, 1.e-12, grad, f_tpd, converg )
1942 
1943  w(1:ncomp) = x_in(1:ncomp)
1944  phi = x_in(ncomp+1)
1945  v = x_in(ncomp+2)
1946  do i = 1, ncomp
1947  rhoi1(i) = xif(i) / v + ( phi - 0.5 ) * w(i)
1948  rhoi2(i) = xif(i) / v + ( phi + 0.5 ) * w(i)
1949  if ( rhoi1(i) < 0.0 ) rhoi1(i) = 0.0
1950  if ( rhoi2(i) < 0.0 ) rhoi2(i) = 0.0
1951  end do
1952 
1953  if ( maxval(abs(grad(:))) < 1.e-7 ) then
1954 
1955  converg = 1
1956  if ( outp >= 2 ) write (*,*) ' '
1957  if ( outp >= 2 ) write (*,*) '========================================================'
1958  if ( outp >= 1 ) write (*,*) 'Newton iteration (subsequent to Rachford-Rice) converged'
1959  if ( outp >= 2 ) write (*,*) '========================================================'
1960  if ( outp >= 2 ) write (*,*) ' x p1 = ',rhoi1( 1:ncomp ) / sum( rhoi1( 1:ncomp ) )
1961  if ( outp >= 2 ) write (*,*) ' x p2 = ',rhoi2( 1:ncomp ) / sum( rhoi2( 1:ncomp ) )
1962  if ( outp >= 2 ) write (*,*) 'eta 1,2 = ',dense(1), dense(2)
1963  if ( outp >= 2 ) write (*,*) ' '
1964  alpha = v * ( 0.5 - phi )*sum( rhoi2(1:ncomp) )
1965 
1966  end if
1967 
1968  deallocate( w, x_in, grad, hessian )
1969 
1970  end if
1971 
1972 
1973 end subroutine rachford_rice
1974 
1975 
1976 
1977 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1978 ! subroutine bubble_point_rachford_rice
1979 !
1980 ! This subroutine performs a bubble- or dew-point calculation for a defined
1981 ! composition (xiF). A bubble point calculation for a given liquid composition
1982 ! is done, when dense(1) is a liquid density (starting value, needs to be
1983 ! possible to converge). A dew point calculation is done, when dense(1) (low
1984 ! starting value) can converge to a vapor density.
1985 !
1986 ! Either, the temperature or the pressure is iterated. The algorith is based on
1987 ! the Rachford-Rice equation.
1988 !
1989 ! input: xi(1,:), t or p
1990 ! And starting values: dense(1), dense(2), p or t
1991 !
1992 ! output: p or t, xi(2,:), rhoi1(:), rhoi1(:)
1993 ! values of: dense(:), lnphi(:,:) are also converged values, if converg=1
1994 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1995 
1996 subroutine bubble_point_rachford_rice ( iterate_t, converg, rhoi1, rhoi2 )
1997 
1998  use basic_variables
1999  implicit none
2000 
2001  !-----------------------------------------------------------------------------
2002  integer, intent(in) :: iterate_t
2003  integer, intent(in out) :: converg
2004  real, dimension(nc), intent(out) :: rhoi1
2005  real, dimension(nc), intent(out) :: rhoi2
2006 
2007  !-----------------------------------------------------------------------------
2008  integer :: k1, k2
2009  integer :: nr_inside_steps, nr_outside_steps
2010  real, parameter :: tol_out = 1.e-11
2011  real, parameter :: tol_in_p = 1.e-10
2012  real, parameter :: tol_in_t = 1.e-8
2013  real :: tol_in
2014  real :: p0, t0, dp, dt, deltap
2015  real :: p_sav, t_sav
2016  real :: damping, acceleration
2017  real :: lnphi(np,nc), Ki(nc)
2018  real :: f0, f1, dfdp
2019  real :: error1, error2 = 0.0
2020  real :: xi_compare(nc)
2021  logical :: SR_exit
2022  !-----------------------------------------------------------------------------
2023 
2024  t_sav = t
2025  p_sav = p
2026 
2027  converg = 0
2028  nphas = 2 ! number of phases
2029 
2030  nr_outside_steps = 1000
2031  nr_inside_steps = 10
2032  if ( iterate_t == 1 ) tol_in = tol_in_t
2033  if ( iterate_t == 0 ) tol_in = tol_in_p
2034  sr_exit = .false.
2035 
2036  xi_compare(:) = xi(2,:)
2037 
2038  ensemble_flag = 'tp'
2039  densta(1) = dense(1) ! Index 1 is for liquid density (here: packing fraction eta)
2040  densta(2) = dense(2) ! Index 2 is for vapour density (here: packing fraction eta)
2041 
2042  dp = 1.0 ! pressure step in unit [Pa]
2043  dt = 1.e-3 ! temperature step in unit [K]
2044  if ( outp >= 2 ) write (*,'(a,4G20.12)') 'bbp-RR, N_out, N_inner, xi(1,1), xi(2,1), p, t, error1'
2045 
2046  !-----------------------------------------------------------------------------
2047  ! start iteration
2048  !-----------------------------------------------------------------------------
2049 
2050  k1 = 0
2051  error1 = tol_out + 1.0
2052  do while ( error1 > tol_out .AND. k1 < nr_outside_steps )
2053 
2054  !--------------------------------------------------------------------------
2055  ! outer loop (converging the mole fractions)
2056  !--------------------------------------------------------------------------
2057 
2058  k1 = k1 + 1
2059 
2060  k2 = 0
2061  error2 = tol_in + 1.0
2062  do while ( error2 > tol_in .AND. k2 < nr_inside_steps )
2063 
2064  !-----------------------------------------------------------------------
2065  ! inner loop (converging the pressure)
2066  !-----------------------------------------------------------------------
2067 
2068  k2 = k2 + 1
2069 
2070  p0 = p
2071  t0 = t
2072  if ( iterate_t == 1 ) then
2073  t = t0 - dt
2074  else
2075  p = p0 - dp
2076  end if
2077 
2078  CALL fugacity ( lnphi )
2079 
2080  if ( abs(1.0-dense(1)/dense(2)) < 1.e-6 .AND. sum(abs(ki(1:ncomp)-1.0)) < 1.e-8 ) exit
2081  ki(1:ncomp) = exp( lnphi(1,1:ncomp) - lnphi(2,1:ncomp) )
2082  f0 = sum( xi(1,1:ncomp) * ki(1:ncomp) ) - 1.0
2083 
2084  p = p0
2085  t = t0
2086  CALL fugacity ( lnphi )
2087 
2088  ki(1:ncomp) = exp( lnphi(1,1:ncomp) - lnphi(2,1:ncomp) )
2089  f1 = sum( xi(1,1:ncomp) * ki(1:ncomp) ) - 1.0
2090 
2091  if ( iterate_t == 1 ) then
2092 
2093  dfdp = ( f1 - f0 ) / dt
2094  if ( dfdp == 0.0 ) dfdp = 0.0000001
2095 
2096  !--- Newton ---------------------------------------------------------
2097  t = t0 - f1 / dfdp
2098 
2099  !--- error ----------------------------------------------------------
2100  error2 = abs( t - t0 )
2101  if ( outp >= 3 ) write (*,'(a,3G20.12)') 'inner loop error', t, t0, error2
2102  if ( error2 > 2.0 .OR. t < 0.0 ) sr_exit = .true.
2103  if ( error2 > 2.0 .OR. t < 0.0 ) exit
2104 
2105  else
2106 
2107  dfdp = ( f1 - f0 ) / ( log(p0) - log(p0 - dp) )
2108  if ( dfdp == 0.0 ) dfdp = 0.0000001
2109 
2110  !--- Newton ---------------------------------------------------------
2111  ! p = EXP( LOG(p0) - f1 / dfdp )
2112  deltap = f1 / dfdp
2113  if ( deltap > 2.0 ) deltap = 2.0
2114  if ( deltap < -2.0 ) deltap = -2.0
2115  p = log( p0) - deltap
2116  p = exp( p )
2117  p = min( p, 1.e9 )
2118 
2119  !--- damping --------------------------------------------------------
2120  damping = 0.6
2121  if ( abs( p / p0 - 1.0 ) > 0.3 ) p = damping * p + ( 1.0 - damping ) * p0
2122 
2123  !--- error ----------------------------------------------------------
2124  error2 = abs( p / p0 - 1.0 )
2125 
2126  if ( outp >= 3 ) write (*,'(a,3G20.12)') 'inner loop error', p, p0, error2
2127  if ( error2 > 10.0 ) sr_exit = .true.
2128  if ( error2 > 10.0 ) exit
2129 
2130  end if
2131 
2132  end do
2133 
2134  if ( sr_exit ) exit
2135 
2136  !--------------------------------------------------------------------------
2137  ! determine x and error of outer loop
2138  !--------------------------------------------------------------------------
2139  xi(2,1:ncomp) = ki(1:ncomp)* xi(1,1:ncomp) / sum( xi(1,1:ncomp) * ki(1:ncomp) )
2140  xi(2,1:ncomp) = max( xi(2,1:ncomp), 1.e-50 )
2141 
2142  if ( sum( xi(2,1:ncomp) ) < 0.4 ) xi(2,1:ncomp) = xi(2,1:ncomp) / sum( xi(2,1:ncomp) )
2143 
2144  error1 = sum( abs( xi_compare(1:ncomp) - xi(2,1:ncomp) ) )
2145 
2146  if ( k1 > 100 ) then
2147  acceleration = real(k1) / 100.0 + 1.0
2148  if ( k1 > 900 ) acceleration = acceleration * 2.0
2149  xi(2,1:ncomp) = xi(2,1:ncomp) + acceleration * ( xi(2,1:ncomp) - xi_compare(1:ncomp) ) ! acceleration
2150  xi(2,1:ncomp) = max( xi(2,1:ncomp), 1.e-50 )
2151  ! write (*,'(a,4g20.12)') 'accel',acceleration,xi(2,1), xi(2,2), sum( xi(2,1:ncomp) )
2152  xi(2,1:ncomp) = xi(2,1:ncomp) / sum( xi(2,1:ncomp) )
2153  end if
2154 
2155  xi_compare(:) = xi(2,:)
2156 
2157  if ( error2 > tol_in ) error1 = error1 + tol_out
2158 
2159  if ( k1 >= 500 ) error1 = error1 / 10.0
2160  if ( k1 >= 1000 ) error1 = error1 / 100.0
2161 
2162  if ( outp >= 2 ) write (*,'(a,2i5,5G18.10)') 'bbp-RR',k1,k2,xi(1,1), xi(2,1), p, t, error1
2163 
2164  if ( k1 == 100 .AND. outp > 0 ) write (*,*) 'bbp-Rachford-Rice: outside loop slow convergence'
2165 
2166  end do
2167 
2168 
2169  !-----------------------------------------------------------------------------
2170  ! if convergence: accept solution
2171  !-----------------------------------------------------------------------------
2172 
2173  if ( error1 < tol_out .AND. error2 < tol_in ) then
2174 
2175  rhoi1( 1:ncomp ) = rhoi_cal( 1, 1:ncomp )
2176  rhoi2( 1:ncomp ) = rhoi_cal( 2, 1:ncomp )
2177  converg = 1
2178  if ( outp >= 2 ) write (*,*) ' '
2179  if ( outp >= 2 ) write (*,*) '========================================================'
2180  if ( outp >= 1 ) write (*,*) 'bubble point Rachford-Rice converged'
2181  if ( outp >= 2 ) write (*,*) '========================================================'
2182  if ( outp >= 2 ) write (*,*) ' x p1 = ',rhoi1( 1:ncomp ) / sum( rhoi1( 1:ncomp ) )
2183  if ( outp >= 2 ) write (*,*) ' x p2 = ',rhoi2( 1:ncomp ) / sum( rhoi2( 1:ncomp ) )
2184  if ( outp >= 2 ) write (*,*) 'eta 1,2 = ',dense(1), dense(2)
2185  if ( outp >= 2 ) write (*,*) ' '
2186 
2187  else
2188 
2189  t = t_sav
2190  p = p_sav
2191  if ( outp >= 1) write (*,*) 'bbp-Rachf.-R.: not converged'
2192 
2193  end if
2194 
2195 end subroutine bubble_point_rachford_rice
2196 
2197 
2198 
2199 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2200 ! subroutine T-x rachford_rice
2201 !
2202 ! This subroutine solves the a T-x flash problem, i.e. for specified T and x of
2203 ! species index 2 in the phase 1 (usually liquid).
2204 !
2205 ! input: xiF(1:ncomp), xi(1,1:ncomp), xi(2,1:ncomp), dense(1), dense(2), t, p
2206 !
2207 ! output: rhoi1( 1:ncomp ), rhoi2( 1:ncomp )
2208 ! xi( 1, 1:ncomp ), xi( 2, 1:ncomp )
2209 ! values of: dense(:), lnphi(:,:) are also converged values, if converg=1
2210 ! ph_frac (phase fraction) could be defined as output, but is not yet.
2211 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2212 
2213 subroutine t_x_rachford_rice ( x12, converg, rhoi1, rhoi2 )
2214 
2215  use basic_variables
2216  implicit none
2217 
2218  !-----------------------------------------------------------------------------
2219  real, intent(in) :: x12
2220  integer, intent(in out) :: converg
2221  real, dimension(nc), intent(out) :: rhoi1
2222  real, dimension(nc), intent(out) :: rhoi2
2223 
2224  !-----------------------------------------------------------------------------
2225  integer :: count
2226  real :: f, fr, dfdx
2227  real :: delta
2228  real :: p0
2229  !-----------------------------------------------------------------------------
2230 
2231  delta = 0.00001
2232 
2233 
2234  count = 0
2235  f = 1.0
2236  do while ( abs( f ) > 1.e-6 .AND. count < 10 )
2237 
2238  count = count + 1
2239  p0 = p
2240 
2241  p = p0 * ( 1.0 + delta )
2242  call rachford_rice ( converg, rhoi1, rhoi2 )
2243  fr = xi(1,2) - x12
2244 
2245  p = p0
2246  call rachford_rice ( converg, rhoi1, rhoi2 )
2247  f = xi(1,2) - x12
2248 
2249  dfdx = ( fr - f ) / ( p0 * delta )
2250 
2251  p = p0 - f / dfdx
2252  !write (*,*) f, p0
2253 
2254  end do
2255 
2256 end subroutine t_x_rachford_rice
2257 
2258 
2259 
2260 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2261 !
2262 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2263 
2264 subroutine occupy_val_init ( rhoi1, rhoi2 )
2265 
2266  use basic_variables
2267  use eos_variables, only: pi, dhs, mseg
2268  implicit none
2269 
2270  !-----------------------------------------------------------------------------
2271  real, dimension(nc), intent(in) :: rhoi1
2272  real, dimension(nc), intent(in) :: rhoi2
2273 
2274  !-----------------------------------------------------------------------------
2275  integer :: i, ph
2276  !-----------------------------------------------------------------------------
2277 
2278  dense(1) = pi/6.0 * sum( rhoi1(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
2279  dense(2) = pi/6.0 * sum( rhoi2(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
2280 
2281  xi( 1, 1:ncomp ) = rhoi1( 1:ncomp ) / sum( rhoi1(1:ncomp) )
2282  xi( 2, 1:ncomp ) = rhoi2( 1:ncomp ) / sum( rhoi2(1:ncomp) )
2283  lnx( 1, 1:ncomp ) = log( xi( 1, 1:ncomp ) )
2284  lnx( 2, 1:ncomp ) = log( xi( 2, 1:ncomp ) )
2285 
2286  val_init(1) = dense(1)
2287  val_init(2) = dense(2)
2288  val_init(3) = t
2289  val_init(4) = p
2290  do ph = 1, 2
2291  do i = 1, ncomp
2292  val_init(4+i+(ph-1)*ncomp) = lnx( ph, i )
2293  end do
2294  end do
2295 
2296 end subroutine occupy_val_init
2297 
2298 
2299 
2300 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2301 !
2302 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2303 
2304 subroutine read_mode_staring_value ( scan_index, outp )
2305 
2306  use utilities, only: file_open
2307  implicit none
2308 
2309  !-----------------------------------------------------------------------------
2310  integer, intent(out) :: scan_index
2311  integer, intent(in out) :: outp
2312 
2313  !-----------------------------------------------------------------------------
2314  character (LEN=50) :: textstring
2315  character (LEN=50) :: start_value_file
2316  !-----------------------------------------------------------------------------
2317 
2318  first_call = .false.
2319  start_value_file = './input_file/starting_value_default.inp'
2320  call file_open ( start_value_file, 84 )
2321  read (84,*) textstring
2322  read (84,*) scan_index, outp
2323  close (84)
2324  if ( outp == 3 ) then
2325  if ( scan_index == 1 ) outp = 0
2326  if ( scan_index == 2 ) outp = 1
2327  end if
2328  if ( outp >= 1 ) write (*,*) 'level of info given to terminal:',outp
2329  write (*,*) 'calculation option composition-scan:',scan_index
2330 
2331 end subroutine read_mode_staring_value
2332 
2333 
2334 
2335 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2336 !
2337 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2338 
2339 subroutine output_rho_info ( rhoi )
2340 
2341  use eos_variables, only: ncomp, nc, pi, dhs, mseg
2342  implicit none
2343  real, dimension(nc), intent(IN) :: rhoi
2344  real :: dens
2345  !-----------------------------------------------------------------------------
2346 
2347  dens = pi/6.0 * sum( rhoi(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
2348 
2349  write (*,*) ' eta = ', dens
2350  write (*,*) ' x = ', rhoi( 1:ncomp ) / sum( rhoi(1:ncomp) )
2351 
2352 end subroutine output_rho_info
2353 
2354 
2355 
2356 end module starting_values
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
starting_values
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module STARTING_VALUES This m...
Definition: modules.f90:251
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29
eos_numerical_derivatives
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:220
starting_values::start_var
subroutine, public start_var(converg)
subroutine start_var
Definition: starting_value.f90:65