fixed circular dependency problem by making the calculation of satura…

…tion pressure of pure compounds an `ArModel` method

src/equilibria/equilibria.f90

@@ -7,9 +7,6 @@ module yaeos__equilibria
    ! Equilibrium State definitions
    use yaeos__equilibria_equilibrium_state, only: EquilibriumState
 
-   ! Pure component saturation pressure
-   use yaeos__equilibria_pure_psat, only: Psat
-
    ! Phase split calculations
    use yaeos__equilibria_flash, only: flash
 

src/models/residual_helmholtz/ar_models.f90

@@ -59,6 +59,7 @@ module yaeos__models_ar
       procedure :: entropy_residual_vt
       procedure :: Cv_residual_vt
       procedure :: Cp_residual_vt
+      procedure :: Psat_pure
    end type ArModel
 
    interface size
@@ -429,7 +430,7 @@ subroutine lnfug_vt(eos, &
       RT = R*T
 
       if (present(lnf) .and. .not. (&
-              present(dlnfdn) &
+         present(dlnfdn) &
          .or. present(dlnfdV) &
          .or. present(dlnfdT) &
          )) then
@@ -440,7 +441,7 @@ subroutine lnfug_vt(eos, &
          where (n /= 0)
             lnf = log(n/totn) + Arn/RT - log(V/(totn*RT))
          endwhere
-         
+
          if (present(P)) P = P_in
 
          return
@@ -471,7 +472,7 @@ subroutine lnfug_vt(eos, &
       end if
 
       if (present(dlnfdV)) then
-         dlnfdV = -dPdn_in/RT         
+         dlnfdV = -dPdn_in/RT
       end if
 
       if (present(dlnfdT)) then
@@ -605,4 +606,42 @@ subroutine Cp_residual_vt(eos, n, V, T, Cp)
 
       Cp = Cv - T*dPdT**2/dPdV - totn*R
    end subroutine Cp_residual_vt
+
+   real(pr) function Psat_pure(eos, ncomp, T)
+      !! Calculation of saturation pressure of a pure component using the
+      !! secant method.
+      class(ArModel), intent(in) :: eos !! Model that will be used
+      integer, intent(in) :: ncomp
+      !! Number of component in the mixture from which the saturation pressure
+      !! will be calculated
+      real(pr), intent(in) :: T !! Temperature [K]
+
+      real(pr) :: P1, P2
+      real(pr) :: f1, f2
+
+      real(pr) :: n(size(eos))
+
+      n = 0
+      n(ncomp) = 1
+
+      P1 = 0.5
+      P2 = 1
+
+      do while(abs(diff(P2)) > 1e-5)
+         f1 = diff(P1)
+         f2 = diff(P2)
+         Psat_pure = (P1 * f2 - P2 * f1)/(f2 - f1)
+         P1 = P2
+         P2 = Psat_pure
+      end do
+   contains
+      real(pr) function diff(P)
+         real(pr), intent(in) :: P
+         real(pr) :: V_l, V_v
+         real(pr) :: phi_v(size(eos)), phi_l(size(eos))
+         call eos%lnphi_pt(n, P=P, T=T, V=V_v, lnPhi=phi_v, root_type="vapor")
+         call eos%lnphi_pt(n, P=P, T=T, V=V_l, lnPhi=phi_l, root_type="liquid")
+         diff = phi_v(ncomp) - phi_l(ncomp)
+      end function diff
+   end function Psat_pure
 end module yaeos__models_ar

src/models/residual_helmholtz/cubic/implementations/implementations.f90

@@ -295,7 +295,6 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model)
       !!
       !! \[\delta_1 = d_1 + d_2 (d_3 - Z_c)^d_4 + d_5 (d_3 - Z_c) ^ d_6\]
       !! \[k = (A_1  Z_c + A_0)\omega^2 + (B_1 Z_c + B_0)\omega + (C_1 Z_c + C_0)\]
-      use yaeos__equilibria, only: Psat
       use yaeos__models_ar_cubic_quadratic_mixing, only: QMR_RKPR
       use yaeos__models_ar_cubic_alphas, only: AlphaRKPR
       real(pr), intent(in) :: tc(:) !! Critical Temperature [K]
@@ -378,7 +377,7 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model)
          do i=1,nc
             diff = 1
             do while (abs(diff) > 1e-6)
-               Psat_i = Psat(model, i, 0.7*Tc(i))
+               Psat_i = model%Psat_pure(i, 0.7*Tc(i))
                diff = (w(i) - (-1 - log10(Psat_i/Pc(i))))
                alpha%k(i) = alpha%k(i) + 0.1*diff
                model%alpha = alpha

src/models/solvers/saturation_point.f90

@@ -6,13 +6,14 @@ module yaeos__m_s_sp
 
 contains
 
-   subroutine saturation_F(model, z, X, ns, S, F, dF)
+   subroutine saturation_F(model, z, X, ns, S, F, dF, dPdVz, dPdVy)
       class(ArModel), intent(in) :: model
       real(pr), intent(in) :: X(:)
       integer, intent(in) :: ns
       real(pr), intent(in) :: S
       real(pr), intent(out) :: F(:)
       real(pr), optional, intent(out) :: dF(:, :)
+      real(pr), intent(out) :: dPdVz, dPdVy
 
       ! Variables
       real(pr) :: T, Vz, Vy
@@ -22,14 +23,14 @@ subroutine saturation_F(model, z, X, ns, S, F, dF)
       real(pr) :: lnfug_z(size(model)), dlnfug_dn_z(size(model), size(model))
       real(pr) :: dlnfug_dT_z(size(model)), dlnfug_dV_z(size(model))
       real(pr) :: dlnfug_dP_z(size(model))
-      real(pr) :: Pz, dPdTz, dPdVz, dPdn_z(size(z))
+      real(pr) :: Pz, dPdTz, dPdn_z(size(z))
 
       ! incipient phase variables
       real(pr) :: y(size(z))
       real(pr) :: lnfug_y(size(model)), dlnfug_dn_y(size(model), size(model))
       real(pr) :: dlnfug_dT_y(size(model)), dlnfug_dV_y(size(model))
       real(pr) :: dlnfug_dP_y(size(model))
-      real(pr) :: Py, dPdTy, dPdVy, dPdn_y(size(z))
+      real(pr) :: Py, dPdTy, dPdn_y(size(z))
 
       real(pr) :: lnPspec
 
@@ -196,23 +197,28 @@ subroutine solve_TP(model, kind, z, X, ns, S, tol, max_iterations, its)
       integer, intent(in) :: max_iterations
       integer, intent(out) :: its
 
-      integer :: innits
+      integer :: nc
 
       real(pr) :: F(size(X))
       real(pr) :: dF(size(X), size(X))
       real(pr) :: dFdS(size(X))
       real(pr) :: dx(size(X))
 
+      nc = size(X) - 2
 
-
-      its= 0
+      its = 0
+      dX = 1
       do while (its < max_iterations)
+         its = its + 1
          call saturation_TP(model=model, z=z, kind=kind, X=X, ns=ns, S=S, F=F, dF=dF, dFdS=dFdS)
 
          dX = solve_system(dF, -F)
+         do while (maxval(abs(dX)) > 1)
+            dX = dX / 2
+         end do
+
          X = X + dX
          if (all(abs(F) < tol)) exit
-         its = its + 1
       end do
 
    end subroutine solve_TP
@@ -227,6 +233,8 @@ subroutine solve_VxVyT(model, z, X, ns, S, tol, max_iterations, its)
       real(pr), intent(in) :: tol
       integer, intent(in) :: max_iterations
       integer, intent(out) :: its
+
+      real(pr) :: dPdVz, dPdVy
 
       integer :: nc
       real(pr) :: F(size(X))
@@ -238,7 +246,7 @@ subroutine solve_VxVyT(model, z, X, ns, S, tol, max_iterations, its)
 
       its = 0
       do while (its < max_iterations)
-         call saturation_F(model, z, X, ns, S, F, dF)
+         call saturation_F(model, z, X, ns, S, F, dF, dPdVz, dPdVy)
          if (all(abs(F) < tol)) exit
 
          dX = solve_system(dF, -F)

test/test_implementations/ar_models/cubics/test_rkpr.f90

@@ -78,10 +78,10 @@ subroutine test_rkpr_cons_mixture(error)
 
       ! Consistency tests
       call check(error, abs(eq31) <= 1e-13)
-      call check(error, maxval(abs(eq33)) < 1e-14)
-      call check(error, maxval(abs(eq34)) < 1e-14)
-      call check(error, abs(eq36) <= 1e-14)
-      call check(error, abs(eq37) <= 1e-14)
+      call check(error, maxval(abs(eq33)) < 1e-13)
+      call check(error, maxval(abs(eq34)) < 1e-13)
+      call check(error, abs(eq36) <= 1e-13)
+      call check(error, abs(eq37) <= 1e-13)
 
       ! ========================================================================
       ! Model with kij and lij

test/test_saturation.f90

@@ -183,7 +183,7 @@ subroutine test_px2_envelope(error)
    end subroutine test_px2_envelope
 
    subroutine test_pure_psat(error)
-      use yaeos, only: pr, ArModel, Psat
+      use yaeos, only: pr, ArModel
       use fixtures_models, only: binary_PR76
       type(error_type), allocatable, intent(out) :: error
       class(ArModel), allocatable :: model
@@ -196,7 +196,7 @@ subroutine test_pure_psat(error)
       Psats_val = [260.37450286310201, 30.028551527997834]
 
          do i=1,2
-            Psats(i) = Psat(model, i, T)
+            Psats(i) = model%Psat_pure(i, T)
          end do
    ! call check(error, maxval(abs(Psats-Psats_val)) < abs_tolerance)
    end subroutine test_pure_psat