APIs for external libraries

app/fit.f90

@@ -1,107 +1,9 @@
-module yaeos__fitting_fit_nrtl
-   use yaeos__constants, only: pr
-   use yaeos__fitting, only: FittingProblem
-   use yaeos__models, only: ArModel, NRTL, CubicEoS, MHV
-   use forsus, only: Substance
-   implicit none
-
-   integer, parameter :: nc = 2
-
-   type, extends(FittingProblem) :: FitMHVNRTL
-   contains
-      procedure :: get_model_from_X => model_from_X
-   end type
-
-contains
-
-   subroutine init_model(problem, sus)
-      use yaeos, only: R, ArModel, CubicEoS, PengRobinson78, RKPR, SoaveRedlichKwong
-      class(FitMHVNRTL), intent(in out) :: problem
-      type(Substance), intent(in) :: sus(2)
-      type(MHV) :: mixrule
-      type(NRTL) :: ge
-      real(pr) :: tc(nc), pc(nc), w(nc), vc(nc), zc(nc)
-      real(pr) :: a(nc, nc), b(nc, nc), c(nc, nc), bi(nc)
-
-      a=0; b=0; c=0
-
-      tc = sus%critical%critical_temperature%value
-      pc = sus%critical%critical_pressure%value/1e5
-      w = sus%critical%acentric_factor%value
-      vc = sus%critical%critical_volume%value
-      zc = pc*vc/(R*tc)
-
-      ge = NRTL(a, b, c)
-
-      ! problem%model = RKPR(tc, pc, w, zc)
-      allocate(CubicEoS :: problem%model)
-      problem%model = SoaveRedlichKwong(tc, pc, w)
-
-      associate(m => problem%model)
-         select type(m)
-         type is (CubicEoS)
-            bi = m%b
-            mixrule = MHV(ge=ge, q=-0.593_pr, b=bi)
-            deallocate(m%mixrule)
-            m%mixrule = mixrule
-         end select
-      end associate
-   end subroutine init_model
-
-   function model_from_X(problem, X) result(model)
-      use yaeos, only: R, RKPR, PengRobinson78, ArModel, QMR, CubicEoS
-      use yaeos__models_ar_cubic_quadratic_mixing, only: RKPR_D1mix
-      real(pr), intent(in) :: X(:)
-      class(FitMHVNRTL), intent(in) :: problem
-      class(ArModel), allocatable :: model
-      type(NRTL) :: ge
-
-      real(pr) :: a(nc, nc), b(nc, nc), c(nc, nc)
-
-      a=0; b=0; c=0
-
-      a(1, 2) = x(1)
-      a(2, 1) = x(2)
-
-      b(1, 2) = x(3)
-      b(2, 1) = x(4)
-
-      c(1, 2) = x(5)
-      c(2, 1) = x(6)
-
-      ge = NRTL(a, b, c)
-
-      associate (pm => problem%model)
-         select type(pm)
-         type is (CubicEoS)
-            model = pm
-         end select
-      end associate
-
-      ! model = problem%model
-
-      select type(model)
-      class is (CubicEoS)
-         associate(mr => model%mixrule)
-            select type (mr)
-            class is (MHV)
-               mr%l(1, 2) = x(7)
-               mr%l(2, 1) = x(7)
-               mr%ge = ge
-               model%del1 = x(8:)
-               model%del2 = (1._pr - model%del1)/(1._pr + model%del1)
-            end select
-         end associate
-      end select
-   end function model_from_X
-end module
-
 program main
    !! Binary system parameter optimization
    use yaeos, only: EquilibriaState, pr, ArModel, PengRobinson78, CubicEoS, saturation_pressure
    use forsus, only: Substance, forsus_dir
    use yaeos__fitting, only: FittingProblem, fobj, optimize
-   use yaeos__fitting_fit_nrtl, only: FitMHVNRTL, init_model
+   use yaeos__fitting_fit_nrtl_mhv, only: FitMHVNRTL, init_model
    integer, parameter :: nc = 2, np=7 + nc
    integer :: i, infile, iostat
 
@@ -161,7 +63,8 @@ program main
    print *, "Xf:", X
 
    if (allocated(model)) deallocate (model)
-   model = prob%get_model_from_X(X)
+   call prob%get_model_from_X(X)
+   model = prob%model
 
    ! ===========================================================================
    ! Write out results and experimental values

example/extra/adiff_pr76.f90

@@ -5,10 +5,8 @@ module hyperdual_pr76
     implicit none
 
     type, extends(ArModelAdiff) :: PR76
-        type(Substances) :: composition
         real(pr), allocatable :: kij(:, :), lij(:, :)
         real(pr), allocatable :: ac(:), b(:), k(:)
-        real(pr), allocatable :: tc(:), pc(:), w(:)
     contains
         procedure :: Ar => arfun
         procedure :: get_v0 => v0
@@ -19,24 +17,24 @@ module hyperdual_pr76
 
 contains
 
-    type(PR76) function setup(tc_in, pc_in, w_in, kij_in, lij_in) result(self)
+    type(PR76) function setup(tc, pc, w, kij, lij) result(self)
         !! Seup an Autodiff_PR76 model
