Merge pull request #163 from thermotools/cpp

C++ wrapper and cmake based build system

.gitignore

@@ -2,6 +2,7 @@ syntax: glob
 
 objects
 bin
+build/
 run*
 doc/dox.out
 doc/doxygen/html

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "external/pFUnit"]
+	path = external/pFUnit
+	url = https://github.com/thermotools/pFUnit.git
+[submodule "external/lapack"]
+	path = external/lapack
+	url = https://github.com/Reference-LAPACK/lapack.git

CMakeLists.txt

@@ -0,0 +1,74 @@
+string(ASCII 27 Esc)
+set(ColorRed "${Esc}[31m")
+set(ColorBlue "${Esc}[34m")
+
+cmake_minimum_required(VERSION 3.5)
+project(ThermoPack LANGUAGES Fortran)
+
+if(NOT MSVC)
+    execute_process(COMMAND bash -c "arch" OUTPUT_VARIABLE PROC)
+    set(tp_flags_common "-cpp -fPIC -fdefault-real-8 -fdefault-double-8 -frecursive -fopenmp -Wno-unused-function -Wno-unused-variable")
+
+    include(CheckCompilerFlag)
+    check_compiler_flag(Fortran "-arch arm64" arm64_supported)
+    if(arm64_supported)
+        set(gf_proc "-arch arm64")
+        set(gf_march "-arch arm64 -fno-expensive-optimizations")
+    else()
+        set(gf_proc "-mieee-fp")
+        set(gf_march "-march=x86-64 -msse2")
+    endif()
+
+    set(tp_debug_flags "${gf_proc} ${tp_flags_common} -g -fbounds-check -fbacktrace -ffpe-trap=invalid,zero,overflow -Wno-unused-dummy-argument -Wall")
+    set(tp_profile_flags "${tp_flags_common} -g -pg")
+    set(tp_optim_flags "${tp_flags_common} -O3 ${gf_march} -funroll-loops")
+
+    if(APPLE)
+        set(CMAKE_FORTRAN_COMPILER /opt/homebrew/bin/gfortran)
+        set(CMAKE_OSX_ARCHITECTURES "arm64" CACHE STRING "The OSX architecture")
+    endif()
+else()
+    # Build using: cmake --build . --config Release
+    set(tp_flags_common "/nologo /fpp /fpe:0 /names:lowercase /assume:underscore /fp:precise /extend-source:132 /iface:cref")
+    set(tp_optim_flags "${tp_flags_common}")
+    set(tp_debug_flags "${tp_flags_common} /check:bounds /traceback")
+endif()
+
+# Check if the build type is not set and set it to Release
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE)
+endif()
+
+# Set FORTRAN compile flags for thermopack and lapack
+set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} ${tp_optim_flags}")
+set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} ${tp_debug_flags}")
+
+if(MSVC)
+    if(POLICY CMP0077)
+        set(CMAKE_POLICY_DEFAULT_CMP0077 NEW)
+    endif()
+    # Disable single/complex/complex16 as it is not used by thermopack
+    set(BUILD_SINGLE OFF)
+    set(BUILD_COMPLEX OFF)
+    set(BUILD_COMPLEX16 OFF)
+    add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/lapack)
+    # Add flags not comatible with lapack
+    set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} /real-size:64 ")
+    set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} /real-size:64 /warn:all,noexternal /check:all,noarg_temp_created,nopointers")
+endif(MSVC)
+
+add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/src "thermopack")
+
+option(test "Build test suite" OFF)
+if(test)
+    cmake_policy(SET CMP0074 NEW)
+    find_package(PFUNIT QUIET)
+    if (NOT PFUNIT_FOUND)
+        message("${ColorRed}Did not find pFUnit. If you have installed pFUnit, make sure that PFUNIT_DIR is set. Falling back to building pFUnit from source...")
+        add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/pFUnit)
+    else()
+        message(STATUS "Found pFUnit: ${PFUNIT_DIR}")
+    endif()
+    message(STATUS "Creating run_unittests target")
+    add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/unittests)
+endif()

