v0.2.0 (#25)

* [Donwhill] migrate to Downhill.jl

* [vtflash] add state to vtflash result

* [test vtflash] just work test

* [vtflash] add directory for vtflash

* [vtflash] nvt gradient and hessian with test

* [vtflash] nvthessian! and tests

Create a function for calculating hessian only.
Add tests for it on 4-component mixture.

Update documentation of nvt.jl functions.

* [vtflash] moved and renamed pressure gradient function to nvt.jl

* [vtflash] add helmholtz energy calculation

* [vtflash] Draft of general interface

For now, PhysicalState variables are implemented and tested against previous code.

Both implementations (general for variables and original) work for now and may be compared.

Soon, original version will be deprecated and granted name of FractionState or smth like that.

* [vtflash] upd vtflashresult interface

* Layout change of mixture to StructVector

* [vtflash] rm some commented code

* [vtflash] mv documentation for variables constructor from concentration and saturation

* [vtflash ratiostate] constructors for RatioState

* [flash physical] copying of x to state wrapper in constrain step

* [vtflash] kwargs for vt_flash to control optimizator

* [vtflash] add errors for constrain step

* [vtflash physical] throwing custom errors

* [vtflash] update for optimizers tolerances and single phase result

* [vtflash physical] upd docstring for script

* [vtflash ratio] add gradient, hessian and closures. It works!

* [vtflash idealidentity] constructors

* [vtflash idealidentity] add gradient, hessian and closures

* [test vtflash] add comparing of vtflash for 3 possible internal variables

* [vtflash newton] mv

* [vtflash ratio] typo

* [vtflash newton] Add support of new interface to vt_flash_newton

Flash by Newton was tested against bfgs for RatioState for four-component mixture.

I faced possible drawbacks of copying in closures like
state1x .= x
for Newton algo the internal x = convert(..., x_), which gives x === x_ for vector of floats.
So, in @debug call there was overwriting in gradient (gprev) call.
Somehow it produced mistakes.

* [vtflash idealidentity] typo-bug

* [vtflash nvt] upd comment

* [vtflash newton] force copy to float of init x

* [test vt_flash] add vt_flash_newton tests

* [test vtflash] tests of states for bfgs and newton based flashes

* [vtflash] parametric constructor for state variables

* [vtflash] occasionally, this seems type-stable

* [vtflash] fix abstractype in field of BrusilovskyEoSMixture

* #22 Rm unnecessary type alias definition

* Convergence condition on chemical potentials and pressure (#24)

* #23 Add convcond to newton optimization
* #23 Convergence condition for newton version of vt-flash
* #23 Update of convergence criteria
* #23 Update conv criteria, now tolerances are same by default, chematol is renamed to chemtol
* #23 Change signature of convcond to (x, xpre, y, ypre, g) -> Bool
* #23 Deprecate unused simplified version of linesearch
* #23 Add convcond for quasinewton vt-flash / Can be used currently only from dev of Downhill
* #23 Move convergence closure to vt_flash.jl source
* #23 Documentation for vt_flash, vt_flash! and vt_flash_newton

* Set version to v0.2.0

* Upd Readme

Co-authored-by: Vasily <[email protected]>

Project.toml

@@ -1,12 +1,13 @@
 name = "CubicEoS"
 uuid = "0772a1fa-d924-41f2-9f07-889679823dc5"
 authors = ["Vasily Pisarev <[email protected]>", "Stepan Zakharov <[email protected]>"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 CubicEoSDatabase = "8742253c-1987-415f-b31d-1ec9670d9ed7"
-DescentMethods = "f4becde8-b16e-4b5a-8f91-16ef0c22c8bc"
+Downhill = "a4c28711-7027-4a57-8564-74545b4697a4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 
 [compat]
 julia = "1.4"

README.md

@@ -1,54 +1,101 @@
 # CubicEoS.jl
-The package implements functions to work with substances and mixtures described by cubic equations of state.
+
+[![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)
+
+CubicEoS.jl implements functions to work with substances and mixtures described by cubic equations of state. This includes
+
+- Basic thermodynamics: pressure, z-factor, chemical potential etc.;
+- Check of NVT stability of single-phase state;
+- Phase equilibrium calculation for NVT variables.
 
 So far, the general cubic equation of state [[Brusilovsky, SPE Reservoir Engineering, February 1992](https://doi.org/10.2118/20180-PA)] is implemented.
 
+Most of implemented functions
+
+- are designed in a zero-allocating way, such functions has optional `buf` keyword;
+- has a desctructive option `func!`.
+
+## Installation
+
+The package registry is **[LammpsToolsRegistry](https://github.com/stepanzh-test-org/LammpsToolsRegistry)**. So, you need to add the registry and then install CubicEoS.jl in a usual Julia-way.
+
+```julia
+]pkg> registry add https://github.com/stepanzh-test-org/LammpsToolsRegistry
+pkg> add CubicEoS
+```
+
 ## Brusilovsky equation of state
 
-To construct an set of parameters for a component, use
+Components and mixtures may be constructed either explicitly or by loading.
+
+A **component** can be defined explicitly
 ```julia
-BrusilovskyEoSComponent(;name="No Name", critical_pressure, critical_temperature, acentric_factor, Omegac, Zc, Psi, molar_mass, carbon_number::Integer)
+BrusilovskyEoSComponent(; name="No Name", critical_pressure, critical_temperature, acentric_factor, Omegac, Zc, Psi, molar_mass, carbon_number::Integer)
 ```
-The values of parameters `Omegac`, `Zc` and `Psi` for hydrocarbons may be found in Brusilovsky's paper (https://doi.org/10.2118/20180-PA).
+where the values of parameters `Omegac`, `Zc` and `Psi` may be found in Brusilovsky's paper (https://doi.org/10.2118/20180-PA).
 
 The temperatures must be in absolute scale (e.g., in Kelvins).
 
-The parameters can be loaded from a file (using **CubicEoSDatabase.jl**):
+Or, the component can be loaded from a file (using **[CubicEoSDatabase.jl](https://github.com/stepanzh/CubicEoSDatabase.jl)**)
 ```julia
-methane = load(BrusilovskyEoSComponent, name = "methane", physics_db = "martinez.csv", eos_db = "brusilovsky.csv")
+methane = CubicEoS.load(BrusilovskyEoSComponent; name="methane"[, custom_databases...])
 ```
 
-Mixtures are constructed via
+A **mixture** is constructed via
+
 ```julia
 BrusilovskyEoSMixture(; components::AbstractVector{<:BrusilovskyEoSComponent}, constant, linear, quadratic)
 ```
 where `constant`, `linear` and `quadratic` are matrices of constant, linear and quadratic in temperature terms for Zudkevitch-Joffe corrections.
 
-The parameters can be loaded from a file (using **CubicEoSDatabase.jl**):
+Or can be loaded from a file (using **[CubicEoSDatabase.jl](https://github.com/stepanzh/CubicEoSDatabase.jl)**):
 ```julia
-c1c5 = load(BrusilovskyEoSMixture, names = ["methane", "n-pentane", comp_physics_db = "martinez.csv", comp_eos_db = "brusilovsky.csv", mix_eos_db = "brusilovsky_mix.csv")
+c1c5 = CubicEoS.load(BrusilovskyEoSMixture; names=("methane", "n-pentane")[, custom_databases...])
 ```
 
 ## Basic thermodynamics
 
-To get the pressure of a pure component:
+Basic thermodynamics includes pressure, Wilson saturation pressure and z-factor (`compressibility`).
+
 ```julia
-pressure(component; nmol, volume, temperature)
+pressure(component or mixture, nmol, volume, temperature)
 ```
 
-To get the estimate of the saturation pressure at a given temperature:
 ```julia
 wilson_saturation_pressure(component, RT)
 ```
 
-To get compressibility at given pressure and temperature:
 ```julia
-compressibility(mixture, molar_composition, pressure, RT, phase = 'g'
+compressibility(mixture, molar_composition, pressure, RT, phase='g')
 ```
 
-## Phase stability
+## Chemical potential
+
+The packages includes functions for calculating activity coefficient and its Jacobian matrix for a mixture defined by Brusilovsky EoS.
+
+```julia
+log_ca = log_c_activity(mixture, nmol, volume, RT)
+log_ca, jac = log_c_activity_wj(mixture, nmol, volume, RT)
+```
+
+In case of a component you may use a mixture of one component.
+
+## NVT phase equilibrium
+
+**Phase stability.** To check if a single-phase state is stable, defined in NVT variables, use
 
-To check if a single-phase state is stable:
 ```julia
-vt_stability(mix::BrusilovskyEoSMixture, nmol::AbstractVector, volume, RT)
+isstable, vtstab_results = vt_stability(mix, nmol, volume, RT)
 ```
+
+**Flash.** To calculate NVT phase equilibrium use
+
+```julia
+# quasi-Newton phase split
+flash_result = vt_flash(mix, nmol, volume, RT, StateVariables[; tol, maxiter])
+# Newton phase split
+flash_result = vt_flash_newton(mix, nmol, volume, RT, StateVariables[; tol, maxiter])
+```
+
+where type `StateVariables` defines an internal variables used at phase split stage in optimization solver (e.g. `CubicEoS.PhysicalState`). Quasi-Newton solver is implemented in **[Downhill.jl](https://github.com/vvpisarev/Downhill.jl)**.
+

src/CubicEoS.jl

@@ -6,7 +6,10 @@ export molar_mass, carbon_number, name, describe, load, components # Do we need
 export ncomponents, pressure, wilson_saturation_pressure, compressibility
 export log_c_activity!, log_c_activity, log_c_activity_wj!, log_c_activity_wj
 export vt_stability, vt_stability_buffer
-export vt_flash
+export vt_flash, vt_flash_newton
+export converged
+
+using Downhill
 
 include("constants.jl")
 include("types.jl")
@@ -15,7 +18,6 @@ include("dbload.jl")
 include("basic_thermo.jl")
 include("chempotential.jl")
 include("vt_stability.jl")
-include("vt_flash.jl")
-include("newton.jl")
+include("vt_flash/vt_flash.jl")
 
 end # module

src/basic_thermo.jl

@@ -157,18 +157,18 @@ function __eos_parameters_impl__(
     ai, aij = auxv, auxm
 
     comp = mixture.components
-    ai .= a_coef.(comp, RT)
+    psi = comp.Psi
+    ac = comp.ac
+    @. ai = ac * (1 + psi * (1 - sqrt(RT / comp.RTc)) )^2
 
-    Bm = Cm = Dm = zero(T)
-    @inbounds for i in eachindex(nmol, comp)
-        Bm += nmol[i] * comp[i].b
-        Cm += nmol[i] * comp[i].c
-        Dm += nmol[i] * comp[i].d
-    end
+    bi, ci, di = comp.b, comp.c, comp.d
+    Bm = dot(nmol, bi)
+    Cm = dot(nmol, ci)
+    Dm = dot(nmol, di)
 
-    temp = RT / GAS_CONSTANT_SI - 273.16
+    tempC = RT / GAS_CONSTANT_SI - 273.16
     eij, gij, hij = mixture.eij, mixture.gij, mixture.hij
-    aij .= (one(T) .- (eij .+ temp .* (gij .+ temp .* hij))) .* sqrt.(ai .* ai')
+    aij .= (one(T) .- (eij .+ tempC .* (gij .+ tempC .* hij))) .* sqrt.(ai .* ai')
     Am = dot(nmol, aij, nmol)  # Am = nmolᵀ aij nmol
     return Am, Bm, Cm, Dm, aij
 end

src/chempotential.jl

@@ -87,6 +87,7 @@ function log_c_activity!(
     nc = ncomponents(mix)
     ∑mol = sum(nmol)
 
+    comp = components(mix)
     # коэффициенты вычислены для состава в [моль]
     A, B, C, D, aij = eos_parameters(mix, nmol, RT; buf = buf)
 
@@ -100,12 +101,7 @@ function log_c_activity!(
     AbyRTCmD = A / (RT * CmD)
     log_VpCbyVpD = log(VpC / VpD)
 
-    @inbounds for i in 1:nc
-        subst = mix.components[i]
-        b = subst.b
-        c = subst.c
-        d = subst.d
-
+    @inbounds for (i, b, c, d) in zip(1:nc, comp.b, comp.c, comp.d)
         # ∂A_i = 2 ∑_j (N_j * a_ij)
         # прежний код использовал sum( генератор с for in eachindex ) - давал 2 аллокации
         ∂A = 2 * dot(nmol, @view aij[:,i])
@@ -154,8 +150,8 @@ function log_c_activity_wj(
     buf::BrusilovskyThermoBuffer = thermo_buffer(mix),
 ) where {T}
     nc = ncomponents(mix)
-    log_ca = Vector{T}(undef, nc)
-    jacobian = similar(log_ca, (nc, nc))
+    log_ca = similar(nmol, T)
+    jacobian = Matrix{T}(undef, nc, nc)
     return log_c_activity_wj!(log_ca, jacobian, mix, nmol, volume, RT; buf = buf)
 end
 
@@ -216,6 +212,7 @@ function __log_c_activity_wj_impl__(
 
     aux1 = buf.vec1
     aux2 = buf.vec2
+    comp = components(mix)
     nc = ncomponents(mix)
     V = volume
     ntotal = sum(nmol)
@@ -228,33 +225,24 @@ function __log_c_activity_wj_impl__(
     log_VmBbyV = log(VmB/volume)
     AbyRTCmD = A / (RT * CmD)
     log_VpCbyVpD = log(VpC / VpD)
-    log2 = log1p(CmD / (V + D)) / (RT * CmD)
+    log2 = 2 * log1p(CmD / (V + D)) / (RT * CmD)
+
+    aux1 .= comp.b ./ VmB
 
-    @inbounds map!(aux1, 1:nc) do i
-        mix[i].b / VmB
-    end
     log_ca .= aux1 .* ntotal .- log_VmBbyV
     jacobian .= aux1 .+ aux1' .+ ntotal .* aux1 .* aux1'
     mul!(aux1, aij, nmol)
-    @inbounds map!(aux2, 1:nc) do i
-        mix[i].c - mix[i].d
-    end
+    aux2 .= comp.c .- comp.d
     log_ca .-= AbyRTCmD .* log_VpCbyVpD .* (2 .* aux1 ./ A .- aux2 ./ CmD)
-    jacobian .-= (2 * log2) .* ((A / CmD^2) .* aux2 .* aux2' .+ aij)
-    jacobian .+= (2 * log2 / CmD) .* (aux2 .* aux1' .+ aux1 .* aux2')
+    jacobian .-= log2 .* ((A / CmD^2) .* aux2 .* aux2' .+ aij)
+    jacobian .+= (log2 / CmD) .* (aux2 .* aux1' .+ aux1 .* aux2')
 
     aux1 .= 2 .* aux1 ./ (RT * CmD) .- aux2 .* (A / (RT * CmD^2))
-    @inbounds map!(aux2, 1:nc) do i
-        mix[i].c / VpC - mix[i].d / VpD
-    end
+    aux2 .= comp.c ./ VpC .- comp.d ./ VpD
     log_ca .-= AbyRTCmD .* aux2
     jacobian .-= aux1 .* aux2'.+ aux2 .* aux1'
-    @inbounds map!(aux1, 1:nc) do i
-        mix[i].d / VpD
-    end
-    @inbounds map!(aux2, 1:nc) do i
-        mix[i].c / VpC
-    end
+    aux1 .= comp.d ./ VpD
+    aux2 .= comp.c ./ VpC
 
     jacobian .-= AbyRTCmD .* (aux1 .* aux1' .- aux2 .* aux2')
     return log_ca, jacobian