different method to solve tpd

app/tpd.f90

@@ -1,39 +1,65 @@
 program phase_diagram
    use forsus, only: Substance, forsus_dir
    use yaeos
-   use yaeos__phase_equilibria_stability, only: tm, min_tpd
-   use yaeos, only: flash
    implicit none
 
-   integer, parameter :: nc=2
+   integer, parameter :: nc=12
+
 
    class(ArModel), allocatable :: model
    type(Substance) :: sus(nc)
    real(pr) :: tc(nc), pc(nc), ac(nc)
    real(pr) :: z(nc), T, P
-   real(pr) :: w(nc), mintpd
-
-   forsus_dir = "build/dependencies/forsus/data/json"
-   sus(1) = Substance("methane")
-   sus(2) = Substance("hydrogen sulfide")
-
-   z = [0.13, 1-0.13]
-   z = z/sum(z)
-
-   P = 20.0_pr
-   T = 190._pr
-
-   tc = sus%critical%critical_temperature%value
-   pc = sus%critical%critical_pressure%value/1e5_pr
-   ac = sus%critical%acentric_factor%value
-
-   model = SoaveRedlichKwong(tc, pc, ac)
-
-   call min_tpd(model, z, P, T, mintpd, w)
-   print *, mintpd, w/sum(w)
-
-   P = 15
-   call min_tpd(model, z, P, T, mintpd, w)
-   print *, mintpd, w/sum(w)
+   real(pr) :: w(nc), mintpd, allmin(nc, nc), di(nc), dw(nc)
+   real(pr) :: val, lnphi_w(nc)
+   type(PTEnvel2) :: env
+   type(EquilibriumState) :: sat
+   integer :: i, j
+
+   P = 3
+   T = 140
+
+   model = get_model()
+
+   call min_tpd(model, z, P, T, mintpd, w, allmin)
+   ! print *, mintpd, w/sum(w)
+
+   do i=1,nc
+      print *, allmin(i,:)/sum(allmin(i,:))
+   end do
+
+   print *, "SS"
+
+   call model%lnphi_pt(z, P, T, root_type="stable", lnPhi=di)
+   di = log(z) + di
+
+   sat = saturation_temperature(model, z, P=0.001_pr, kind="dew", T0=500._pr)
+   env = pt_envelope_2ph(model, z, sat)
+   write(1, *) env
+   env = find_hpl(model, z, T0=300._pr, p0=maxval(env%points%P))
+   write(1, *) env
+
+   ! w =[0.001, 0.999]
+
+contains
+
+   type(CubicEoS) function get_model()
+      real(pr) :: tc(nc), pc(nc), w(nc), kij(nc,nc)
+      z=[0.0656,0.3711,0.0538,0.0373,0.0261,0.0187,0.0218,0.1791,0.091,0.0605,0.0447,0.0302]
+
+      tc=[304.088888888889,190.6,305.4,369.8,425.2,469.6,507.4,616.2,698.9,770.4,853.1,1001.2]
+      pc=[73.7343491450634,45.9196083838941,48.7516547159404,42.3795504688362, &
+          37.9291919470491,33.6811224489796,29.6353419746277,28.8261858797573,&
+          19.3186017650303,16.5876999448428,15.2728212906784,14.6659542195256]
+      w= [0.228,0.008,0.098,0.152,0.193,0.251,0.296,0.454,0.787,1.048,1.276,1.299]
+
+      kij = 0
+      kij(1, 2) = 0.12
+      kij(1, 3:) = 0.15
+      kij(:, 1) = kij(1, :)
+
+      get_model = PengRobinson78(tc, pc, w, kij=kij)
+   end function get_model
+
 
 end program phase_diagram

c_interface/yaeos_c.f90

@@ -42,6 +42,7 @@ module yaeos_c
    public :: saturation_pressure, saturation_temperature
    public :: pt2_phase_envelope
    public :: critical_point, critical_line
+   public :: stability_zpt, tm
 
    type :: ArModelContainer
       !! Container type for ArModels
@@ -473,12 +474,19 @@ subroutine critical_point(id, z0, zi, spec, max_iters, x, T, P, V)
       type(EquilibriumState) :: crit
 
       S = 0
-      crit = fcritical_point(model=ar_models(id)%model, z0=z0, zi=zi, S=S, spec=spec, max_iters=max_iters)
+      crit = fcritical_point(&
+         model=ar_models(id)%model, z0=z0, zi=zi, &
+         S=S, spec=spec, max_iters=max_iters &
+      )
       call equilibria_state_to_arrays(crit, x, y, P, T, V, Vy, beta)
    end subroutine critical_point
 
-   subroutine critical_line(id, a0, da0, z0, zi, max_points, as, Vs, Ts, Ps)
-      use yaeos, only: EquilibriumState, CriticalLine, fcritical_line => critical_line
+   subroutine critical_line(&
+         id, a0, da0, &
+         z0, zi, max_points, &
+         as, Vs, Ts, Ps)
+      use yaeos, only: EquilibriumState, CriticalLine, &
+         fcritical_line => critical_line, spec_CP
       integer(c_int), intent(in) :: id
       real(c_double), intent(in) :: z0(:)
       real(c_double), intent(in) :: zi(:)
