Merge pull request #72 from ipqa-research/dev

Road to v0.3

.gitignore

@@ -5,6 +5,7 @@ tools/tapeout/pr_diff.f90
 .gdb_history
 doc/ford_site
 *.mod
+*.smod
 *mypy*
 *.so
 tools/uml
@@ -22,3 +23,24 @@ tools/notebooks/database/molarmass.csv
 tools/notebooks/database/ogUNIFAC/ogUNIFAC_groups.csv
 app/plot
 app/plot.gnu
+test/data/case1.nml
+10K
+D
+env0.dat
+env1.dat
+env2.dat
+env3.dat
+env4.dat
+env5.dat
+env6.dat
+envel
+fit_case
+fort.1
+fort.2
+hv
+kijvar
+log_srk
+log_srk-hv
+srk
+srk_hv
+fort.2

README.md

@@ -72,6 +72,16 @@ call pressure(model, z, V, T, P, dPdN=dPdN)
 print *, dPdN
 ```
 
+Examples of code with simple applications showing the capabilities of `yaeos`
+can be found at [example/tutorials](example/tutorials). Each example can be run
+with:
+
+```bash
+ fpm run --example <example name here>
+```
+
+Not providing any example will show all the possible examples that can be run.
+
 # How to install/run it
 `yaeos` is intended to use as a [`fpm`](fpm.fortran-lang.org)
 
@@ -98,8 +108,32 @@ cd yaeos
 fpm run
 ```
 
+## Developing with vscode
+If your intention is either to develop for `yaeos` or to explore in more detail
+the library with debugging. We provide some predefined defuaults to work with
+`vscode`. You can add them to the cloned repository by running:
+
+```bash
+git clone https://github.com/ipqa-research/vscode-fortran .vscode
+```
+
+From the project main directory 
+
+## Available examples
+In this repository we provide a series of examples of the different things that
+can be calculated with `yaeos`. The source codes for the examples can be seen
+at the [example/tutorials](example/tutorials) directory.
+
+All the examples can be run with
+
+```bash
+fpm run --example <example_name_here>
+```
+
 # Including new models with Automatic Differentiation.
