Merge pull request #19 from ipqa-research/main

Update from main

.gitmodules

@@ -0,0 +1,3 @@
+[submodule ".vscode"]
+	path = .vscode
+	url = git@github.com:ipqa-research/vscode-fortran.git

README.md

@@ -90,6 +90,6 @@ call fugacity_vt(eos, n, v, t, lnfug=lnfug, dlnphidn=dlnphidn)
 
 ## Documentation
 The latest API documentation for the `main` branch can be found
-[here](https://fedebenelli.github.io/yaeos). This was generated from the source
+[here](https://ipqa-research.github.io/yaeos). This was generated from the source
 code using [FORD](https://github.com/Fortran-FOSS-Programmers/ford). We're
 working in extending it more.

app/flash.f90

@@ -0,0 +1,34 @@
+program flasher
+    use yaeos, only: pr, EquilibriaState, flash, PengRobinson76, ArModel, fugacity_tp
+    implicit none
+
+    class(ArModel), allocatable :: model
+    type(EquilibriaState) :: flash_result
+
+    real(pr) :: tc(2), pc(2), w(2)
+
+    real(pr) :: n(2), p, v, t, k0(2)
+    integer :: iters, i, j
+
+    print *, "FLASH EXAMPLE:"
+
+    n = [0.4, 0.6]
+    tc = [190.564, 425.12]
+    pc = [45.99, 37.96]
+    w = [0.0115478, 0.200164]
+    model = PengRobinson76(tc, pc, w)
+
+    P = 60
+    t = 294
+    k0 = (PC/P)*exp(5.373*(1 + w)*(1 - TC/T))
+    print *, k0
+
+    flash_result = flash(model, n, t=t, p_spec=p, k0=k0, iters=iters)
+
+    print *, "X:", flash_result%x, sum(flash_result%x)
+    print *, "Y:", flash_result%y, sum(flash_result%y)
+    print *, "Vx: ", flash_result%Vx
+    print *, "Vy: ", flash_result%Vy
+    print *, "P: ", flash_result%p
+    print *, "T: ", flash_result%T
+end program

example/benchmark.f90

@@ -46,7 +46,7 @@ subroutine yaeos_run(n, dn, model_name)
         if (dn) then
             call fugacity_vt(model, z, V, T, P, lnfug, dlnPhidP, dlnphidT, dlnphidn)
         else
-            call fugacity_vt(model, z, V, T, P, lnfug, dlnPhidP, dlnphidT)
+            call fugacity_vt(model, z, V, T, lnfug=lnfug)
         end if
     end subroutine
 
@@ -57,13 +57,13 @@ subroutine benchmarks
         real(8) :: et, st
 
 
-        do n=1,5
+        do n=1,30
             time = 0
             std = 0
             mean = 0
             do i=1,nevals
                 call cpu_time(st)
-                    call yaeos_run(n, .true., "Adiff PR76")
+                    call yaeos_run(n, .true., "Analytic PR76")
                 call cpu_time(et)
 
                 time = (et-st)*1e6

example/demo.f90

@@ -1,7 +1,9 @@
 program examples
     use bench, only: benchmarks
     use hyperdual_pr76, only: adiff_pr76
+    use flashing, only: run_flashes
 
-    call benchmarks
-    call adiff_pr76
+    ! call benchmarks
+    ! call adiff_pr76
+    call run_flashes
 end program

example/flash.f90

@@ -0,0 +1,33 @@
+module flashing
+    use yaeos, only: pr, EquilibriaState, flash, PengRobinson76, ArModel, fugacity_tp
+    implicit none
+contains
+
+    subroutine run_flashes
+        class(ArModel), allocatable :: model
+        type(EquilibriaState) :: flash_result
+
+        real(pr) :: tc(2), pc(2), w(2)
+
+        real(pr) :: n(2), p, v, t, k0(2)
+        integer :: iters, i, j
+
+        print *, "FLASH EXAMPLE:"
+
+        n = [0.4, 0.6]
+        tc = [190.564, 425.12]
+        pc = [45.99, 37.96]
+        w = [0.0115478, 0.200164]
+        model = PengRobinson76(tc, pc, w)
+
+        P = 60
+        t = 294
+        k0 = (PC/P)*exp(5.373*(1 + w)*(1 - TC/T))
+        print *, k0
+
+        flash_result = flash(model, n, t=t, p_spec=p, k0=k0, iters=iters)
+        print *, "X:", flash_result%x, sum(flash_result%x)
+        print *, "Y:", flash_result%y, sum(flash_result%y)
+
+    end subroutine
+end module

example/models/models.f90

@@ -0,0 +1,21 @@
+module models
+    use yaeos, only: pr, R, CubicEoS
+
+contains
+
+    type(CubicEoS) function binary_PR76() result(eos)
+        use yaeos, only: PengRobinson76
+        integer, parameter :: n=2
+        real(pr) :: tc(n), pc(n), w(n)
+        real(pr) :: kij(n, n), lij(n, n)
+        tc = [190._pr, 310._pr]
+        pc = [14._pr, 30._pr]
+        w = [0.001_pr, 0.03_pr]
+
+        kij = reshape([0., 0.0, 0.0, 0.], [n,n]) 
+        lij = kij / 2 
+
+        eos = PengRobinson76(tc, pc, w, kij, lij)
+    end function
+
+end module

src/constants.f90

@@ -4,7 +4,7 @@ module yaeos_constants
 
    implicit none
 
-   integer, parameter :: pr = real64 !! Machine Precision
+   integer, parameter :: pr = real64 !! Used precision
    real(pr), parameter :: R = 0.08314472_pr !! Ideal Gas constant
    character(len=254) :: database_path = "database" !! Path to find database
    character(len=1) :: path_sep = "/" !! File separator (to preprocess on Win or Mac/linux)

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

@@ -58,13 +58,13 @@ type(CubicEoS) function PengRobinson76(tc, pc, w, kij, lij) result(model)
         if (present(kij)) then
             mixrule%k = kij
         else
-            mixrule%k = reshape([(i, i=1,nc**2)], [nc, nc])
+            mixrule%k = reshape([(0, i=1,nc**2)], [nc, nc])
         endif
 
         if (present(lij)) then
             mixrule%l = lij
         else
-            mixrule%l = reshape([(i, i=1,nc**2)], [nc, nc])
+            mixrule%l = reshape([(0, i=1,nc**2)], [nc, nc])
         endif
 
         model%components = composition

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

@@ -52,13 +52,13 @@ type(CubicEoS) function SoaveRedlichKwong(tc, pc, w, kij, lij) result(model)
         if (present(kij)) then
             mixrule%k = kij
         else
-            mixrule%k = reshape([(i, i=1,nc**2)], [nc, nc])
+            mixrule%k = reshape([(0, i=1,nc**2)], [nc, nc])
         endif
 
         if (present(lij)) then
             mixrule%l = lij
         else
-            mixrule%l = reshape([(i, i=1,nc**2)], [nc, nc])
+            mixrule%l = reshape([(0, i=1,nc**2)], [nc, nc])
         endif
 
         model%components = composition

src/phase_equilibria/equilibria_state.f90

@@ -0,0 +1,14 @@
+module yaeos_equilibria_equilibria_state
+   use yaeos_constants, only: pr
+   implicit none
+
+   type :: EquilibriaState
+      integer :: iters !! Iterations needed to reach the state
+      real(pr), allocatable :: y(:) !! Vapour molar fractions
+      real(pr), allocatable :: x(:) !! Liquid molar fractions
+      real(pr) :: Vx !! Liquid volume [L/mol]
+      real(pr) :: Vy !! Vapor volume [L/mol]
+      real(pr) :: t !! Temperature [K]
+      real(pr) :: p !! Pressure [bar]
+   end type
+end module