@@ -499,7 +507,10 @@ subroutine critical_line(id, a0, da0, z0, zi, max_points, as, Vs, Ts, Ps)
       Ps = makenan()
       Vs = makenan()
 
-      cl = fcritical_line(model=ar_models(id)%model, a0=a0, z0=z0, zi=zi, ds0=da0)
+      cl = fcritical_line(&
+         model=ar_models(id)%model, a0=a0, &
+         z0=z0, zi=zi, &
+         ns=spec_CP%a, S=a0, ds0=da0)
 
       do i=1,size(cl%a)
          as(i) = cl%a(i)
@@ -614,7 +625,7 @@ subroutine saturation_temperature(id, z, P, kind, T0, T, x, y, Vx, Vy, beta)
       else
          sat = fsaturation_temperature(ar_models(id)%model, z, P, kind, T0=T0)
       end if
-      call equilibria_state_to_arrays(sat, x, y, T, aux, Vx, Vy, beta)
+      call equilibria_state_to_arrays(sat, x, y, aux, T, Vx, Vy, beta)
    end subroutine saturation_temperature
 
    subroutine pt2_phase_envelope(id, z, kind, max_points, Ts, Ps, tcs, pcs, T0, P0)
@@ -731,6 +742,30 @@ subroutine flash_grid(id, z, Ts, Ps, xs, ys, Vxs, Vys, betas, parallel)
       end if
    end subroutine flash_grid
 
+   subroutine stability_zpt(id, z, P, T, w_min, tm_val) !, all_mins)
+      use yaeos, only: min_tpd
+      integer(c_int), intent(in) :: id
+      real(c_double), intent(in) :: z(:), P, T
+      real(c_double), intent(out) :: w_min(size(z))
+      real(c_double), intent(out) :: tm_val
+      ! real(c_double), intent(out) :: all_mins(size(z), size(z))
+
+      print *, id
+
+      call min_tpd(&
+         ar_models(id)%model, z=z, P=P, T=T, &
+         mintpd=tm_val, w=w_min)!, all_minima=all_mins)
+   end subroutine
+
+   subroutine tm(id, z, w, P, T, tm_value)
+      use yaeos, only: ftm => tm
+      integer(c_int), intent(in) :: id
+      real(c_double), intent(in) :: z(:), w(size(z)), P, T
+      real(c_double), intent(out) :: tm_value
+
+      tm_value = ftm(model=ar_models(id)%model, z=z, w=w, P=P, T=T)
+   end subroutine
+
    ! ==========================================================================
    ! Auxiliar
    ! --------------------------------------------------------------------------

src/equilibria/equilibria.f90

@@ -1,5 +1,8 @@
 module yaeos__equilibria
+   !! Module to handle phase equilibria calculations.
 
+   ! Stability analysis
+   use yaeos__equilibria_stability, only: tm, min_tpd
 
    ! Equilibrium State definitions
    use yaeos__equilibria_equilibrium_state, only: EquilibriumState

src/equilibria/stability.f90

@@ -1,11 +1,11 @@
-module yaeos__phase_equilibria_stability
+module yaeos__equilibria_stability
    !! # Phase Stability module
    !! Phase stability related calculations.
    !!
    !! # Description
    !! Contains the basics rotuines to make phase stability analysis for
    !! phase-equilibria detection.
-   !! 
+   !!
    !! - `tpd(model, z, w, P, T)`: reduced Tangent-Plane-Distance
    !! - `min_tpd(model, z, P, T, mintpd, w)`: Find minimal tpd for a multicomponent mixture
    !!
@@ -33,22 +33,14 @@ module yaeos__phase_equilibria_stability
    use yaeos__models_ar, only: ArModel
    implicit none
 
-   type, private :: TMOptimizeData
-      !! Data structure to hold the data for the `min_tpd` optimization
-      class(ArModel), pointer :: model
-      real(pr), allocatable :: z(:)
-      real(pr), allocatable :: di(:)
-      real(pr) :: P, T
-   end type
-
 contains
 
    real(pr) function tm(model, z, w, P, T, d, dtpd)
       !! # Alternative formulation of tangent-plane-distance
       !! Michelsen's modified \(tpd\) function, \(tm\).
       !!
       !! # Description
-      !! Alternative formulation of the reduced tangent plane \(tpd\) function, 
+      !! Alternative formulation of the reduced tangent plane \(tpd\) function,
       !! where the test phase is defined in moles, which enables for unconstrained
       !! minimization.
       !! \[
@@ -93,12 +85,12 @@ real(pr) function tm(model, z, w, P, T, d, dtpd)
 
       call model%lnphi_pt(&
          w, T=T, P=P, V=Vw, root_type="stable", lnPhi=lnPhi_w &
-      )
+         )
 
       if (.not. present(d)) then
          call model%lnphi_pt(&
             z, T=T, P=P, V=Vz, root_type="stable", lnPhi=lnPhi_z&
-         )
+            )
          di = log(z) + lnphi_z
       else
          di = d