addon/cppExamples/.gitignore

@@ -0,0 +1 @@
+build/

addon/cppExamples/CMakeLists.txt

@@ -0,0 +1,53 @@
+# Example for Compiling and linking C++ executables using thermopack.
+# This example assumes that libthermopack.dylib has already been built and installed with `make install`
+# The environment variable THERMOPACK_DIR should point to the root directory of ThermoPack (where thermopack-config.cmake is found)
+# this can be set with `export THERMOPACK_DIR=<path/to/thermopack>`.
+#
+# Note: On macOS (arm64), you might need to
+#   export CC=/opt/homebrew/bin/gcc-13
+#   export CXX=/opt/homebrew/bin/g++-13
+# before running cmake, in order to force compilation for arm64.
+
+project(cppExamples)
+cmake_minimum_required(VERSION 3.2)
+
+string(ASCII 27 Esc)
+set(ColourRed "${Esc}[31m")
+set(ColourBlue "${Esc}[34m")
+
+set(THERMOPACK_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../..")
+
+# Searches for thermopack-config.cmake in the directory ${THERMOPACK_DIR}, which can be set with `export THERMOPACK_DIR=<path/to/thermopack>`
+# Once thermopack is found, it provides some convenience variables shown below, as well as the exported target `thermopack` which can be used for linking directly.
+find_package(THERMOPACK REQUIRED)
+if (NOT THERMOPACK_INSTALLED)
+    message("${ColourRed}Thermopack was found but is not installed. Exiting ...")
+    return()
+endif()
+
+message("${ColourBlue}Thermopack was found at : ${THERMOPACK_ROOT}")
+message("${ColourBlue}Using libthermopack at : ${THERMOPACK_LIB}")
+message("${ColourBlue}Using cppThermopack headers at : ${THERMOPACK_INCLUDE}")
+message("${ColourBlue}Linking thermopack using exported target thermopack")
+
+set(CMAKE_BUILD_TYPE "Release")
+set(CMAKE_CXX_FLAGS "-Wall -Wextra -Wfatal-errors -std=c++17")
+
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -O0")
+
+
+add_executable(basic basic.cpp)
+target_link_libraries(basic PUBLIC thermopack)
+
+add_executable(tp_properties tp_properties.cpp)
+target_link_libraries(tp_properties PUBLIC thermopack)
+
+add_executable(tv_properties tv_properties.cpp)
+target_link_libraries(tv_properties PUBLIC thermopack)
+
+add_executable(flashes flashes.cpp)
+target_link_libraries(flashes PUBLIC thermopack)
+
+add_executable(saturation saturation.cpp)
+target_link_libraries(saturation PUBLIC thermopack)

addon/cppExamples/basic.cpp

@@ -0,0 +1,44 @@
+#include <iostream>
+#include <vector>
+#include <cppThermopack/cubic.h>
+#include <cppThermopack/saftvrmie.h>
+
+// Note: The structs `Property` and `VectorProperty` are found in utils.h, which is implicitly included in all headers
+// through thermo.h
+// Also note: Order of bools for computing differentials follows the same convention as pycThermoPack (i.e. order of bools is equal to order of input-parameters)
+// In general: The method signatures are kept consistent with pycThermopack, such that the pycThermopack docs can serve as a coarse lookup for the methods in the cpp wrapper
+//             methods that return a single value (i.e. pressure_tv, specific_volume, etc.) will return a `Property`, which has an implicit conversion to `double`, while
+//             methods that return several values (i.e. chemical_potential_tv) will return a `VectorProperty`, which has an implicit conversion to std::vector<double>
+
+int main(){
+    PengRobinson pr("N2,O2"); // Initialise Peng-Robinson EoS
+    double T = 300.;
+    double V = 1e3;
+    vector1d n = {1., 2.};
+
+    double pres = pr.pressure_tv(T, V, n); // Compute pressure (implicit conversion to double)
+    Property p = pr.pressure_tv(T, V, n, false, true, true); // Computing property and differentials
+    std::cout << "P : \n" << p << std::endl;
+    std::cout << "Pressures are equal : " << ((pres == p.value()) ? "true" : "false") << std::endl;
+    double dpdv = p.dv(); // Accessing differentials
+    std::cout << "Accessed dpdv : " << dpdv << std::endl;
+
+    std::vector<double> chem_pot = pr.chemical_potential_tv(T, V, n); // Compute chemical potential (implicit conversion to std::vector)
+    VectorProperty mu = pr.chemical_potential_tv(T, V, n, true, false, true); // Computing vector-valued property and differentials
+    std::cout << "mu :\n" << mu;
+    std::vector<double> dmudt = mu.dt(); // Accessing the differentials
+    std::vector<std::vector<double>> dmudn = mu.dn();
+
+    Saftvrmie svrm("AR,KR"); // Initialise SAFT-VR Mie EoS
+    V = 1e-3;
+    std::cout << "SVRM chemical potential : \n" << svrm.chemical_potential_tv(T, V, n, false, true) << std::endl;
+
+
+    Saftvrmie svrm_c1("C1");
+    p = 1e6;
+    T = 120.;
+    double vl = svrm_c1.specific_volume(T, p, {1}, Phase::liq);
+    double vg = svrm_c1.specific_volume(T, p, {1}, Phase::vap);
+    std::cout << "At T = " << T << " K, p = " << p / 1e5 << " bar, specific volume of methane is : " << vl << " / " << vg << std::endl;
+    return 0;
+}

addon/cppExamples/flashes.cpp

@@ -0,0 +1,40 @@
+#include <cppThermopack/cubic.h>
+#include <vector>
+#include <iostream>
+
+int main(){
+    SoaveRedlichKwong srk("PROP1OL,NC6");
+    std::cout << "Demonstrating flashes with SRK ...\n";
+
+    double T1{348.}, p1{1e5};
+    std::vector<double> z = {0.4, 0.6};
+    FlashResult flsh = srk.two_phase_tpflash(T1, p1, z);
+    std::cout << "Initial state is \n" << flsh;
+
+    double h = srk.enthalpy(T1, p1, flsh.x, srk.LIQPH) * flsh.betaL 
+                + srk.enthalpy(T1, p1, flsh.y, srk.VAPPH) * flsh.betaV;
+    std::cout << "Molar enthalpy at (T, p) = (" << T1 << ", " << p1 << ") : " << h << " J / mol\n";
+
+    double p2 = 5e3;
+    std::cout << "\nIsenthalpic decompression to " << p2 / 1e5 << " bar ...\n";
+    flsh = srk.two_phase_phflash(p2, h, z);
+    std::cout << "After isenthalpic decompression: \n" << flsh;
+
+    double s = srk.entropy(flsh.T, flsh.p, flsh.x, srk.LIQPH) * flsh.betaL
+                + srk.entropy(flsh.T, flsh.p, flsh.y, srk.VAPPH) * flsh.betaV;
+    std::cout << "Molar entropy after decompression : " << s << " J / mol K\n";
+    std::cout << "\nIsentropic compression to 2 bar ...\n";
+    flsh = srk.two_phase_psflash(2e5, s, z);
+    std::cout << "After isentropic compression : \n" << flsh;
+
+    double v = srk.specific_volume(flsh.T, flsh.p, flsh.x, srk.LIQPH) * flsh.betaL
+                + srk.specific_volume(flsh.T, flsh.p, flsh.y, srk.VAPPH) * flsh.betaV;
+    double u = srk.internal_energy_tv(flsh.T, v, flsh.z);
+    std::cout << "Internal energy after isentropic compression : " << u << " J / mol\n";
+
+    double q = 30.; // J / mol, Heat to add
+    std::cout << "\nIsochoric heating with " << q << " J / mol heat ...\n";
+    flsh = srk.two_phase_uvflash(u + q, v, z);
+    std::cout << "After isochoric heating : \n" << flsh;
+    return 0;
+}

addon/cppExamples/saturation.cpp

@@ -0,0 +1,27 @@
+#include "cppThermopack/thermo.h"
+#include "cppThermopack/pcsaft.h"
+#include <vector>
+#include <iostream>
+
+int main(){
+    Pcsaft pcs("C3,NC6,ETOH"); // Initialize pcsaft for propane/hexane/ethanol mixture
+
+    std::vector<double> z = {0.2, 0.5, 0.3}; // Total composition
+
+    std::vector<double> crit_tvp = pcs.critical(z);
+
+    double T{0.6 * crit_tvp[0]}; // Ensuring that we compute at some sub-critical temperature (0.6 * Tc)
+    std::cout << "Computing dew- and bubble pressure at " << T << " K:\n";
+    double p_dew = pcs.dew_pressure(T, z).p;
+    double p_bub = pcs.bubble_pressure(T, z).p;
+
+    std::cout << pcs.dew_pressure(T, z) << "\n";
+    std::cout << pcs.bubble_pressure(T, z) << "\n";
+
+    double p = 0.5 * (p_dew + p_bub);
+    std::cout << "\nComputing dew- and bubble temperature at " << p / 1e5 << " bar\n";
+    std::cout << pcs.dew_temperature(p, z) << "\n";
+    std::cout << pcs.bubble_temperature(p, z) << "\n";
+
+    return 0;
+}

addon/cppExamples/tp_properties.cpp

@@ -0,0 +1,23 @@
+#include <cppThermopack/saftvrmie.h>
+#include <iostream>
+#include <vector>
+
+int main(){
+    std::string comps = "C3,C1";
+    Saftvrmie svrm(comps);
+    std::cout << "Initialized SAFT-VR Mie for " << comps << "\n";
+    double T = 300;
+    double p = 1e6;
+    std::vector<double> z = {0.7, 0.3};
+
+    double h = svrm.enthalpy(T, p, z, svrm.VAPPH);
+    std::cout << "Vapour phase specific enthalpy : " << h << " J / mol\n";
+    std::cout << "Liquid phase specific enthalpy : " << static_cast<double>(svrm.enthalpy(T, p, z, svrm.LIQPH)) << " J / mol\n";
+
+    Property s_vap = svrm.entropy(T, p, z, svrm.VAPPH, false, false, true);
+    std::cout << "Vapour phase partial molar entropies : " << s_vap.dn()[0] << ", " << s_vap.dn()[0] << " J / mol^2 K\n";
+    Property s_liq = svrm.entropy(T, p, z, svrm.LIQPH, false, false, true);
+    std::cout << "Liquid phase partial molar entropies : " << s_liq.dn()[0] << ", " << s_liq.dn()[1] << " J / mol^2 K\n";
+
+    return 0;
+}

addon/cppExamples/tv_properties.cpp

@@ -0,0 +1,59 @@
+#include <cppThermopack/pcsaft.h>
+#include <iostream>
+#include <vector>
+
+// A couple convenience methods
+inline std::vector<double> operator*(double d, const std::vector<double>& v){
+    std::vector<double> v_out{v};
+    for (double& vi : v_out) vi *= d;
+    return v_out;
+}
+
+inline std::vector<double> operator*(const std::vector<double>& v, double d){
+    return d * v;
+}
+
+inline std::ostream& operator<<(std::ostream& os, const std::vector<double> v){
+    for (auto d : v) os << d << ", ";
+    return os;
+}
+
+// This is the example
+int main(){
+    std::string comps = "C1,ETOH,PROP1OL";
+    Pcsaft pcs(comps);
+    std::cout << "Initialized PC-SAFT for components " << comps << "\n";
+
+    double T{400.}, V{0.3};
+    double n_tot = 15;
+    std::vector<double> z = {0.2, 0.5, 0.5};
+    std::vector<double> n = n_tot * z;
+
+    std::cout << "Computing TV properties at T = " << T << " K, V = " << V << " m^3, n_tot = " << n_tot << " mol ... \n\n";
+
+    std::cout << "Pressure : " << pcs.pressure_tv(T, V, n) / 1e5 << " bar\n";
+    std::cout << "Internal energy : " << pcs.internal_energy_tv(T, V, n) / 1e3 << " kJ\n";
+    std::cout << "Entropy : " << pcs.entropy_tv(T, V, n) << " J / K\n";
+    std::cout << "Enthalpy : " << pcs.enthalpy_tv(T, V, n) / 1e3 << " kJ\n";
+    std::cout << "Helmholtz energy : " << pcs.helmholtz_tv(T, V, n) / 1e3 << " kJ\n";
+    std::cout << "Chemical potential : " << pcs.chemical_potential_tv(T, V, n) << " J / mol\n";
+    std::cout << "Logarithm of fugacity coeff. : " << pcs.fugacity_tv(T, V, n) << "\n";
+
+    std::cout << "\nComputing derivatives ...\n";
+    // Bools mark what derivatives to compute. 
+    // Order is dt, dv, dn
+    Property H = pcs.enthalpy_tv(T, V, n, true, true, true);
+    std::cout << "dHdT : " << H.dt() << " J / K, dHdV : " << H.dv() << " J / m^3, dHdn_i : " << H.dn() << "\n";
+
+    Property A = pcs.helmholtz_tv(T, V, n, false, false, true);
+    std::cout << "dAdn : " << A.dn() << "\n";
+    VectorProperty mu = pcs.chemical_potential_tv(T, V, n, true, false, true);
+    std::cout << "mu   : " << mu.value() << "\n";
+    std::cout << "Chemical potential with derivatives : \n" << mu << "\n";
+
+    VectorProperty lnphi = pcs.fugacity_tv(T, V, n, false, true, true);
+    std::cout << "Logarithm of fugacity coefficients : \n" << lnphi;
+
+    std::vector<std::vector<double>> dmudn = mu.dn(); // Extracting derivatives
+    return 0;
+}

addon/cppThermopack/cppThermopack/cubic.h

@@ -0,0 +1,43 @@
+#pragma once
+#include <string>
+#include "dllmacros.h"
+#include "thermo.h"
+
+extern "C" {
+    void get_export_name(eoslibinit, init_cubic)(char* comps, char* eos, char* mixing, char* alpha, char* ref,
+						int* volume_shift, size_t comp_strlen, size_t eos_strlen, size_t mixing_strlen,
+						size_t alpha_strlen, size_t ref_strlen);
+
+}
+
+class Cubic : public Thermo{
+    public:
+    Cubic(std::string comps, std::string eos,
+	    std::string mixing="vdW", std::string alpha="Classic",
+	    std::string ref="Default", bool volume_shift=false)
+	: Thermo(comps)
+    {
+	activate();
+	int vol_shift_int = volume_shift ? true_int : 0;
+	get_export_name(eoslibinit, init_cubic)(comps.data(), eos.data(), mixing.data(), alpha.data(), ref.data(),
+	    &vol_shift_int, comps.size(), eos.size(), mixing.size(), alpha.size(), ref.size());
+
+    }
+
+};
+
+class PengRobinson : public Cubic{
+    public:
+    PengRobinson(std::string comps, std::string mixing="vdW", std::string alpha="Classic",
+	    std::string ref="Default", bool volume_shift=false)
+	: Cubic(comps, "PR", mixing, alpha, ref, volume_shift)
+	{}
+};
+
+class SoaveRedlichKwong : public Cubic{
+    public:
+    SoaveRedlichKwong(std::string comps, std::string mixing="vdW", std::string alpha="Classic",
+	    std::string ref="Default", bool volume_shift=false)
+	: Cubic(comps, "SRK", mixing, alpha, ref, volume_shift)
+	{}
+};