improved solvers a bit

src/equilibria/boundaries/phase_envelopes_pt.f90

@@ -107,7 +107,7 @@ function pt_envelope_2ph(&
 
       where(z == 0)
          X(:nc) = 0
-      end where   
+      end where
 
       X(nc+1) = log(first_point%T)
       X(nc+2) = log(first_point%P)
@@ -197,7 +197,7 @@ subroutine foo(X, ns, S, F, dF, dFdS)
          df(:nc, nc + 2) = P * (dlnphi_dp_y - dlnphi_dp_z)
 
          df(nc + 1, :nc) = y
-         
+
          do i=1,nc
             if (z(i) == 0) then
                F(i) = 0
@@ -254,18 +254,23 @@ subroutine update_spec(X, ns, S, dS, dXdS, step_iters)
             ] &
             )
 
-         do while(abs(dXdS(nc+1)*dS) < 0.01 .and. abs(dXdS(nc+2)*dS) < 0.01)
+         ! Avoid small steps on T or P
+         do while(&
+            abs(dXdS(nc+1)*dS) < 0.05 &
+            .and. abs(dXdS(nc+2)*dS) < 0.05 &
+            .and. dS /= 0)
             dS = dS * 1.1
          end do
 
+         ! Dont make big steps in compositions
          do while(maxval(abs(dXdS(:nc)*dS)) > 0.1 * maxval(abs(X(:nc))))
             dS = 0.7*dS
          end do
 
          if (present(maximum_pressure)) then
             if (X(nc+2) > log(maximum_pressure)) dS = 0
          end if
-         
+
          call save_point(X, step_iters)
          call detect_critical(X, dXdS, ns, S, dS)
       end subroutine update_spec
@@ -340,7 +345,7 @@ subroutine detect_critical(X, dXdS, ns, S, dS)
 
          inner = 0
          do while (&
-                  maxval(abs(X(:nc))) < 0.01 &
+            maxval(abs(X(:nc))) < 0.07 &
             .and. inner < 5000)
             ! If near a critical point, jump over it
             inner = inner + 1
@@ -427,13 +432,13 @@ type(PTEnvel2) function find_hpl(model, z, T0, P0)
       !! of the same component in a pure phase. This is done for each component
       !! in the mixture. The component with the highest temperature is selected
       !! as it should be the first one appearing. If all components have a
-      !! negative difference then the mixture is probably stable at all 
+      !! negative difference then the mixture is probably stable at all
       !! temperatures.
       class(ArModel), intent(in) :: model !! Equation of state model
       real(pr), intent(in) :: z(:) !! Mole fractions
       real(pr), intent(in) :: T0 !! Initial temperature [K]
       real(pr), intent(in) :: P0 !! Search pressure [bar]
-      
+
       integer :: i
       real(pr) :: y(size(z))
       real(pr) :: lnphi_y(size(z)), lnphi_z(size(z))
@@ -447,13 +452,13 @@ type(PTEnvel2) function find_hpl(model, z, T0, P0)
       do ncomp=1,nc
          y = 0
          y(ncomp) = 1
-         
+
          do i=int(T0), 1, -10
             T = real(i, pr)
             call model%lnphi_pt(z, P, T, root_type="liquid", lnPhi=lnphi_z)
             call model%lnphi_pt(y, P, T, root_type="liquid", lnPhi=lnphi_y)
 
-            ! Fugacity of the component ncomp 
+            ! Fugacity of the component ncomp
             ! z * phi_i_mixture / phi_i_pure
             ! if eq > 1 then the fugacity in the mixture is above the pure,
             ! so the component is more stable on another phase
@@ -483,6 +488,6 @@ type(PTEnvel2) function find_hpl(model, z, T0, P0)
       find_hpl = pt_envelope_2ph( &
          model, z, fr, &
          specified_variable_0=nc+2, delta_0=-5.0_pr, iterations=1000)
-   end function
+   end function find_hpl
 
 end module yaeos__equilibria_boundaries_phase_envelopes_pt

src/equilibria/critical.f90

@@ -159,7 +159,7 @@ subroutine foo(X, ns, S, F, dF, dFdS)
 
          real(pr) :: z(size(u)), Xsol(3)
 
-         if (X(1) > 1) then
+         if (X(CPSpec%a) > 1) then
             return
          end if
 
@@ -179,12 +179,10 @@ subroutine update_specification(X, ns, S, dS, dXdS, iterations)
          real(pr), intent(in out) ::  dXdS(:) !! \(\frac{dX}{dS}\)
          integer, intent(in) :: iterations !! Iterations needed to converge point
 
-
-
          ns = maxloc(abs(dXdS), dim=1)
          dS = dXdS(ns)*dS
          dXdS = dXdS/dXdS(ns)
-         if (exp(X(4)) > max_P) then
+         if (exp(X(CPSpec%P)) > max_P) then
             dS = 0
          end if
 

src/equilibria/saturations_points.f90