-        real(pr) :: tc_in(:)
-        real(pr) :: pc_in(:)
-        real(pr) :: w_in(:)
-        real(pr) :: kij_in(:, :)
-        real(pr) :: lij_in(:, :)
-
-        self%tc = tc_in
-        self%pc = pc_in
-        self%w = w_in
-
-        self%ac = 0.45723553_pr * R**2 * self%tc**2 / self%pc
-        self%b = 0.07779607_pr * R * self%tc/self%pc
-        self%k = 0.37464_pr + 1.54226_pr * self%w - 0.26993_pr * self%w**2
-
-        self%kij = kij_in
-        self%lij = lij_in
+        real(pr) :: tc(:)
+        real(pr) :: pc(:)
+        real(pr) :: w(:)
+        real(pr) :: kij(:, :)
+        real(pr) :: lij(:, :)
+
+        self%components%tc = tc
+        self%components%pc = pc
+        self%components%w = w
+
+        self%ac = 0.45723553_pr * R**2 * tc**2 / pc
+        self%b = 0.07779607_pr * R * tc/pc
+        self%k = 0.37464_pr + 1.54226_pr * w - 0.26993_pr * w**2
+
+        self%kij = kij
+        self%lij = lij
     end function
 
     function arfun(self, n, v, t) result(ar)
@@ -50,8 +48,8 @@ function arfun(self, n, v, t) result(ar)
 
         integer :: i, j
 
