diff --git a/.gitignore b/.gitignore
index 8ad801ef9d08adbc64a2a35fd431d083ffd262a4..b0d1bc9682a3e9bec185b874dfd5b6fbfd5f6f8d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
/bld
+/bldgcc
/build
/interface/teqpversion.hpp
/pypirc
diff --git a/dev/docker/boost_bcp_docker/Dockerfile b/dev/docker/boost_bcp_docker/Dockerfile
index fa45df86246522b769b89201292ba85aea3f3161..93e2f008038fbf92801bf72a2d4aa0a5b0dad930 100644
--- a/dev/docker/boost_bcp_docker/Dockerfile
+++ b/dev/docker/boost_bcp_docker/Dockerfile
@@ -5,15 +5,15 @@ FROM ubuntu:20.04
RUN apt-get -y -m update && DEBIAN_FRONTEND=noninteractive apt-get install -y cmake g++ git zip wget xz-utils
RUN mkdir /boost && \
- wget -c https://boostorg.jfrog.io/artifactory/main/release/1.77.0/source/boost_1_77_0.tar.gz -O - | tar -xz -C /boost && \
- cd /boost/boost_1_77_0/ && \
+ wget -c --no-check-certificate https://boostorg.jfrog.io/artifactory/main/release/1.82.0/source/boost_1_82_0.tar.gz -O - | tar -xz -C /boost && \
+ cd /boost/boost_1_82_0/ && \
./bootstrap.sh && \
./b2 tools/bcp
-WORKDIR /boost/boost_1_77_0
+WORKDIR /boost/boost_1_82_0
RUN mkdir /boost_teqp && \
bin.v2/tools/bcp/gcc-9/release/link-static/bcp multiprecision/cpp_bin_float.hpp multiprecision/cpp_complex.hpp multiprecision/eigen.hpp functional/hash.hpp numeric/odeint.hpp typeof/incr_registration_group.hpp mp11.hpp algorithm/string/join.hpp asio/thread_pool.hpp asio/post.hpp /boost_teqp && \
zip -r /boost_teqp.zip /boost_teqp && \
tar cJf /boost_teqp.tar.xz /boost_teqp
-CMD cp /*.tar.xz /output
\ No newline at end of file
+CMD cp /*.tar.xz /output
diff --git a/dev/docker/boost_bcp_docker/boost_teqp.tar.xz b/dev/docker/boost_bcp_docker/boost_teqp.tar.xz
index a73fe375fb8a8245e6893bc29fcbb7f3df9c8c99..64568da7a176b86508b46ecba1d666154d2df57c 100644
Binary files a/dev/docker/boost_bcp_docker/boost_teqp.tar.xz and b/dev/docker/boost_bcp_docker/boost_teqp.tar.xz differ
diff --git a/dev/docker/clangdev/Dockerfile b/dev/docker/clangdev/Dockerfile
new file mode 100644
index 0000000000000000000000000000000000000000..6c43c4517dc212290ef888b00ac0539ed130197e
--- /dev/null
+++ b/dev/docker/clangdev/Dockerfile
@@ -0,0 +1,13 @@
+## Just use docker-compose up to run
+
+# Or, pick a different tag (change the tag) to use a different version of clang
+FROM silkeh/clang:16
+
+RUN apt-get -y -m update && DEBIAN_FRONTEND=noninteractive apt-get install -y cmake git zip nano gcovr ninja-build
+
+RUN DEBIAN_FRONTEND=noninteractive apt-get install -y python3 libpython3-dev python3.9-distutils
+
+# Run the catch exe, generating gcov output
+CMD mkdir build && cd build && \
+ cmake /teqp -GNinja -DCMAKE_BUILD_TYPE=Release && \
+ cmake --build . --target teqpcpp --parallel 4
diff --git a/dev/docker/clangdev/docker-compose.yml b/dev/docker/clangdev/docker-compose.yml
new file mode 100644
index 0000000000000000000000000000000000000000..149152336301ab69ebec02b5142d9b983314fc48
--- /dev/null
+++ b/dev/docker/clangdev/docker-compose.yml
@@ -0,0 +1,15 @@
+version: '3.2'
+
+services:
+ app:
+ build:
+ context: ./
+ dockerfile: Dockerfile
+ environment:
+ DEBIAN_FRONTEND: noninteractive
+ mem_limit: 10000000000
+ volumes:
+ - type: bind
+ source: ../../..
+ target: /teqp
+ read_only: true
diff --git a/doc/source/models/SAFT-VR-Mie-Polar.ipynb b/doc/source/models/SAFT-VR-Mie-Polar.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..0570ca1fe5ed12fbafe7bfe6ceddabf8e2c14f1e
--- /dev/null
+++ b/doc/source/models/SAFT-VR-Mie-Polar.ipynb
@@ -0,0 +1,121 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d1842386",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import teqp\n",
+ "teqp.__version__"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "16f05ec8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas\n",
+ "import matplotlib.pyplot as plt\n",
+ "import CoolProp.CoolProp as CP\n",
+ "import scipy.integrate\n",
+ "import math"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "be8268a7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ek = 100\n",
+ "sigma_m = 3e-10\n",
+ "factor = 1.0/(4*math.pi*8.8541878128e-12*1.380649e-23)\n",
+ " \n",
+ "N_A = 6.022e23\n",
+ "fig, (ax1, ax2) = plt.subplots(2, 1)\n",
+ "\n",
+ "# # From https://arxiv.org/pdf/mtrl-th/9501001.pdf which pulled from M. van Leeuwen and B. Smit, Phys. Rev. Lett. 71, 3991 (1993)\n",
+ "# mustar2 = [2.5, 3.0, 3.5, 4.0]\n",
+ "# T = [2.63, 3.35, 4.20, 5.07]\n",
+ "# rho = [0.29, 0.25, 0.24, 0.24]\n",
+ "# ax1.plot(mustar2, T, 'd')\n",
+ "# ax2.plot(mustar2, rho, 'd')\n",
+ "\n",
+ "mustar2 = [1,2,3,4]\n",
+ "T = [1.41, 1.60, 1.82, 2.06]\n",
+ "rho = [0.30, 0.31, 0.312, 0.289]\n",
+ "ax1.plot(mustar2, T, 's')\n",
+ "ax2.plot(mustar2, rho, 's')\n",
+ "\n",
+ "for polar_model in ['GrossVrabec','GubbinsTwu+GubbinsTwu','GubbinsTwu+Luckas']:\n",
+ " # Comparing with Hentschke, DOI: https://doi.org/10.1103/physreve.75.011506\n",
+ " x = []; y = []; TT = []; DD = []\n",
+ " rhostar_guess = 0.27\n",
+ " Tstar_guess = 1.5\n",
+ " for mustar2 in np.arange(0.001, 15, 0.1):\n",
+ " z = np.array([1.0])\n",
+ " model = teqp.make_model({\n",
+ " \"kind\": 'SAFT-VR-Mie',\n",
+ " \"model\": {\n",
+ " \"polar_model\": polar_model,\n",
+ " \"coeffs\": [{\n",
+ " \"name\": \"Stockmayer\",\n",
+ " \"BibTeXKey\": \"me\",\n",
+ " \"m\": 1.0,\n",
+ " \"epsilon_over_k\": ek, # [K]\n",
+ " \"sigma_m\": sigma_m,\n",
+ " \"lambda_r\": 12.0,\n",
+ " \"lambda_a\": 6.0,\n",
+ " \"(mu^*)^2\": mustar2,\n",
+ " \"nmu\": 1\n",
+ " }]\n",
+ " }\n",
+ " })\n",
+ "\n",
+ " T, rho = model.solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*sigma_m**3))\n",
+ " # Store the values\n",
+ " x.append(mustar2)\n",
+ " TT.append(T/ek)\n",
+ " DD.append(rho*N_A*sigma_m**3)\n",
+ " # Update the guess for the next calculation\n",
+ " Tstar_guess = TT[-1]\n",
+ " rhostar_guess = DD[-1]\n",
+ "\n",
+ " ax1.plot(x, TT, label=polar_model)\n",
+ " ax2.plot(x, DD)\n",
+ " \n",
+ "ax1.legend(loc='best')\n",
+ "ax1.set(ylabel=r'$T^*$')\n",
+ "ax2.set(ylabel=r'$\\rho^*$', xlabel=r'$(\\mu^*)^2$')\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 61aff76f4bc05caca67806fb4d907d25b62678b7..ec495583c966792cc97a437635cff50175c4ce3d 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -6,6 +6,7 @@
#include "teqp/algorithms/critical_tracing.hpp"
#include "teqp/algorithms/critical_pure.hpp"
#include "teqp/algorithms/VLE_types.hpp"
+#include "teqp/algorithms/VLE_pure.hpp"
#include <Eigen/Dense>
// Imports from boost for numerical integration
@@ -55,223 +56,6 @@ auto linsolve(const A& a, const B& b) {
return a.matrix().colPivHouseholderQr().solve(b.matrix()).array().eval();
}
-template<typename Model, typename TYPE = double, ADBackends backend = ADBackends::autodiff>
-class IsothermPureVLEResiduals {
-public:
- using EigenArray = Eigen::Array<TYPE, 2, 1>;
- using EigenArray1 = Eigen::Array<TYPE, 1, 1>;
- using EigenMatrix = Eigen::Array<TYPE, 2, 2>;
-private:
- const Model& m_model;
- const TYPE m_T;
- const Eigen::ArrayXd molefracs;
- EigenMatrix J;
- EigenArray y;
-public:
- std::size_t icall = 0;
- double Rr, R0;
-
- IsothermPureVLEResiduals(const Model& model, const TYPE& T, const std::optional<Eigen::ArrayXd>& molefracs_ = std::nullopt) : m_model(model), m_T(T),
- molefracs( (molefracs_) ? molefracs_.value() : Eigen::ArrayXd::Ones(1,1)) {
- Rr = m_model.R(molefracs);
- R0 = m_model.R(molefracs);
- };
-
- const auto& get_errors() { return y; };
-
- auto call(const EigenArray& rhovec) {
- assert(rhovec.size() == 2);
-
- const EigenArray1 rhovecL = rhovec.head(1);
- const EigenArray1 rhovecV = rhovec.tail(1);
- const auto rhomolarL = rhovecL.sum(), rhomolarV = rhovecV.sum();
-
- using tdx = TDXDerivatives<Model,TYPE,Eigen::ArrayXd>;
-
- const TYPE &T = m_T;
- //const TYPE R = m_model.R(molefracs);
- double R0_over_Rr = R0 / Rr;
-
- auto derL = tdx::template get_Ar0n<2, backend>(m_model, T, rhomolarL, molefracs);
- auto pRTL = rhomolarL*(R0_over_Rr + derL[1]); // p/(R*T)
- auto dpRTLdrhoL = R0_over_Rr + 2*derL[1] + derL[2];
- auto hatmurL = derL[1] + derL[0] + R0_over_Rr*log(rhomolarL);
- auto dhatmurLdrho = (2*derL[1] + derL[2])/rhomolarL + R0_over_Rr/rhomolarL;
-
- auto derV = tdx::template get_Ar0n<2, backend>(m_model, T, rhomolarV, molefracs);
- auto pRTV = rhomolarV*(R0_over_Rr + derV[1]); // p/(R*T)
- auto dpRTVdrhoV = R0_over_Rr + 2*derV[1] + derV[2];
- auto hatmurV = derV[1] + derV[0] + R0_over_Rr *log(rhomolarV);
- auto dhatmurVdrho = (2*derV[1] + derV[2])/rhomolarV + R0_over_Rr/rhomolarV;
-
- y(0) = pRTL - pRTV;
- J(0, 0) = dpRTLdrhoL;
- J(0, 1) = -dpRTVdrhoV;
-
- y(1) = hatmurL - hatmurV;
- J(1, 0) = dhatmurLdrho;
- J(1, 1) = -dhatmurVdrho;
-
- icall++;
- return y;
- }
- auto Jacobian(const EigenArray& rhovec){
- return J;
- }
- //auto numJacobian(const EigenArray& rhovec) {
- // EigenArray plus0 = rhovec, plus1 = rhovec;
- // double dr = 1e-6 * rhovec[0];
- // plus0[0] += dr; plus1[1] += dr;
- // EigenMatrix J;
- // J.col(0) = (call(plus0) - call(rhovec))/dr;
- // J.col(1) = (call(plus1) - call(rhovec))/dr;
- // return J;
- //}
-};
-
-template<typename Residual, typename Scalar>
-auto do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
- using EArray = Eigen::Array<Scalar, 2, 1>;
- auto rhovec = (EArray() << rhoL, rhoV).finished();
- auto r0 = resid.call(rhovec);
- auto J = resid.Jacobian(rhovec);
- for (int iter = 0; iter < maxiter; ++iter){
- if (iter > 0) {
- r0 = resid.call(rhovec);
- J = resid.Jacobian(rhovec);
- }
- auto v = J.matrix().colPivHouseholderQr().solve(-r0.matrix()).array().eval();
- auto rhovecnew = (rhovec + v).eval();
- //double r00 = static_cast<double>(r0[0]);
- //double r01 = static_cast<double>(r0[1]);
-
- // If the solution has stopped improving, stop. The change in rhovec is equal to v in infinite precision, but
- // not when finite precision is involved, use the minimum non-denormal float as the determination of whether
- // the values are done changing
- auto minval = std::numeric_limits<Scalar>::epsilon();
- //double minvaldbl = static_cast<double>(minval);
- if (((rhovecnew - rhovec).cwiseAbs() < minval).all()) {
- break;
- }
- if ((r0.cwiseAbs() < minval).all()) {
- break;
- }
- rhovec = rhovecnew;
- }
- return rhovec;
-}
-
-template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
-auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter, const std::optional<Eigen::ArrayXd>& molefracs = std::nullopt) {
- Eigen::ArrayXd molefracs_{Eigen::ArrayXd::Ones(1,1)};
- if (molefracs){ molefracs_ = molefracs.value(); }
- auto res = IsothermPureVLEResiduals<Model, Scalar, backend>(model, T, molefracs_);
- return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
-}
-
-/***
- * \brief Calculate the derivative of vapor pressure with respect to temperature
- * \param model The model to operate on
- * \param T Temperature
- * \param rhoL Liquid density
- * \param rhoV Vapor density
- *
- * Based upon
- * \f[
- * \frac{dp_{\sigma}}{dT} = \frac{h''-h'}{T(v''-v')} = \frac{s''-s'}{v''-v'}
- * \f]
- * where the \f$h''-h'\f$ is given by the difference in residual enthalpy \f$h''-h' = h^r''-h^r'\f$ because the ideal-gas parts cancel
- */
-template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
-auto dpsatdT_pure(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, const std::optional<Eigen::ArrayXd>& molefracs = std::nullopt) {
-
- auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
- if (molefracs){ molefrac = molefracs.value(); }
-
- using tdx = TDXDerivatives<decltype(model), double, decltype(molefrac)>;
- using iso = IsochoricDerivatives<decltype(model), double, decltype(molefrac)>;
-
- auto R = model.R(molefrac);
-
- auto hrVLERTV = tdx::get_Ar01(model, T, rhoV, molefrac) + tdx::get_Ar10(model, T, rhoV, molefrac);
- auto hrVLERTL = tdx::get_Ar01(model, T, rhoL, molefrac) + tdx::get_Ar10(model, T, rhoL, molefrac);
- auto deltahr_over_T = R*(hrVLERTV-hrVLERTL);
- auto dpsatdT = deltahr_over_T/(1/rhoV-1/rhoL); // From Clausius-Clapeyron; dp/dT = Deltas/Deltav = Deltah/(T*Deltav); Delta=V-L
- return dpsatdT;
-}
-
-/***
- \brief Starting at the critical point, trace the VLE down to a temperature of interest
-
- \note This method only works for well-behaved EOS, notably absent from that category are EOS in the multiparameter category with orthobaric scaling exponent not equal to 0.5 at the critical point. Most other analytical EOS work fine
-
- The JSON data structure defines the variables that need to be specified.
-
- In the current implementation, there are a few steps:
- 1. Solve for the true critical point satisfying \f$(\partial p/\partial \rho)_{T}=(\partial^2p/\partial\rho^2)_{T}=0\f$
- 2. Take a small step away from the critical point (this is where the beta=0.5 assumption is invoked)
- 3. Integrate from the near critical temperature to the temperature of interest
- */
-template<typename Model>
-auto pure_trace_VLE(const Model& model, const double T, const nlohmann::json &spec){
- // Start at the true critical point, from the specified guess value
- nlohmann::json pure_spec;
- Eigen::ArrayXd z{Eigen::ArrayXd::Ones(1,1)};
- if (spec.contains("pure_spec")){
- pure_spec = spec.at("pure_spec");
- z = Eigen::ArrayXd(pure_spec.at("alternative_length").get<int>()); z.setZero();
- z(pure_spec.at("alternative_pure_index").get<int>()) = 1;
- }
-
- auto [Tc, rhoc] = solve_pure_critical(model, spec.at("Tcguess").get<double>(), spec.at("rhocguess").get<double>(), pure_spec);
-
- // Small step towards lower temperature close to critical temperature
- double Tclose = spec.at("Tred").get<double>()*Tc;
- auto rhoLrhoV = extrapolate_from_critical(model, Tc, rhoc, Tclose, z);
- auto rhoLrhoVpolished = pure_VLE_T(model, Tclose, rhoLrhoV[0], rhoLrhoV[1], spec.value("NVLE", 10), z);
- if (rhoLrhoVpolished[0] == rhoLrhoVpolished[1]){
- throw teqp::IterationError("Converged to trivial solution");
- }
-
- // "Integrate" down to temperature of interest
- int Nstep = spec.at("Nstep");
- double R = model.R(z);
- bool with_deriv = spec.at("with_deriv");
- double dT = -(Tclose-T)/(Nstep-1);
- using tdx = TDXDerivatives<decltype(model)>;
- for (auto T_: Eigen::ArrayXd::LinSpaced(Nstep, Tclose, T)){
- rhoLrhoVpolished = pure_VLE_T(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], spec.value("NVLE", 10), z);
-
- //auto pL = rhoLrhoVpolished[0]*R*T_*(1 + tdx::get_Ar01(model, T_, rhoLrhoVpolished[0], z));
- //auto pV = rhoLrhoVpolished[1]*R*T_*(1 + tdx::get_Ar01(model, T_, rhoLrhoVpolished[1], z));
- //std::cout << pL << " " << pV << " " << pL/pV-1 << std::endl;
- if (with_deriv){
- // Get drho/dT for both phases
- auto dpsatdT = dpsatdT_pure(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], z);
- auto get_drhodT = [&z, R, dpsatdT](const Model& model, double T, double rho){
- auto dpdrho = R*T*(1 + 2*tdx::get_Ar01(model, T, rho, z) + tdx::get_Ar02(model, T, rho, z));
- auto dpdT = R*rho*(1 + tdx::get_Ar01(model, T, rho, z) - tdx::get_Ar11(model, T, rho, z));
- return -dpdT/dpdrho + dpsatdT/dpdrho;
- };
- auto drhodTL = get_drhodT(model, T_, rhoLrhoVpolished[0]);
- auto drhodTV = get_drhodT(model, T_, rhoLrhoVpolished[1]);
- // Use the obtained derivative to calculate the step in rho from deltarho = (drhodT)*dT
- auto DeltarhoL = dT*drhodTL, DeltarhoV = dT*drhodTV;
- rhoLrhoVpolished[0] += DeltarhoL;
- rhoLrhoVpolished[1] += DeltarhoV;
- }
-
- // Updated values for densities at new T
- if (!std::isfinite(rhoLrhoVpolished[0])){
- throw teqp::IterationError("The density is no longer valid; try increasing Nstep");
- }
- if (rhoLrhoVpolished[0] == rhoLrhoVpolished[1]){
- throw teqp::IterationError("Converged to trivial solution; try increasing Nstep");
- }
- }
- return rhoLrhoVpolished;
-}
-
/***
* \brief Do a vapor-liquid phase equilibrium problem for a mixture (binary only for now) with mole fractions specified in the liquid phase
* \param model The model to operate on
@@ -289,8 +73,8 @@ auto pure_trace_VLE(const Model& model, const double T, const nlohmann::json &sp
* this component will not be allowed to change (they will stay zero, avoiding the possibility that
* they go to a negative value, which can cause trouble for some EOS)
*/
-template<typename Model, typename Scalar, typename Vector>
-auto mix_VLE_Tx(const Model& model, Scalar T, const Vector& rhovecL0, const Vector& rhovecV0, const Vector& xspec, double atol, double reltol, double axtol, double relxtol, int maxiter) {
+inline auto mix_VLE_Tx(const AbstractModel& model, double T, const Eigen::ArrayXd& rhovecL0, const Eigen::ArrayXd& rhovecV0, const Eigen::ArrayXd& xspec, double atol, double reltol, double axtol, double relxtol, int maxiter) {
+ using Scalar = double;
const Eigen::Index N = rhovecL0.size();
auto lengths = (Eigen::ArrayX<Eigen::Index>(3) << rhovecL0.size(), rhovecV0.size(), xspec.size()).finished();
@@ -300,18 +84,17 @@ auto mix_VLE_Tx(const Model& model, Scalar T, const Vector& rhovecL0, const Vect
Eigen::MatrixXd J(2 * N, 2 * N), r(2 * N, 1), x(2 * N, 1);
x.col(0).array().head(N) = rhovecL0;
x.col(0).array().tail(N) = rhovecV0;
- using isochoric = IsochoricDerivatives<Model, Scalar, Vector>;
Eigen::Map<Eigen::ArrayXd> rhovecL(&(x(0)), N);
Eigen::Map<Eigen::ArrayXd> rhovecV(&(x(0 + N)), N);
- auto RT = model.R(xspec) * T;
+ auto RT = model.get_R(xspec) * T;
VLE_return_code return_code = VLE_return_code::unset;
for (int iter = 0; iter < maxiter; ++iter) {
- auto [PsirL, PsirgradL, hessianL] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecL);
- auto [PsirV, PsirgradV, hessianV] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecV);
+ auto [PsirL, PsirgradL, hessianL] = model.build_Psir_fgradHessian_autodiff(T, rhovecL);
+ auto [PsirV, PsirgradV, hessianV] = model.build_Psir_fgradHessian_autodiff(T, rhovecV);
auto rhoL = rhovecL.sum();
auto rhoV = rhovecV.sum();
Scalar pL = rhoL * RT - PsirL + (rhovecL.array() * PsirgradL.array()).sum(); // The (array*array).sum is a dot product
@@ -422,10 +205,9 @@ struct hybrj_functor__mix_VLE_Tp : Functor<double>
const VectorXd::Index n = x.size() / 2;
Eigen::Map<const Eigen::ArrayXd> rhovecL(&(x(0)), n);
Eigen::Map<const Eigen::ArrayXd> rhovecV(&(x(0 + n)), n);
- auto RT = model.R((rhovecL / rhovecL.sum()).eval()) * T;
- using isochoric = IsochoricDerivatives<Model, Scalar, VectorXd>;
- auto [PsirL, PsirgradL, hessianL] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecL);
- auto [PsirV, PsirgradV, hessianV] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecV);
+ auto RT = model.get_R((rhovecL / rhovecL.sum()).eval()) * T;
+ auto [PsirL, PsirgradL, hessianL] = model.build_Psir_fgradHessian_autodiff(T, rhovecL);
+ auto [PsirV, PsirgradV, hessianV] = model.build_Psir_fgradHessian_autodiff(T, rhovecV);
auto rhoL = rhovecL.sum();
auto rhoV = rhovecV.sum();
Scalar pL = rhoL * RT - PsirL + (rhovecL.array() * PsirgradL.array()).sum(); // The (array*array).sum is a dot product
@@ -458,10 +240,9 @@ struct hybrj_functor__mix_VLE_Tp : Functor<double>
assert(J.rows() == 2*n);
assert(J.cols() == 2*n);
- auto RT = model.R((rhovecL / rhovecL.sum()).eval()) * T;
- using isochoric = IsochoricDerivatives<Model, Scalar, VectorXd>;
- auto [PsirL, PsirgradL, hessianL] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecL);
- auto [PsirV, PsirgradV, hessianV] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecV);
+ auto RT = model.get_R((rhovecL / rhovecL.sum()).eval()) * T;
+ auto [PsirL, PsirgradL, hessianL] = model.build_Psir_fgradHessian_autodiff(T, rhovecL);
+ auto [PsirV, PsirgradV, hessianV] = model.build_Psir_fgradHessian_autodiff(T, rhovecV);
auto dpdrhovecL = RT + (hessianL * rhovecL.matrix()).array();
auto dpdrhovecV = RT + (hessianV * rhovecV.matrix()).array();
@@ -502,8 +283,8 @@ struct hybrj_functor__mix_VLE_Tp : Functor<double>
* \param rhovecV0 Initial values for vapor mole concentrations
* \param flags Flags controlling the iteration and stopping conditions
*/
-template<typename Model, typename Scalar, typename Vector>
-auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhovecL0, const Vector& rhovecV0, const std::optional<MixVLETpFlags>& flags_ = std::nullopt) {
+
+inline auto mix_VLE_Tp(const AbstractModel& model, double T, double pgiven, const Eigen::ArrayXd& rhovecL0, const Eigen::ArrayXd& rhovecV0, const std::optional<MixVLETpFlags>& flags_ = std::nullopt) {
auto flags = flags_.value_or(MixVLETpFlags{});
@@ -519,7 +300,7 @@ auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhove
VLE_return_code return_code = VLE_return_code::unset;
std::string message = "";
- using FunctorType = hybrj_functor__mix_VLE_Tp<Model>;
+ using FunctorType = hybrj_functor__mix_VLE_Tp<AbstractModel>;
FunctorType functor(model, T, pgiven);
Eigen::VectorXd initial_r(2 * N); initial_r.setZero();
functor(x, initial_r);
@@ -615,8 +396,8 @@ auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhove
* \param flags Additional flags
*/
-template<typename Model, typename Scalar, typename Vector>
-auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec, Scalar T0, const Vector& rhovecL0, const Vector& rhovecV0, const std::optional<MixVLEpxFlags>& flags_ = std::nullopt) {
+inline auto mixture_VLE_px(const AbstractModel& model, double p_spec, const Eigen::ArrayXd& xmolar_spec, double T0, const Eigen::ArrayXd& rhovecL0, const Eigen::ArrayXd& rhovecV0, const std::optional<MixVLEpxFlags>& flags_ = std::nullopt) {
+ using Scalar = double;
auto flags = flags_.value_or(MixVLEpxFlags{});
@@ -636,7 +417,6 @@ auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec
x(0) = T0;
x.segment(1, N) = rhovecL0;
x.tail(N) = rhovecV0;
- using isochoric = IsochoricDerivatives<Model, Scalar>;
Eigen::Map<Eigen::ArrayXd> rhovecL(&(x(1)), N);
Eigen::Map<Eigen::ArrayXd> rhovecV(&(x(1 + N)), N);
@@ -647,15 +427,15 @@ auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec
for (int iter = 0; iter < flags.maxiter; ++iter) {
- auto RL = model.R(xmolar_spec);
+ auto RL = model.get_R(xmolar_spec);
auto RLT = RL * T;
auto RVT = RLT; // Note: this should not be exactly the same if you use mole-fraction-weighted gas constants
// calculations from the EOS in the isochoric thermodynamics formalism
- auto [PsirL, PsirgradL, hessianL] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecL);
- auto [PsirV, PsirgradV, hessianV] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecV);
- auto DELTAdmu_dT_res = (isochoric::build_d2PsirdTdrhoi_autodiff(model, T, rhovecL.eval())
- - isochoric::build_d2PsirdTdrhoi_autodiff(model, T, rhovecV.eval())).eval();
+ auto [PsirL, PsirgradL, hessianL] = model.build_Psir_fgradHessian_autodiff(T, rhovecL);
+ auto [PsirV, PsirgradV, hessianV] = model.build_Psir_fgradHessian_autodiff(T, rhovecV);
+ auto DELTAdmu_dT_res = (model.build_d2PsirdTdrhoi_autodiff(T, rhovecL.eval())
+ - model.build_d2PsirdTdrhoi_autodiff(T, rhovecV.eval())).eval();
auto make_diag = [](const Eigen::ArrayXd& v) -> Eigen::ArrayXXd {
Eigen::MatrixXd A = Eigen::MatrixXd::Identity(v.size(), v.size());
@@ -690,10 +470,10 @@ auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec
J.block(0, 1, N, N) = HtotL; // These are the concentration derivatives
J.block(0, N+1, N, N) = -HtotV; // These are the concentration derivatives
// Pressure contributions in Jacobian
- J(N, 0) = isochoric::get_dpdT_constrhovec(model, T, rhovecL)/p_spec;
+ J(N, 0) = model.get_dpdT_constrhovec(T, rhovecL)/p_spec;
J.block(N, 1, 1, N) = dpdrhovecL.transpose()/p_spec;
// No vapor concentration derivatives
- J(N+1, 0) = isochoric::get_dpdT_constrhovec(model, T, rhovecV)/p_spec;
+ J(N+1, 0) = model.get_dpdT_constrhovec(T, rhovecV)/p_spec;
// No liquid concentration derivatives
J.block(N+1, N+1, 1, N) = dpdrhovecV.transpose()/p_spec;
// Mole fraction contributions in Jacobian
@@ -741,13 +521,11 @@ auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec
return std::make_tuple(return_code, T, rhovecLfinal, rhovecVfinal);
}
-
-template<class Model, class Scalar, class VecType>
-auto get_drhovecdp_Tsat(const Model& model, const Scalar &T, const VecType& rhovecL, const VecType& rhovecV) {
+inline auto get_drhovecdp_Tsat(const AbstractModel& model, const double &T, const Eigen::ArrayXd& rhovecL, const Eigen::ArrayXd& rhovecV) {
//tic = timeit.default_timer();
- using id = IsochoricDerivatives<Model, Scalar, VecType>;
- Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hliq = id::build_Psi_Hessian_autodiff(model, T, rhovecL).eval();
- Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hvap = id::build_Psi_Hessian_autodiff(model, T, rhovecV).eval();
+ using Scalar = double;
+ Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hliq = model.build_Psi_Hessian_autodiff(T, rhovecL).eval();
+ Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hvap = model.build_Psi_Hessian_autodiff(T, rhovecV).eval();
//Hvap[~np.isfinite(Hvap)] = 1e20;
//Hliq[~np.isfinite(Hliq)] = 1e20;
@@ -768,10 +546,10 @@ auto get_drhovecdp_Tsat(const Model& model, const Scalar &T, const VecType& rhov
}
else{
// Special treatment for infinite dilution
- auto murL = id::build_Psir_gradient_autodiff(model, T, rhovecL);
- auto murV = id::build_Psir_gradient_autodiff(model, T, rhovecV);
- auto RL = model.R(rhovecL / rhovecL.sum());
- auto RV = model.R(rhovecV / rhovecV.sum());
+ auto murL = model.build_Psir_gradient_autodiff(T, rhovecL);
+ auto murV = model.build_Psir_gradient_autodiff(T, rhovecV);
+ auto RL = model.get_R(rhovecL / rhovecL.sum());
+ auto RV = model.get_R(rhovecV / rhovecV.sum());
// First, for the liquid part
for (auto i = 0; i < N; ++i) {
@@ -817,15 +595,13 @@ auto get_drhovecdp_Tsat(const Model& model, const Scalar &T, const VecType& rhov
/**
* Derivative of molar concentration vectors w.r.t. p along an isobar of the phase envelope for binary mixtures
*/
-template<class Model, class Scalar, class VecType>
-auto get_drhovecdT_psat(const Model& model, const Scalar &T, const VecType& rhovecL, const VecType& rhovecV) {
-
- using id = IsochoricDerivatives<Model, Scalar, VecType>;
+inline auto get_drhovecdT_psat(const AbstractModel& model, const double &T, const Eigen::ArrayXd& rhovecL, const Eigen::ArrayXd& rhovecV) {
+ using Scalar = double;
if (rhovecL.size() != 2) { throw std::invalid_argument("Binary mixtures only"); }
assert(rhovecL.size() == rhovecV.size());
- Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hliq = id::build_Psi_Hessian_autodiff(model, T, rhovecL).eval();
- Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hvap = id::build_Psi_Hessian_autodiff(model, T, rhovecV).eval();
+ Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hliq = model.build_Psi_Hessian_autodiff(T, rhovecL).eval();
+ Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hvap = model.build_Psi_Hessian_autodiff(T, rhovecV).eval();
auto N = rhovecL.size();
Eigen::MatrixXd A = decltype(Hliq)::Zero(N, N);
@@ -840,9 +616,9 @@ auto get_drhovecdT_psat(const Model& model, const Scalar &T, const VecType& rhov
A(1, 0) = Hliq.row(0).dot(rhovecL.matrix());
A(1, 1) = Hliq.row(1).dot(rhovecL.matrix());
- auto DELTAdmu_dT = (id::get_dchempotdT_autodiff(model, T, rhovecV) - id::get_dchempotdT_autodiff(model, T, rhovecL)).eval();
- b(0) = DELTAdmu_dT.matrix().dot(rhovecV.matrix()) - id::get_dpdT_constrhovec(model, T, rhovecV);
- b(1) = -id::get_dpdT_constrhovec(model, T, rhovecL);
+ auto DELTAdmu_dT = (model.get_dchempotdT_autodiff(T, rhovecV) - model.get_dchempotdT_autodiff(T, rhovecL)).eval();
+ b(0) = DELTAdmu_dT.matrix().dot(rhovecV.matrix()) - model.get_dpdT_constrhovec(T, rhovecV);
+ b(1) = -model.get_dpdT_constrhovec(T, rhovecL);
// Calculate the derivatives of the liquid phase
drhovecdT_liq = linsolve(A, b);
// Calculate the derivatives of the vapor phase
@@ -850,17 +626,17 @@ auto get_drhovecdT_psat(const Model& model, const Scalar &T, const VecType& rhov
}
else{
// Special treatment for infinite dilution
- auto murL = id::build_Psir_gradient_autodiff(model, T, rhovecL);
- auto murV = id::build_Psir_gradient_autodiff(model, T, rhovecV);
- auto RL = model.R(rhovecL / rhovecL.sum());
- auto RV = model.R(rhovecV / rhovecV.sum());
+ auto murL = model.build_Psir_gradient_autodiff(T, rhovecL);
+ auto murV = model.build_Psir_gradient_autodiff(T, rhovecV);
+ auto RL = model.get_R(rhovecL / rhovecL.sum());
+ auto RV = model.get_R(rhovecV / rhovecV.sum());
// The dot product contains terms of the type:
// rho'_i (R ln(rho"_i /rho'_i) + d mu ^ r"_i/d T - d mu^r'_i/d T)
// Residual contribution to the difference in temperature derivative of chemical potential
// It should be fine to evaluate with zero densities:
- auto DELTAdmu_dT_res = (id::build_d2PsirdTdrhoi_autodiff(model, T, rhovecV) - id::build_d2PsirdTdrhoi_autodiff(model, T, rhovecL)).eval();
+ auto DELTAdmu_dT_res = (model.build_d2PsirdTdrhoi_autodiff(T, rhovecV) - model.build_d2PsirdTdrhoi_autodiff(T, rhovecL)).eval();
// Now the ideal-gas part causes trouble, so multiply by the rhovec, once with liquid, another with vapor
// Start off with the assumption that the rhovec is all positive (fix elements later)
auto DELTAdmu_dT_rhoV_ideal = (rhovecV*(RV*log(rhovecV/rhovecL))).eval();
@@ -874,8 +650,8 @@ auto get_drhovecdT_psat(const Model& model, const Scalar &T, const VecType& rhov
}
double DELTAdmu_dT_rhoV = rhovecV.matrix().dot(DELTAdmu_dT_res.matrix()) + DELTAdmu_dT_rhoV_ideal.sum();
- b(0) = DELTAdmu_dT_rhoV - id::get_dpdT_constrhovec(model, T, rhovecV);
- b(1) = -id::get_dpdT_constrhovec(model, T, rhovecL);
+ b(0) = DELTAdmu_dT_rhoV - model.get_dpdT_constrhovec(T, rhovecV);
+ b(1) = -model.get_dpdT_constrhovec(T, rhovecL);
// First, for the liquid part
for (auto i = 0; i < N; ++i) {
@@ -933,20 +709,18 @@ auto get_drhovecdT_psat(const Model& model, const Scalar &T, const VecType& rhov
* To keep the vapor mole fraction constant, just swap the input molar concentrations to this function, the first concentration
* vector is always the one with fixed mole fractions
*/
-template<class Model, class Scalar, class VecType>
-auto get_drhovecdT_xsat(const Model& model, const Scalar& T, const VecType& rhovecL, const VecType& rhovecV) {
- using id = IsochoricDerivatives<Model, Scalar, VecType>;
-
+inline auto get_drhovecdT_xsat(const AbstractModel& model, const double& T, const Eigen::ArrayXd& rhovecL, const Eigen::ArrayXd& rhovecV) {
+ using Scalar = double;
if (rhovecL.size() != 2) { throw std::invalid_argument("Binary mixtures only"); }
assert(rhovecL.size() == rhovecV.size());
- VecType molefracL = rhovecL / rhovecL.sum();
- VecType deltas = (id::get_dchempotdT_autodiff(model, T, rhovecV) - id::get_dchempotdT_autodiff(model, T, rhovecL)).eval();
- Scalar deltabeta = (id::get_dpdT_constrhovec(model, T, rhovecV)- id::get_dpdT_constrhovec(model, T, rhovecL));
- VecType deltarho = (rhovecV - rhovecL).eval();
+ Eigen::ArrayXd molefracL = rhovecL / rhovecL.sum();
+ Eigen::ArrayXd deltas = (model.get_dchempotdT_autodiff(T, rhovecV) - model.get_dchempotdT_autodiff(T, rhovecL)).eval();
+ Scalar deltabeta = (model.get_dpdT_constrhovec(T, rhovecV)- model.get_dpdT_constrhovec(T, rhovecL));
+ Eigen::ArrayXd deltarho = (rhovecV - rhovecL).eval();
- Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hliq = id::build_Psi_Hessian_autodiff(model, T, rhovecL).eval();
- Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hvap = id::build_Psi_Hessian_autodiff(model, T, rhovecV).eval();
+ Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hliq = model.build_Psi_Hessian_autodiff(T, rhovecL).eval();
+ Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Hvap = model.build_Psi_Hessian_autodiff(T, rhovecV).eval();
Eigen::MatrixXd drhodT_liq, drhodT_vap;
if ((rhovecL != 0).all() && (rhovecV != 0).all()) {
@@ -986,21 +760,20 @@ auto get_drhovecdT_xsat(const Model& model, const Scalar& T, const VecType& rhov
* \left(\frac{dp}{dT}\right)_{x, \sigma} = \left(\frac{\partial p}{\partial T}\right)_{\vec\rho}\frac{dT}{dT} + \sum_k \left(\frac{\partial p}{\partial \rho_k}\right)_{T,\rho_{j\neq k}} \left(\frac{\partial \rho_k}{\partial T}\right)_{x,\sigma}
* \f]
*/
-template<class Model, class Scalar, class VecType>
-auto get_dpsat_dTsat_isopleth(const Model& model, const Scalar& T, const VecType& rhovecL, const VecType& rhovecV) {
+template<typename Model = AbstractModel>
+auto get_dpsat_dTsat_isopleth(const Model& model, const double& T, const Eigen::ArrayXd& rhovecL, const Eigen::ArrayXd& rhovecV) {
// Derivative along phase envelope at constant composition (correct, tested)
auto [drhovecLdT_xsat, drhovecVdT_xsat] = get_drhovecdT_xsat(model, T, rhovecL, rhovecV);
// And the derivative of the total density
auto drhoLdT_sat = drhovecLdT_xsat.sum();
- using tdx = TDXDerivatives<Model, Scalar, VecType>;
double rhoL = rhovecL.sum();
auto molefracL = rhovecL / rhoL;
- auto RT = model.R(molefracL) * T;
- auto derivs = tdx::template get_Ar0n<2>(model, T, rhoL, molefracL);
+ auto RT = model.get_R(molefracL) * T;
+ auto derivs = model.get_Ar02n(T, rhoL, molefracL);
auto dpdrho = RT*(1 + 2 * derivs[1] + derivs[2]);
- Scalar dpdT = model.R(molefracL) * rhoL * (1 + derivs[1] - tdx::get_Ar11(model, T, rhoL, molefracL));
+ double dpdT = model.get_R(molefracL) * rhoL * (1 + derivs[1] - model.get_Ar11(T, rhoL, molefracL));
auto der = dpdT + dpdrho * drhoLdT_sat;
return der;
@@ -1016,8 +789,7 @@ auto get_dpsat_dTsat_isopleth(const Model& model, const Scalar& T, const VecType
* \brief Trace an isotherm with parametric tracing
* \ note If options.revision is 2, the data will be returned in the "data" field, otherwise the data will be returned as root array
*/
-template<typename Model, typename Scalar, typename VecType>
-auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, VecType rhovecV0, const std::optional<TVLEOptions>& options = std::nullopt)
+inline auto trace_VLE_isotherm_binary(const AbstractModel &model, double T, const Eigen::ArrayXd& rhovecL0, const Eigen::ArrayXd& rhovecV0, const std::optional<TVLEOptions>& options = std::nullopt)
{
// Get the options, or the default values if not provided
TVLEOptions opt = options.value_or(TVLEOptions{});
@@ -1117,9 +889,8 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]), N);
auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]) + N, N);
auto rhototL = rhovecL.sum(), rhototV = rhovecV.sum();
- using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
- double pL = rhototL * model.R(rhovecL / rhovecL.sum())*T + id::get_pr(model, T, rhovecL);
- double pV = rhototV * model.R(rhovecV / rhovecV.sum())*T + id::get_pr(model, T, rhovecV);
+ double pL = rhototL * model.get_R(rhovecL / rhovecL.sum())*T + model.get_pr(T, rhovecL);
+ double pV = rhototV * model.get_R(rhovecV / rhovecV.sum())*T + model.get_pr(T, rhovecV);
// Store the derivative
try {
@@ -1144,9 +915,8 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
{"drho/dt", last_drhodt}
};
if (opt.calc_criticality) {
- using ct = CriticalTracing<Model, Scalar, VecType>;
- point["crit. conditions L"] = ct::get_criticality_conditions(model, T, rhovecL);
- point["crit. conditions V"] = ct::get_criticality_conditions(model, T, rhovecV);
+ point["crit. conditions L"] = model.get_criticality_conditions(T, rhovecL);
+ point["crit. conditions V"] = model.get_criticality_conditions(T, rhovecV);
}
JSONdata.push_back(point);
//std::cout << JSONdata.back().dump() << std::endl;
@@ -1198,14 +968,12 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]) + N, N);
auto x = rhovecL / rhovecL.sum();
auto y = rhovecV / rhovecV.sum();
- using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
- double p = rhovecL.sum()*model.R(x)*T + id::get_pr(model, T, rhovecL);
+ double p = rhovecL.sum()*model.get_R(x)*T + model.get_pr(T, rhovecL);
// Check if the solution has gone mechanically unstable
if (opt.calc_criticality) {
- using ct = CriticalTracing<Model, Scalar, VecType>;
- auto condsL = ct::get_criticality_conditions(model, T, rhovecL);
- auto condsV = ct::get_criticality_conditions(model, T, rhovecV);
+ auto condsL = model.get_criticality_conditions(T, rhovecL);
+ auto condsV = model.get_criticality_conditions(T, rhovecV);
if (condsL[0] < opt.crit_termination || condsV[0] < opt.crit_termination){
return true;
}
@@ -1228,7 +996,7 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]), N).eval();
auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[0 + N]), N).eval();
auto x = (Eigen::ArrayXd(2) << rhovecL(0) / rhovecL.sum(), rhovecL(1) / rhovecL.sum()).finished(); // Mole fractions in the liquid phase (to be kept constant)
- auto [return_code, rhovecLnew, rhovecVnew] = mix_VLE_Tx(model, T, rhovecL, rhovecV, x, 1e-10, 1e-8, 1e-10, 1e-8, 10);
+ auto [return_code, rhovecLnew, rhovecVnew] = model.mix_VLE_Tx(T, rhovecL, rhovecV, x, 1e-10, 1e-8, 1e-10, 1e-8, 10);
// If the step is accepted, copy into x again ...
auto rhovecLview = Eigen::Map<Eigen::ArrayXd>(&(x0[0]), N);
@@ -1263,8 +1031,8 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
/***
* \brief Trace an isobar with parametric tracing
*/
-template<typename Model, typename Scalar, typename VecType>
-auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rhovecL0, VecType rhovecV0, const std::optional<PVLEOptions>& options = std::nullopt)
+template<typename Model = AbstractModel>
+auto trace_VLE_isobar_binary(const Model& model, double p, double T0, const Eigen::ArrayXd& rhovecL0, const Eigen::ArrayXd& rhovecV0, const std::optional<PVLEOptions>& options = std::nullopt)
{
// Get the options, or the default values if not provided
PVLEOptions opt = options.value_or(PVLEOptions{});
@@ -1368,9 +1136,8 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(x0[1]), N);
auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[1]) + N, N);
auto rhototL = rhovecL.sum(), rhototV = rhovecV.sum();
- using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
- double pL = rhototL * model.R(rhovecL / rhovecL.sum()) * T + id::get_pr(model, T, rhovecL);
- double pV = rhototV * model.R(rhovecV / rhovecV.sum()) * T + id::get_pr(model, T, rhovecV);
+ double pL = rhototL * model.R(rhovecL / rhovecL.sum()) * T + model.get_pr(T, rhovecL);
+ double pV = rhototV * model.R(rhovecV / rhovecV.sum()) * T + model.get_pr(T, rhovecV);
// Store the derivative
try {
@@ -1395,9 +1162,8 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
{"drho/dt", last_drhodt}
};
if (opt.calc_criticality) {
- using ct = CriticalTracing<Model, Scalar, VecType>;
- point["crit. conditions L"] = ct::get_criticality_conditions(model, T, rhovecL);
- point["crit. conditions V"] = ct::get_criticality_conditions(model, T, rhovecV);
+ point["crit. conditions L"] = model.get_criticality_conditions(T, rhovecL);
+ point["crit. conditions V"] = model.get_criticality_conditions(T, rhovecV);
}
JSONdata.push_back(point);
//std::cout << JSONdata.back().dump() << std::endl;
@@ -1452,9 +1218,8 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
auto y = rhovecV / rhovecV.sum();
// Check if the solution has gone mechanically unstable
if (opt.calc_criticality) {
- using ct = CriticalTracing<Model, Scalar, VecType>;
- auto condsL = ct::get_criticality_conditions(model, T, rhovecL);
- auto condsV = ct::get_criticality_conditions(model, T, rhovecV);
+ auto condsL = model.get_criticality_conditions(T, rhovecL);
+ auto condsV = model.get_criticality_conditions(T, rhovecV);
if (condsL[0] < 1e-12 || condsV[0] < 1e-12) {
return true;
}
@@ -1475,7 +1240,7 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(x0[1]), N).eval();
auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[1 + N]), N).eval();
auto x = (Eigen::ArrayXd(2) << rhovecL(0) / rhovecL.sum(), rhovecL(1) / rhovecL.sum()).finished(); // Mole fractions in the liquid phase (to be kept constant)
- auto [return_code, Tnew, rhovecLnew, rhovecVnew] = mixture_VLE_px(model, p, x, T, rhovecL, rhovecV);
+ auto [return_code, Tnew, rhovecLnew, rhovecVnew] = model.mixture_VLE_px(p, x, T, rhovecL, rhovecV);
// If the step is accepted, copy into x again ...
x0[0] = Tnew;
@@ -1493,4 +1258,33 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
return JSONdata;
}
+#define VLE_FUNCTIONS_TO_WRAP \
+ X(trace_VLE_isobar_binary) \
+ X(trace_VLE_isotherm_binary) \
+ X(get_dpsat_dTsat_isopleth) \
+ X(get_drhovecdT_xsat) \
+ X(get_drhovecdT_psat) \
+ X(get_drhovecdp_Tsat) \
+ X(trace_critical_arclength_binary) \
+ X(mixture_VLE_px) \
+ X(mix_VLE_Tp) \
+ X(mix_VLE_Tx)
+
+#define X(f) template <typename TemplatedModel, typename ...Params, \
+typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, TemplatedModel>::value>::type> \
+inline auto f(const TemplatedModel& model, Params&&... params){ \
+ auto view = teqp::cppinterface::adapter::make_cview(model); \
+ const AbstractModel& am = *view.get(); \
+ return f(am, std::forward<Params>(params)...); \
+}
+ VLE_FUNCTIONS_TO_WRAP
+#undef X
+#undef VLE_FUNCTIONS_TO_WRAP
+
+//
+//template <typename TemplatedModel, typename ...Params, typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, TemplatedModel>::value>::type>
+//auto get_drhovecdT_psat(const TemplatedModel& model, Params&&... params){
+// return get_drhovecdT_psat(teqp::cppinterface::adapter::make_cview(model), std::forward<Params>(params)...);
+//}
+
}; /* namespace teqp*/
diff --git a/include/teqp/algorithms/VLE_pure.hpp b/include/teqp/algorithms/VLE_pure.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dacb80692f5dfb5e295e822e9c8f1b40e174aa9c
--- /dev/null
+++ b/include/teqp/algorithms/VLE_pure.hpp
@@ -0,0 +1,259 @@
+//
+// VLE_pure.hpp
+// teqp
+//
+// Created by Bell, Ian H. (Fed) on 5/3/23.
+//
+#pragma once
+#include <type_traits>
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/derivs.hpp"
+
+namespace teqp{
+
+template<typename Model> using is_AbstractModel = typename std::is_base_of<teqp::cppinterface::AbstractModel, Model>;
+template<typename Model> using is_not_AbstractModel = std::negation<is_AbstractModel<Model>>;
+
+template<typename Model, typename TYPE=double, ADBackends backend = ADBackends::autodiff>
+class IsothermPureVLEResiduals {
+public:
+ using EigenArray = Eigen::Array<TYPE, 2, 1>;
+ using EigenArray1 = Eigen::Array<TYPE, 1, 1>;
+ using EigenMatrix = Eigen::Array<TYPE, 2, 2>;
+private:
+ const Model& m_model;
+ const TYPE m_T;
+ const Eigen::ArrayXd molefracs;
+ EigenMatrix J;
+ EigenArray y;
+public:
+ std::size_t icall = 0;
+ double Rr, R0;
+
+ IsothermPureVLEResiduals(const Model& model, const TYPE& T, const std::optional<Eigen::ArrayXd>& molefracs_ = std::nullopt) : m_model(model), m_T(T),
+ molefracs( (molefracs_) ? molefracs_.value() : Eigen::ArrayXd::Ones(1,1)) {
+ if constexpr(is_not_AbstractModel<Model>::value){
+ Rr = m_model.R(molefracs);
+ R0 = m_model.R(molefracs);
+ }
+ else{
+ Rr = m_model.get_R(molefracs);
+ R0 = m_model.get_R(molefracs);
+ }
+ };
+
+ const auto& get_errors() { return y; };
+
+ template<typename Rho>
+ auto get_der(const Rho& rho){
+ if constexpr(is_not_AbstractModel<Model>::value){
+ using tdx = TDXDerivatives<Model,TYPE,Eigen::ArrayXd>;
+ return tdx::template get_Ar0n<2, backend>(m_model, m_T, rho, molefracs);
+ }
+ else{
+ return m_model.get_Ar02n(m_T, rho, molefracs);
+ }
+ }
+
+ auto call(const EigenArray& rhovec) {
+ assert(rhovec.size() == 2);
+
+ const EigenArray1 rhovecL = rhovec.head(1);
+ const EigenArray1 rhovecV = rhovec.tail(1);
+ const auto rhomolarL = rhovecL.sum(), rhomolarV = rhovecV.sum();
+
+ const TYPE &T = m_T;
+ //const TYPE R = m_model.R(molefracs);
+ double R0_over_Rr = R0 / Rr;
+
+ auto derL = get_der(rhomolarL);
+ auto pRTL = rhomolarL*(R0_over_Rr + derL[1]); // p/(R*T)
+ auto dpRTLdrhoL = R0_over_Rr + 2*derL[1] + derL[2];
+ auto hatmurL = derL[1] + derL[0] + R0_over_Rr*log(rhomolarL);
+ auto dhatmurLdrho = (2*derL[1] + derL[2])/rhomolarL + R0_over_Rr/rhomolarL;
+
+ auto derV = get_der(rhomolarV);
+ auto pRTV = rhomolarV*(R0_over_Rr + derV[1]); // p/(R*T)
+ auto dpRTVdrhoV = R0_over_Rr + 2*derV[1] + derV[2];
+ auto hatmurV = derV[1] + derV[0] + R0_over_Rr *log(rhomolarV);
+ auto dhatmurVdrho = (2*derV[1] + derV[2])/rhomolarV + R0_over_Rr/rhomolarV;
+
+ y(0) = pRTL - pRTV;
+ J(0, 0) = dpRTLdrhoL;
+ J(0, 1) = -dpRTVdrhoV;
+
+ y(1) = hatmurL - hatmurV;
+ J(1, 0) = dhatmurLdrho;
+ J(1, 1) = -dhatmurVdrho;
+
+ icall++;
+ return y;
+ }
+ auto Jacobian(const EigenArray& rhovec){
+ return J;
+ }
+ //auto numJacobian(const EigenArray& rhovec) {
+ // EigenArray plus0 = rhovec, plus1 = rhovec;
+ // double dr = 1e-6 * rhovec[0];
+ // plus0[0] += dr; plus1[1] += dr;
+ // EigenMatrix J;
+ // J.col(0) = (call(plus0) - call(rhovec))/dr;
+ // J.col(1) = (call(plus1) - call(rhovec))/dr;
+ // return J;
+ //}
+};
+
+template<typename Residual, typename Scalar=double>
+auto do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
+ using EArray = Eigen::Array<Scalar, 2, 1>;
+ auto rhovec = (EArray() << rhoL, rhoV).finished();
+ auto r0 = resid.call(rhovec);
+ auto J = resid.Jacobian(rhovec);
+ for (int iter = 0; iter < maxiter; ++iter){
+ if (iter > 0) {
+ r0 = resid.call(rhovec);
+ J = resid.Jacobian(rhovec);
+ }
+ auto v = J.matrix().colPivHouseholderQr().solve(-r0.matrix()).array().eval();
+ auto rhovecnew = (rhovec + v).eval();
+ //double r00 = static_cast<double>(r0[0]);
+ //double r01 = static_cast<double>(r0[1]);
+
+ // If the solution has stopped improving, stop. The change in rhovec is equal to v in infinite precision, but
+ // not when finite precision is involved, use the minimum non-denormal float as the determination of whether
+ // the values are done changing
+ auto minval = std::numeric_limits<Scalar>::epsilon();
+ //double minvaldbl = static_cast<double>(minval);
+ if (((rhovecnew - rhovec).cwiseAbs() < minval).all()) {
+ break;
+ }
+ if ((r0.cwiseAbs() < minval).all()) {
+ break;
+ }
+ rhovec = rhovecnew;
+ }
+ return rhovec;
+}
+
+inline auto pure_VLE_T(const teqp::cppinterface::AbstractModel& model, double T, double rhoL, double rhoV, int maxiter, const std::optional<Eigen::ArrayXd>& molefracs = std::nullopt) {
+ Eigen::ArrayXd molefracs_{Eigen::ArrayXd::Ones(1,1)};
+ if (molefracs){ molefracs_ = molefracs.value(); }
+ auto res = IsothermPureVLEResiduals<teqp::cppinterface::AbstractModel>(model, T, molefracs_);
+ return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
+}
+
+/***
+ * \brief Calculate the derivative of vapor pressure with respect to temperature
+ * \param model The model to operate on
+ * \param T Temperature
+ * \param rhoL Liquid density
+ * \param rhoV Vapor density
+ *
+ * Based upon
+ * \f[
+ * \frac{dp_{\sigma}}{dT} = \frac{h''-h'}{T(v''-v')} = \frac{s''-s'}{v''-v'}
+ * \f]
+ * where the \f$h''-h'\f$ is given by the difference in residual enthalpy \f$h''-h' = h^r''-h^r'\f$ because the ideal-gas parts cancel
+ */
+inline auto dpsatdT_pure(const teqp::cppinterface::AbstractModel& model, double T, double rhoL, double rhoV, const std::optional<Eigen::ArrayXd>& molefracs = std::nullopt) {
+
+ auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+ if (molefracs){ molefrac = molefracs.value(); }
+
+ auto R = model.get_R(molefrac);
+
+ auto hrVLERTV = model.get_Ar01(T, rhoV, molefrac) + model.get_Ar10(T, rhoV, molefrac);
+ auto hrVLERTL = model.get_Ar01(T, rhoL, molefrac) + model.get_Ar10(T, rhoL, molefrac);
+ auto deltahr_over_T = R*(hrVLERTV-hrVLERTL);
+ auto dpsatdT = deltahr_over_T/(1/rhoV-1/rhoL); // From Clapeyron; dp/dT = Deltas/Deltav = Deltah/(T*Deltav); Delta=V-L
+ return dpsatdT;
+}
+
+/***
+ \brief Starting at the critical point, trace the VLE down to a temperature of interest
+
+ \note This method only works for well-behaved EOS, notably absent from that category are EOS in the multiparameter category with orthobaric scaling exponent not equal to 0.5 at the critical point. Most other analytical EOS work fine
+
+ The JSON data structure defines the variables that need to be specified.
+
+ In the current implementation, there are a few steps:
+ 1. Solve for the true critical point satisfying \f$(\partial p/\partial \rho)_{T}=(\partial^2p/\partial\rho^2)_{T}=0\f$
+ 2. Take a small step away from the critical point (this is where the beta=0.5 assumption is invoked)
+ 3. Integrate from the near critical temperature to the temperature of interest
+ */
+inline auto pure_trace_VLE(const AbstractModel& model, const double T, const nlohmann::json &spec){
+ // Start at the true critical point, from the specified guess value
+ nlohmann::json pure_spec;
+ Eigen::ArrayXd z{Eigen::ArrayXd::Ones(1,1)};
+ if (spec.contains("pure_spec")){
+ pure_spec = spec.at("pure_spec");
+ z = Eigen::ArrayXd(pure_spec.at("alternative_length").get<int>()); z.setZero();
+ z(pure_spec.at("alternative_pure_index").get<int>()) = 1;
+ }
+
+ auto [Tc, rhoc] = solve_pure_critical(model, spec.at("Tcguess").get<double>(), spec.at("rhocguess").get<double>(), pure_spec);
+
+ // Small step towards lower temperature close to critical temperature
+ double Tclose = spec.at("Tred").get<double>()*Tc;
+ auto rhoLrhoV = extrapolate_from_critical(model, Tc, rhoc, Tclose, z);
+ auto rhoLrhoVpolished = pure_VLE_T(model, Tclose, rhoLrhoV[0], rhoLrhoV[1], spec.value("NVLE", 10), z);
+ if (rhoLrhoVpolished[0] == rhoLrhoVpolished[1]){
+ throw teqp::IterationError("Converged to trivial solution");
+ }
+
+ // "Integrate" down to temperature of interest
+ int Nstep = spec.at("Nstep");
+ double R = model.R(z);
+ bool with_deriv = spec.at("with_deriv");
+ double dT = -(Tclose-T)/(Nstep-1);
+ using tdx = TDXDerivatives<decltype(model)>;
+ for (auto T_: Eigen::ArrayXd::LinSpaced(Nstep, Tclose, T)){
+ rhoLrhoVpolished = pure_VLE_T(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], spec.value("NVLE", 10), z);
+
+ //auto pL = rhoLrhoVpolished[0]*R*T_*(1 + tdx::get_Ar01(model, T_, rhoLrhoVpolished[0], z));
+ //auto pV = rhoLrhoVpolished[1]*R*T_*(1 + tdx::get_Ar01(model, T_, rhoLrhoVpolished[1], z));
+ //std::cout << pL << " " << pV << " " << pL/pV-1 << std::endl;
+ if (with_deriv){
+ // Get drho/dT for both phases
+ auto dpsatdT = dpsatdT_pure(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], z);
+ auto get_drhodT = [&z, R, dpsatdT](const AbstractModel& model, double T, double rho){
+ auto dpdrho = R*T*(1 + 2*model.get_Ar01(T, rho, z) + model.get_Ar02(T, rho, z));
+ auto dpdT = R*rho*(1 + model.get_Ar01(T, rho, z) - model.get_Ar11(T, rho, z));
+ return -dpdT/dpdrho + dpsatdT/dpdrho;
+ };
+ auto drhodTL = get_drhodT(model, T_, rhoLrhoVpolished[0]);
+ auto drhodTV = get_drhodT(model, T_, rhoLrhoVpolished[1]);
+ // Use the obtained derivative to calculate the step in rho from deltarho = (drhodT)*dT
+ auto DeltarhoL = dT*drhodTL, DeltarhoV = dT*drhodTV;
+ rhoLrhoVpolished[0] += DeltarhoL;
+ rhoLrhoVpolished[1] += DeltarhoV;
+ }
+
+ // Updated values for densities at new T
+ if (!std::isfinite(rhoLrhoVpolished[0])){
+ throw teqp::IterationError("The density is no longer valid; try increasing Nstep");
+ }
+ if (rhoLrhoVpolished[0] == rhoLrhoVpolished[1]){
+ throw teqp::IterationError("Converged to trivial solution; try increasing Nstep");
+ }
+ }
+ return rhoLrhoVpolished;
+}
+
+#define VLE_PURE_FUNCTIONS_TO_WRAP \
+ X(dpsatdT_pure) \
+ X(pure_VLE_T) \
+ X(pure_trace_VLE)
+
+#define X(f) template <typename TemplatedModel, typename ...Params, \
+typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, TemplatedModel>::value>::type> \
+inline auto f(const TemplatedModel& model, Params&&... params){ \
+ auto view = teqp::cppinterface::adapter::make_cview(model); \
+ const AbstractModel& am = *view.get(); \
+ return f(am, std::forward<Params>(params)...); \
+}
+VLE_PURE_FUNCTIONS_TO_WRAP
+#undef X
+#undef VLE_PURE_FUNCTIONS_TO_WRAP
+
+}
diff --git a/include/teqp/algorithms/VLLE.hpp b/include/teqp/algorithms/VLLE.hpp
index 61b4ddce0724c73d61b9cf613dc21514a4e8eedd..e8c18665abcaa41b8f62b1c53fd0d2969d0efe0d 100644
--- a/include/teqp/algorithms/VLLE.hpp
+++ b/include/teqp/algorithms/VLLE.hpp
@@ -3,9 +3,12 @@
#include "teqp/derivs.hpp"
#include "teqp/exceptions.hpp"
#include "teqp/algorithms/VLLE_types.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
namespace teqp {
namespace VLLE {
+
+ using namespace teqp::cppinterface;
/***
* \brief Do a vapor-liquid-liquid phase equilibrium problem for a mixture (binary only for now)
* \param model The model to operate on
@@ -19,8 +22,8 @@ namespace VLLE {
* \param relxtol Relative tolerance on steps in independent variables
* \param maxiter Maximum number of iterations permitted
*/
- template<typename Model, typename Scalar, typename Vector>
- auto mix_VLLE_T(const Model& model, Scalar T, const Vector& rhovecVinit, const Vector& rhovecL1init, const Vector& rhovecL2init, double atol, double reltol, double axtol, double relxtol, int maxiter) {
+
+ auto mix_VLLE_T(const AbstractModel& model, double T, const EArrayd& rhovecVinit, const EArrayd& rhovecL1init, const EArrayd& rhovecL2init, double atol, double reltol, double axtol, double relxtol, int maxiter) {
const Eigen::Index N = rhovecVinit.size();
Eigen::MatrixXd J(3 * N, 3 * N); J.setZero();
@@ -29,8 +32,6 @@ namespace VLLE {
x.head(N) = rhovecVinit;
x.segment(N, N) = rhovecL1init;
x.tail(N) = rhovecL2init;
-
- using isochoric = IsochoricDerivatives<Model, Scalar, Vector>;
Eigen::Map<Eigen::ArrayXd> rhovecV (&(x(0)), N);
Eigen::Map<Eigen::ArrayXd> rhovecL1(&(x(0+N)), N);
@@ -40,23 +41,23 @@ namespace VLLE {
for (int iter = 0; iter < maxiter; ++iter) {
- auto [PsirV, PsirgradV, hessianV] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecV);
- auto [PsirL1, PsirgradL1, hessianL1] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecL1);
- auto [PsirL2, PsirgradL2, hessianL2] = isochoric::build_Psir_fgradHessian_autodiff(model, T, rhovecL2);
+ auto [PsirV, PsirgradV, hessianV] = model.build_Psir_fgradHessian_autodiff(T, rhovecV);
+ auto [PsirL1, PsirgradL1, hessianL1] = model.build_Psir_fgradHessian_autodiff(T, rhovecL1);
+ auto [PsirL2, PsirgradL2, hessianL2] = model.build_Psir_fgradHessian_autodiff(T, rhovecL2);
- auto HtotV = isochoric::build_Psi_Hessian_autodiff(model, T, rhovecV);
- auto HtotL1 = isochoric::build_Psi_Hessian_autodiff(model, T, rhovecL1);
- auto HtotL2 = isochoric::build_Psi_Hessian_autodiff(model, T, rhovecL2);
+ auto HtotV = model.build_Psi_Hessian_autodiff(T, rhovecV);
+ auto HtotL1 = model.build_Psi_Hessian_autodiff(T, rhovecL1);
+ auto HtotL2 = model.build_Psi_Hessian_autodiff(T, rhovecL2);
auto zV = rhovecV/rhovecV.sum(), zL1 = rhovecL1 / rhovecL1.sum(), zL2 = rhovecL2 / rhovecL2.sum();
- double RTL1 = model.R(zL1)*T, RTL2 = model.R(zL2)*T, RTV = model.R(zV)*T;
+ double RTL1 = model.get_R(zL1)*T, RTL2 = model.get_R(zL2)*T, RTV = model.get_R(zV)*T;
auto rhoL1 = rhovecL1.sum();
auto rhoL2 = rhovecL2.sum();
auto rhoV = rhovecV.sum();
- Scalar pL1 = rhoL1 * RTL1 - PsirL1 + (rhovecL1.array() * PsirgradL1.array()).sum(); // The (array*array).sum is a dot product
- Scalar pL2 = rhoL2 * RTL2 - PsirL2 + (rhovecL2.array() * PsirgradL2.array()).sum(); // The (array*array).sum is a dot product
- Scalar pV = rhoV * RTV - PsirV + (rhovecV.array() * PsirgradV.array()).sum();
+ double pL1 = rhoL1 * RTL1 - PsirL1 + (rhovecL1.array() * PsirgradL1.array()).sum(); // The (array*array).sum is a dot product
+ double pL2 = rhoL2 * RTL2 - PsirL2 + (rhovecL2.array() * PsirgradL2.array()).sum(); // The (array*array).sum is a dot product
+ double pV = rhoV * RTV - PsirV + (rhovecV.array() * PsirgradV.array()).sum();
auto dpdrhovecL1 = RTL1 + (hessianL1 * rhovecL1.matrix()).array();
auto dpdrhovecL2 = RTL2 + (hessianL2 * rhovecL2.matrix()).array();
auto dpdrhovecV = RTV + (hessianV * rhovecV.matrix()).array();
@@ -100,7 +101,7 @@ namespace VLLE {
// If the solution has stopped improving, stop. The change in x is equal to dx in infinite precision, but
// not when finite precision is involved, use the minimum non-denormal float as the determination of whether
// the values are done changing
- if (((x.array() - dx.array()).cwiseAbs() < std::numeric_limits<Scalar>::min()).all()) {
+ if (((x.array() - dx.array()).cwiseAbs() < std::numeric_limits<double>::min()).all()) {
return_code = VLLE_return_code::xtol_satisfied;
break;
}
@@ -142,8 +143,8 @@ namespace VLLE {
* \param model The Model to be used for the thermodynamics
* \param traces The nlohmann::json formatted information from the traces, perhaps obtained from trace_VLE_isotherm_binary
*/
- template<typename Model>
- auto find_VLLE_T_binary(const Model& model, const std::vector<nlohmann::json>& traces, const std::optional<VLLEFinderOptions> options = std::nullopt) {
+
+ auto find_VLLE_T_binary(const AbstractModel& model, const std::vector<nlohmann::json>& traces, const std::optional<VLLEFinderOptions> options = std::nullopt) {
std::vector<double> x, y;
auto opt = options.value_or(VLLEFinderOptions{});
diff --git a/include/teqp/algorithms/critical_pure.hpp b/include/teqp/algorithms/critical_pure.hpp
index 04f9055fc1f8ec0712d5842d6bfae4aaf383eeef..869e433f670749d9841db3acb1a21389f027846d 100644
--- a/include/teqp/algorithms/critical_pure.hpp
+++ b/include/teqp/algorithms/critical_pure.hpp
@@ -6,20 +6,23 @@
#include "teqp/derivs.hpp"
#include "teqp/algorithms/rootfinding.hpp"
#include "teqp/exceptions.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
#include <optional>
+using namespace teqp::cppinterface;
+
namespace teqp {
+
/**
* Calculate the criticality conditions for a pure fluid and its Jacobian w.r.t. the temperature and density
* for additional fine tuning with multi-variate rootfinding
*
*/
- template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
- auto get_pure_critical_conditions_Jacobian(const Model& model, const Scalar T, const Scalar rho,
+ inline auto get_pure_critical_conditions_Jacobian(const AbstractModel& model, const double T, const double rho,
const std::optional<std::size_t>& alternative_pure_index = std::nullopt, const std::optional<std::size_t>& alternative_length = std::nullopt) {
- using tdx = TDXDerivatives<Model>;
Eigen::ArrayXd z;
if (!alternative_pure_index) {
z = (Eigen::ArrayXd(1) << 1.0).finished();
@@ -34,9 +37,9 @@ namespace teqp {
throw teqp::InvalidArgument("The provided alternative index of " + std::to_string(index) + " is out of range");
}
}
- auto R = model.R(z);
+ auto R = model.get_R(z);
- auto ders = tdx::template get_Ar0n<4, backend>(model, T, rho, z);
+ auto ders = model.get_Ar04n(T, rho, z);
auto dpdrho = R * T * (1 + 2 * ders[1] + ders[2]); // Should be zero at critical point
auto d2pdrho2 = R * T / rho * (2 * ders[1] + 4 * ders[2] + ders[3]); // Should be zero at critical point
@@ -58,9 +61,9 @@ namespace teqp {
*/
// Note: these derivatives are expressed in terms of 1/T and rho as independent variables
- auto Ar11 = tdx::template get_Arxy<1, 1, backend>(model, T, rho, z);
- auto Ar12 = tdx::template get_Arxy<1, 2, backend>(model, T, rho, z);
- auto Ar13 = tdx::template get_Arxy<1, 3, backend>(model, T, rho, z);
+ auto Ar11 = model.get_Ar11(T, rho, z);
+ auto Ar12 = model.get_Ar12(T, rho, z);
+ auto Ar13 = model.get_Ar13(T, rho, z);
auto d3pdrho3 = R * T / (rho * rho) * (6 * ders[2] + 6 * ders[3] + ders[4]);
auto d_dpdrho_dT = R * (-(Ar12 + 2 * Ar11) + ders[2] + 2 * ders[1] + 1);
@@ -76,7 +79,17 @@ namespace teqp {
return std::make_tuple(resids, J);
}
- template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
+ template <typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff,
+ typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, Model>::value>::type>
+ auto get_pure_critical_conditions_Jacobian(const Model& model, const Scalar T, const Scalar rho,
+ const std::optional<std::size_t>& alternative_pure_index = std::nullopt, const std::optional<std::size_t>& alternative_length = std::nullopt) {
+ using namespace teqp::cppinterface::adapter;
+ auto view_ = std::unique_ptr<AbstractModel>(view(model));
+// static_assert(std::is_base_of<teqp::cppinterface::AbstractModel, std::decay_t<decltype(*view_)> >::value);
+ return get_pure_critical_conditions_Jacobian(*(view_), T, rho, alternative_pure_index, alternative_length);
+ }
+
+ template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff, typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, Model>::value>::type>
auto solve_pure_critical(const Model& model, const Scalar T0, const Scalar rho0, const std::optional<nlohmann::json>& flags = std::nullopt) {
auto x = (Eigen::ArrayXd(2) << T0, rho0).finished();
int maxsteps = 10;
@@ -103,14 +116,47 @@ namespace teqp {
return a.matrix().colPivHouseholderQr().solve(b.matrix()).array().eval();
};
for (auto counter = 0; counter < maxsteps; ++counter) {
- auto [resids, Jacobian] = get_pure_critical_conditions_Jacobian<Model, Scalar, backend>(model, x[0], x[1], alternative_pure_index, alternative_length);
+ auto [resids, Jacobian] = get_pure_critical_conditions_Jacobian(model, x[0], x[1], alternative_pure_index, alternative_length);
+ auto v = linsolve(Jacobian, -resids);
+ x += v;
+ }
+ return std::make_tuple(x[0], x[1]);
+ }
+
+ inline auto solve_pure_critical(const AbstractModel& model, const double T0, const double rho0, const std::optional<nlohmann::json>& flags = std::nullopt) {
+ auto x = (Eigen::ArrayXd(2) << T0, rho0).finished();
+ int maxsteps = 10;
+ std::optional<std::size_t> alternative_pure_index;
+ std::optional<std::size_t> alternative_length;
+ if (flags){
+ if (flags.value().contains("maxsteps")){
+ maxsteps = flags.value().at("maxsteps");
+ }
+ if (flags.value().contains("alternative_pure_index")){
+ auto i = flags.value().at("alternative_pure_index").get<int>();
+ if (i < 0){ throw teqp::InvalidArgument("alternative_pure_index cannot be less than 0"); }
+ alternative_pure_index = i;
+ }
+ if (flags.value().contains("alternative_length")){
+ auto i = flags.value().at("alternative_length").get<int>();
+ if (i < 2){ throw teqp::InvalidArgument("alternative_length cannot be less than 2"); }
+ alternative_length = i;
+ }
+ }
+ // A convenience method to make linear system solving more concise with Eigen datatypes
+ // All arguments are converted to matrices, the solve is done, and an array is returned
+ auto linsolve = [](const auto& a, const auto& b) {
+ return a.matrix().colPivHouseholderQr().solve(b.matrix()).array().eval();
+ };
+ for (auto counter = 0; counter < maxsteps; ++counter) {
+ auto [resids, Jacobian] = get_pure_critical_conditions_Jacobian(model, x[0], x[1], alternative_pure_index, alternative_length);
auto v = linsolve(Jacobian, -resids);
x += v;
}
return std::make_tuple(x[0], x[1]);
}
- template<typename Model, typename Scalar>
+ template <typename Model, typename Scalar, typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, Model>::value>::type>
Scalar get_Brho_critical_extrap(const Model& model, const Scalar& Tc, const Scalar& rhoc, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
using tdx = TDXDerivatives<Model, Scalar>;
@@ -130,8 +176,41 @@ namespace teqp {
return Brho;
}
- template<typename Model, typename Scalar>
-Eigen::Array<double, 2, 1> extrapolate_from_critical(const Model& model, const Scalar& Tc, const Scalar& rhoc, const Scalar& T, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
+
+ template <typename Model, typename Scalar, typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, Model>::value>::type>
+ Eigen::Array<double, 2, 1> extrapolate_from_critical(const Model& model, const Scalar& Tc, const Scalar& rhoc, const Scalar& T, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
+ auto Brho = get_Brho_critical_extrap(model, Tc, rhoc, z);
+
+ auto drhohat_dT = Brho / Tc;
+ auto dT = T - Tc;
+
+ auto drhohat = dT * drhohat_dT;
+ auto rholiq = -drhohat / sqrt(1 - T / Tc) + rhoc;
+ auto rhovap = drhohat / sqrt(1 - T / Tc) + rhoc;
+ return (Eigen::ArrayXd(2) << rholiq, rhovap).finished();
+ }
+
+
+
+ inline double get_Brho_critical_extrap(const AbstractModel& model, const double& Tc, const double& rhoc, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
+
+ auto z_ = (Eigen::ArrayXd(1) << 1.0).finished();
+ if (z){
+ z_ = z.value();
+ }
+ auto R = model.get_R(z_);
+ auto ders = model.get_Ar04n(Tc, rhoc, z_);
+ //auto dpdrho = R*Tc*(1 + 2 * ders[1] + ders[2]); // Should be zero
+ //auto d2pdrho2 = R*Tc/rhoc*(2 * ders[1] + 4 * ders[2] + ders[3]); // Should be zero
+ auto d3pdrho3 = R * Tc / (rhoc * rhoc) * (6 * ders[2] + 6 * ders[3] + ders[4]);
+ auto Ar11 = model.get_Ar11(Tc, rhoc, z_);
+ auto Ar12 = model.get_Ar12(Tc, rhoc, z_);
+ auto d2pdrhodT = R * (1 + 2 * ders[1] + ders[2] - 2 * Ar11 - Ar12);
+ auto Brho = sqrt(6 * d2pdrhodT * Tc / d3pdrho3);
+ return Brho;
+ }
+
+ inline auto extrapolate_from_critical(const AbstractModel& model, const double& Tc, const double& rhoc, const double& T, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
auto Brho = get_Brho_critical_extrap(model, Tc, rhoc, z);
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 3b73fbc3bb302923f49812341da42255ee59235e..731117f25fc4063168e7f2670c97ea399ff25215 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -36,19 +36,18 @@ struct CriticalTracing {
return std::make_tuple(es.eigenvalues(), es.eigenvectors());
}
- static auto eigen_problem(const Model& model, const Scalar T, const VecType& rhovec, const std::optional<VecType>& alignment_v0 = std::nullopt) {
+ static auto eigen_problem(const AbstractModel& model, const Scalar T, const VecType& rhovec, const std::optional<VecType>& alignment_v0 = std::nullopt) {
EigenData ed;
auto N = rhovec.size();
Eigen::ArrayX<bool> mask = (rhovec != 0).eval();
- using id = IsochoricDerivatives<decltype(model)>;
-
// Build the Hessian for the residual part;
#if defined(USE_AUTODIFF)
- auto H = id::build_Psir_Hessian_autodiff(model, T, rhovec);
+ auto H = model.build_Psir_Hessian_autodiff(T, rhovec);
#else
+ using id = IsochoricDerivatives<decltype(model)>;
auto H = id::build_Psir_Hessian_mcx(model, T, rhovec);
#endif
// ... and add ideal-gas terms to H
@@ -123,11 +122,11 @@ struct CriticalTracing {
EigenData ei;
};
- static auto get_minimum_eigenvalue_Psi_Hessian(const Model& model, const Scalar T, const VecType& rhovec) {
+ static auto get_minimum_eigenvalue_Psi_Hessian(const AbstractModel& model, const Scalar T, const VecType& rhovec) {
return eigen_problem(model, T, rhovec).eigenvalues[0];
}
- static auto get_derivs(const Model& model, const Scalar T, const VecType& rhovec, const std::optional<VecType>& alignment_v0 = std::nullopt) {
+ static auto get_derivs(const AbstractModel& model, const double T, const VecType& rhovec, const std::optional<VecType>& alignment_v0 = std::nullopt) {
auto molefrac = rhovec / rhovec.sum();
auto R = model.R(molefrac);
@@ -147,19 +146,20 @@ struct CriticalTracing {
}
#if defined(USE_AUTODIFF)
- // Calculate the first through fourth derivative of Psi^r w.r.t. sigma_1
- ArrayXdual4th v0(ei.v0.size()); for (auto i = 0; i < ei.v0.size(); ++i) { v0[i] = ei.v0[i]; }
- ArrayXdual4th rhovecad(rhovec.size()); for (auto i = 0; i < rhovec.size(); ++i) { rhovecad[i] = rhovec[i]; }
- dual4th varsigma{ 0.0 };
- auto wrapper = [&rhovecad, &v0, &T, &model](const auto& sigma_1) {
- auto rhovecused = (rhovecad + sigma_1 * v0).eval();
- auto rhotot = rhovecused.sum();
- auto molefrac = (rhovecused / rhotot).eval();
- return eval(model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot);
- };
- auto psir_derivs_ = derivatives(wrapper, wrt(varsigma), at(varsigma));
- Eigen::ArrayXd psir_derivs; psir_derivs.resize(5);
- for (auto i = 0; i < 5; ++i) { psir_derivs[i] = psir_derivs_[i]; }
+ auto psir_derivs = model.get_Psir_sigma_derivs(T, rhovec, ei.v0);
+// // Calculate the first through fourth derivative of Psi^r w.r.t. sigma_1
+// ArrayXdual4th v0(ei.v0.size()); for (auto i = 0; i < ei.v0.size(); ++i) { v0[i] = ei.v0[i]; }
+// ArrayXdual4th rhovecad(rhovec.size()); for (auto i = 0; i < rhovec.size(); ++i) { rhovecad[i] = rhovec[i]; }
+// dual4th varsigma{ 0.0 };
+// auto wrapper = [&rhovecad, &v0, &T, &model](const auto& sigma_1) {
+// auto rhovecused = (rhovecad + sigma_1 * v0).eval();
+// auto rhotot = rhovecused.sum();
+// auto molefrac = (rhovecused / rhotot).eval();
+// return eval(model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot);
+// };
+// auto psir_derivs_ = derivatives(wrapper, wrt(varsigma), at(varsigma));
+// Eigen::ArrayXd psir_derivs; psir_derivs.resize(5);
+// for (auto i = 0; i < 5; ++i) { psir_derivs[i] = psir_derivs_[i]; }
#else
using namespace mcx;
@@ -202,11 +202,11 @@ struct CriticalTracing {
return std::any_of(std::begin(foo), std::end(foo), [](const auto x) { return x; });
}
- static auto get_drhovec_dT_crit(const Model& model, const Scalar T, const VecType& rhovec) {
+ static auto get_drhovec_dT_crit(const AbstractModel& model, const Scalar& T, const VecType& rhovec) {
// The derivatives of total Psi w.r.t.sigma_1 (numerical for residual, analytic for ideal)
// Returns a tuple, with residual, ideal, total dicts with of number of derivatives, value of derivative
- auto all_derivs = get_derivs(model, T, rhovec, Eigen::ArrayXd());
+ auto all_derivs = get_derivs(model, T, rhovec, std::nullopt);
auto derivs = all_derivs.tot;
// The temperature derivative of total Psi w.r.t.T from a centered finite difference in T
@@ -278,7 +278,7 @@ struct CriticalTracing {
return drhovec_dT;
}
- static auto get_criticality_conditions(const Model& model, const Scalar T, const VecType& rhovec) {
+ static auto get_criticality_conditions(const AbstractModel& model, const Scalar T, const VecType& rhovec) {
auto derivs = get_derivs(model, T, rhovec, Eigen::ArrayXd());
return (Eigen::ArrayXd(2) << derivs.tot[2], derivs.tot[3]).finished();
}
@@ -288,7 +288,7 @@ struct CriticalTracing {
*
* There is a typo in the paper: the last step should be rho + 2*drho*dir
*/
- static auto EigenVectorPath(const Model& model, const Scalar T, const VecType& rhovec, const VecType &u1, Scalar drho) {
+ static auto EigenVectorPath(const AbstractModel& model, const Scalar T, const VecType& rhovec, const VecType &u1, Scalar drho) {
VecType dir = u1/6.0;
VecType rhoshift = rhovec + drho*u1;
auto e1 = eigen_problem(model, T, rhoshift, u1);
@@ -313,7 +313,7 @@ struct CriticalTracing {
/**
* \brief The LocalStability function of Algorithm 2 from Deiters and Bell: https://dx.doi.org/10.1021/acs.iecr.0c03667
*/
- static auto is_locally_stable(const Model& model, const Scalar T, const VecType& rhovec, const Scalar stability_rel_drho) {
+ static auto is_locally_stable(const AbstractModel& model, const Scalar T, const VecType& rhovec, const Scalar stability_rel_drho) {
// At the critical point
auto ezero = eigen_problem(model, T, rhovec);
Scalar lambda1 = ezero.eigenvalues(0);
@@ -331,7 +331,7 @@ struct CriticalTracing {
}
/// Polish a critical point while keeping the overall composition constant and iterating for temperature and overall density
- static auto critical_polish_fixedmolefrac(const Model& model, const Scalar T, const VecType& rhovec, const Scalar z0) {
+ static auto critical_polish_fixedmolefrac(const AbstractModel& model, const Scalar T, const VecType& rhovec, const Scalar z0) {
auto polish_x_resid = [&model, &z0](const auto& x) {
auto T = x[0];
Eigen::ArrayXd rhovec(2); rhovec << z0*x[1], (1-z0)*x[1];
@@ -351,7 +351,7 @@ struct CriticalTracing {
Eigen::ArrayXd rhovecsoln(2); rhovecsoln << x[1]*z0, x[1] * (1 - z0);
return std::make_tuple(x[0], rhovecsoln);
}
- static auto critical_polish_fixedrho(const Model& model, const Scalar T, const VecType& rhovec, const int i) {
+ static auto critical_polish_fixedrho(const AbstractModel& model, const Scalar T, const VecType& rhovec, const int i) {
Scalar rhoval = rhovec[i];
auto polish_x_resid = [&model, &i, &rhoval](const auto& x) {
auto T = x[0];
@@ -372,7 +372,7 @@ struct CriticalTracing {
Eigen::ArrayXd rho = x.tail(x.size() - 1);
return std::make_tuple(x[0], rho);
}
- static auto critical_polish_fixedT(const Model& model, const Scalar T, const VecType& rhovec) {
+ static auto critical_polish_fixedT(const AbstractModel& model, const Scalar T, const VecType& rhovec) {
auto polish_T_resid = [&model, &T](const auto& x) {
auto derivs = get_derivs(model, T, x);
return (Eigen::ArrayXd(2) << derivs.tot[2], derivs.tot[3]).finished();
@@ -388,7 +388,7 @@ struct CriticalTracing {
return x;
}
- static auto trace_critical_arclength_binary(const Model& model, const Scalar& T0, const VecType& rhovec0, const std::optional<std::string>& filename_ = std::nullopt, const std::optional<TCABOptions> &options_ = std::nullopt) -> nlohmann::json {
+ static auto trace_critical_arclength_binary(const AbstractModel& model, const Scalar& T0, const VecType& rhovec0, const std::optional<std::string>& filename_ = std::nullopt, const std::optional<TCABOptions> &options_ = std::nullopt) -> nlohmann::json {
std::string filename = filename_.value_or("");
TCABOptions options = options_.value_or(TCABOptions{});
@@ -474,10 +474,9 @@ struct CriticalTracing {
// Calculate some other parameters, for debugging, or scientific interest
auto rhotot = rhovec.sum();
- using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
- double p = rhotot * model.R(rhovec / rhovec.sum()) * T + id::get_pr(model, T, rhovec);
+ double p = rhotot * model.R(rhovec / rhovec.sum()) * T + model.get_pr(T, rhovec);
auto conditions = get_criticality_conditions(model, T, rhovec);
- double splus = id::get_splus(model, T, rhovec);
+ double splus = model.get_splus(T, rhovec);
auto dxdt = x0;
xprime(x0, dxdt, -1.0);
@@ -507,9 +506,8 @@ struct CriticalTracing {
std::stringstream out;
auto rhotot = rhovec.sum();
double z0 = rhovec[0] / rhotot;
- using id = IsochoricDerivatives<decltype(model)>;
auto conditions = get_criticality_conditions(model, T, rhovec);
- out << z0 << "," << rhovec[0] << "," << rhovec[1] << "," << T << "," << rhotot * model.R(rhovec / rhovec.sum()) * T + id::get_pr(model, T, rhovec) << "," << c << "," << dt << "," << conditions(0) << "," << conditions(1) << std::endl;
+ out << z0 << "," << rhovec[0] << "," << rhovec[1] << "," << T << "," << rhotot * model.R(rhovec / rhovec.sum()) * T + model.get_pr(T, rhovec) << "," << c << "," << dt << "," << conditions(0) << "," << conditions(1) << std::endl;
std::string sout(out.str());
std::cout << sout;
if (ofs.is_open()) {
@@ -595,22 +593,22 @@ struct CriticalTracing {
}
// The polishers to be tried, in order, to polish the critical point
- using PolisherType = std::function<std::tuple<double, VecType>(const Model, const Scalar, const VecType&)>;
+ using PolisherType = std::function<std::tuple<double, VecType>(const AbstractModel&, const Scalar, const VecType&)>;
std::vector<PolisherType> polishers = {
- [&](const Model& model, const Scalar T, const VecType& rhovec) {
+ [&](const AbstractModel& model, const Scalar T, const VecType& rhovec) {
auto [Tnew, rhovecnew] = critical_polish_fixedmolefrac(model, T, rhovec, z0);
return std::make_tuple(Tnew, rhovecnew);
},
- [&](const Model& model, const Scalar T, const VecType& rhovec) {
+ [&](const AbstractModel& model, const Scalar T, const VecType& rhovec) {
auto rhovecnew = critical_polish_fixedT(model, T, rhovec);
return std::make_tuple(T, rhovecnew);
},
- [&](const Model& model, const Scalar T, const VecType& rhovec) {
+ [&](const AbstractModel& model, const Scalar T, const VecType& rhovec) {
int i = 0;
auto [Tnew, rhovecnew] = critical_polish_fixedrho(model, T, rhovec, i);
return std::make_tuple(Tnew, rhovecnew);
},
- [&](const Model& model, const Scalar T, const VecType& rhovec) {
+ [&](const AbstractModel& model, const Scalar T, const VecType& rhovec) {
int i = 1;
auto [Tnew, rhovecnew] = critical_polish_fixedrho(model, T, rhovec, i);
return std::make_tuple(Tnew, rhovecnew);
@@ -746,13 +744,39 @@ struct CriticalTracing {
* \f]
* where the derivatives on the RHS without the subscript $\CRL$ are homogeneous derivatives taken at the mixture statepoint, and are NOT along the critical curve.
*/
- static auto get_dp_dT_crit(const Model& model, const Scalar& T, const VecType& rhovec) {
- using id = IsochoricDerivatives<Model, Scalar, VecType>;
- auto dpdrhovec = id::get_dpdrhovec_constT(model, T, rhovec);
- Scalar v = id::get_dpdT_constrhovec(model, T, rhovec) + (dpdrhovec.array()*get_drhovec_dT_crit(model, T, rhovec).array()).sum();
+ static auto get_dp_dT_crit(const AbstractModel& model, const Scalar& T, const VecType& rhovec) {
+ auto dpdrhovec = model.get_dpdrhovec_constT(T, rhovec);
+ Scalar v = model.get_dpdT_constrhovec(T, rhovec) + (dpdrhovec.array()*get_drhovec_dT_crit(model, T, rhovec).array()).sum();
return v;
}
-
-}; // namespace VecType
+
+#define CRIT_FUNCTIONS_TO_WRAP \
+ X(get_dp_dT_crit) \
+ X(trace_critical_arclength_binary) \
+ X(critical_polish_fixedmolefrac) \
+ X(get_drhovec_dT_crit) \
+ X(get_derivs) \
+ X(eigen_problem)
+
+#define X(f) template <typename TemplatedModel, typename ...Params, \
+typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, TemplatedModel>::value>::type> \
+static auto f(const TemplatedModel& model, Params&&... params){ \
+ auto view = teqp::cppinterface::adapter::make_cview(model); \
+ const AbstractModel& am = *view.get(); \
+ return f(am, std::forward<Params>(params)...); \
+}
+ CRIT_FUNCTIONS_TO_WRAP
+#undef X
+#undef CRIT_FUNCTIONS_TO_WRAP
+
+//template <typename TemplatedModel, typename ...Params,
+//typename = typename std::enable_if<not std::is_base_of<teqp::cppinterface::AbstractModel, TemplatedModel>::value>::type>
+//static auto critical_polish_fixedmolefrac(const TemplatedModel& model, Params&&... params){
+// auto view = teqp::cppinterface::adapter::make_cview(model);
+// const AbstractModel& am = *view.get();
+// return critical_polish_fixedmolefrac(am, std::forward<Params>(params)...);
+//}
+
+}; // namespace Critical
}; // namespace teqp
diff --git a/include/teqp/algorithms/rootfinding.hpp b/include/teqp/algorithms/rootfinding.hpp
index 4f2e7921a5bf68bdd4bcc060ac6a0a7b83dd3044..e9e54f02979be3f1dfc25289d3ef5a70d2a8d4cf 100644
--- a/include/teqp/algorithms/rootfinding.hpp
+++ b/include/teqp/algorithms/rootfinding.hpp
@@ -18,6 +18,9 @@ auto NewtonRaphson(Callable f, const Inputs& args, double tol) {
Eigen::ArrayXd v = J.colPivHouseholderQr().solve(-r0.matrix());
x += v;
auto err = r0.matrix().norm();
+ if (!std::isfinite(err)){
+ throw std::invalid_argument("err is now NaN");
+ }
if (err < tol) {
break;
}
diff --git a/include/teqp/cpp/deriv_adapter.hpp b/include/teqp/cpp/deriv_adapter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..11a2272ea6bb255731a7a1b6ac20691ab6aeaf25
--- /dev/null
+++ b/include/teqp/cpp/deriv_adapter.hpp
@@ -0,0 +1,183 @@
+#pragma once
+
+#include "teqp/derivs.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/exceptions.hpp"
+
+namespace teqp{
+namespace cppinterface{
+namespace adapter{
+
+// The ownership
+template<typename ModelType>
+struct Owner{
+private:
+ ModelType model;
+public:
+ auto& get_ref(){ return model; };
+ const auto& get_cref() const { return model; };
+ const std::type_index index;
+ Owner(ModelType&& m) : model(m), index(std::type_index(typeid(ModelType))) {};
+};
+
+template<typename ModelType>
+struct ConstViewer{
+private:
+ const ModelType& model;
+public:
+ auto& get_ref(){ return model; };
+ const auto& get_cref() const { return model; };
+ const std::type_index index;
+ ConstViewer(ModelType& m) : model(m), index(std::type_index(typeid(ModelType))) {};
+};
+
+namespace internal{
+ template<class T>struct tag{using type=T;};
+}
+
+/**
+ This class holds a const reference to a class, and exposes an interface that matches that used in AbstractModel
+
+ The exposed methods cover all the derivative methods that are obtained by derivatives of the model
+ */
+template<typename ModelPack>
+class DerivativeAdapter : public teqp::cppinterface::AbstractModel{
+private:
+ ModelPack mp;
+public:
+ auto& get_ModelPack_ref(){ return mp; }
+ const auto& get_ModelPack_cref() const { return mp; }
+
+ template<typename T>
+ DerivativeAdapter(internal::tag<T> tag_, const T&& mp): mp(mp) {} ;
+
+ const std::type_index& get_type_index() const override {
+ return mp.index;
+ };
+
+// template<typename T>
+// DerivativeAdapter(const Owner<T>&& mp): mp(mp) {} ;
+//
+// template<typename T>
+// DerivativeAdapter(const ConstViewer<T>&& mp): mp(mp) {} ;
+
+ const AllowedModels& get_model() const override { throw teqp::NotImplementedError(""); };
+ AllowedModels& get_mutable_model() override { throw teqp::NotImplementedError(""); };
+
+ virtual double get_R(const EArrayd& molefrac) const override {
+ return mp.get_cref().R(molefrac);
+ };
+
+ virtual double get_Arxy(const int NT, const int ND, const double T, const double rhomolar, const EArrayd& molefrac) const override{
+ return TDXDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_Ar(NT, ND, mp.get_cref(), T, rhomolar, molefrac);
+ };
+
+ // Here X-Macros are used to create functions like get_Ar00, get_Ar01, ....
+#define X(i,j) virtual double get_Ar ## i ## j(const double T, const double rho, const REArrayd& molefrac) const override { return TDXDerivatives<decltype(mp.get_cref()), double, EArrayd>::template get_Arxy<i,j>(mp.get_cref(), T, rho, molefrac); };
+ ARXY_args
+#undef X
+ // And like get_Ar01n, get_Ar02n, ....
+#define X(i) virtual EArrayd get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefrac) const override { auto vals = TDXDerivatives<decltype(mp.get_cref()), double, EArrayd>::template get_Ar0n<i>(mp.get_cref(), T, rho, molefrac); return Eigen::Map<Eigen::ArrayXd>(&(vals[0]), vals.size()); };
+ AR0N_args
+#undef X
+
+ // Virial derivatives
+ virtual double get_B2vir(const double T, const EArrayd& z) const override {
+ return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_B2vir(mp.get_cref(), T, z);
+ };
+ virtual std::map<int, double> get_Bnvir(const int Nderiv, const double T, const EArrayd& z) const override {
+ return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_Bnvir_runtime(Nderiv, mp.get_cref(), T, z);
+ };
+ virtual double get_B12vir(const double T, const EArrayd& z) const override {
+ return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_B12vir(mp.get_cref(), T, z);
+ };
+ virtual double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const override {
+ return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_dmBnvirdTm_runtime(Nderiv, NTderiv, mp.get_cref(), T, molefrac);
+ };
+
+ // Derivatives from isochoric thermodynamics (all have the same signature within each block), and they differ by their output argument
+#define X(f) virtual double f(const double T, const EArrayd& rhovec) const override { return IsochoricDerivatives<decltype(mp.get_cref()), double, EArrayd>::f(mp.get_cref(), T, rhovec); };
+ ISOCHORIC_double_args
+#undef X
+#define X(f) virtual EArrayd f(const double T, const EArrayd& rhovec) const override { return IsochoricDerivatives<decltype(mp.get_cref()), double, EArrayd>::f(mp.get_cref(), T, rhovec); };
+ ISOCHORIC_array_args
+#undef X
+#define X(f) virtual EMatrixd f(const double T, const EArrayd& rhovec) const override { return IsochoricDerivatives<decltype(mp.get_cref()), double, EArrayd>::f(mp.get_cref(), T, rhovec); };
+ ISOCHORIC_matrix_args
+#undef X
+#define X(f) virtual std::tuple<double, Eigen::ArrayXd, Eigen::MatrixXd> f(const double T, const EArrayd& rhovec) const override { return IsochoricDerivatives<decltype(mp.get_cref()), double, EArrayd>::f(mp.get_cref(), T, rhovec); };
+ ISOCHORIC_multimatrix_args
+#undef X
+ virtual Eigen::ArrayXd get_Psir_sigma_derivs(const double T, const EArrayd& rhovec, const EArrayd& v) const override{
+ return IsochoricDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_Psir_sigma_derivs(mp.get_cref(), T, rhovec, v);
+ };
+
+ virtual EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z ) const override { throw teqp::NotImplementedError("Not available"); };
+};
+
+template<typename TemplatedModel> auto view(const TemplatedModel& tp){
+ ConstViewer cv{tp};
+ return new DerivativeAdapter<decltype(cv)>(internal::tag<decltype(cv)>{}, std::move(cv));
+}
+template<typename TemplatedModel> auto own(const TemplatedModel&& tp){
+ Owner o(std::move(tp));
+ return new DerivativeAdapter<decltype(o)>(internal::tag<decltype(o)>{}, std::move(o));
+}
+
+template<typename TemplatedModel> auto make_owned(const TemplatedModel& tmodel){
+ using namespace teqp::cppinterface;
+ return std::unique_ptr<AbstractModel>(own(std::move(tmodel)));
+};
+
+template<typename TemplatedModel> auto make_cview(const TemplatedModel& tmodel){
+ using namespace teqp::cppinterface;
+ return std::unique_ptr<AbstractModel>(view(tmodel));
+};
+
+/**
+ Get a const reference to the model
+
+ Available for both ownership and const viewer holder types
+ */
+template<typename ModelType>
+const ModelType& get_model_cref(const AbstractModel *am)
+{
+ if (am == nullptr){
+ throw teqp::InvalidArgument("Argument to get_model_cref is a nullptr");
+ }
+ const auto* mptr = dynamic_cast<const DerivativeAdapter<ConstViewer<const ModelType>>*>(am);
+ const auto* mptr2 = dynamic_cast<const DerivativeAdapter<Owner<const ModelType>>*>(am);
+ if (mptr != nullptr){
+ return mptr->get_ModelPack_cref().get_cref();
+ }
+ else if (mptr2 != nullptr){
+ return mptr2->get_ModelPack_cref().get_cref();
+ }
+ else{
+ throw teqp::InvalidArgument("Unable to cast model to desired type");
+ }
+}
+
+/**
+ Get a mutable reference to the model
+
+ Only available when the holder type is ownership (not available for const viewer holder type)
+ */
+template<typename ModelType>
+ModelType& get_model_ref(AbstractModel *am)
+{
+ if (am == nullptr){
+ throw teqp::InvalidArgument("Argument to get_model_ref is a nullptr");
+ }
+ auto* mptr2 = dynamic_cast<DerivativeAdapter<Owner<ModelType>>*>(am);
+ if (mptr2 != nullptr){
+ return mptr2->get_ModelPack_ref().get_ref();
+ }
+ else{
+ throw teqp::InvalidArgument("Unable to cast model to desired type; only the Owner ownership model is allowed");
+ }
+}
+
+}
+}
+}
diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index 7ce5bee4476ad924dbcc4d83ede38cf0d4b4d5bb..6413f989993a9f88dcb65fdb90ab05345d23b737 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -1,5 +1,6 @@
#pragma once
#include <memory>
+#include <typeindex>
#include <Eigen/Dense>
#include "nlohmann/json.hpp"
@@ -47,7 +48,8 @@ using REMatrixd = Eigen::Ref<const Eigen::Array<double, Eigen::Dynamic, Eigen::D
// Functions that return a double, take T and rhovec as arguments
#define ISOCHORIC_double_args \
X(get_pr) \
- X(get_splus)
+ X(get_splus) \
+ X(get_dpdT_constrhovec)
#define ISOCHORIC_array_args \
X(build_Psir_gradient_autodiff) \
@@ -55,12 +57,16 @@ using REMatrixd = Eigen::Ref<const Eigen::Array<double, Eigen::Dynamic, Eigen::D
X(get_dchempotdT_autodiff) \
X(get_fugacity_coefficients) \
X(get_partial_molar_volumes) \
- X(build_d2PsirdTdrhoi_autodiff)
+ X(build_d2PsirdTdrhoi_autodiff) \
+ X(get_dpdrhovec_constT)
#define ISOCHORIC_matrix_args \
X(build_Psir_Hessian_autodiff) \
X(build_Psi_Hessian_autodiff)
+#define ISOCHORIC_multimatrix_args \
+ X(build_Psir_fgradHessian_autodiff)
+
namespace teqp {
namespace cppinterface {
@@ -88,7 +94,11 @@ namespace teqp {
virtual AllowedModels& get_mutable_model() = 0;
virtual ~AbstractModel() = default;
+ virtual const std::type_index& get_type_index() const = 0;
+
+
virtual double get_R(const EArrayd&) const = 0;
+ double R(const EArrayd& x) const { return get_R(x); };
virtual double get_Arxy(const int, const int, const double, const double, const EArrayd&) const = 0;
@@ -100,7 +110,7 @@ namespace teqp {
#define X(i) virtual EArrayd get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefrac) const = 0;
AR0N_args
#undef X
- virtual double get_neff(const double, const double, const EArrayd&) const = 0;
+
// Virial derivatives
virtual double get_B2vir(const double T, const EArrayd& z) const = 0;
@@ -118,44 +128,51 @@ namespace teqp {
#define X(f) virtual EMatrixd f(const double T, const EArrayd& rhovec) const = 0;
ISOCHORIC_matrix_args
#undef X
+ #define X(f) virtual std::tuple<double, Eigen::ArrayXd, Eigen::MatrixXd> f(const double T, const EArrayd& rhovec) const = 0;
+ ISOCHORIC_multimatrix_args
+ #undef X
+ virtual Eigen::ArrayXd get_Psir_sigma_derivs(const double T, const EArrayd& rhovec, const EArrayd& v) const = 0;
+
+ double get_neff(const double, const double, const EArrayd&) const;
virtual EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z ) const = 0;
- virtual std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& = std::nullopt) const = 0;
- virtual std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tgiven) const = 0;
- virtual std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index=-1, int alternative_length=2) const = 0;
- virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
- virtual double dpsatdT_pure(const double T, const double rhoL, const double rhoV) const = 0;
-
- virtual std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
- virtual std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
- virtual double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
- virtual nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovec0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &) const = 0;
- virtual nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &) const = 0;
- virtual std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const = 0;
- virtual MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const = 0;
- virtual std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags) const = 0;
-
- virtual std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const = 0;
- virtual std::vector<nlohmann::json> find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options = std::nullopt) const = 0;
-
- virtual nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>&, const std::optional<TCABOptions> & = std::nullopt) const = 0;
- virtual EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const = 0;
- virtual double get_dp_dT_crit(const double T, const REArrayd& rhovec) const = 0;
- virtual EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const = 0;
- virtual EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>&) const = 0;
- virtual double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const = 0;
+ std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& = std::nullopt) const ;
+ EArray2 extrapolate_from_critical(const double Tc, const double rhoc, const double Tgiven) const;
+ std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index=-1, int alternative_length=2) const;
+
+ EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const;
+ double dpsatdT_pure(const double T, const double rhoL, const double rhoV) const;
+
+ virtual std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const;
+ virtual std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const;
+ virtual double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const;
+ virtual nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovec0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> & = std::nullopt) const;
+ virtual nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> & = std::nullopt) const;
+ virtual std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const;
+ virtual MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags = std::nullopt) const;
+ virtual std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags = std::nullopt) const;
+
+ std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const;
+ std::vector<nlohmann::json> find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options = std::nullopt) const;
+
+ virtual nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& = std::nullopt, const std::optional<TCABOptions> & = std::nullopt) const;
+ virtual EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const;
+ virtual double get_dp_dT_crit(const double T, const REArrayd& rhovec) const;
+ virtual EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const;
+ virtual EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& = std::nullopt) const;
+ virtual double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const;
};
// Generic JSON-based interface where the model description is encoded as JSON
- std::shared_ptr<AbstractModel> make_model(const nlohmann::json &);
+ std::unique_ptr<AbstractModel> make_model(const nlohmann::json &);
// Expose specialized factory functions for different models
// Mostly these are just adapter functions that prepare some
// JSON and pass it to the make_model function
// ....
- std::shared_ptr<AbstractModel> make_multifluid_model(
+ std::unique_ptr<AbstractModel> make_multifluid_model(
const std::vector<std::string>& components,
const std::string& coolprop_root,
const std::string& BIPcollectionpath = {},
@@ -163,6 +180,6 @@ namespace teqp {
const std::string& departurepath = {}
);
- std::shared_ptr<AbstractModel> emplace_model(AllowedModels&& model);
+ std::unique_ptr<AbstractModel> build_model_ptr(const nlohmann::json& json);
}
}
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index c6df0701b41e14999a5cf2aa93751a4ce9a69acf..bfc9432fb8ec4c4c83d8696d29b8aa50710c03a1 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -3,6 +3,7 @@
#include <complex>
#include <map>
#include <tuple>
+#include <numeric>
#include "teqp/types.hpp"
#include "teqp/exceptions.hpp"
@@ -31,6 +32,7 @@ typename ContainerType::value_type derivT(const FuncType& f, TType T, const Cont
return f(std::complex<TType>(T, h), rho).imag() / h;
}
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
/***
* \brief Given a function, use multicomplex derivatives to calculate the derivative with
* respect to the first variable which here is temperature
@@ -42,6 +44,7 @@ typename ContainerType::value_type derivTmcx(const FuncType& f, TType T, const C
auto ders = diff_mcx1(wrapper, T, 1);
return ders[0];
}
+#endif
/***
* \brief Given a function, use complex step derivatives to calculate the derivative with respect
@@ -163,11 +166,11 @@ struct TDXDerivatives {
}
}
else if constexpr (iT > 0 && iD == 0) {
- auto Trecip = 1.0 / T;
+ Scalar Trecip = 1.0 / T;
if constexpr (be == ADBackends::autodiff) {
// If a pure derivative, then we can use autodiff::Real for that variable and Scalar for other variable
autodiff::Real<iT, Scalar> Trecipad = Trecip;
- auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(1.0/Trecip__, rho, molefrac); };
+ auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(forceeval(1.0/Trecip__), rho, molefrac); };
return powi(Trecip, iT)*derivatives(f, along(1), at(Trecipad))[iT];
}
else if constexpr (iT == 1 && be == ADBackends::complex_step) {
@@ -196,7 +199,7 @@ struct TDXDerivatives {
return eval(w.alpha(T_, rho_, molefrac)); };
auto wrts = std::tuple_cat(build_duplicated_tuple<iT>(std::ref(Trecipad)), build_duplicated_tuple<iD>(std::ref(rhoad)));
auto der = derivatives(f, std::apply(wrt_helper(), wrts), at(Trecipad, rhoad));
- return powi(1.0 / T, iT) * powi(rho, iD) * der[der.size() - 1];
+ return powi(forceeval(1.0 / T), iT) * powi(rho, iD) * der[der.size() - 1];
}
#if defined(TEQP_MULTICOMPLEX_ENABLED)
else if constexpr (be == ADBackends::multicomplex) {
@@ -337,11 +340,11 @@ struct TDXDerivatives {
template<int Nderiv, ADBackends be = ADBackends::autodiff, class AlphaWrapper>
static auto get_Agenn0(const AlphaWrapper& w, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
std::valarray<Scalar> o(Nderiv+1);
- auto Trecip = 1.0 / T;
+ Scalar Trecip = 1.0 / T;
if constexpr (be == ADBackends::autodiff) {
// If a pure derivative, then we can use autodiff::Real for that variable and Scalar for other variable
autodiff::Real<Nderiv, Scalar> Trecipad = Trecip;
- auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(1.0/Trecip__, rho, molefrac); };
+ auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(forceeval(1.0/Trecip__), rho, molefrac); };
auto ders = derivatives(f, along(1), at(Trecipad));
for (auto n = 0; n <= Nderiv; ++n) {
o[n] = powi(Trecip, n) * ders[n];
@@ -483,9 +486,10 @@ struct VirialDerivatives {
{
std::map<int, double> dnalphardrhon;
if constexpr(be == ADBackends::autodiff){
- autodiff::HigherOrderDual<Nderiv, double> rhodual = 0.0;
auto f = [&model, &T, &molefrac](const auto& rho_) { return model.alphar(T, rho_, molefrac); };
- auto derivs = derivatives(f, wrt(rhodual), at(rhodual));
+ autodiff::Real<Nderiv, Scalar> rhoreal = 0.0;
+ auto derivs = derivatives(f, along(1), at(rhoreal));
+
for (auto n = 1; n < Nderiv; ++n){
dnalphardrhon[n] = derivs[n];
}
@@ -519,7 +523,7 @@ struct VirialDerivatives {
}
/// This version of the get_Bnvir takes the maximum number of derivatives as a runtime argument
- /// and then forwards all arguments to the templated function
+ /// and then forwards all arguments to the corresponding templated function
template <ADBackends be = ADBackends::autodiff>
static auto get_Bnvir_runtime(const int Nderiv, const Model& model, const Scalar &T, const VectorType& molefrac) {
switch(Nderiv){
@@ -698,7 +702,7 @@ struct IsochoricDerivatives{
*/
static auto get_dPsirdT_constrhovec(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot_ = rhovec.sum();
- auto molefrac = rhovec / rhotot_;
+ auto molefrac = (rhovec / rhotot_).eval();
autodiff::Real<1, Scalar> Tad = T;
auto f = [&model, &rhotot_, &molefrac](const auto& T_) {return rhotot_*model.R(molefrac)*T_*model.alphar(T_, rhotot_, molefrac); };
return derivatives(f, along(1), at(Tad))[1];
@@ -764,6 +768,7 @@ struct IsochoricDerivatives{
return H;
}
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
/***
* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations (residual contribution only)
*
@@ -785,6 +790,7 @@ struct IsochoricDerivatives{
auto H = get_Hessian<mattype, fcn_t, VectorType, HessianMethods::Multiple>(func, rho);
return H;
}
+#endif
/***
* \brief Gradient of Psir = ar*rho w.r.t. the molar concentrations
@@ -802,6 +808,7 @@ struct IsochoricDerivatives{
return val;
}
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
/***
* \brief Gradient of Psir = ar*rho w.r.t. the molar concentrations
*
@@ -821,6 +828,7 @@ struct IsochoricDerivatives{
}
return out;
}
+#endif
/***
* \brief Gradient of Psir = ar*rho w.r.t. the molar concentrations
*
@@ -963,6 +971,21 @@ struct IsochoricDerivatives{
auto denominator = -pow2(rhotot)*RT*(1 + 2*ders[1] + ders[2]);
return (numerator/denominator).eval();
}
+
+ static VectorType get_Psir_sigma_derivs(const Model& model, const Scalar& T, const VectorType& rhovec, const VectorType& v) {
+ autodiff::Real<4, double> sigma = 0.0;
+ auto rhovecad = rhovec.template cast<decltype(sigma)>(), vad = v.template cast<decltype(sigma)>();
+ auto wrapper = [&rhovecad, &vad, &T, &model](const auto& sigma_1) {
+ auto rhovecused = (rhovecad + sigma_1 * vad).eval();
+ auto rhotot = rhovecused.sum();
+ auto molefrac = (rhovecused / rhotot).eval();
+ return forceeval(model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot);
+ };
+ auto der = derivatives(wrapper, along(1), at(sigma));
+ VectorType ret(der.size());
+ for (auto i = 0; i < ret.size(); ++i){ ret[i] = der[i];}
+ return ret;
+ }
};
template<int Nderivsmax, AlphaWrapperOption opt>
diff --git a/include/teqp/exceptions.hpp b/include/teqp/exceptions.hpp
index 90080a39376a2a7a5a5ed12cb969225fd108ded2..4b6d1e4323fc87f50dd5e3e6cbda6edf5447a65d 100644
--- a/include/teqp/exceptions.hpp
+++ b/include/teqp/exceptions.hpp
@@ -38,4 +38,9 @@ namespace teqp {
};
using IterationError = IterationFailure;
+ class NotImplementedError : public teqpException {
+ public:
+ NotImplementedError(const std::string& msg) : teqpException(200, msg) {};
+ };
+
}; // namespace teqp
diff --git a/include/teqp/json_builder.hpp b/include/teqp/json_builder.hpp
index 4b795d2f76e74adc4379d9ac75cd736d499be4a1..0b2d73f3cda6509d94a26f846d33719d738980fa 100644
--- a/include/teqp/json_builder.hpp
+++ b/include/teqp/json_builder.hpp
@@ -5,8 +5,10 @@
#include "teqp/models/mie/mie.hpp"
#include "teqp/exceptions.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
#include "nlohmann/json.hpp"
+#include <memory>
namespace teqp {
@@ -42,46 +44,7 @@ namespace teqp {
return CPA::CPAfactory(spec);
}
else if (kind == "PCSAFT") {
- using namespace PCSAFT;
- std::optional<Eigen::ArrayXXd> kmat;
- if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
- kmat = build_square_matrix(spec["kmat"]);
- }
-
- if (spec.contains("names")){
- std::vector<std::string> names = spec["names"];
- if (kmat && kmat.value().rows() != names.size()){
- throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
- }
- return PCSAFTMixture(names, kmat.value_or(Eigen::ArrayXXd{}));
- }
- else if (spec.contains("coeffs")){
- std::vector<SAFTCoeffs> coeffs;
- for (auto j : spec["coeffs"]) {
- SAFTCoeffs c;
- c.name = j.at("name");
- c.m = j.at("m");
- c.sigma_Angstrom = j.at("sigma_Angstrom");
- c.epsilon_over_k = j.at("epsilon_over_k");
- c.BibTeXKey = j.at("BibTeXKey");
- if (j.contains("(mu^*)^2") && j.contains("nmu")){
- c.mustar2 = j.at("(mu^*)^2");
- c.nmu = j.at("nmu");
- }
- if (j.contains("(Q^*)^2") && j.contains("nQ")){
- c.Qstar2 = j.at("(Q^*)^2");
- c.nQ = j.at("nQ");
- }
- coeffs.push_back(c);
- }
- if (kmat && kmat.value().rows() != coeffs.size()){
- throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
- }
- return PCSAFTMixture(coeffs, kmat.value_or(Eigen::ArrayXXd{}));
- }
- else{
- throw std::invalid_argument("you must provide names or coeffs, but not both");
- }
+ return PCSAFT::PCSAFTfactory(spec);
}
else if (kind == "SAFT-VR-Mie") {
return SAFTVRMie::SAFTVRMiefactory(spec);
diff --git a/include/teqp/math/pow_templates.hpp b/include/teqp/math/pow_templates.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..aebbc96ab17e7588d980c05fd74186c8d2268f95
--- /dev/null
+++ b/include/teqp/math/pow_templates.hpp
@@ -0,0 +1,13 @@
+#pragma once
+#include "teqp/types.hpp"
+
+namespace teqp{
+template<typename A> auto POW2(const A& x) { return forceeval(x*x); }
+template<typename A> auto POW3(const A& x) { return forceeval(POW2(x)*x); }
+template<typename A> auto POW4(const A& x) { return forceeval(POW2(x)*POW2(x)); }
+template<typename A> auto POW5(const A& x) { return forceeval(POW2(x)*POW3(x)); }
+template<typename A> auto POW7(const A& x) { return forceeval(POW2(x)*POW5(x)); }
+template<typename A> auto POW8(const A& x) { return forceeval(POW4(x)*POW4(x)); }
+template<typename A> auto POW10(const A& x) { return forceeval(POW2(x)*POW8(x)); }
+template<typename A> auto POW12(const A& x) { return forceeval(POW4(x)*POW8(x)); }
+}
diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp
index 6d8dd256b940fb823097e4c681e4fc2b1724f56a..cd1e16d398af4ed5c3b3214e5cc8c743aac7c91c 100644
--- a/include/teqp/models/cubics.hpp
+++ b/include/teqp/models/cubics.hpp
@@ -12,6 +12,7 @@ Implementations of the canonical cubic equations of state
#include "teqp/constants.hpp"
#include "teqp/exceptions.hpp"
#include "cubicsuperancillary.hpp"
+#include "teqp/json_tools.hpp"
#include "nlohmann/json.hpp"
@@ -20,16 +21,16 @@ Implementations of the canonical cubic equations of state
namespace teqp {
/**
-* \brief The standard alpha function used by Peng-Robinson and SRK
-*/
+ * \brief The standard alpha function used by Peng-Robinson and SRK
+ */
template<typename NumType>
class BasicAlphaFunction {
private:
NumType Tci, ///< The critical temperature
- mi; ///< The "m" parameter
+ mi; ///< The "m" parameter
public:
BasicAlphaFunction(NumType Tci, NumType mi) : Tci(Tci), mi(mi) {};
-
+
template<typename TType>
auto operator () (const TType& T) const {
return forceeval(pow2(forceeval(1.0 + mi * (1.0 - sqrt(T / Tci)))));
@@ -44,18 +45,18 @@ class GenericCubic {
protected:
std::valarray<NumType> ai, bi;
const NumType Delta1, Delta2, OmegaA, OmegaB;
- int superanc_index;
+ int superanc_index;
const AlphaFunctions alphas;
Eigen::ArrayXXd kmat;
-
+
nlohmann::json meta;
-
+
template<typename TType, typename IndexType>
auto get_ai(TType T, IndexType i) const { return ai[i]; }
-
+
template<typename TType, typename IndexType>
auto get_bi(TType T, IndexType i) const { return bi[i]; }
-
+
template<typename IndexType>
void check_kmat(IndexType N) {
if (kmat.cols() != kmat.rows()) {
@@ -68,10 +69,10 @@ protected:
throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components [" + std::to_string(N) + "]");
}
};
-
+
public:
GenericCubic(NumType Delta1, NumType Delta2, NumType OmegaA, NumType OmegaB, int superanc_index, const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa, const AlphaFunctions& alphas, const Eigen::ArrayXXd& kmat)
- : Delta1(Delta1), Delta2(Delta2), OmegaA(OmegaA), OmegaB(OmegaB), superanc_index(superanc_index), alphas(alphas), kmat(kmat)
+ : Delta1(Delta1), Delta2(Delta2), OmegaA(OmegaA), OmegaB(OmegaB), superanc_index(superanc_index), alphas(alphas), kmat(kmat)
{
ai.resize(Tc_K.size());
bi.resize(Tc_K.size());
@@ -81,13 +82,13 @@ public:
}
check_kmat(ai.size());
};
-
+
void set_meta(const nlohmann::json& j) { meta = j; }
auto get_meta() const { return meta; }
auto get_kmat() const { return kmat; }
-
+
/// Return a tuple of saturated liquid and vapor densities for the EOS given the temperature
- /// Uses the superancillary equations from Bell and Deiters:
+ /// Uses the superancillary equations from Bell and Deiters:
auto superanc_rhoLV(double T) const {
if (ai.size() != 1) {
throw std::invalid_argument("function only available for pure species");
@@ -96,18 +97,18 @@ public:
auto b = get_b(T, z);
auto Ttilde = R(z)*T*b/get_a(T,z);
return std::make_tuple(
- CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)/b,
- CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)/b
- );
+ CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)/b,
+ CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)/b
+ );
}
-
+
const NumType Ru = get_R_gas<double>(); /// Universal gas constant, exact number
-
+
template<class VecType>
auto R(const VecType& molefrac) const {
return Ru;
}
-
+
template<typename TType, typename CompType>
auto get_a(TType T, const CompType& molefracs) const {
std::common_type_t<TType, decltype(molefracs[0])> a_ = 0.0;
@@ -124,7 +125,7 @@ public:
}
return forceeval(a_);
}
-
+
template<typename TType, typename CompType>
auto get_b(TType /*T*/, const CompType& molefracs) const {
std::common_type_t<TType, decltype(molefracs[0])> b_ = 0.0;
@@ -133,11 +134,11 @@ public:
}
return forceeval(b_);
}
-
+
template<typename TType, typename RhoType, typename MoleFracType>
auto alphar(const TType& T,
- const RhoType& rho,
- const MoleFracType& molefrac) const
+ const RhoType& rho,
+ const MoleFracType& molefrac) const
{
if (molefrac.size() != alphas.size()) {
throw std::invalid_argument("Sizes do not match");
@@ -155,16 +156,16 @@ auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen:
double Delta1 = 1;
double Delta2 = 0;
AcentricType m = 0.48 + 1.574 * acentric - 0.176 * acentric * acentric;
-
+
std::vector<AlphaFunctionOptions> alphas;
for (auto i = 0; i < Tc_K.size(); ++i) {
alphas.emplace_back(BasicAlphaFunction(Tc_K[i], m[i]));
}
-
+
// See https://doi.org/10.1021/acs.iecr.1c00847
double OmegaA = 1.0 / (9.0 * (cbrt(2) - 1));
double OmegaB = (cbrt(2) - 1) / 3;
-
+
nlohmann::json meta = {
{"Delta1", Delta1},
{"Delta2", Delta2},
@@ -172,12 +173,22 @@ auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen:
{"OmegaB", OmegaB},
{"kind", "Soave-Redlich-Kwong"}
};
-
+
auto cub = GenericCubic(Delta1, Delta2, OmegaA, OmegaB, CubicSuperAncillary::SRK_CODE, Tc_K, pc_K, std::move(alphas), kmat);
cub.set_meta(meta);
return cub;
}
+/// A JSON-based factory function for the canonical SRK model
+inline auto make_canonicalSRK(const nlohmann::json& spec){
+ std::valarray<double> Tc_K = spec.at("Tcrit / K"), pc_Pa = spec.at("pcrit / Pa"), acentric = spec.at("acentric");
+ Eigen::ArrayXXd kmat(0, 0);
+ if (spec.contains("kmat")){
+ kmat = build_square_matrix(spec.at("kmat"));
+ }
+ return canonical_SRK(Tc_K, pc_Pa, acentric, kmat);
+}
+
template <typename TCType, typename PCType, typename AcentricType>
auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::ArrayXXd& kmat = {}) {
double Delta1 = 1+sqrt(2.0);
@@ -211,4 +222,14 @@ auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::
return cub;
}
+/// A JSON-based factory function for the canonical SRK model
+inline auto make_canonicalPR(const nlohmann::json& spec){
+ std::valarray<double> Tc_K = spec.at("Tcrit / K"), pc_Pa = spec.at("pcrit / Pa"), acentric = spec.at("acentric");
+ Eigen::ArrayXXd kmat(0, 0);
+ if (spec.contains("kmat")){
+ kmat = build_square_matrix(spec.at("kmat"));
+ }
+ return canonical_PR(Tc_K, pc_Pa, acentric, kmat);
+}
+
}; // namespace teqp
diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index f2c0ae06914156f52b2be53f87cc537b2e896570..f9fc2964796561910be79ab3ae78d86492c851e7 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -14,13 +14,16 @@
#include "teqp/json_tools.hpp"
#include "teqp/exceptions.hpp"
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
#include "MultiComplex/MultiComplex.hpp"
+#endif
#include "multifluid_eosterms.hpp"
#include "multifluid_reducing.hpp"
#include <boost/algorithm/string/join.hpp>
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
// See https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
namespace Eigen {
template<typename TN> struct NumTraits<mcx::MultiComplex<TN>> : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
@@ -36,6 +39,7 @@ namespace Eigen {
};
};
}
+#endif
namespace teqp{
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index 324c134cea6ef2ba8e7a866bcb7f52674f50a789..04b534b0476045b4273655086b2e64f28b4cc49a 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -410,19 +410,48 @@ public:
}
};
-inline auto PCSAFTfactory(const nlohmann::json& json) {
- std::vector<SAFTCoeffs> coeffs;
- for (auto j : json) {
- SAFTCoeffs c;
- c.name = j.at("name");
- c.m = j.at("m");
- c.sigma_Angstrom = j.at("sigma_Angstrom");
- c.epsilon_over_k = j.at("epsilon_over_k");
- c.BibTeXKey = j.at("BibTeXKey");
- coeffs.push_back(c);
+/// A JSON-based factory function for the PC-SAFT model
+inline auto PCSAFTfactory(const nlohmann::json& spec) {
+ std::optional<Eigen::ArrayXXd> kmat;
+ if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
+ kmat = build_square_matrix(spec["kmat"]);
}
- return PCSAFTMixture(coeffs);
-};
+
+ if (spec.contains("names")){
+ std::vector<std::string> names = spec["names"];
+ if (kmat && kmat.value().rows() != names.size()){
+ throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
+ }
+ return PCSAFTMixture(names, kmat.value_or(Eigen::ArrayXXd{}));
+ }
+ else if (spec.contains("coeffs")){
+ std::vector<SAFTCoeffs> coeffs;
+ for (auto j : spec["coeffs"]) {
+ SAFTCoeffs c;
+ c.name = j.at("name");
+ c.m = j.at("m");
+ c.sigma_Angstrom = j.at("sigma_Angstrom");
+ c.epsilon_over_k = j.at("epsilon_over_k");
+ c.BibTeXKey = j.at("BibTeXKey");
+ if (j.contains("(mu^*)^2") && j.contains("nmu")){
+ c.mustar2 = j.at("(mu^*)^2");
+ c.nmu = j.at("nmu");
+ }
+ if (j.contains("(Q^*)^2") && j.contains("nQ")){
+ c.Qstar2 = j.at("(Q^*)^2");
+ c.nQ = j.at("nQ");
+ }
+ coeffs.push_back(c);
+ }
+ if (kmat && kmat.value().rows() != coeffs.size()){
+ throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
+ }
+ return PCSAFTMixture(coeffs, kmat.value_or(Eigen::ArrayXXd{}));
+ }
+ else{
+ throw std::invalid_argument("you must provide names or coeffs, but not both");
+ }
+}
} /* namespace PCSAFT */
}; // namespace teqp
diff --git a/include/teqp/models/saft/correlation_integrals.hpp b/include/teqp/models/saft/correlation_integrals.hpp
index 6b1084aa75e8396a2bca13a3cddf13299cbe5978..2a2e0416cd8b6cd718161cc8f652fa0f4848b905 100644
--- a/include/teqp/models/saft/correlation_integrals.hpp
+++ b/include/teqp/models/saft/correlation_integrals.hpp
@@ -7,6 +7,12 @@
namespace teqp {
namespace SAFTpolar{
+// |x| = sqrt(x^2), the latter is well-suited to differentiation
+template<typename X>
+inline auto differentiable_abs(const X& x){
+ return forceeval(sqrt(x*x));
+};
+
static const std::map<int, std::array<double, 12>> Luckas_J_coeffs = {
{4, {-1.38410152e00, -7.05792933e-01, 2.60947023e00, 1.96828333e01, 1.13619510e01, -2.98510490e01, -3.15686398e01, -2.00943290e01, 5.11029320e01, 1.44194150e01, 9.40061069e00, -2.36844608e01}},
{5, {-6.89702637e-01, -1.62382602e-01, 1.16302441e00, 1.42067443e01, 4.59642681e00, -1.81421003e01, -2.45012804e01, -8.42839734e00, 3.25579587e01, 1.16339969e01, 4.00080085e00, -1.54419815e01}},
@@ -17,7 +23,7 @@ static const std::map<int, std::array<double, 12>> Luckas_J_coeffs = {
{10, {-2.30585563e-01, 1.86890174e-02, 3.31660880e-01, 8.88275358e00, 5.90996555e-01, -9.19786291e00, -1.79530632e01, -1.12502567e00, 1.87717642e01, 9.25235661e00, 5.55909350e-01, -9.47729033e00}},
{11, {-2.19221482e-01, 2.05982865e-02, 3.05823216e-01, 8.58399185e00, 4.60555635e-01, -8.78667564e00, -1.78552884e01, -8.80186074e-01, 1.84146650e01, 9.31691675e00, 4.43027169e-01, -9.40608814e00}},
{12, {-2.13591301e-01, 2.15809377e-02, 2.89113218e-01, 8.37291362e00, 3.68864477e-01, -8.49747359e00, -1.79188678e01, -7.07751483e-01, 1.82906357e01, 9.45598216e00, 3.64636215e-01, -9.44579903e00}},
- {13, {-2.11696363e-01, 2.20767096e-02, 2.78453228e-01, 8.23384704e00, 3.02478832e-01, -8.30238277e00, 1.81183432e01, -5.82685843e-01, 1.83497132e01, 9.65764061e00, 3.08765033e-01, -9.57226243e00}},
+ {13, {-2.11696363e-01, 2.20767096e-02, 2.78453228e-01, 8.23384704e00, 3.02478832e-01, -8.30238277e00, -1.81183432e01, -5.82685843e-01, 1.83497132e01, 9.65764061e00, 3.08765033e-01, -9.57226243e00}},
{14, {-2.12238116e-01, 2.23157164e-02, 2.71883622e-01, 8.15208365e00, 2.53073339e-01, -8.17919021e00, -1.84273649e01, -4.89407317e-01, 1.85502740e01, 9.90906951e00, 2.67965960e-01, -9.76495200e00}},
{15, {-2.14335965e-01, 2.24240261e-02, 2.68094773e-01, 8.11899188e00, 2.15506735e-01, -8.11465705e00, -1.88310645e01, -4.18309476e-01, 1.88679367e01, 1.02033085e01, 2.37674032e-01, -1.00120648e01}},
};
@@ -38,14 +44,12 @@ public:
template<typename TType, typename RhoType>
auto get_J(const TType& Tstar, const RhoType& rhostar) const{
- auto Z_1 = 0.3 + 0.05*n;
- auto Z_2 = 1.0/n;
- auto A_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
- auto A_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
- auto A_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
- // |x| = sqrt(x^2), the latter is well-suited to differentiation
- auto differentiable_abs = [](const auto& x){ return sqrt(x*x); };
- std::common_type_t<TType, RhoType> out = (A_0 + A_1*pow(Tstar, Z_1) + A_2*pow(Tstar, Z_2))*exp(1.0/(Tstar + 4.0/pow(differentiable_abs(log(rhostar/sqrt(2.0))), 3.0)));
+ double Z_1 = 0.3 + 0.05*n;
+ double Z_2 = 1.0/n;
+ RhoType A_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
+ RhoType A_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
+ RhoType A_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
+ std::common_type_t<TType, RhoType> out = (A_0 + A_1*pow(Tstar, Z_1) + A_2*pow(Tstar, Z_2))*exp(1.0/(Tstar + 4.0/pow(differentiable_abs(log(forceeval(rhostar/sqrt(2.0)))), 3.0)));
return out;
}
};
@@ -73,13 +77,13 @@ public:
template<typename TType, typename RhoType>
auto get_K(const TType& Tstar, const RhoType& rhostar) const{
- auto Z_1 = 2.0;
- auto Z_2 = 3.0;
- auto Z_3 = 4.0;
- auto b_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
- auto b_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
- auto b_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
- auto b_3 = a03 + a13*rhostar + a23*rhostar*rhostar + a33*rhostar*rhostar*rhostar;
+ double Z_1 = 2.0;
+ double Z_2 = 3.0;
+ double Z_3 = 4.0;
+ RhoType b_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
+ RhoType b_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
+ RhoType b_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
+ RhoType b_3 = a03 + a13*rhostar + a23*rhostar*rhostar + a33*rhostar*rhostar*rhostar;
std::common_type_t<TType, RhoType> out = b_0 + b_1*Tstar + b_2*pow(exp(pow(1.0-rhostar/sqrt(2.0),Z_3)), Z_1) + b_3*pow(exp(pow(1.0-rhostar/sqrt(2.0),Z_3)), Z_2);
return out;
}
diff --git a/include/teqp/models/saft/polar_terms.hpp b/include/teqp/models/saft/polar_terms.hpp
index 96389a5f5d76e51821b5c29f53b737649053918a..25ddbcf25a50db5743eb601cb9bc8973b752f177 100644
--- a/include/teqp/models/saft/polar_terms.hpp
+++ b/include/teqp/models/saft/polar_terms.hpp
@@ -12,20 +12,13 @@
#include "correlation_integrals.hpp"
#include <optional>
#include <Eigen/Dense>
+#include "teqp/math/pow_templates.hpp"
+#include <variant>
namespace teqp{
namespace SAFTpolar{
-template<typename A> auto POW2(const A& x) { return forceeval(x*x); }
-template<typename A> auto POW3(const A& x) { return forceeval(POW2(x)*x); }
-template<typename A> auto POW4(const A& x) { return forceeval(POW2(x)*POW2(x)); }
-template<typename A> auto POW5(const A& x) { return forceeval(POW2(x)*POW3(x)); }
-template<typename A> auto POW7(const A& x) { return forceeval(POW2(x)*POW5(x)); }
-template<typename A> auto POW8(const A& x) { return forceeval(POW4(x)*POW4(x)); }
-template<typename A> auto POW10(const A& x) { return forceeval(POW2(x)*POW8(x)); }
-template<typename A> auto POW12(const A& x) { return forceeval(POW4(x)*POW8(x)); }
-
/// Eq. 10 from Gross and Vrabec
template <typename Eta, typename MType, typename TType>
auto get_JDD_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
@@ -37,7 +30,7 @@ auto get_JDD_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(5) << -0.5873164, 1.2489132, -0.5085280, 0, 0).finished();
static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(5) << 3.4869576, -14.915974, 15.372022, 0, 0).finished();
- std::common_type_t<Eta, MType> summer = 0.0;
+ std::common_type_t<Eta, MType, TType> summer = 0.0;
for (auto n = 0; n < 5; ++n){
auto anij = a_0[n] + (mij-1)/mij*a_1[n] + (mij-1)/mij*(mij-2)/mij*a_2[n]; // Eq. 12
auto bnij = b_0[n] + (mij-1)/mij*b_1[n] + (mij-1)/mij*(mij-2)/mij*b_2[n]; // Eq. 13
@@ -71,7 +64,7 @@ auto get_JQQ_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(5) << -0.8137340, 10.064030, -10.876631, 0.0, 0.0).finished();
static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(5) << 6.8682675, -5.1732238, -17.240207, 0.0, 0.0).finished();
- std::common_type_t<Eta, MType> summer = 0.0;
+ std::common_type_t<Eta, MType, TType> summer = 0.0;
for (auto n = 0; n < 5; ++n){
auto anij = a_0[n] + (mij-1)/mij*a_1[n] + (mij-1)/mij*(mij-2)/mij*a_2[n]; // Eq. 12
auto bnij = b_0[n] + (mij-1)/mij*b_1[n] + (mij-1)/mij*(mij-2)/mij*b_2[n]; // Eq. 13
@@ -127,7 +120,7 @@ public:
auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
auto sigmaij = (sigma[i]+sigma[j])/2;
- auto Tstarij = T/epskij;
+ auto Tstarij = forceeval(T/epskij);
auto mij = std::min(sqrt(m[i]*m[j]), 2.0);
summer += x[i]*x[j]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*POW3(sigma[i]*sigma[j]/sigmaij)*ninj*mustar2[i]*mustar2[j]*get_JDD_2ij(eta, mij, Tstarij);
}
@@ -219,7 +212,7 @@ public:
auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
auto sigmaij = (sigma[i]+sigma[j])/2;
- auto Tstarij = T/epskij;
+ auto Tstarij = forceeval(T/epskij);
auto mij = std::min(sqrt(m[i]*m[j]), 2.0);
summer += x[i]*x[j]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*ninj*Qstar2[i]*Qstar2[j]*get_JQQ_2ij(eta, mij, Tstarij);
}
@@ -233,11 +226,11 @@ public:
auto get_alpha3QQ(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
const auto& x = mole_fractions; // concision
const auto& sigma = sigma_Angstrom; // concision
- const auto N = mole_fractions.size();
+ const std::size_t N = mole_fractions.size();
std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
- for (auto k = 0; k < N; ++k){
+ for (std::size_t i = 0; i < N; ++i){
+ for (std::size_t j = 0; j < N; ++j){
+ for (std::size_t k = 0; k < N; ++k){
auto ninjnk = nQ[i]*nQ[j]*nQ[k];
if (ninjnk > 0){
// Lorentz-Berthelot mixing rules for sigma
@@ -280,6 +273,17 @@ enum class multipolar_argument_spec {
TK_rhoNm3_molefractions
};
+template<typename type>
+struct MultipolarContributionGrossVrabecTerms{
+ type alpha2DD;
+ type alpha3DD;
+ type alphaDD;
+ type alpha2QQ;
+ type alpha3QQ;
+ type alphaQQ;
+ type alpha;
+};
+
class MultipolarContributionGrossVrabec{
public:
static constexpr multipolar_argument_spec arg_spec = multipolar_argument_spec::TK_rhoNA3_packingfraction_molefractions;
@@ -302,7 +306,7 @@ public:
template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
auto eval(const TTYPE& T, const RhoType& rho_A3, const EtaType& eta, const VecType& mole_fractions) const {
- using type = std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])>;
+ using type = std::common_type_t<TTYPE, RhoType, EtaType, decltype(mole_fractions[0])>;
type alpha2DD = 0.0, alpha3DD = 0.0, alphaDD = 0.0;
if (di && di.value().has_a_polar){
alpha2DD = di.value().get_alpha2DD(T, rho_A3, eta, mole_fractions);
@@ -319,20 +323,31 @@ public:
auto alpha = forceeval(alphaDD + alphaQQ);
- struct Terms{
- type alpha2DD;
- type alpha3DD;
- type alphaDD;
- type alpha2QQ;
- type alpha3QQ;
- type alphaQQ;
- type alpha;
- };
- return Terms{alpha2DD, alpha3DD, alphaDD, alpha2QQ, alpha3QQ, alphaQQ, alpha};
+ return MultipolarContributionGrossVrabecTerms<type>{alpha2DD, alpha3DD, alphaDD, alpha2QQ, alpha3QQ, alphaQQ, alpha};
}
};
+template<typename type>
+struct MultipolarContributionGubbinsTwuTermsGT{
+ type alpha2;
+ type alpha3;
+ type alpha;
+};
+
+template<typename KType, typename RhoType, typename Txy>
+auto get_Kijk(const KType& Kint, const RhoType& rhostar, const Txy &Tstarij, const Txy &Tstarik, const Txy &Tstarjk){
+ return forceeval(pow(forceeval(Kint.get_K(Tstarij, rhostar)*Kint.get_K(Tstarik, rhostar)*Kint.get_K(Tstarjk, rhostar)), 1.0/3.0));
+};
+
+// Special treatment needed here because the 334,445 term is negative, so negative*negative*negative is negative, and negative^{1/3} is undefined
+// First flip the sign on the triple product, do the evaluation, the flip it back. Not documented in Gubbins&Twu, but this seem reasonable,
+// in the spirit of the others.
+template<typename KType, typename RhoType, typename Txy>
+auto get_Kijk_334445(const KType& Kint, const RhoType& rhostar, const Txy &Tstarij, const Txy &Tstarik, const Txy &Tstarjk){
+ return forceeval(-pow(-forceeval(Kint.get_K(Tstarij, rhostar)*Kint.get_K(Tstarik, rhostar)*Kint.get_K(Tstarjk, rhostar)), 1.0/3.0));
+};
+
/**
\tparam JIntegral A type that can be indexed with a single integer n to give the J^{(n)} integral
\tparam KIntegral A type that can be indexed with a two integers a and b to give the K(a,b) integral
@@ -346,6 +361,7 @@ public:
private:
const Eigen::ArrayXd sigma_m, epsilon_over_k, mubar2, Qbar2;
const bool has_a_polar;
+ const Eigen::ArrayXd sigma_m3, sigma_m5;
const JIntegral J6{6};
const JIntegral J8{8};
@@ -359,9 +375,10 @@ private:
const KIntegral K444_555{444, 555};
const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2
const double PI_ = static_cast<double>(EIGEN_PI);
+ Eigen::MatrixXd SIGMAIJ, EPSKIJ;
public:
- MultipolarContributionGubbinsTwu(const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mubar2, const Eigen::ArrayX<double> &Qbar2) : sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2), has_a_polar(mubar2.cwiseAbs().sum() > 0 || Qbar2.cwiseAbs().sum() > 0) {
+ MultipolarContributionGubbinsTwu(const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mubar2, const Eigen::ArrayX<double> &Qbar2) : sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2), has_a_polar(mubar2.cwiseAbs().sum() > 0 || Qbar2.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)) {
// Check lengths match
if (sigma_m.size() != mubar2.size()){
throw teqp::InvalidArgument("bad size of mubar2");
@@ -369,121 +386,146 @@ public:
if (sigma_m.size() != Qbar2.size()){
throw teqp::InvalidArgument("bad size of Qbar2");
}
+ // Pre-calculate the mixing terms;
+ const std::size_t N = sigma_m.size();
+ SIGMAIJ.resize(N, N); EPSKIJ.resize(N, N);
+ for (std::size_t i = 0; i < N; ++i){
+ for (std::size_t j = 0; j < N; ++j){
+ // Lorentz-Berthelot mixing rules
+ double epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
+ double sigmaij = (sigma_m[i]+sigma_m[j])/2;
+ SIGMAIJ(i, j) = sigmaij;
+ EPSKIJ(i, j) = epskij;
+ }
+ }
}
MultipolarContributionGubbinsTwu& operator=( const MultipolarContributionGubbinsTwu& ) = delete; // non copyable
template<typename TTYPE, typename RhoType, typename RhoStarType, typename VecType>
auto get_alpha2(const TTYPE& T, const RhoType& rhoN, const RhoStarType& rhostar, const VecType& mole_fractions) const{
const auto& x = mole_fractions; // concision
- const auto& sigma = sigma_m; // concision
- const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> alpha2_112 = 0.0, alpha2_123 = 0.0, alpha2_224 = 0.0;
+ const std::size_t N = mole_fractions.size();
+ using XTtype = std::common_type_t<TTYPE, decltype(mole_fractions[0])>;
+ std::common_type_t<TTYPE, RhoType, RhoStarType, decltype(mole_fractions[0])> alpha2_112 = 0.0, alpha2_123 = 0.0, alpha2_224 = 0.0;
- const auto factor_112 = forceeval(-2.0*PI_*rhoN/3.0); //*POW2(4*PI_*epsilon_0)
- const auto factor_123 = forceeval(-PI_*rhoN/3.0);
- const auto factor_224 = forceeval(-14.0*PI_*rhoN/5.0);
-
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
- // Lorentz-Berthelot mixing rules
- auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
- auto sigmaij = (sigma[i]+sigma[j])/2;
-
- auto Tstari = T/epsilon_over_k[i], Tstarj = T/epsilon_over_k[j];
- auto leading = x[i]*x[j]/(Tstari*Tstarj); // common for all alpha_2 terms
- auto Tstarij = forceeval(T/epskij);
-
- alpha2_112 += factor_112*leading*POW3(sigma[i]*sigma[j])/POW3(sigmaij)*mubar2[i]*mubar2[j]*J6.get_J(Tstarij, rhostar);
- alpha2_123 += factor_123*leading*POW3(sigma[i])*POW5(sigma[j])/POW5(sigmaij)*mubar2[i]*Qbar2[j]*J8.get_J(Tstarij, rhostar);
- alpha2_224 += factor_224*leading*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*Qbar2[i]*Qbar2[j]*J10.get_J(Tstarij, rhostar);
+ const RhoType factor_112 = -2.0*PI_*rhoN/3.0;
+ const RhoType factor_123 = -PI_*rhoN/3.0;
+ const RhoType factor_224 = -14.0*PI_*rhoN/5.0;
+
+ for (std::size_t i = 0; i < N; ++i){
+ for (std::size_t j = 0; j < N; ++j){
+
+ TTYPE Tstari = forceeval(T/EPSKIJ(i, i)), Tstarj = forceeval(T/EPSKIJ(j, j));
+ XTtype leading = forceeval(x[i]*x[j]/(Tstari*Tstarj)); // common for all alpha_2 terms
+ TTYPE Tstarij = forceeval(T/EPSKIJ(i, j));
+ double sigmaij = SIGMAIJ(i,j);
+ {
+ double dbl = sigma_m3[i]*sigma_m3[j]/powi(sigmaij,3)*mubar2[i]*mubar2[j];
+ alpha2_112 += leading*dbl*J6.get_J(Tstarij, rhostar);
+ }
+ {
+ double dbl = sigma_m3[i]*sigma_m5[j]/powi(sigmaij,5)*mubar2[i]*Qbar2[j];
+ alpha2_123 += leading*dbl*J8.get_J(Tstarij, rhostar);
+ }
+ {
+ double dbl = sigma_m5[i]*sigma_m5[j]/powi(sigmaij,7)*Qbar2[i]*Qbar2[j];
+ alpha2_224 += leading*dbl*J10.get_J(Tstarij, rhostar);
+ }
}
}
-
- return forceeval(alpha2_112 + 2.0*alpha2_123 + alpha2_224);
+ return forceeval(factor_112*alpha2_112 + 2.0*factor_123*alpha2_123 + factor_224*alpha2_224);
}
+
template<typename TTYPE, typename RhoType, typename RhoStarType, typename VecType>
auto get_alpha3(const TTYPE& T, const RhoType& rhoN, const RhoStarType& rhostar, const VecType& mole_fractions) const{
- const auto& x = mole_fractions; // concision
- const auto& sigma = sigma_m; // concision
- const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summerA_112_112_224 = 0.0, summerA_112_123_213 = 0.0, summerA_123_123_224 = 0.0, summerA_224_224_224 = 0.0;
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summerB_112_112_112 = 0.0, summerB_112_123_123 = 0.0, summerB_123_123_224 = 0.0, summerB_224_224_224 = 0.0;
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
-
- // Lorentz-Berthelot mixing rules
- auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
- auto sigmaij = (sigma[i]+sigma[j])/2;
-
- auto Tstari = T/epsilon_over_k[i], Tstarj = T/epsilon_over_k[j];
- auto Tstarij = forceeval(T/epskij);
-
- auto leading = x[i]*x[j]/pow(Tstari*Tstarj, 3.0/2.0); // common for all alpha_3A terms
-
- summerA_112_112_224 += leading*pow(sigma[i]*sigma[j], 11.0/2.0)/POW8(sigmaij)*mubar2[i]*mubar2[j]*sqrt(Qbar2[i]*Qbar2[j])*J11.get_J(Tstarij, rhostar);
- summerA_112_123_213 += leading*pow(sigma[i]*sigma[j], 11.0/2.0)/POW8(sigmaij)*mubar2[i]*mubar2[j]*sqrt(Qbar2[i]*Qbar2[j])*J11.get_J(Tstarij, rhostar);
- summerA_123_123_224 += leading*pow(sigma[i], 11.0/2.0)*pow(sigma[j], 15.0/2.0)/POW10(sigmaij)*mubar2[i]*sqrt(Qbar2[i])*pow(Qbar2[j], 3.0/2.0)*J13.get_J(Tstarij, rhostar);
- summerA_224_224_224 += leading*pow(sigma[i]*sigma[j], 15.0/2.0)/POW12(sigmaij)*Qbar2[i]*Qbar2[j]*J15.get_J(Tstarij, rhostar);
+ const VecType& x = mole_fractions; // concision
+ const std::size_t N = mole_fractions.size();
+ using type = std::common_type_t<TTYPE, RhoType, RhoStarType, decltype(mole_fractions[0])>;
+ using XTtype = std::common_type_t<TTYPE, decltype(mole_fractions[0])>;
+ type summerA_112_112_224 = 0.0, summerA_112_123_213 = 0.0, summerA_123_123_224 = 0.0, summerA_224_224_224 = 0.0;
+ type summerB_112_112_112 = 0.0, summerB_112_123_123 = 0.0, summerB_123_123_224 = 0.0, summerB_224_224_224 = 0.0;
+
+ for (std::size_t i = 0; i < N; ++i){
+ for (std::size_t j = 0; j < N; ++j){
+
+ TTYPE Tstari = forceeval(T/EPSKIJ(i,i)), Tstarj = forceeval(T/EPSKIJ(j,j));
+ TTYPE Tstarij = forceeval(T/EPSKIJ(i,j));
+
+ XTtype leading = forceeval(x[i]*x[j]/pow(forceeval(Tstari*Tstarj), 3.0/2.0)); // common for all alpha_3A terms
+ double sigmaij = SIGMAIJ(i,j);
+ double POW4sigmaij = powi(sigmaij, 4);
+ double POW8sigmaij = POW4sigmaij*POW4sigmaij;
+ double POW10sigmaij = powi(sigmaij, 10);
+ double POW12sigmaij = POW4sigmaij*POW8sigmaij;
- for (auto k = 0; k < N; ++k){
- auto Tstark = T/epsilon_over_k[k];
- auto epskik = sqrt(epsilon_over_k[i]*epsilon_over_k[k]);
- auto epskjk = sqrt(epsilon_over_k[j]*epsilon_over_k[k]);
- auto Tstarik = T/epskik;
- auto Tstarjk = T/epskjk;
-
+ {
+ double dbl = pow(sigma_m[i]*sigma_m[j], 11.0/2.0)/POW8sigmaij*mubar2[i]*mubar2[j]*sqrt(Qbar2[i]*Qbar2[j]);
+ summerA_112_112_224 += leading*dbl*J11.get_J(Tstarij, rhostar);
+ }
+ {
+ double dbl = pow(sigma_m[i]*sigma_m[j], 11.0/2.0)/POW8sigmaij*mubar2[i]*mubar2[j]*sqrt(Qbar2[i]*Qbar2[j]);
+ summerA_112_123_213 += leading*dbl*J11.get_J(Tstarij, rhostar);
+ }
+ {
+ double dbl = pow(sigma_m[i], 11.0/2.0)*pow(sigma_m[j], 15.0/2.0)/POW10sigmaij*mubar2[i]*sqrt(Qbar2[i])*pow(Qbar2[j], 3.0/2.0);
+ summerA_123_123_224 += leading*dbl*J13.get_J(Tstarij, rhostar);
+ }
+ {
+ double dbl = pow(sigma_m[i]*sigma_m[j], 15.0/2.0)/POW12sigmaij*Qbar2[i]*Qbar2[j];
+ summerA_224_224_224 += leading*dbl*J15.get_J(Tstarij, rhostar);
+ }
+
+ for (std::size_t k = 0; k < N; ++k){
+ TTYPE Tstark = forceeval(T/EPSKIJ(k,k));
+ TTYPE Tstarik = forceeval(T/EPSKIJ(i,k));
+ TTYPE Tstarjk = forceeval(T/EPSKIJ(j,k));
+ double sigmaik = SIGMAIJ(i,k), sigmajk = SIGMAIJ(j,k);
+
// Lorentz-Berthelot mixing rules for sigma
- auto sigmaij = (sigma[i]+sigma[j])/2;
- auto sigmaik = (sigma[i]+sigma[k])/2;
- auto sigmajk = (sigma[j]+sigma[k])/2;
- auto leadingijk = x[i]*x[j]*x[k]/(Tstari*Tstarj*Tstark);
-
- auto get_Kijk = [&](const auto& Kint){
- return forceeval(pow(Kint.get_K(Tstarij, rhostar)*Kint.get_K(Tstarik, rhostar)*Kint.get_K(Tstarjk, rhostar), 1.0/3.0));
- };
-
- // Special treatment needed here because the 334,445 term is negative, so negative*negative*negative is negative, and negative^{1/3} is undefined
- // First flip the sign on the triple product, do the evaluation, the flip it back. Not documented in Gubbins&Twu, but this seem reasonable,
- // in the spirit of the others.
- auto get_Kijk_334445 = [&](const auto& Kint){
- return forceeval(-pow(-Kint.get_K(Tstarij, rhostar)*Kint.get_K(Tstarik, rhostar)*Kint.get_K(Tstarjk, rhostar), 1.0/3.0));
- };
-
+ XTtype leadingijk = forceeval(x[i]*x[j]*x[k]/(Tstari*Tstarj*Tstark));
+
if (std::abs(mubar2[i]*mubar2[j]*mubar2[k]) > 0){
- auto K222333 = get_Kijk(K222_333);
- summerB_112_112_112 += leadingijk*POW3(sigma[i]*sigma[j]*sigma[k])/(sigmaij*sigmaik*sigmajk)*mubar2[i]*mubar2[j]*mubar2[k]*K222333;
+ auto K222333 = get_Kijk(K222_333, rhostar, Tstarij, Tstarik, Tstarjk);
+ double dbl = sigma_m3[i]*sigma_m3[j]*sigma_m3[k]/(sigmaij*sigmaik*sigmajk)*mubar2[i]*mubar2[j]*mubar2[k];
+ summerB_112_112_112 += forceeval(leadingijk*dbl*K222333);
}
if (std::abs(mubar2[i]*mubar2[j]*Qbar2[k]) > 0){
- auto K233344 = get_Kijk(K233_344);
- summerB_112_123_123 += leadingijk*POW3(sigma[i]*sigma[j])*POW5(sigma[k])/(sigmaij*POW2(sigmaik*sigmajk))*mubar2[i]*mubar2[j]*Qbar2[k]*K233344;
+ auto K233344 = get_Kijk(K233_344, rhostar, Tstarij, Tstarik, Tstarjk);
+ double dbl = sigma_m3[i]*sigma_m3[j]*sigma_m5[k]/(sigmaij*POW2(sigmaik*sigmajk))*mubar2[i]*mubar2[j]*Qbar2[k];
+ summerB_112_123_123 += leadingijk*dbl*K233344;
}
if (std::abs(mubar2[i]*Qbar2[j]*Qbar2[k]) > 0){
- auto K334445 = get_Kijk_334445(K334_445);
- summerB_123_123_224 += leadingijk*POW3(sigma[i])*POW5(sigma[j]*sigma[k])/(POW2(sigmaij*sigmaik)*POW3(sigmajk))*mubar2[i]*Qbar2[j]*Qbar2[k]*K334445;
+ auto K334445 = get_Kijk_334445(K334_445, rhostar, Tstarij, Tstarik, Tstarjk);
+ double dbl = sigma_m3[i]*sigma_m5[j]*sigma_m5[k]/(POW2(sigmaij*sigmaik)*POW3(sigmajk))*mubar2[i]*Qbar2[j]*Qbar2[k];
+ summerB_123_123_224 += leadingijk*dbl*K334445;
}
if (std::abs(Qbar2[i]*Qbar2[j]*Qbar2[k]) > 0){
- auto K444555 = get_Kijk(K444_555);
- summerB_224_224_224 += leadingijk*POW5(sigma[i]*sigma[j]*sigma[k])/(POW3(sigmaij*sigmaik*sigmajk))*Qbar2[i]*Qbar2[j]*Qbar2[k]*K444555;
+ auto K444555 = get_Kijk(K444_555, rhostar, Tstarij, Tstarik, Tstarjk);
+ double dbl = POW5(sigma_m[i]*sigma_m[j]*sigma_m[k])/(POW3(sigmaij*sigmaik*sigmajk))*Qbar2[i]*Qbar2[j]*Qbar2[k];
+ summerB_224_224_224 += leadingijk*dbl*K444555;
}
}
}
}
- auto alpha3A_112_112_224 = 8.0*PI_*rhoN/25.0*summerA_112_112_224;
- auto alpha3A_112_123_213 = 8.0*PI_*rhoN/75.0*summerA_112_123_213;
- auto alpha3A_123_123_224 = 8.0*PI_*rhoN/35.0*summerA_123_123_224;
- auto alpha3A_224_224_224 = 144.0*PI_*rhoN/245.0*summerA_224_224_224;
-
- auto alpha3A = forceeval(3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224);
-
- auto alpha3B_112_112_112 = 32.0*POW3(PI_)*POW2(rhoN)/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112;
- auto alpha3B_112_123_123 = 64.0*POW3(PI_)*POW2(rhoN)/315.0*sqrt(3.0*PI_)*summerB_112_123_123;
- auto alpha3B_123_123_224 = -32.0*POW3(PI_)*POW2(rhoN)/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224;
- auto alpha3B_224_224_224 = 32.0*POW3(PI_)*POW2(rhoN)/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224;
-
- auto alpha3B = forceeval(alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224);
+ type alpha3A_112_112_224 = 8.0*PI_*rhoN/25.0*summerA_112_112_224;
+ type alpha3A_112_123_213 = 8.0*PI_*rhoN/75.0*summerA_112_123_213;
+ type alpha3A_123_123_224 = 8.0*PI_*rhoN/35.0*summerA_123_123_224;
+ type alpha3A_224_224_224 = 144.0*PI_*rhoN/245.0*summerA_224_224_224;
+
+ type alpha3A = 3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224;
+
+ RhoType rhoN2 = rhoN*rhoN;
+
+ type alpha3B_112_112_112 = 32.0*POW3(PI_)*rhoN2/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112;
+ type alpha3B_112_123_123 = 64.0*POW3(PI_)*rhoN2/315.0*sqrt(3.0*PI_)*summerB_112_123_123;
+ type alpha3B_123_123_224 = -32.0*POW3(PI_)*rhoN2/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224;
+ type alpha3B_224_224_224 = 32.0*POW3(PI_)*rhoN2/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224;
+
+ type alpha3B = alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224;
+
return forceeval(alpha3A + alpha3B);
}
@@ -492,11 +534,14 @@ public:
*/
template<typename TTYPE, typename RhoType, typename VecType>
auto eval(const TTYPE& T, const RhoType& rhoN, const VecType& mole_fractions) const {
+ error_if_expr(T); error_if_expr(rhoN); error_if_expr(mole_fractions[0]);
using type = std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])>;
// Calculate the effective reduced diameter (cubed) to be used for evaluation
// Eq. 24 from Gubbins
type sigma_x3 = 0.0;
+
+ error_if_expr(sigma_m3);
auto N = mole_fractions.size();
for (auto i = 0; i < N; ++i){
for (auto j = 0; j < N; ++j){
@@ -512,17 +557,9 @@ public:
alpha3 = get_alpha3(T, rhoN, rhostar, mole_fractions);
alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
}
-
- using alpha2_t = decltype(alpha2);
- using alpha3_t = decltype(alpha3);
- using alpha_t = decltype(alpha);
- struct Terms{
- alpha2_t alpha2;
- alpha3_t alpha3;
- alpha_t alpha;
- };
- return Terms{alpha2, alpha3, alpha};
+ return MultipolarContributionGubbinsTwuTermsGT<type>{alpha2, alpha3, alpha};
}
+
};
/// The variant containing the multipolar types that can be provided
diff --git a/include/teqp/models/saftvrmie.hpp b/include/teqp/models/saftvrmie.hpp
index d69e8b5529945851b78fa29276a80bf080f1ee43..8f98a182960f398b3436036829fcbbb0ed65e252 100644
--- a/include/teqp/models/saftvrmie.hpp
+++ b/include/teqp/models/saftvrmie.hpp
@@ -12,8 +12,9 @@
#include "teqp/exceptions.hpp"
#include "teqp/constants.hpp"
#include "teqp/math/quadrature.hpp"
-#include "teqp/models/pcsaft.hpp"
+#include "teqp/models/saft/polar_terms.hpp"
#include <optional>
+#include <variant>
namespace teqp {
namespace SAFTVRMie {
@@ -207,19 +208,6 @@ struct SAFTVRMieChainContributionTerms{
return get_cij(lambda_r_ij + lambda_a_ij);
}
- template<typename X>
- auto POW2(const X& x) const{
- return forceeval(x*x);
- };
- template<typename X>
- auto POW3(const X& x) const{
- return forceeval(x*POW2(x));
- };
- template<typename X>
- auto POW4(const X& x) const{
- return forceeval(POW2(x)*POW2(x));
- };
-
public:
// One entry per component
@@ -847,27 +835,30 @@ public:
// First values for the Mie chain with dispersion (always included)
error_if_expr(T); error_if_expr(rhomolar);
auto vals = terms.get_core_calcs(T, rhomolar, mole_fractions);
- auto alphar = forceeval(vals.alphar_mono + vals.alphar_chain);
+ using type = std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])>;
+ type alphar = vals.alphar_mono + vals.alphar_chain;
+ type packing_fraction = vals.zeta[3];
- if (polar){ // polar term is present
- using mas = SAFTpolar::multipolar_argument_spec;
- auto visitor = [&](const auto& contrib){
- if constexpr(std::decay_t<decltype(contrib)>::arg_spec == mas::TK_rhoNA3_packingfraction_molefractions){
- auto rho_A3 = forceeval(rhomolar*N_A*1e-30);
- auto packing_fraction = vals.zeta[3];
- auto alpha = contrib.eval(T, rho_A3, packing_fraction, mole_fractions).alpha;
- alphar += alpha;
- }
- else if constexpr(std::decay_t<decltype(contrib)>::arg_spec == mas::TK_rhoNm3_molefractions){
- auto rhoN_m3 = forceeval(rhomolar*N_A);
- alphar += contrib.eval(T, rhoN_m3, mole_fractions).alpha;
- }
- else{
- throw teqp::InvalidArgument("Don't know how to handle this kind of arguments in polar term");
- }
- };
- std::visit(visitor, polar.value());
- }
+ if (polar){ // polar term is present
+ using mas = SAFTpolar::multipolar_argument_spec;
+ auto visitor = [&T, &rhomolar, &mole_fractions, &packing_fraction](const auto& contrib) -> type {
+
+ constexpr mas arg_spec = std::decay_t<decltype(contrib)>::arg_spec;
+ if constexpr(arg_spec == mas::TK_rhoNA3_packingfraction_molefractions){
+ RhoType rho_A3 = rhomolar*N_A*1e-30;
+ type alpha = contrib.eval(T, rho_A3, packing_fraction, mole_fractions).alpha;
+ return alpha;
+ }
+ else if constexpr(arg_spec == mas::TK_rhoNm3_molefractions){
+ RhoType rhoN_m3 = rhomolar*N_A;
+ return contrib.eval(T, rhoN_m3, mole_fractions).alpha;
+ }
+ else{
+ throw teqp::InvalidArgument("Don't know how to handle this kind of arguments in polar term");
+ }
+ };
+ alphar += std::visit(visitor, polar.value());
+ }
return forceeval(alphar);
}
@@ -988,23 +979,21 @@ inline auto SAFTVRMiefactory(const nlohmann::json & spec){
auto polar = MultipolarContributionGrossVrabec(ms, sigma_ms*1e10, epsks, mustar2, nmu, Qstar2, nQ);
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
}
- else if (polar_model == "GubbinsTwu+Luckas"){
+ if (polar_model == "GubbinsTwu+Luckas"){
using MCGTL = MultipolarContributionGubbinsTwu<LuckasJIntegral, LuckasKIntegral>;
auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval();
auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval();
auto polar = MCGTL(sigma_ms, epsks, mubar2, Qbar2);
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
}
- else if (polar_model == "GubbinsTwu+GubbinsTwu"){
+ if (polar_model == "GubbinsTwu+GubbinsTwu"){
using MCGG = MultipolarContributionGubbinsTwu<GubbinsTwuJIntegral, GubbinsTwuKIntegral>;
auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval();
auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval();
auto polar = MCGG(sigma_ms, epsks, mubar2, Qbar2);
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
}
- else{
- throw teqp::InvalidArgument("didn't understand this polar_model:"+polar_model);
- }
+ throw teqp::InvalidArgument("didn't understand this polar_model:"+polar_model);
}
}
else{
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index 4f66d0ddccd958e5e04d963c754f31a4e0e485d3..bcebf2193a70615156cec65ed36eae6c3a5e5e3b 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -15,7 +15,9 @@ namespace teqp {
// compilation speed hit invoked by these headers
#if !defined(NO_TYPES_HEADER)
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
#include "MultiComplex/MultiComplex.hpp"
+#endif
#include <valarray>
#include <chrono>
@@ -72,7 +74,9 @@ namespace teqp {
// See https://stackoverflow.com/a/41438758
template<typename T> struct is_mcx_t : public std::false_type {};
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
template<typename T> struct is_mcx_t<mcx::MultiComplex<T>> : public std::true_type {};
+#endif
// Extract the underlying value from more complicated numerical types, like complex step types with
// a tiny increment in the imaginary direction
@@ -198,7 +202,7 @@ namespace teqp {
}
return o;
}
-
+#if defined(TEQP_MULTICOMPLEX_ENABLED)
template<typename T>
auto pow(const mcx::MultiComplex<T>& x, const Eigen::ArrayXd& e) {
Eigen::Array<mcx::MultiComplex<T>, Eigen::Dynamic, 1> o(e.size());
@@ -207,6 +211,7 @@ namespace teqp {
}
return o;
}
+#endif
}; // namespace teqp
diff --git a/interface/CPP/AbstractModel_implementation.cpp b/interface/CPP/AbstractModel_implementation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f97939bc3a962fe35cb5f3e486e196ea7af9406
--- /dev/null
+++ b/interface/CPP/AbstractModel_implementation.cpp
@@ -0,0 +1,94 @@
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/algorithms/critical_pure.hpp"
+#include "teqp/algorithms/VLE_pure.hpp"
+#include "teqp/algorithms/VLE.hpp"
+#include "teqp/algorithms/VLLE.hpp"
+
+namespace teqp{
+ namespace cppinterface{
+
+ double AbstractModel::get_neff(const double T, const double rho, const EArrayd& molefracs) const {
+ return -3.0*(this->get_Ar01(T, rho, molefracs) - this->get_Ar11(T, rho, molefracs) )/this->get_Ar20(T,rho,molefracs);
+ };
+
+ std::tuple<double, double> AbstractModel::solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& flags) const {
+ return teqp::solve_pure_critical(*this, T, rho, flags.value_or(nlohmann::json{}));
+ }
+ std::tuple<EArrayd, EMatrixd> AbstractModel::get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const {
+ return teqp::get_pure_critical_conditions_Jacobian(*this, T, rho, alternative_pure_index, alternative_length);
+ }
+ EArray2 AbstractModel::extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const {
+ return teqp::extrapolate_from_critical(*this, Tc, rhoc, Tnew);
+ }
+
+ EArray2 AbstractModel::pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const {
+ return teqp::pure_VLE_T(*this, T, rhoL, rhoV, maxiter);
+ }
+
+ double AbstractModel::dpsatdT_pure(const double T, const double rhoL, const double rhoV) const {
+ return teqp::dpsatdT_pure(*this, T, rhoL, rhoV);
+ }
+
+ std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> AbstractModel::mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const{
+
+ return VLLE::mix_VLLE_T(*this, T, rhovecVinit, rhovecL1init, rhovecL2init, atol, reltol, axtol, relxtol, maxiter);
+ }
+
+ std::vector<nlohmann::json> AbstractModel::find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options) const{
+ return VLLE::find_VLLE_T_binary(*this, traces, options);;
+ }
+
+ std::tuple<VLE_return_code,EArrayd,EArrayd> AbstractModel::mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const{
+ return teqp::mix_VLE_Tx(*this, T, rhovecL0, rhovecV0, xspec, atol, reltol, axtol, relxtol, maxiter);
+
+ }
+ MixVLEReturn AbstractModel::mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const{
+ return teqp::mix_VLE_Tp(*this, T, pgiven, rhovecL0, rhovecV0, flags);
+ }
+ std::tuple<VLE_return_code,double,EArrayd,EArrayd> AbstractModel::mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags) const{
+ return teqp::mixture_VLE_px(*this, p_spec, xmolar_spec, T0, rhovecL0, rhovecV0, flags);
+ }
+
+ std::tuple<EArrayd, EArrayd> AbstractModel::get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
+ return teqp::get_drhovecdp_Tsat(*this, T, rhovecL, rhovecV);
+ }
+ std::tuple<EArrayd, EArrayd> AbstractModel::get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
+ return teqp::get_drhovecdT_psat(*this, T, rhovecL, rhovecV);
+ }
+ double AbstractModel::get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
+ return teqp::get_dpsat_dTsat_isopleth(*this, T, rhovecL, rhovecV);
+ }
+ nlohmann::json AbstractModel::trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const{
+ return teqp::trace_VLE_isotherm_binary(*this, T0, rhovecL0, rhovecV0, options);
+ }
+ nlohmann::json AbstractModel::trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const{
+ return teqp::trace_VLE_isobar_binary(*this, p, T0, rhovecL0, rhovecV0, options);
+ }
+
+ nlohmann::json AbstractModel::trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename, const std::optional<TCABOptions> &options) const {
+ using crit = teqp::CriticalTracing<decltype(*this), double, std::decay_t<decltype(rhovec0)>>;
+ return crit::trace_critical_arclength_binary(*this, T0, rhovec0, filename , options);
+ }
+ EArrayd AbstractModel::get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const {
+ using crit = teqp::CriticalTracing<decltype(*this), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_drhovec_dT_crit(*this, T, rhovec);
+ }
+ double AbstractModel::get_dp_dT_crit(const double T, const REArrayd& rhovec) const {
+ using crit = teqp::CriticalTracing<decltype(*this), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_dp_dT_crit(*this, T, rhovec);
+ }
+ EArray2 AbstractModel::get_criticality_conditions(const double T, const REArrayd& rhovec) const {
+ using crit = teqp::CriticalTracing<decltype(*this), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_criticality_conditions(*this, T, rhovec);
+ }
+ EigenData AbstractModel::eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const {
+ using crit = teqp::CriticalTracing<decltype(*this), double, std::decay_t<decltype(rhovec)>>;
+ return crit::eigen_problem(*this, T, rhovec, alignment_v0.value_or(Eigen::ArrayXd()));
+ }
+ double AbstractModel::get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const {
+ using crit = teqp::CriticalTracing<decltype(*this), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_minimum_eigenvalue_Psi_Hessian(*this, T, rhovec);
+ }
+
+ }
+}
diff --git a/interface/CPP/model_factories.cpp b/interface/CPP/model_factories.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..edcfa57068aefedf55074ca58b0cfae5a5ee49fe
--- /dev/null
+++ b/interface/CPP/model_factories.cpp
@@ -0,0 +1,11 @@
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
+#include "teqp/models/saftvrmie.hpp"
+
+namespace teqp{
+ namespace cppinterface{
+ std::unique_ptr<teqp::cppinterface::AbstractModel> make_SAFTVRMie(const nlohmann::json &j){
+ return teqp::cppinterface::adapter::make_owned(SAFTVRMie::SAFTVRMiefactory(j));
+ }
+ }
+}
diff --git a/interface/CPP/teqp_impl_Arxy.cpp b/interface/CPP/teqp_impl_Arxy.cpp
deleted file mode 100644
index 2326d6a6527eb1ce423aa8d4ed579edf8a59d966..0000000000000000000000000000000000000000
--- a/interface/CPP/teqp_impl_Arxy.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-#include "teqp/derivs.hpp"
-
-#include "teqpcpp.cpp"
-using MI = teqp::cppinterface::ModelImplementer;
-
-// Derivatives from isochoric thermodynamics (all have the same signature)
-#define X(f) \
-double MI::f(const double T, const EArrayd& rhovec) const { \
- return std::visit([&](const auto& model) { \
- using id = IsochoricDerivatives<decltype(model), double, EArrayd>; \
- return id::f(model, T, rhovec); \
- }, m_model); \
-}
-ISOCHORIC_double_args
-#undef X
-
-#define X(f) \
-EArrayd MI::f(const double T, const EArrayd& rhovec) const { \
- return std::visit([&](const auto& model) { \
- using id = IsochoricDerivatives<decltype(model), double, EArrayd>; \
- return id::f(model, T, rhovec); \
- }, m_model); \
-}
-ISOCHORIC_array_args
-#undef X
-
-#define X(f) \
-EMatrixd MI::f(const double T, const EArrayd& rhovec) const { \
- return std::visit([&](const auto& model) { \
- using id = IsochoricDerivatives<decltype(model), double, EArrayd>; \
- return id::f(model, T, rhovec); \
- }, m_model); \
-}
-ISOCHORIC_matrix_args
-#undef X
-
-double MI::get_Arxy(const int NT, const int ND, const double T, const double rho, const EArrayd& molefracs) const {
- return std::visit([&](const auto& model) {
- using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
- return tdx::template get_Ar(NT, ND, model, T, rho, molefracs);
- }, m_model);
-}
-
-double MI::get_neff(const double T, const double rho, const EArrayd& molefracs) const {
- return std::visit([&](const auto& model) {
- using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
- return tdx::template get_neff(model, T, rho, molefracs);
- }, m_model);
-}
-
-EArray33d MI::get_deriv_mat2(const double T, double rho, const EArrayd& z) const {
- return std::visit([&](const auto& model) {
- // Although the template argument suggests that only residual terms
- // are returned, also the ideal-gas ones are returned because the
- // ideal-gas term is required to implement alphar which just redirects
- // to alphaig
- return DerivativeHolderSquare<2, AlphaWrapperOption::residual>(model, T, rho, z).derivs;
- }, m_model);
-}
diff --git a/interface/CPP/teqp_impl_VLE.cpp b/interface/CPP/teqp_impl_VLE.cpp
deleted file mode 100644
index 91be45db837653c24104a13b6c0af5adf376606e..0000000000000000000000000000000000000000
--- a/interface/CPP/teqp_impl_VLE.cpp
+++ /dev/null
@@ -1,58 +0,0 @@
-#include "teqpcpp.cpp"
-#include "teqp/algorithms/VLE.hpp"
-
-using namespace teqp;
-using MI = teqp::cppinterface::ModelImplementer;
-
-std::tuple<EArrayd, EArrayd> MI::get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
- return std::visit([&](const auto& model) {
- return teqp::get_drhovecdp_Tsat(model, T, rhovecL, rhovecV);
- }, m_model);
-}
-std::tuple<EArrayd, EArrayd> MI::get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
- return std::visit([&](const auto& model) {
- return teqp::get_drhovecdT_psat(model, T, rhovecL, rhovecV);
- }, m_model);
-}
-double MI::get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
- return std::visit([&](const auto& model) {
- return teqp::get_dpsat_dTsat_isopleth(model, T, rhovecL, rhovecV);
- }, m_model);
-}
-
-nlohmann::json MI::trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const{
- return std::visit([&](const auto& model) {
- return teqp::trace_VLE_isotherm_binary(model, T0, rhovecL0, rhovecV0, options);
- }, m_model);
-}
-nlohmann::json MI::trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const{
- return std::visit([&](const auto& model) {
- return teqp::trace_VLE_isobar_binary(model, p, T0, rhovecL0, rhovecV0, options);
- }, m_model);
-}
-std::tuple<VLE_return_code,EArrayd,EArrayd> MI::mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const{
- return std::visit([&](const auto& model) {
- return teqp::mix_VLE_Tx(model, T, rhovecL0, rhovecV0, xspec, atol, reltol, axtol, relxtol, maxiter);
- }, m_model);
-}
-MixVLEReturn MI::mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const{
- return std::visit([&](const auto& model) {
- return teqp::mix_VLE_Tp(model, T, pgiven, rhovecL0, rhovecV0, flags);
- }, m_model);
-}
-std::tuple<VLE_return_code,double,EArrayd,EArrayd> MI::mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags ) const{
- return std::visit([&](const auto& model) {
- return teqp::mixture_VLE_px(model, p_spec, xmolar_spec, T0, rhovecL0, rhovecV0, flags);
- }, m_model);
-}
-EArray2 MI::pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const {
- return std::visit([&](const auto& model) {
- return teqp::pure_VLE_T(model, T, rhoL, rhoV, maxiter);
- }, m_model);
-}
-
-double MI::dpsatdT_pure(const double T, const double rhoL, const double rhoV) const {
- return std::visit([&](const auto& model) {
- return teqp::dpsatdT_pure(model, T, rhoL, rhoV);
- }, m_model);
-}
diff --git a/interface/CPP/teqp_impl_VLLE.cpp b/interface/CPP/teqp_impl_VLLE.cpp
deleted file mode 100644
index eae56964194423ca62aeb5f955378a341e62e59b..0000000000000000000000000000000000000000
--- a/interface/CPP/teqp_impl_VLLE.cpp
+++ /dev/null
@@ -1,17 +0,0 @@
-#include "teqpcpp.cpp"
-#include "teqp/algorithms/VLLE.hpp"
-
-using namespace teqp;
-using MI = teqp::cppinterface::ModelImplementer;
-
-std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> MI::mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const{
- return std::visit([&](const auto& model) {
- return VLLE::mix_VLLE_T(model, T, rhovecVinit, rhovecL1init, rhovecL2init, atol, reltol, axtol, relxtol, maxiter);
- }, m_model);
-}
-
-std::vector<nlohmann::json> MI::find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options) const{
- return std::visit([&](const auto& model) {
- return VLLE::find_VLLE_T_binary(model, traces, options);
- }, m_model);
-}
diff --git a/interface/CPP/teqp_impl_crit.cpp b/interface/CPP/teqp_impl_crit.cpp
deleted file mode 100644
index 8f43d964d02e1e1668210f19e3cdd8c7209f66ec..0000000000000000000000000000000000000000
--- a/interface/CPP/teqp_impl_crit.cpp
+++ /dev/null
@@ -1,61 +0,0 @@
-#include "teqp/algorithms/critical_tracing.hpp"
-
-#include "teqpcpp.cpp"
-
-using namespace teqp;
-using MI = teqp::cppinterface::ModelImplementer;
-
-nlohmann::json MI::trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename, const std::optional<TCABOptions> &options) const {
- 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, filename , options);
- }, m_model);
-}
-EArrayd MI::get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const {
- return std::visit([&](const auto& model) {
- using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
- return crit::get_drhovec_dT_crit(model, T, rhovec);
- }, m_model);
-}
-double MI::get_dp_dT_crit(const double T, const REArrayd& rhovec) const {
- return std::visit([&](const auto& model) {
- using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
- return crit::get_dp_dT_crit(model, T, rhovec);
- }, m_model);
-}
-EArray2 MI::get_criticality_conditions(const double T, const REArrayd& rhovec) const {
- return std::visit([&](const auto& model) {
- using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
- return crit::get_criticality_conditions(model, T, rhovec);
- }, m_model);
-}
-EigenData MI::eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const {
- return std::visit([&](const auto& model) {
- using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
- return crit::eigen_problem(model, T, rhovec, alignment_v0.value_or(Eigen::ArrayXd()));
- }, m_model);
-}
-double MI::get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const {
- return std::visit([&](const auto& model) {
- using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
- return crit::get_minimum_eigenvalue_Psi_Hessian(model, T, rhovec);
- }, m_model);
-}
-
-
-std::tuple<double, double> MI::solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& flags) const {
- return std::visit([&](const auto& model) {
- return teqp::solve_pure_critical(model, T, rho, flags.value_or(nlohmann::json{}));
- }, m_model);
-}
-std::tuple<EArrayd, EMatrixd> MI::get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const {
- return std::visit([&](const auto& model) {
- return teqp::get_pure_critical_conditions_Jacobian(model, T, rho, alternative_pure_index, alternative_length);
- }, m_model);
-}
-std::tuple<double, double> MI::extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const {
- return std::visit([&](const auto& model) {
- auto mat = teqp::extrapolate_from_critical(model, Tc, rhoc, Tnew);
- return std::make_tuple(mat[0], mat[1]);
- }, m_model);
-}
diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp
index ea2de8e4bed0ce713927eb557753c2ddd0dd0135..3e6b7e32fb5275029e6384e4783c8e0c2c03af23 100644
--- a/interface/CPP/teqp_impl_factory.cpp
+++ b/interface/CPP/teqp_impl_factory.cpp
@@ -1,15 +1,61 @@
-#include "teqpcpp.cpp"
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/json_builder.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
-using namespace teqp;
-using namespace teqp::cppinterface;
+namespace teqp {
+ namespace cppinterface {
-std::shared_ptr<AbstractModel> teqp::cppinterface::make_model(const nlohmann::json& j) {
- return std::make_shared<ModelImplementer>(build_model(j));
-}
-std::shared_ptr<AbstractModel> teqp::cppinterface::make_multifluid_model(const std::vector<std::string>& components, const std::string& coolprop_root, const std::string& BIPcollectionpath, const nlohmann::json& flags, const std::string& departurepath) {
- return std::make_shared<ModelImplementer>(build_multifluid_model(components, coolprop_root, BIPcollectionpath, flags, departurepath));
-}
+ std::unique_ptr<teqp::cppinterface::AbstractModel> make_SAFTVRMie(const nlohmann::json &j);
+
+ using makefunc = std::function<std::unique_ptr<teqp::cppinterface::AbstractModel>(const nlohmann::json &j)>;
+ using namespace teqp::cppinterface::adapter;
+
+ static std::unordered_map<std::string, makefunc> pointer_factory = {
+ {"vdW1", [](const nlohmann::json& spec){ return make_owned(vdWEOS1(spec.at("a"), spec.at("b"))); }},
+ {"vdW", [](const nlohmann::json& spec){ return make_owned(vdWEOS<double>(spec.at("Tcrit / K"), spec.at("pcrit / Pa"))); }},
+
+ {"PR", [](const nlohmann::json& spec){ return make_owned(make_canonicalPR(spec));}},
+ {"SRK", [](const nlohmann::json& spec){ return make_owned(make_canonicalSRK(spec));}},
+ {"CPA", [](const nlohmann::json& spec){ return make_owned(CPA::CPAfactory(spec));}},
+ {"PCSAFT", [](const nlohmann::json& spec){ return make_owned(PCSAFT::PCSAFTfactory(spec));}},
+
+ {"multifluid", [](const nlohmann::json& spec){ return make_owned(multifluidfactory(spec));}},
+ {"SW_EspindolaHeredia2009", [](const nlohmann::json& spec){ return make_owned(squarewell::EspindolaHeredia2009(spec.at("lambda")));}},
+ {"EXP6_Kataoka1992", [](const nlohmann::json& spec){ return make_owned(exp6::Kataoka1992(spec.at("alpha")));}},
+ {"AmmoniaWaterTillnerRoth", [](const nlohmann::json& spec){ return make_owned(AmmoniaWaterTillnerRoth());}},
+ {"LJ126_TholJPCRD2016", [](const nlohmann::json& spec){ return make_owned(build_LJ126_TholJPCRD2016());}},
+ {"LJ126_KolafaNezbeda1994", [](const nlohmann::json& spec){ return make_owned(LJ126KolafaNezbeda1994());}},
+ {"LJ126_Johnson1993", [](const nlohmann::json& spec){ return make_owned(LJ126Johnson1993());}},
+ {"Mie_Pohl2023", [](const nlohmann::json& spec){ return make_owned(Mie::Mie6Pohl2023(spec.at("lambda_a")));}},
+ {"2CLJF-Dipole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_dipole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));}},
+ {"2CLJF-Quadrupole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_quadrupole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));}},
+ {"IdealHelmholtz", [](const nlohmann::json& spec){ return make_owned(IdealHelmholtz(spec));}},
+
+ // Implemented in its own compilation unit to help with compilation time
+ {"SAFT-VR-Mie", [](const nlohmann::json& spec){ return make_SAFTVRMie(spec); }}
+ };
-std::shared_ptr<AbstractModel> teqp::cppinterface::emplace_model(AllowedModels&& model){
- return std::make_shared<ModelImplementer>(std::move(model));
+ std::unique_ptr<teqp::cppinterface::AbstractModel> build_model_ptr(const nlohmann::json& json) {
+
+ // Extract the name of the model and the model parameters
+ std::string kind = json.at("kind");
+ auto spec = json.at("model");
+
+ auto itr = pointer_factory.find(kind);
+ if (itr != pointer_factory.end()){
+ return (itr->second)(spec);
+ }
+ else{
+ throw std::invalid_argument("Don't understand \"kind\" of: " + kind);
+ }
+ }
+
+ std::unique_ptr<AbstractModel> make_multifluid_model(const std::vector<std::string>& components, const std::string& coolprop_root, const std::string& BIPcollectionpath, const nlohmann::json& flags, const std::string& departurepath) {
+ return make_owned(build_multifluid_model(components, coolprop_root, BIPcollectionpath, flags, departurepath));
+ }
+
+ std::unique_ptr<AbstractModel> make_model(const nlohmann::json& j) {
+ return build_model_ptr(j);
+ }
+ }
}
diff --git a/interface/CPP/teqp_impl_isochoric.cpp b/interface/CPP/teqp_impl_isochoric.cpp
deleted file mode 100644
index e3d7fa9982e3bfba20e10e54b233cc0379123070..0000000000000000000000000000000000000000
--- a/interface/CPP/teqp_impl_isochoric.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-#include "teqp/derivs.hpp"
-
-#include "teqpcpp.cpp"
-using MI = teqp::cppinterface::ModelImplementer;
-
-// Here XMacros are used to create functions like get_Ar00, get_Ar01, ....
-#define X(i,j) \
-double MI::get_Ar ## i ## j(const double T, const double rho, const REArrayd& molefracs) const { \
- return std::visit([&](const auto& model) { \
- using tdx = teqp::TDXDerivatives<decltype(model), double, REArrayd>; \
- return tdx::template get_Arxy<i,j>(model, T, rho, molefracs); \
- }, m_model); \
-}
-ARXY_args
-#undef X
-
-// Here XMacros are used to create functions like get_Ar01n, get_Ar02n, ....
-#define X(i) \
-EArrayd MI::get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefracs) const { \
- return std::visit([&](const auto& model) -> EArrayd { \
- using tdx = teqp::TDXDerivatives<decltype(model), double, REArrayd>; \
- auto vals = tdx::template get_Ar0n<i>(model, T, rho, molefracs); \
- return Eigen::Map<Eigen::ArrayXd>(&(vals[0]), vals.size()); \
- }, m_model); \
-}
-AR0N_args
-#undef X
diff --git a/interface/CPP/teqp_impl_virials.cpp b/interface/CPP/teqp_impl_virials.cpp
deleted file mode 100644
index 35aaab070085da4bb98104baad6ff07646509141..0000000000000000000000000000000000000000
--- a/interface/CPP/teqp_impl_virials.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#include "teqp/derivs.hpp"
-#include "teqpcpp.cpp"
-
-using MI = teqp::cppinterface::ModelImplementer;
-
-double MI::get_B2vir(const double T, const EArrayd& molefrac) const {
- return std::visit([&](const auto& model) {
- using vd = VirialDerivatives<decltype(model), double, RAX>;
- return vd::get_B2vir(model, T, molefrac);
- }, m_model);
-}
-double MI::get_B12vir(const double T, const EArrayd& molefrac) const {
- return std::visit([&](const auto& model) {
- using vd = VirialDerivatives<decltype(model), double, RAX>;
- return vd::get_B12vir(model, T, molefrac);
- }, m_model);
-}
-std::map<int, double> MI::get_Bnvir(const int Nderiv, const double T, const EArrayd& molefrac) const {
- return std::visit([&](const auto& model) {
- using vd = VirialDerivatives<decltype(model), double, RAX>;
- return vd::get_Bnvir_runtime(Nderiv, model, T, molefrac);
- }, m_model);
-}
-double MI::get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const {
- return std::visit([&](const auto& model) {
- using vd = VirialDerivatives<decltype(model), double, RAX>;
- return vd::get_dmBnvirdTm_runtime(Nderiv, NTderiv, model, T, molefrac);
- }, m_model);
-}
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
deleted file mode 100644
index 8bcd676637595c5c889d7dd19a164c03446538ad..0000000000000000000000000000000000000000
--- a/interface/CPP/teqpcpp.cpp
+++ /dev/null
@@ -1,105 +0,0 @@
-#include "teqp/cpp/teqpcpp.hpp"
-#include "teqp/json_builder.hpp"
-
-namespace teqp {
- namespace cppinterface {
-
- class ModelImplementer : public AbstractModel {
- private:
- template<typename cls>
- const cls& get_or_fail(const std::string& typestr) const{
- if (std::holds_alternative<cls>(m_model)){
- return std::get<cls>(m_model);
- }
- else{
- throw std::invalid_argument("This method is only available for models of the type " + std::string(typestr));
- }
- }
- protected:
- AllowedModels m_model;
- using RAX = Eigen::Ref<const Eigen::ArrayXd>;
-
- public:
- ModelImplementer(AllowedModels&& model) : m_model(model) {};
-
- const AllowedModels& get_model() const override{
- return m_model;
- }
- AllowedModels& get_mutable_model() override{
- return m_model;
- }
-
- double get_R(const EArrayd& molefracs) const override {
- return std::visit([&](const auto& model) {
- return model.R(molefracs);
- }, m_model);
- }
-
- nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename, const std::optional<TCABOptions> &options) const override;
- EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const override;
- double get_dp_dT_crit(const double T, const REArrayd& rhovec) const override;
- EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const override;
- EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const override;
- double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const override ;
-
- std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& flags) const override ;
- std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const override ;
- std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const override ;
-
- // Derivatives from isochoric thermodynamics (all have the same signature)
- #define X(f) \
- double f(const double T, const EArrayd& rhovec) const override ;
- ISOCHORIC_double_args
- #undef X
-
- #define X(f) \
- EArrayd f(const double T, const EArrayd& rhovec) const override ;
- ISOCHORIC_array_args
- #undef X
-
- #define X(f) \
- EMatrixd f(const double T, const EArrayd& rhovec) const override ;
- ISOCHORIC_matrix_args
- #undef X
-
- EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z) const override ;
-
- double get_B2vir(const double T, const EArrayd& molefrac) const override;
- double get_B12vir(const double T, const EArrayd& molefrac) const override;
- std::map<int, double> get_Bnvir(const int Nderiv, const double T, const EArrayd& molefrac) const override;
- double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const override;
-
- double get_Arxy(const int NT, const int ND, const double T, const double rho, const EArrayd& molefracs) const override;
- // Here XMacros are used to create functions like get_Ar00, get_Ar01, ....
- #define X(i,j) \
- double get_Ar ## i ## j(const double T, const double rho, const REArrayd& molefracs) const override;
- ARXY_args
- #undef X
- // Here XMacros are used to create functions like get_Ar01n, get_Ar02n, ....
- #define X(i) \
- EArrayd get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefracs) const override;
- AR0N_args
- #undef X
-
- double get_neff(const double T, const double rho, const EArrayd& molefracs) const override;
-
- EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const override;
- double dpsatdT_pure(const double T, const double rhoL, const double rhoV) const override;
- std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
- std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
- double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
-
- nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const override;
- nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const override;
- std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override;
- MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const override;
- std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags ) const override;
-
- std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override;
-
- std::vector<nlohmann::json> find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options) const override;
- };
-
-
- }
-}
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 668ce120a345c5549a9c208cc8abbeb548e4a737..d0d9450d15b743e1e9292b71935f1e031f8be920 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -12,6 +12,7 @@
#include "teqp/cpp/teqpcpp.hpp"
#include "teqp/models/multifluid_ancillaries.hpp"
#include "teqp/algorithms/iteration.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
namespace py = pybind11;
using namespace py::literals;
@@ -45,33 +46,41 @@ void add_multifluid(py::module& m){
m.def("get_departure_json", &get_departure_json, py::arg("name"), py::arg("root"));
}
+template<typename TYPE>
+const TYPE& get_typed(const py::object& o){
+ using namespace teqp::cppinterface;
+ using namespace teqp::cppinterface::adapter;
+ // Cast Python-wrapped AbstractModel to the C++ AbstractModel
+ const AbstractModel* am = o.cast<const AbstractModel *>();
+ // Cast the C++ AbstractModel to the derived adapter class
+ return get_model_cref<TYPE>(am);
+}
+
+template<typename TYPE>
+TYPE& get_mutable_typed(py::object& o){
+ using namespace teqp::cppinterface;
+ using namespace teqp::cppinterface::adapter;
+ // Cast Python-wrapped AbstractModel to the C++ AbstractModel
+ AbstractModel* am = o.cast<AbstractModel *>();
+ // Cast the C++ AbstractModel to the derived adapter class
+ return get_model_ref<TYPE>(am);
+}
+
void add_multifluid_mutant(py::module& m) {
using namespace teqp;
using namespace teqp::cppinterface;
// A typedef for the base model
using MultiFluid = decltype(build_multifluid_model(std::vector<std::string>{"", ""}, "", ""));
-
+
// Wrap the function for generating a multifluid mutant
m.def("_build_multifluid_mutant", [](const py::object& o, const nlohmann::json &j){
- const auto& model = std::get<MultiFluid>(o.cast<const AbstractModel *>()->get_model());
- AllowedModels mutant{build_multifluid_mutant(model, j)};
- return emplace_model(std::move(mutant));
+ const MultiFluid& model = get_typed<MultiFluid>(o);
+ auto mutant{build_multifluid_mutant(model, j)};
+ return teqp::cppinterface::adapter::make_owned(mutant);
});
}
-template<typename TYPE>
-const TYPE& get_typed(py::object& o){
- using namespace teqp::cppinterface;
- return std::get<TYPE>(o.cast<const AbstractModel *>()->get_model());
-}
-
-template<typename TYPE>
-TYPE& get_mutable_typed(py::object& o){
- using namespace teqp::cppinterface;
- return std::get<TYPE>(o.cast<AbstractModel *>()->get_mutable_model());
-}
-
template<typename TYPE>
void attach_multifluid_methods(py::object&obj){
auto setattr = py::getattr(obj, "__setattr__");
@@ -85,26 +94,52 @@ void attach_multifluid_methods(py::object&obj){
setattr("get_alpharij", MethodType(py::cpp_function([](py::object& o, const int i, const int j, const double tau, const double delta){ return get_typed<TYPE>(o).dep.get_alpharij(i,j,tau,delta); }, "self"_a, "i"_a, "j"_a, "tau"_a, "delta"_a), obj));
}
+// Type index variables matching the model types, used for runtime attachment of model-specific methods
+const std::type_index vdWEOS1_i{std::type_index(typeid(vdWEOS1))};
+const std::type_index PCSAFT_i{std::type_index(typeid(PCSAFT_t))};
+const std::type_index SAFTVRMie_i{std::type_index(typeid(SAFTVRMie_t))};
+const std::type_index canonical_cubic_i{std::type_index(typeid(canonical_cubic_t))};
+const std::type_index AmmoniaWaterTillnerRoth_i{std::type_index(typeid(AmmoniaWaterTillnerRoth))};
+const std::type_index idealgas_i{std::type_index(typeid(idealgas_t))};
+const std::type_index multifluid_i{std::type_index(typeid(multifluid_t))};
+const std::type_index multifluidmutant_i{std::type_index(typeid(multifluidmutant_t))};
+const std::type_index SW_EspindolaHeredia2009_i{std::type_index(typeid(SW_EspindolaHeredia2009_t))};
+const std::type_index EXP6_Kataoka1992_i{std::type_index(typeid(EXP6_Kataoka1992_t))};
+const std::type_index twocenterLJF_i{std::type_index(typeid(twocenterLJF_t))};
+
+/**
+ At runtime we can add additional model-specific methods that only apply for a particular model. We take in a Python-wrapped
+ object and use runtime instrospection to figure out what kind of model it is. Then, we attach methods that are evaluated
+ ad *runtime* to call methods of the instance. This is tricky, although it works just fine.
+
+ */
// You cannot know at runtime what is contained in the model so you must iterate
// over possible model types and attach methods accordingly
void attach_model_specific_methods(py::object& obj){
using namespace teqp::cppinterface;
- const auto& model = obj.cast<AbstractModel *>()->get_model();
- auto setattr = py::getattr(obj, "__setattr__");
+
+ // Get things from Python (these are just for convenience to save typing)
auto MethodType = py::module_::import("types").attr("MethodType");
+ auto setattr = py::getattr(obj, "__setattr__");
+
+ AbstractModel* am = obj.cast<AbstractModel *>();
+ if (am == nullptr){
+ throw std::invalid_argument("Bad cast of argument to C++ AbstractModel type");
+ }
+ const auto& index = am->get_type_index();
- if (std::holds_alternative<vdWEOS1>(model)){
+ if (index == vdWEOS1_i){
setattr("get_a", MethodType(py::cpp_function([](py::object& o){ return get_typed<vdWEOS1>(o).get_a(); }), obj));
setattr("get_b", MethodType(py::cpp_function([](py::object& o){ return get_typed<vdWEOS1>(o).get_b(); }), obj));
}
- else if (std::holds_alternative<PCSAFT_t>(model)){
+ else if (index == PCSAFT_i){
setattr("get_m", MethodType(py::cpp_function([](py::object& o){ return get_typed<PCSAFT_t>(o).get_m(); }), obj));
setattr("get_sigma_Angstrom", MethodType(py::cpp_function([](py::object& o){ return get_typed<PCSAFT_t>(o).get_sigma_Angstrom(); }), obj));
setattr("get_epsilon_over_k_K", MethodType(py::cpp_function([](py::object& o){ return get_typed<PCSAFT_t>(o).get_epsilon_over_k_K(); }), obj));
setattr("max_rhoN", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<PCSAFT_t>(o).max_rhoN(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<PCSAFT_t>(o).get_kmat(); }), obj));
}
- else if (std::holds_alternative<SAFTVRMie_t>(model)){
+ else if (index == SAFTVRMie_i){
setattr("get_m", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_m(); }), obj));
setattr("get_sigma_Angstrom", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_sigma_Angstrom(); }), obj));
setattr("get_sigma_m", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_sigma_m(); }), obj));
@@ -115,16 +150,15 @@ void attach_model_specific_methods(py::object& obj){
setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_kmat(); }), obj));
setattr("get_core_calcs", MethodType(py::cpp_function([](py::object& o, double T, double rhomolar, REArrayd& molefrac){ return get_typed<SAFTVRMie_t>(o).get_core_calcs(T, rhomolar, molefrac); }, "self"_a, "T"_a, "rhomolar"_a, "molefrac"_a), obj));
}
- else if (std::holds_alternative<canonical_cubic_t>(model)){
+ else if (index == canonical_cubic_i){
setattr("get_a", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<canonical_cubic_t>(o).get_a(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
setattr("get_b", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<canonical_cubic_t>(o).get_b(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
setattr("superanc_rhoLV", MethodType(py::cpp_function([](py::object& o, double T){ return get_typed<canonical_cubic_t>(o).superanc_rhoLV(T); }, "self"_a, "T"_a), obj));
setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<canonical_cubic_t>(o).get_kmat(); }), obj));
setattr("get_meta", MethodType(py::cpp_function([](py::object& o){ return get_typed<canonical_cubic_t>(o).get_meta(); }), obj));
setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed<canonical_cubic_t>(o).set_meta(s); }, "self"_a, "s"_a), obj));
-
}
- else if (std::holds_alternative<AmmoniaWaterTillnerRoth>(model)){
+ else if (index == AmmoniaWaterTillnerRoth_i){
setattr("TcNH3", get_typed<AmmoniaWaterTillnerRoth>(obj).TcNH3);
setattr("vcNH3", get_typed<AmmoniaWaterTillnerRoth>(obj).vcNH3);
setattr("get_Tr", MethodType(py::cpp_function([](py::object& o, REArrayd& molefrac){ return get_typed<AmmoniaWaterTillnerRoth>(o).get_Treducing(molefrac); }, "self"_a, "molefrac"_a), obj));
@@ -138,7 +172,7 @@ void attach_model_specific_methods(py::object& obj){
return m.alphar_departure(tau, delta_, xNH3).imag()/h;
}, "self"_a, "tau"_a, "delta"_a, "xNH3"_a), obj));
}
- else if (std::holds_alternative<multifluid_t>(model)){
+ else if (index == multifluid_i){
attach_multifluid_methods<multifluid_t>(obj);
setattr("build_ancillaries", MethodType(py::cpp_function([](py::object& o, std::optional<int> i = std::nullopt){
const auto& c = get_typed<multifluid_t>(o);
@@ -154,7 +188,7 @@ void attach_model_specific_methods(py::object& obj){
return teqp::MultiFluidVLEAncillaries(jancillaries);
}, "self"_a, py::arg_v("i", std::nullopt, "None")), obj));
}
- else if (std::holds_alternative<idealgas_t>(model)){
+ else if (index == idealgas_i){
// Here X-Macros are used to create functions like get_Aig00, get_Aig01, ....
#define X(i,j) setattr(stringify(get_Aig ## i ## j), MethodType(py::cpp_function([](py::object& o, double T, double rho, REArrayd& molefrac){ \
using tdx = teqp::TDXDerivatives<idealgas_t, double, REArrayd>; \
@@ -163,17 +197,17 @@ void attach_model_specific_methods(py::object& obj){
ARXY_args
#undef X
}
- else if (std::holds_alternative<multifluidmutant_t>(model)){
+ else if (index == multifluidmutant_i){
attach_multifluid_methods<multifluidmutant_t>(obj);
}
- else if (std::holds_alternative<SW_EspindolaHeredia2009_t>(model)){
+ else if (index == SW_EspindolaHeredia2009_i){
// Have to use a method because lambda is a reserved word in Python
setattr("get_lambda", MethodType(py::cpp_function([](py::object& o){ return get_typed<SW_EspindolaHeredia2009_t>(o).get_lambda(); }, "self"_a), obj));
}
- else if (std::holds_alternative<EXP6_Kataoka1992_t>(model)){
+ else if (index == EXP6_Kataoka1992_i){
setattr("alpha", get_typed<EXP6_Kataoka1992_t>(obj).get_alpha());
}
- else if (std::holds_alternative<twocenterLJF_t>(model)){
+ else if (index == twocenterLJF_i){
setattr("Lstar", get_typed<twocenterLJF_t>(obj).L);
setattr("mustar_sq", get_typed<twocenterLJF_t>(obj).mu_sq);
}
@@ -314,7 +348,7 @@ void init_teqp(py::module& m) {
add_multifluid_mutant(m);
using am = teqp::cppinterface::AbstractModel;
- py::class_<AbstractModel, std::shared_ptr<AbstractModel>>(m, "AbstractModel", py::dynamic_attr())
+ py::class_<AbstractModel, std::unique_ptr<AbstractModel>>(m, "AbstractModel", py::dynamic_attr())
.def("get_R", &am::get_R, "molefrac"_a.noconvert())
diff --git a/setup.py b/setup.py
index b48177b92efca46dbeb1d2e901e09212bf27ae5e..3b710e756ad0981e582f157aaf299dda5ebe66ce 100644
--- a/setup.py
+++ b/setup.py
@@ -247,7 +247,7 @@ try:
description='Templated EQuation of state Package',
long_description='',
ext_modules=[CMakeExtension('teqp.teqp')], # teqp.teqp is the extension module that lives inside the teqp package
- packages=['teqp'],
+ packages=['teqp','teqp.fluiddata.dev.fluids','teqp.fluiddata.dev.mixtures'],
include_package_data=True,
cmdclass=dict(build_ext=CMakeBuild),
zip_safe=False,
diff --git a/src/bench_AbstractModel.cpp b/src/bench_AbstractModel.cpp
index 96bd5e59683288ca0599048407b89fdfccbd9930..34acbaab26e54f9b1cb79a0969ddb62d0b3a4cf2 100644
--- a/src/bench_AbstractModel.cpp
+++ b/src/bench_AbstractModel.cpp
@@ -3,6 +3,8 @@
#include "teqp/cpp/teqpcpp.hpp"
#include "teqp/cpp/derivs.hpp"
+#include "teqp/derivs.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
using namespace teqp;
@@ -18,7 +20,6 @@ TEST_CASE("multifluid derivatives", "[mf]")
}};
//std::cout << j.dump(2);
auto am = teqp::cppinterface::make_model(j);
- auto am2 = teqp::cppinterface::make_vdW1(2, 3);
auto z = (Eigen::ArrayXd(1) << 1.0).finished();
auto rhovec = 300.0* z;
@@ -51,3 +52,33 @@ TEST_CASE("multifluid derivatives", "[mf]")
return teqp::cppinterface::build_iteration_Jv(vars, mat, mat2, 8.3144, 300.0, 300.0, z);
};
}
+
+TEST_CASE("multifluid derivatives via DerivativeAdapter", "[mf]")
+{
+ nlohmann::json j = {
+ {"kind", "multifluid"},
+ {"model", {
+ {"components", {"../mycp/dev/fluids/Methane.json"}},
+ {"BIP", "../mycp/dev/mixtures/mixture_binary_pairs.json"},
+ {"departure", "../mycp/dev/mixtures/mixture_departure_functions.json"}
+ }
+ }};
+ auto am = teqp::cppinterface::make_model(j);
+ auto model = multifluidfactory(j["model"]);
+
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ using namespace cppinterface;
+ using vd = VirialDerivatives<decltype(model), double, decltype(z)>;
+
+ BENCHMARK("B4 natively") {
+ return vd::get_Bnvir<4>(model, 300.0, z);
+ };
+ BENCHMARK("B4 via AbstractModel") {
+ return am->get_Bnvir(4, 300, z);
+ };
+ BENCHMARK("B4 via DerivativeAdapter") {
+ using namespace teqp::cppinterface::adapter;
+ return view(model)->get_Bnvir(4, 300, z);
+ };
+}
diff --git a/src/tests/catch_test_SAFTpolar.cxx b/src/tests/catch_test_SAFTpolar.cxx
index 1d6772028584d6a5f9d55d52e9f4ce26acd239c2..0fd88a3d7727df7d83f90de08abad55ed697486d 100644
--- a/src/tests/catch_test_SAFTpolar.cxx
+++ b/src/tests/catch_test_SAFTpolar.cxx
@@ -36,11 +36,11 @@ TEST_CASE("Evaluation of J^{(n)}", "[checkJvals]")
{
double Tstar = 2.0, rhostar = 0.4;
SECTION("Luckas v. G&T"){
- for (auto n = 4; n < 37; ++n){
+ for (auto n = 4; n < 16; ++n){
CAPTURE(n);
auto L = LuckasJIntegral(n).get_J(Tstar, rhostar);
auto G = GubbinsTwuJIntegral(n).get_J(Tstar, rhostar);
- CHECK(L == Approx(G));
+ CHECK(L == Approx(G).margin(0.05));
}
}
}
@@ -55,7 +55,7 @@ TEST_CASE("Evaluation of K^{(n,m)}", "[checkKvals]")
auto G = GubbinsTwuKIntegral(n,m).get_K(Tstar, rhostar);
CAPTURE(n);
CAPTURE(m);
- CHECK(L == Approx(G));
+ CHECK(L == Approx(G).margin(std::abs(L)/2));
}
}
}
@@ -183,9 +183,10 @@ TEST_CASE("Check Stockmayer critical points with polarity terms", "[SAFTVRMiepol
// }
SECTION("With multipolar terms"){
- for (std::string polar_model : {"GrossVrabec", "GubbinsTwu+Luckas", "GubbinsTwu+GubbinsTwu"}){
- const bool print = true;
+ for (std::string polar_model : {"GubbinsTwu+Luckas", "GubbinsTwu+GubbinsTwu", "GrossVrabec"}){
+ const bool print = false;
double Tstar_guess = Tstar_guess_init, rhostar_guess = rhostar_guess_init;
+ if (print) std::cout << "===== " << polar_model << " =====" << std::endl;
if (print) std::cout << "(mu^*)^2, T^*, rho^*" << std::endl;
auto critpure = nonpolar_model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)));
@@ -196,7 +197,8 @@ TEST_CASE("Check Stockmayer critical points with polarity terms", "[SAFTVRMiepol
double mustar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23);
double Qstar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23);
j["model"]["polar_model"] = polar_model;
- for (double mustar2 = 0.0; mustar2 < 5; mustar2 += 0.1){
+
+ for (double mustar2 = 0.1; mustar2 < 5; mustar2 += 0.1){
j["model"]["coeffs"][0]["mu_Cm"] = sqrt(mustar2/mustar2factor*(ek*pow(sigma_m, 3)));
j["model"]["coeffs"][0]["nmu"] = 1;
auto model = teqp::cppinterface::make_model(j);
@@ -242,7 +244,7 @@ TEST_CASE("Check Stockmayer critical points with polarity terms", "[SAFTVRMiepol
std::valarray<std::tuple<double, double>> mu2Q22CLJ = {{6,2}, {6,4}, {12,2}, {12,4}};
SECTION("With mu&Q terms"){
- const bool print = true;
+ const bool print = false;
double Tstar_guess = Tstar_guess_init, rhostar_guess = rhostar_guess_init;
if (print) std::cout << "(mu^*)^2, T^*, rho^*" << std::endl;
diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx
index 1fb237d31566b3497dfd07fc0070140c51001d0b..f325e7df147bb73c86f98aa03bbb6f9bafe5cfe7 100644
--- a/src/tests/catch_test_cubics.cxx
+++ b/src/tests/catch_test_cubics.cxx
@@ -129,7 +129,8 @@ TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixtu
std::valarray<double> Tc_K = { 190.564, 154.581},
pc_Pa = { 4599200, 5042800},
acentric = { 0.011, 0.022};
- auto model = canonical_PR(Tc_K, pc_Pa, acentric);
+ const auto modelptr = teqp::cppinterface::adapter::make_owned(canonical_PR(Tc_K, pc_Pa, acentric));
+ const auto& model = teqp::cppinterface::adapter::get_model_cref<canonical_cubic_t>(modelptr.get());
const auto N = Tc_K.size();
using state_type = std::vector<double>;
REQUIRE(N == 2);
@@ -164,8 +165,7 @@ TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixtu
auto rho = rhovecL.sum();
auto molefrac = rhovecL / rhovecL.sum();
- using id = IsochoricDerivatives<decltype(model)>;
- auto pfromderiv = rho * model.R(molefrac) * T + id::get_pr(model, T, rhovecL);
+ auto pfromderiv = rho * modelptr->R(molefrac) * T + modelptr->get_pr(T, rhovecL);
return pfromderiv;
};