From 3d30449c107f046b84a5a3e1fd0bb600bb98608f Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Sun, 27 Mar 2022 14:24:37 -0400
Subject: [PATCH] C++ interface (#10)
* A fully encapsulated library wrapper with C++ interface
* Refactor types.hpp to hide the big headers in a macro conditional section
---
CMakeLists.txt | 21 ++++++++++++--
include/teqp/json_builder.hpp | 18 ++----------
include/teqp/models/cubics.hpp | 8 +++---
include/teqp/models/fwd.hpp | 32 +++++++++++++++++++++
include/teqp/models/multifluid_reducing.hpp | 24 +++++-----------
include/teqp/models/pcsaft.hpp | 7 +++--
include/teqp/types.hpp | 25 ++++++++++++----
interface/CPP/teqpcpp.cpp | 17 +++++++++++
interface/CPP/teqpcpp.hpp | 15 ++++++++++
interface/CPP/test_teqpcpp.cpp | 17 +++++++++++
10 files changed, 135 insertions(+), 49 deletions(-)
create mode 100644 include/teqp/models/fwd.hpp
create mode 100644 interface/CPP/teqpcpp.cpp
create mode 100644 interface/CPP/teqpcpp.hpp
create mode 100644 interface/CPP/test_teqpcpp.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 85fecd4..83a4c8e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -41,16 +41,19 @@ set(AUTODIFF_BUILD_PYTHON FALSE CACHE BOOL "No autodiff python")
set(AUTODIFF_BUILD_EXAMPLES FALSE CACHE BOOL "No autodiff examples")
add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/externals/autodiff")
-# Find eigen library
-# find_package(Eigen3 REQUIRED)
-
# Turn on more useful diagnostic messages in nlohmann::json, for instance if you are accessing a field that doesn't exist
set(JSON_Diagnostics TRUE CACHE BOOL "Turn on more helpful diagnostics in nlohmann::json")
+
+
option (TEQP_NO_PYTHON
"Enable to NOT build the Python interface"
OFF)
+option (TEQP_TESTTEQPCPP
+ "Enable to add a target with a test of the C++ interface"
+ OFF)
+
option (TEQP_TEQPC
"Enable to build the shared library with extern \"C\" interface"
OFF)
@@ -71,6 +74,18 @@ option (TEQP_MULTIPRECISION_ENABLED
"Enable the use of boost::multiprecision"
OFF)
+# Add a static library with the C++ interface that uses only STL
+# types so that recompilation of a library that uses teqp
+# doesn't require a full compile for a single LOC change
+add_library(teqpcpp STATIC "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/teqpcpp.cpp")
+target_link_libraries(teqpcpp PUBLIC teqpinterface PUBLIC autodiff)
+
+
+if (TEQP_TESTTEQPCPP)
+ add_executable(test_teqpcpp "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/test_teqpcpp.cpp")
+ target_link_libraries(test_teqpcpp PUBLIC teqpcpp)
+endif()
+
if (TEQP_JAVASCRIPT_MODULE)
# cmake -DTEQP_JAVASCRIPT_MODULE=ON
# -DCMAKE_TOOLCHAIN_FILE=${EMSCRIPTEN}/cmake/Platform/Emscripten.cmake
diff --git a/include/teqp/json_builder.hpp b/include/teqp/json_builder.hpp
index 60167a3..a7664cf 100644
--- a/include/teqp/json_builder.hpp
+++ b/include/teqp/json_builder.hpp
@@ -1,10 +1,6 @@
#pragma once
-#include "teqp/models/vdW.hpp"
-#include "teqp/models/cubics.hpp"
-#include "teqp/models/CPA.hpp"
-#include "teqp/models/pcsaft.hpp"
-#include "teqp/models/multifluid.hpp"
+#include "teqp/models/fwd.hpp"
#include "teqp/exceptions.hpp"
@@ -12,17 +8,7 @@
namespace teqp {
- using vad = std::valarray<double>;
-
- // Define the EOS types by interrogating the types returned by the respective factory function
- using cub = decltype(canonical_PR(vad{}, vad{}, vad{}));
- using cpatype = decltype(CPA::CPAfactory(nlohmann::json{}));
- using pcsafttype = decltype(PCSAFT::PCSAFTfactory(nlohmann::json{}));
- using multifluidtype = decltype(multifluidfactory(nlohmann::json{}));
-
- using AllowedModels = std::variant<vdWEOS1, cub, cpatype, pcsafttype, multifluidtype>;
-
- AllowedModels build_model(const nlohmann::json& json) {
+ static AllowedModels build_model(const nlohmann::json& json) {
// Extract the name of the model and the model parameters
std::string kind = json.at("kind");
diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp
index 9d47cf3..e8a4371 100644
--- a/include/teqp/models/cubics.hpp
+++ b/include/teqp/models/cubics.hpp
@@ -29,7 +29,7 @@ public:
template<typename TType>
auto operator () (const TType& T) const {
- return forceeval(powi(forceeval(1.0 + mi * (1.0 - sqrt(T / Tci))), 2));
+ return forceeval(pow2(forceeval(1.0 + mi * (1.0 - sqrt(T / Tci)))));
}
};
@@ -60,7 +60,7 @@ public:
ai.resize(Tc_K.size());
bi.resize(Tc_K.size());
for (auto i = 0; i < Tc_K.size(); ++i) {
- ai[i] = OmegaA * powi(Ru * Tc_K[i], 2) / pc_Pa[i];
+ ai[i] = OmegaA * pow2(Ru * Tc_K[i]) / pc_Pa[i];
bi[i] = OmegaB * Ru * Tc_K[i] / pc_Pa[i];
}
k = std::valarray<std::valarray<NumType>>(std::valarray<NumType>(0.0, Tc_K.size()), Tc_K.size());
@@ -169,10 +169,10 @@ auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric) {
std::vector<AlphaFunctionOptions> alphas;
for (auto i = 0; i < Tc_K.size(); ++i) {
if (acentric[i] < 0.491) {
- m[i] = 0.37464 + 1.54226*acentric[i] - 0.26992*powi(acentric[i], 2);
+ m[i] = 0.37464 + 1.54226*acentric[i] - 0.26992*pow2(acentric[i]);
}
else {
- m[i] = 0.379642 + 1.48503*acentric[i] -0.164423*powi(acentric[i], 2) + 0.016666*powi(acentric[i], 3);
+ m[i] = 0.379642 + 1.48503*acentric[i] -0.164423*pow2(acentric[i]) + 0.016666*pow3(acentric[i]);
}
alphas.emplace_back(BasicAlphaFunction(Tc_K[i], m[i]));
}
diff --git a/include/teqp/models/fwd.hpp b/include/teqp/models/fwd.hpp
new file mode 100644
index 0000000..9a65966
--- /dev/null
+++ b/include/teqp/models/fwd.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+/**
+* The name of this file is currently a bit of a misnomer while
+* we think about whether it is possible to forward-declare the models
+*
+* It seems like perhaps that is not possible. For now this header
+* just provides the main variant with the model definitions.
+*/
+
+#include <valarray>
+#include <variant>
+
+#include "teqp/models/vdW.hpp"
+#include "teqp/models/cubics.hpp"
+#include "teqp/models/CPA.hpp"
+#include "teqp/models/pcsaft.hpp"
+#include "teqp/models/multifluid.hpp"
+
+namespace teqp {
+
+ using vad = std::valarray<double>;
+
+ // Define the EOS types by interrogating the types returned by the respective factory function
+ using AllowedModels = std::variant<
+ vdWEOS1,
+ decltype(canonical_PR(vad{}, vad{}, vad{})),
+ decltype(CPA::CPAfactory(nlohmann::json{})),
+ decltype(PCSAFT::PCSAFTfactory(nlohmann::json{})),
+ decltype(multifluidfactory(nlohmann::json{}))
+ >;
+}
\ No newline at end of file
diff --git a/include/teqp/models/multifluid_reducing.hpp b/include/teqp/models/multifluid_reducing.hpp
index e51c571..01c9d95 100644
--- a/include/teqp/models/multifluid_reducing.hpp
+++ b/include/teqp/models/multifluid_reducing.hpp
@@ -157,15 +157,6 @@ namespace teqp {
private:
Eigen::MatrixXd YT, Yv;
- template <typename Num>
- auto cube(Num x) const {
- return forceeval(x * x * x);
- }
- template <typename Num>
- auto square(Num x) const {
- return forceeval(x * x);
- }
-
public:
const Eigen::MatrixXd betaT, gammaT, betaV, gammaV;
const Eigen::ArrayXd Tc, vc;
@@ -185,8 +176,8 @@ namespace teqp {
for (auto j = i + 1; j < N; ++j) {
YT(i, j) = betaT(i, j) * gammaT(i, j) * sqrt(Tc[i] * Tc[j]);
YT(j, i) = betaT(j, i) * gammaT(j, i) * sqrt(Tc[i] * Tc[j]);
- Yv(i, j) = 1.0 / 8.0 * betaV(i, j) * gammaV(i, j) * cube(cbrt(vc[i]) + cbrt(vc[j]));
- Yv(j, i) = 1.0 / 8.0 * betaV(j, i) * gammaV(j, i) * cube(cbrt(vc[i]) + cbrt(vc[j]));
+ Yv(i, j) = 1.0 / 8.0 * betaV(i, j) * gammaV(i, j) * pow3(cbrt(vc[i]) + cbrt(vc[j]));
+ Yv(j, i) = 1.0 / 8.0 * betaV(j, i) * gammaV(j, i) * pow3(cbrt(vc[i]) + cbrt(vc[j]));
}
}
}
@@ -197,13 +188,13 @@ namespace teqp {
auto N = z.size();
typename MoleFractions::value_type sum1 = 0.0;
for (auto i = 0; i < N; ++i) {
- sum1 = sum1 + square(z[i]) * Yc[i];
+ sum1 = sum1 + pow2(z[i]) * Yc[i];
}
typename MoleFractions::value_type sum2 = 0.0;
for (auto i = 0; i < N - 1; ++i) {
for (auto j = i + 1; j < N; ++j) {
- sum2 = sum2 + 2.0 * z[i] * z[j] * (z[i] + z[j]) / (square(beta(i, j)) * z[i] + z[j]) * Yij(i, j);
+ sum2 = sum2 + 2.0 * z[i] * z[j] * (z[i] + z[j]) / (pow2(beta(i, j)) * z[i] + z[j]) * Yij(i, j);
}
}
@@ -217,8 +208,7 @@ namespace teqp {
class MultiFluidInvariantReducingFunction {
private:
Eigen::MatrixXd YT, Yv;
- template <typename Num> auto cube(Num x) const { return x * x * x; }
- template <typename Num> auto square(Num x) const { return x * x; }
+
public:
const Eigen::MatrixXd phiT, lambdaT, phiV, lambdaV;
const Eigen::ArrayXd Tc, vc;
@@ -238,8 +228,8 @@ namespace teqp {
for (auto j = 0; j < N; ++j) {
YT(i, j) = sqrt(Tc[i] * Tc[j]);
YT(j, i) = sqrt(Tc[i] * Tc[j]);
- Yv(i, j) = 1.0 / 8.0 * cube(cbrt(vc[i]) + cbrt(vc[j]));
- Yv(j, i) = 1.0 / 8.0 * cube(cbrt(vc[i]) + cbrt(vc[j]));
+ Yv(i, j) = 1.0 / 8.0 * pow3(cbrt(vc[i]) + cbrt(vc[j]));
+ Yv(j, i) = 1.0 / 8.0 * pow3(cbrt(vc[i]) + cbrt(vc[j]));
}
}
}
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index ee8f4e8..9fd6d3f 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -219,11 +219,12 @@ public:
auto get_sigma_Angstrom() { return sigma_Angstrom; }
auto get_epsilon_over_k_K() { return epsilon_over_k; }
- void print_info() {
- std::cout << "i m sigma / A e/kB / K \n ++++++++++++++" << std::endl;
+ auto print_info() {
+ std::string s = std::string("i m sigma / A e/kB / K \n ++++++++++++++") + "\n";
for (auto i = 0; i < m.size(); ++i) {
- std::cout << i << " " << m[i] << " " << sigma_Angstrom[i] << " " << epsilon_over_k[i] << std::endl;
+ s += std::to_string(i) + " " + std::to_string(m[i]) + " " + std::to_string(sigma_Angstrom[i]) + " " + std::to_string(epsilon_over_k[i]) + "\n";
}
+ return s;
}
template<typename VecType>
double max_rhoN(double T, const VecType& mole_fractions) {
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index 1db3936..c9ab7ba 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -1,10 +1,23 @@
#pragma once
+#include <vector>
+#include "Eigen/Dense"
+
+namespace teqp {
+ template<typename T> inline T pow2(const T& x) { return x * x; }
+ template<typename T> inline T pow3(const T& x) { return x * x * x; }
+
+ inline auto toeig(const std::vector<double>& v) -> Eigen::ArrayXd { return Eigen::Map<const Eigen::ArrayXd>(&(v[0]), v.size()); }
+}
+
+// Everything inside this NO_TYPES_HEADER section involves slow and large headers
+// so if you are just defining types (but not their implementations), you can avoid the
+// compilation speed hit invoked by these headers
+#if !defined(NO_TYPES_HEADER)
+
#include "MultiComplex/MultiComplex.hpp"
-#include <vector>
#include <valarray>
-#include <set>
#include <chrono>
#if defined(TEQP_MULTIPRECISION_ENABLED)
@@ -17,8 +30,6 @@
#include <autodiff/forward/dual/eigen.hpp>
using namespace autodiff;
-#include "Eigen/Dense"
-
namespace teqp {
// Registration of types that are considered to be containers
@@ -77,7 +88,7 @@ namespace teqp {
}
}
- inline auto toeig (const std::vector<double>& v) -> Eigen::ArrayXd { return Eigen::Map<const Eigen::ArrayXd>(&(v[0]), v.size()); }
+
class Timer {
private:
@@ -175,4 +186,6 @@ namespace teqp {
return o;
}
-}; // namespace teqp
\ No newline at end of file
+}; // namespace teqp
+
+#endif
\ No newline at end of file
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
new file mode 100644
index 0000000..45c8356
--- /dev/null
+++ b/interface/CPP/teqpcpp.cpp
@@ -0,0 +1,17 @@
+#include "teqpcpp.hpp"
+#include "teqp/derivs.hpp"
+#include "teqp/algorithms/critical_tracing.hpp"
+
+double teqp::cppinterface::get_Arxy(const teqp::AllowedModels&model, const int NT, const int ND, const double T, const double rho, const Eigen::ArrayXd &molefracs){
+ return std::visit([&](const auto& model) {
+ using tdx = teqp::TDXDerivatives<decltype(model), double, std::decay_t<decltype(molefracs)>>;
+ return tdx::get_Ar(NT, ND, model, T, rho, molefracs);
+ }, model);
+}
+
+nlohmann::json teqp::cppinterface::trace_critical_arclength_binary(const teqp::AllowedModels& model, const double T0, const Eigen::ArrayXd &rhovec0) {
+ return std::visit([&](const auto& model) {
+ using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec0)>>;
+ return crit::trace_critical_arclength_binary(model, T0, rhovec0, "");
+ }, model);
+}
\ No newline at end of file
diff --git a/interface/CPP/teqpcpp.hpp b/interface/CPP/teqpcpp.hpp
new file mode 100644
index 0000000..207f848
--- /dev/null
+++ b/interface/CPP/teqpcpp.hpp
@@ -0,0 +1,15 @@
+#pragma once
+
+#include <valarray>
+
+// Pulls in the AllowedModels variant
+#include "teqp/models/fwd.hpp"
+
+namespace teqp {
+ namespace cppinterface {
+
+ // Wrapper functions
+ double get_Arxy(const AllowedModels& model, const int NT, const int ND, const double T, const double rho, const Eigen::ArrayXd& molefracs);
+ nlohmann::json trace_critical_arclength_binary(const AllowedModels& model, const double T0, const Eigen::ArrayXd& rhovec0);
+ }
+}
\ No newline at end of file
diff --git a/interface/CPP/test_teqpcpp.cpp b/interface/CPP/test_teqpcpp.cpp
new file mode 100644
index 0000000..1bcfa9b
--- /dev/null
+++ b/interface/CPP/test_teqpcpp.cpp
@@ -0,0 +1,17 @@
+#define NO_TYPES_HEADER
+#include "teqpcpp.hpp"
+
+using namespace teqp::cppinterface;
+
+int main(){
+ teqp::AllowedModels model = teqp::build_multifluid_model({ "Methane","Ethane" }, "../mycp");
+
+ //teqp::AllowedModels model = teqp::vdWEOS1(3, 0.1);
+
+ auto z = (Eigen::ArrayXd(2) << 0.5, 0.5).finished( );
+
+ double Ar01 = get_Arxy(model, 0, 1, 300, 3, z);
+
+ double Tc1 = 371.1;
+ auto cr = trace_critical_arclength_binary(model, Tc1, z);
+}
\ No newline at end of file
--
GitLab