-        associate(pc => self%pc, ac => self%ac, b => self%b, k => self%k, &
-                  kij => self%kij, lij => self%lij, tc => self%tc &
+        associate(pc => self%components%pc, ac => self%ac, b => self%b, k => self%k, &
+                  kij => self%kij, lij => self%lij, tc => self%components%tc &
         )
             a = 1.0_pr + k * (1.0_pr - sqrt(t/tc))
             a = ac * a ** 2

example/extra/hard_spheres_mixing.f90

@@ -0,0 +1,209 @@
+module yaeos__HardSpheresCubicEoS
+   use yaeos, only: pr, R, Substances
+   use yaeos__ar_models_hyperdual, only: ArModelAdiff
+   use yaeos__autodiff
+   implicit none
+
+   type, extends(ArModelAdiff) :: HardSpheresCubicEoS
+      real(pr), allocatable :: ac(:), b(:), k(:), kij(:, :), lij(:, :)
+   contains
+      procedure :: Ar => arfun
+      procedure :: get_v0 => v0
+   end type HardSpheresCubicEoS
+
+   real(pr), parameter :: del1 = 1._pr + sqrt(2._pr)
+   real(pr), parameter :: del2 = 1._pr - sqrt(2._pr)
+
+contains
+
+   function arfun(self, n, v, t) result(Ar)
+      class(HardSpheresCubicEoS) :: self !! Model
+      type(hyperdual), intent(in) :: n(:) !! Number of moles vector
+      type(hyperdual), intent(in) :: v !! Volume [L/mol]
+      type(hyperdual), intent(in) :: t !! Temperature [K]
+      type(hyperdual) :: Ar !! Residual Helmholtz energy
+
+      type(hyperdual), dimension(size(n)) :: a
+      type(hyperdual) :: nij
+
+      ! Cubic Parameters
+      type(hyperdual) :: Ar_att
+      type(hyperdual) :: amix, bmix
+      type(hyperdual) :: b_v
+
+      ! HardMixing parameters
+      type(hyperdual) :: Ar_rep
+      type(hyperdual) :: lambda(0:3), eta
+      type(hyperdual) :: l1l2_l3l0
+      type(hyperdual) :: l23_l0l32
+      type(hyperdual) :: logContribution
+      real(pr), parameter    :: xi = 4.0_pr
+
+      integer :: i, j
+      integer :: nc
+
+      nc = size(n)
+
+      ! Alpha function
+      associate(&
+         ac => self%ac, b => self%b, k => self%k, &
+         tc => self%components%tc, &
+         kij => self%kij, lij => self%lij)
+
+         ! Attractive parameter
+         a = self%ac*(1.0_pr + self%k*(1.0_pr - sqrt(t/self%components%tc)))**2
+
+         ! Mixing rule
+         amix = 0.0_pr
+         bmix = 0.0_pr
+         lambda = 0.0_pr
+         do i = 1, nc
+            do j = 1, nc
+               nij = n(i)*n(j)
+               amix = amix + nij * sqrt(a(i)*a(j)) * (1 - kij(i, j))
+               bmix = bmix + nij * 0.5_pr * (b(i) + b(j)) *(1 - lij(i, j))
+            end do
+            lambda(1) = lambda(1) + n(i)*b(i)**(1.0_pr/3.0_pr)
+            lambda(2) = lambda(2) + n(i)*b(i)**(2.0_pr/3.0_pr)
+         end do
+         bmix = bmix/sum(n)
+
+         lambda(0) = sum(n)
+         lambda(3) = bmix
+         eta = bmix/v/xi
+         l1l2_l3l0 = lambda(1)*lambda(2)/lambda(3)/lambda(0)
+         l23_l0l32 = lambda(2)**3/lambda(0)/lambda(3)**2
+         logContribution = log(1._pr - xi*eta)/xi
+
+         Ar_rep = -sum(n)*((1._pr + 3._pr*l1l2_l3l0)*logContribution &
+                          + 3._pr/xi*(l23_l0l32 - 1._pr/2._pr - l1l2_l3l0/2._pr) &
+                          *(eta + logContribution))
+         b_v = bmix/v
+         Ar_att = (- sum(n) * log(1.0_pr - b_v) &
+             - amix / (R*t*bmix)*1.0_pr / (del1 - del2) &
+             * log((1.0_pr + del1 * b_v) / (1.0_pr + del2 * b_v)) &
+             ) * (R * t)
+
+         ar = Ar_rep + Ar_att
+      end associate
+   end function arfun
+
+   function v0(self, n, P, T)
+      class(HardSpheresCubicEoS), intent(in) :: self !! Model
+      real(pr), intent(in) :: n(:) !! Moles vector
+      real(pr), intent(in) :: P !! Pressure [bar]
+      real(pr), intent(in) :: T !! Temperature [K]
+      real(pr) :: v0
+
+      v0 = sum(n * self%b)
+   end function v0
+
+
+   subroutine main
+      use yaeos
+      use forsus, only: Substance, forsus_dir, forsus_default_dir
+      use hyperdual_pr76, only: hPr76 => PR76, set_hpr => setup
+
+      integer, parameter :: nc=2
+
+      real(pr) :: n(nc), V, P, Phs, T
+      real(pr) :: tc(nc), pc(nc), w(nc), kij(nc,nc), lij(nc,nc)
+      type(Substance) :: sus(nc)
+
+      type(CubicEoS) :: pr76
+      type(hPR76) :: hdpr76
+      type(HardSpheresCubicEoS) :: hspr76
+
+      type(EquilibriaState) :: eq
+      type(PTEnvel2) :: env
+
+      real(pr) :: Ar, ArT, ArV, ArTV, ArT2, ArV2, Arn(nc), Artn(nc), ArVn(nc), arn2(nc,nc)
+
+      integer :: i
+
+      forsus_dir = "build/dependencies/forsus/" // forsus_default_dir
+
+      sus(1) = Substance("methane")
+      sus(2) = Substance("n-decane")
+
+      tc = sus%critical%critical_temperature%value
+      pc = sus%critical%critical_pressure%value/1e5
+      w = sus%critical%acentric_factor%value
+
+      tc(2) = 874.0
+      pc(2) = 6.8
+      w(2) = 1.52596
+
+      kij = 0
+      lij = 0
+
+      pr76 = PengRobinson76(tc, pc, w, kij, lij)
+      hdpr76 = set_hpr(tc, pc, w, kij, lij)
+
+      ! Copy PR76 into HSPR76
+      hspr76%components%tc = tc
+      hspr76%components%pc = pc
+      hspr76%components%w = w
+      hspr76%ac = pr76%ac
+      hspr76%b = pr76%b
+
+      associate(alpha => pr76%alpha)
+         select type(alpha)
+          type is (AlphaSoave)
+            hspr76%k = alpha%k
+         end select
+      end associate
+
+      hspr76%kij = kij
+      hspr76%lij = lij
+
+      n = [0.3, 0.7]
+      V = 1
+      T = 400
+
+      block
+         real(pr) :: dPdV, k, khs
+         do i=1,100
+            P = real(i, pr)
+            call volume(pr76, n, P, T, V, root_type="stable")
+            call pressure(pr76, n, V, T, P, dPdV)
+            k = -1/V * 1/dPdV
+
+            call volume(hspr76, n, P, T, V, root_type="stable")
+            call pressure(hspr76, n, V, T, P, dPdV)
+            khs = -1/V * 1/dPdV
+
+            print *, P, k, khs
+         end do
+      end block
+
+      P = 50
+      do i=1,99
+         n(1) = real(i)/100
+         n(2) = 1 - n(1)
+
+         eq = saturation_pressure(pr76, n, T=T, kind="bubble", P0=P)
+         P = eq%p
+         if (eq%iters < 1000) write(1, *) eq%x(1), eq%y(1), eq%p
+
+         eq = saturation_pressure(hspr76, n, T=T, kind="bubble", p0=eq%P)
+         if (eq%iters < 1000) write(2, *) eq%x(1), eq%y(1), eq%p
+      end do
+      call exit
+
+      T = 150
+      eq = saturation_pressure(pr76, n, T=T, kind="bubble")
+      print *, eq%iters, eq
+      env = pt_envelope_2ph(pr76, n, eq)
+
+      print *, size(env%points)
+      write(1, *) env
+
+      eq = saturation_pressure(hspr76, n, T=T, kind="bubble", p0=eq%P)
+      print *, eq%iters, eq
+
+      env = pt_envelope_2ph(hspr76, n, eq)
+      print *, size(env%points)
+      write(2, *) env
+   end subroutine main
+end module yaeos__HardSpheresCubicEoS