-We are using the `hyperdual` module developed by [Philipp Rehner](https://github.com/prehner) and [Gernot Bauer](https://github.com/g-bauer)
+We are using the `hyperdual` module developed by 
+[Philipp Rehner](https://github.com/prehner) 
+and [Gernot Bauer](https://github.com/g-bauer)
 
 > The automatic differentiation API isn't fully optimized yet so performance is
 > much slower than it should be.

STYLE_GUIDE.md

@@ -17,9 +17,9 @@ This style guide is a living document and proposed changes may be adopted after
 
 * Source files should contain at most one `program`, `module`, or `submodule`
 * The filename should match the program or module name and have the file extension `.f90` or `.F90` if preprocessing is required
-* `module` names should include it's subdirectory. Using `yaeos` for the parent
+* `module` names should include it's subdirectory, using `yaeos__` for the parent
   `src` directory. For example the module in `src/phase_equilibria/flash.f90`
-  should be named `yaeos_phase_equilibria_flash`
+  should be named `yaeos__phase_equilibria_flash`.
 * If the interface and implementation is split using submodules the implementation submodule file should have the same name as the
   interface (parent) module but end in `_implementation`
   E.g., `string_class.f90` and `string_class_implementation.f90`

app/fit.f90

@@ -0,0 +1,173 @@
+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)
+
+      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
+   integer, parameter :: nc = 2, np=7 + nc
+   integer :: i, infile, iostat
+
+   type(EquilibriaState), allocatable :: exp_points(:)
+   type(EquilibriaState) :: point
+
+   type(FitMHVNRTL) :: prob
+   type(Substance) :: sus(2)
+
+   class(ArModel), allocatable :: model
+
+   real(pr) :: T, P, x1, y1, kij, X(np), told, error
+   character(len=14) :: kind
+
+   ! ==========================================================================
+   ! Setup components and read data file
+   ! --------------------------------------------------------------------------
+   forsus_dir = "build/dependencies/forsus/data/json"
+   sus(1) = Substance("nitrogen", only=["critical"])
+   sus(2) = Substance("n-octane", only=["critical"])
+
+   allocate (exp_points(0))
+   open (newunit=infile, file="fit_case", iostat=iostat)
+   do
+      read (infile, *, iostat=iostat) kind, t, p, x1, y1
+      if (iostat /= 0) exit
+      select case (kind)
+      case ("bubble", "dew", "liquid-liquid")
+         point = EquilibriaState( &
+                 kind=kind, T=T, P=P, x=[x1, 1 - x1], y=[y1, 1 - y1], &
+                 Vx=0._pr, Vy=0._pr, iters=0, beta=0._pr &
+                 )
+      end select
+      exp_points = [exp_points, point]
+   end do
+   close (infile)
+
+   ! ==========================================================================
+   ! Setup optimization problem and call the optimization function
+   ! --------------------------------------------------------------------------
+   call init_model(prob, sus)
+
+   prob%experimental_points = exp_points
+
+   X = 0
+   X(1:2) = [0.1, 0.3]
+   X(5:6) = [0.1, 0.2]
+   X(8:) = [1, 1]
+
+   prob%parameter_step = [(0.5_pr, i=1,size(x))]
+   prob%solver_tolerance = 1e-7
+   prob%verbose = .true.
+
+   print *, "X0:", X
+   error = optimize(X, prob)
+   print *, "FO:", error
+   print *, "Xf:", X
+
+   if (allocated(model)) deallocate (model)
+   model = prob%get_model_from_X(X)
+
+   ! ===========================================================================
+   ! Write out results and experimental values
+   ! ---------------------------------------------------------------------------
+   told = exp_points(1)%T
+   do i = 1, size(exp_points)
+      point = saturation_pressure( &
+              model, exp_points(i)%x, exp_points(i)%t, kind="bubble", &
+              p0=exp_points(i)%p, y0=exp_points(i)%y &
+              )
+      if (told /= point%t) write (*, "(/)")
+      print *, exp_points(i)%x(1), exp_points(i)%y(1), exp_points(i)%P, &
+         point%x(1), point%y(1), point%P
+      told = point%T
+   end do
+end program main

app/isotherm.f90

@@ -1,25 +1,32 @@
 program main
-    use yaeos, only: pr, ArModel, PengRobinson76, volume
+    use yaeos, only: pr, ArModel, PengRobinson76, volume, pressure, SoaveRedlichKwong
+    use forsus, only: Substance, forsus_dir
     implicit none
     class(ArModel), allocatable :: eos
 
     integer, parameter :: nc = 2
     real(pr), dimension(nc) :: n, tc, pc, w
     real(pr), dimension(nc, nc) :: kij, lij
 
+    type(Substance) :: sus(2)
+
     real(pr) :: v
     real(pr) :: p0, pf, dp, p
     real(pr) :: t0, tf, dt, t
     integer :: i, j, n_p_points=500, n_t_points=5
 
-    n = [0.3_pr, 0.7_pr]
-    tc = [190._pr, 310._pr]
-    pc = [14._pr, 30._pr]
-    w = [0.001_pr, 0.03_pr]
-    kij = reshape([0., 0.1, 0.1, 0.], [nc,nc])
+    forsus_dir = "build/dependencies/forsus/data/json"
+    sus(1) = Substance("propane", only=["critical"])
+    sus(2) = Substance("n-butane", only=["critical"])
+
+    n = [0.7, 0.3]
+    tc = sus%critical%critical_temperature%value
+    pc = sus%critical%critical_pressure%value/1e5
+    w = sus%critical%acentric_factor%value
+    kij = reshape([0., 0.0, 0.0, 0.], [nc,nc])
     lij = kij / 2
 
-    eos = PengRobinson76(tc, pc, w, kij, lij)
+    eos = SoaveRedlichKwong(tc, pc, w, kij, lij)
 
     p0 = 100
     pf = 0.1
@@ -29,15 +36,24 @@ program main
     tf = 350
     dt = (tf-t0)/n_t_points
 
-    do j=0, n_t_points - 1
-        t = t0 + dt * j
-        print *, "# ", t
-        do i=0,n_p_points-1
-            p = p0 + dp * i
-            call volume(eos, n, p, t, v, root_type="stable")
-            print *, v, p
-        end do
-        print *, ""
-        print *, ""
+    T = 300
+    do i=1,1000
+        V = real(i, pr)/1000._pr
+        call pressure(eos, n, V, T, P=P)
+        print *, V, P
     end do
+
+
+
+    ! do j=0, n_t_points - 1
+    !     t = t0 + dt * j
+    !     print *, "# ", t
+    !     do i=0,n_p_points-1
+    !         p = p0 + dp * i
+    !         call volume(eos, n, p, t, v, root_type="stable")
+    !         print *, v, p
+    !     end do
+    !     print *, ""
+    !     print *, ""
+    ! end do
 end program main