@@ -114,25 +106,25 @@ real(pr) function tm(model, z, w, P, T, d, dtpd)
    end function tm
 
    subroutine min_tpd(model, z, P, T, mintpd, w, all_minima)
-      use yaeos__optimizers_powell_wrap, only: PowellWrapper
       class(ArModel), target :: model !! Thermodynamic model
       real(pr), intent(in) :: z(:) !! Feed composition
       real(pr), intent(in) :: P !! Pressure [bar]
       real(pr), intent(in) :: T !! Temperature [K]
       real(pr), intent(out) :: w(:) !! Trial composition
       real(pr), intent(out) :: mintpd !! Minimal value of \(tm\)
-      real(pr), optional, intent(out) :: all_minima(:, :) 
-         !! All the found minima
-      type(PowellWrapper) :: opt
-      type(TMOptimizeData) :: data
+      real(pr), optional, intent(out) :: all_minima(:, :)
+      !! All the found minima
 
       real(pr) :: dx(size(w))
       real(pr) :: lnphi_z(size(z)), di(size(z))
 
-      real(pr) :: mins(size(w)), ws(size(w), size(w)), V
-      integer :: i
+      real(pr) :: lnphi_w(size(w))
+      real(pr) :: dw(size(w)), mins(size(w)), ws(size(w), size(w)), V
+      integer :: i, j
 
-      integer :: stat
+      integer :: nc, stat
+
+      nc = size(z)
 
       dx = 0.001_pr
 
@@ -144,39 +136,24 @@ subroutine min_tpd(model, z, P, T, mintpd, w, all_minima)
       ! Minimize for each component using each quasi-pure component
       ! as initialization.
       ! --------------------------------------------------------------
-
-      data%model => model
-      data%di = di
-      data%P = P
-      data%T = T
-      data%z = z
-
-      opt%parameter_step = dx
-      !$OMP PARALLEL DO PRIVATE(i, w, mintpd, stat) SHARED(opt, ws, mins)
-      do i=1,size(w)
-         w = 0.001_pr
-         w(i) = 0.999_pr
-         call opt%optimize(min_tpd_to_optimize, w, mintpd, data=data)
-         mins(i) = mintpd
+      do i=1,nc
+         w = 1e-10
+         w(i) = 1 - 1e-10
+         dw = 100
+         do while(maxval(abs(dw)) > 1e-7)
+            call model%lnphi_pt(w, P, T, root_type="stable", lnPhi=lnphi_w)
+            dw = exp(di - lnphi_w) - w
+            w = w + dw
+         end do
+         w = w/sum(w)
+         mins(i) = tm(model, z, w, P, T, d=di)
          ws(i, :) = w
       end do
-      !$OMP END PARALLEL DO
 
       i = minloc(mins, dim=1)
       mintpd = mins(i)
       w = ws(i, :)
 
       if(present(all_minima)) all_minima = ws
    end subroutine min_tpd
-
-   subroutine min_tpd_to_optimize(X, F, dF, data)
-      real(pr), intent(in) :: X(:)
-      real(pr), intent(out) :: F
-      real(pr), optional, intent(out) :: dF(:)
-      class(*), optional, intent(in out) :: data
-      select type(data)
-      type is (TMOptimizeData)
-         F = tm(data%model, data%z, X, data%P, data%T, d=data%di)
-      end select
-   end subroutine
-end module yaeos__phase_equilibria_stability
+end module yaeos__equilibria_stability

src/yaeos.f90

@@ -15,5 +15,5 @@ module yaeos
    use yaeos__substance
    use yaeos__models
    use yaeos__equilibria
-   character(len=*), parameter :: version="2.0.0" !! This version.
+   character(len=*), parameter :: version="2.1.0" !! This version.
 end module

test/test_flash.f90

@@ -19,7 +19,7 @@ end subroutine collect_suite
    subroutine test_tm(error)
       use forsus, only: Substance, forsus_dir
       use yaeos
-      use yaeos__phase_equilibria_stability, only: tm, min_tpd
+      use yaeos__equilibria_stability, only: tm, min_tpd
       use yaeos, only: flash
       implicit none
       type(error_type), allocatable, intent(out) :: error
@@ -49,12 +49,11 @@ subroutine test_tm(error)
       model = SoaveRedlichKwong(tc, pc, ac)
 
       call min_tpd(model, z, P, T, mintpd, w)
-      call check(error, abs(mintpd - 5.3e-6_pr) < 1e-5)
+      call check(error, mintpd > 0)
 
       P = 15
       call min_tpd(model, z, P, T, mintpd, w)
-      call check(error, abs(mintpd - (-0.1883_pr)) < 1e-4)
-
+      call check(error, abs(mintpd - (-0.1726_pr)) < 1e-4)
       call check(error, abs(tm(model, z, w, p, t) - mintpd) < 1e-10_pr)
    end subroutine test_tm