@@ -208,7 +208,7 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0
       if (present (t0)) then
          t = t0
       else
-         t = 150._pr
+         t = 250._pr
       end if
 
       if (present(y0)) then
@@ -259,7 +259,7 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0
 
          if (.not. ieee_is_finite(step) .or. ieee_is_nan(step)) exit
 
-         do while (abs(step) > 0.5*T .or. T - step < 0)
+         do while (T - step < 0)
             if (isnan(step)) step = 10
             step = step/2
          end do
@@ -269,8 +269,7 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0
          if (abs(step) < tol .and. abs(f) < tol) exit
       end do
       ! ========================================================================
-      !if (its >= iters_first_step) then
-      if (.true.) then
+      if (its >= iters_first_step) then
          block
             real(pr) :: X(size(n)+2), S
             integer :: ns, nc

src/math/continuation.f90

@@ -165,6 +165,8 @@ function continuation(&
             if (stop(X, ns, S, dS, dXdS, newton_its)) exit
          end if
 
+         if (dS == 0) exit
+
          X = X + dXdS * dS
          S = X(ns)
       end do

test/test_critical.f90

@@ -1,21 +1,29 @@
 program main
    !! Test the calculation of critical lines
    use yaeos
+   use yaeos__math, only: interpol
+   use fortime, only: Timer
    implicit none
+
+   logical :: WRITE_FILES=.true.
 
    integer, parameter :: nc=12
+   integer :: a_nearest
+   integer :: npt=6
 
    type(CubicEoS) :: model
    type(EquilibriumState) :: sat, crit
    type(PTEnvel2) :: env
    type(CriticalLine) :: cl
+   type(Timer) :: tim
 
-   real(pr) ::  V, T, P, a
+   real(pr) ::  V, T, P, a, da
 
    real(pr) :: z(nc)
    real(pr) :: z0(nc)
    real(pr) :: zi(nc)
 
+
    integer :: i
 
    model = get_model()
@@ -25,8 +33,10 @@ program main
    z = a*zi + (1-a)*z0
 
    ! Get the full phase envelope of the fluid
-   sat = saturation_temperature(model, z, P=0.0001_pr, kind="dew")
+   call tim%timer_start()
+   sat = saturation_temperature(model, z, P=0.1_pr, kind="dew")
    env = pt_envelope_2ph(model, z, sat, maximum_pressure=1000._pr)
+   call tim%timer_stop()
 
    ! Calculate the critical point
    T = sum(model%components%Tc * z)
@@ -41,16 +51,33 @@ program main
    end if
 
    ! Now test the critical lines
+   print *, "CL"
+   call tim%timer_start()
    cl = critical_line(model, a0=a, z0=z0, zi=zi, dS0=0.1_pr, max_points=5000)
+   call tim%timer_stop()
 
-
-   do i=1,5
-      a = cl%a(i)
+   if (WRITE_FILES) call write_cl
+
+
+   da = (cl%a(size(cl%a)) - cl%a(1))/(npt+1)
+   do i=1,npt
+      a = cl%a(1) + da*i
       z = a*zi + (1-a)*z0
+      call tim%timer_start()
       sat = saturation_temperature(model, z, P=0.1_pr, kind="dew")
-      env = pt_envelope_2ph(model, z, sat, maximum_pressure=1000._pr)
+      env = pt_envelope_2ph(model, z, sat, maximum_pressure=2000._pr, points=1000, delta_0=1.5_pr)
+      print *, "Envelope", i, size(env%points)
+      call tim%timer_stop()
+      if (WRITE_FILES) call write_env
+
+      a_nearest = minloc(abs(cl%a - a), dim=1)
 
-      if (sum(([cl%T(i), cl%P(i)] - [env%cps(1)%T, env%cps(1)%P]))**2 > 1e-2) then
+      T = interpol(cl%a(a_nearest), cl%a(a_nearest+1), cl%T(a_nearest), cl%T(a_nearest+1), a)
+      P = interpol(cl%a(a_nearest), cl%a(a_nearest+1), cl%P(a_nearest), cl%P(a_nearest+1), a)
+
+      if (maxval(([T, P] - [env%cps(1)%T, env%cps(1)%P])) > 10) then
+         write(*, *) [T, P] 
+         write(*, *) [env%cps(1)%T, env%cps(1)%P]
          error stop "Critical line failed"
       end if
    end do
@@ -63,17 +90,35 @@ type(CubicEoS) function get_model()
       zi=[1.0,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
 
       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]
+         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]
+         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
+
+   subroutine write_cl
+      do i=1,size(cl%a)
+         write(1, *) cl%a(i), cl%T(i), cl%P(i)
+      end do
+      write(1, *)
+      write(1, *)
+   end subroutine
+
+   subroutine write_env
+      integer :: i
+      do i=1,size(env%points)
+         write(1, *) a, env%points(i)%T, env%points(i)%P
+      end do
+      write(1, *)
+      write(1, *)
+      call flush(1)
+   end subroutine write_env
 end program main