diff --git a/doc/source/models/MultifluidPlusAC.ipynb b/doc/source/models/MultifluidPlusAC.ipynb new file mode 100644 index 00000000..393b46d7 --- /dev/null +++ b/doc/source/models/MultifluidPlusAC.ipynb @@ -0,0 +1,254 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "ef2685a7", + "metadata": {}, + "source": [ + "# Multifluid+Activity\n", + "\n", + "Following the derivations in [Jaeger et al.](https://doi.org/10.1016/j.fluid.2018.04.015), the multifluid plus activity coefficient model can be defined as a the conventional corresponding states portion\n", + "$$\n", + "\\alpha^{\\rm r}_{\\rm CS} = \\sum_i \\alpha^{\\rm r}_{oi}(\\tau, \\delta)\n", + "$$\n", + "plus the departure term coming from the activity coefficient model\n", + "$$\n", + "\\alpha^{\\rm Dep}(T,\\rho,\\bar x) = \\frac{\\ln(1+b_{\\rm mix}\\rho)}{\\ln(1+\\frac{1}{u})}\\left[\\frac{g^{\\rm E,R}_{\\rm GE}(T,\n", + "\\bar x)}{RT}-\\sum_{i=1}^Nx_i\\left(\\alpha^{\\rm r}_{oi}(\\delta_{\\rm ref},\\tau)-\\alpha^{\\rm r}_{oi}(\\delta_{i,ref}, \\tau_i)\\right)\\right]\n", + "$$\n", + "with \n", + "$$\n", + "b_{\\rm mix} = \\sum_i b_ix_i\n", + "$$\n", + "$$\n", + "\\frac{g^{\\rm E,R}_{\\rm GE}(T,\n", + "\\bar x)}{RT} = \\sum_ix_i\\ln\\left(\\gamma_i^{\\rm R}\\right)\n", + "$$\n", + "\n", + "$$\n", + "\\delta_{\\rm ref} = \\frac{1}{ub_{\\rm mix}\\rho_r(\\bar x)}\n", + "$$\n", + "$$\n", + "\\delta_{i,{\\rm ref}} = \\frac{1}{ub_i\\rho_{c,i}}\n", + "$$\n", + "and the conventional reducing function is used to define $T_{\\rm r}(\\bar x)$ and $\\rho_{\\rm r}(\\bar x)$.\n", + "\n", + "Any activity coefficient model can be used for $\\gamma_i^{\\rm R}$, but so far (as of version 0.22) only the Wilson model is available due to how well it works with the advanced cubic mixing rules." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "937c2216", + "metadata": {}, + "outputs": [], + "source": [ + "import io, teqp, numpy as np, CoolProp.CoolProp as CP\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4e03c647", + "metadata": {}, + "outputs": [], + "source": [ + "# Four isotherms of experimental data for CO2 + Nitrogen from doi: 10.1016/j.fluid.2016.05.015\n", + "import pandas\n", + "dat = pandas.read_csv(io.StringIO(\"\"\"PointID y1 uy1 x1 ux1 p/bar up T/K\n", + "1 0.0274 0.0007 0.0068 0.0002 59.830 0.053 293.10 \n", + "2 0.0664 0.0014 0.0183 0.0004 64.864 0.080 293.10 \n", + "3 0.0978 0.0020 0.0298 0.0007 69.772 0.080 293.10 \n", + "4 0.1199 0.0024 0.0424 0.0009 74.737 0.080 293.10 \n", + "5 0.1219 0.0028 0.1132 0.0023 89.869 0.080 293.10 \n", + "6 0.1339 0.0024 0.0995 0.0022 89.198 0.080 293.10 \n", + "7 0.1399 0.0026 0.0943 0.0020 88.853 0.080 293.10 \n", + "8 0.1461 0.0027 0.0823 0.0019 86.962 0.080 293.10 \n", + "9 0.1466 0.0028 0.0778 0.0017 85.942 0.080 293.10 \n", + "10 0.1466 0.0028 0.0772 0.0016 85.868 0.080 293.10 \n", + "1 0.1378 0.0027 0.0159 0.0004 42.667 0.051 273.08 \n", + "2 0.2143 0.0038 0.0297 0.0007 49.547 0.051 273.08 \n", + "3 0.2612 0.0043 0.0411 0.0009 55.238 0.051 273.08 \n", + "4 0.3209 0.0049 0.0609 0.0013 65.069 0.088 273.08 \n", + "5 0.3554 0.0051 0.0786 0.0016 73.395 0.088 273.08 \n", + "6 0.3758 0.0052 0.0978 0.0019 81.061 0.088 273.08 \n", + "7 0.3903 0.0053 0.1190 0.0023 90.706 0.088 273.08 \n", + "8 0.3914 0.0053 0.1477 0.0028 100.966 0.088 273.08 \n", + "9 0.3879 0.0053 0.1614 0.0030 104.806 0.088 273.08 \n", + "10 0.3724 0.0052 0.1875 0.0033 110.846 0.088 273.08\n", + "11 0.3550 0.0051 0.2068 0.0036 114.105 0.088 273.08\n", + "12 0.2727 0.0044 0.2531 0.0041 118.020 0.088 273.08\n", + "13 0.3343 0.0049 0.2268 0.0038 116.295 0.088 273.08\n", + "1 0.2048 0.0038 0.0106 0.0003 25.754 0.050 253.05 \n", + "2 0.3019 0.0049 0.0217 0.0005 30.479 0.050 253.05 \n", + "3 0.4638 0.0056 0.0436 0.0010 45.352 0.050 253.05 \n", + "4 0.5319 0.0056 0.0647 0.0014 58.188 0.050 253.05 \n", + "5 0.5854 0.0054 0.1077 0.0021 78.315 0.084 253.05 \n", + "6 0.5979 0.0054 0.1497 0.0028 98.276 0.084 253.05\n", + "7 0.5898 0.0054 0.1801 0.0032 109.241 0.084 253.05\n", + "8 0.5042 0.0057 0.0570 0.0012 51.343 0.084 253.05\n", + "9 0.5644 0.0055 0.0861 0.0017 67.594 0.084 253.05\n", + "10 0.5949 0.0054 0.1267 0.0024 86.883 0.084 253.05\n", + "11 0.5826 0.0054 0.2015 0.0035 116.614 0.084 253.05\n", + "12 0.5537 0.0055 0.2431 0.0040 129.873 0.084 253.05\n", + "13 0.4973 0.0055 0.2971 0.0046 139.161 0.084 253.05\n", + "14 0.4971 0.0055 0.2972 0.0046 139.261 0.084 253.05\n", + "1 0.7076 0.0050 0.0257 0.0006 27.983 0.056 223.10\n", + "2 0.7774 0.0041 0.0522 0.0011 44.918 0.056 223.10\n", + "3 0.8077 0.0036 0.0930 0.0019 64.906 0.081 223.10\n", + "4 0.8131 0.0035 0.1261 0.0024 84.799 0.081 223.10\n", + "5 0.8057 0.0035 0.1584 0.0029 104.410 0.081 223.10\n", + "6 0.7843 0.0038 0.1982 0.0035 125.782 0.081 223.10\n", + "7 0.7533 0.0041 0.2380 0.0040 144.287 0.081 223.10\n", + "8 0.7150 0.0045 0.2813 0.0044 159.015 0.081 223.10\n", + "9 0.6942 0.0047 0.3064 0.0047 165.347 0.081 223.10\n", + "\"\"\"), sep='\\s+', engine='python')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74404b8c", + "metadata": {}, + "outputs": [], + "source": [ + "Tc_K = np.array([304.21, 126.19])\n", + "pc_Pa = np.array([7.383e6, 3395800.0])\n", + "\n", + "def get_bi():\n", + " OmegaB = 0.08664 \n", + " R = 8.31446261815324\n", + " return (OmegaB*R*Tc_K/pc_Pa).tolist()\n", + "\n", + "donor = teqp.build_multifluid_model([\"CO2\", \"Nitrogen\"], teqp.get_datapath())\n", + "vc_m3mol = donor.get_vcvec()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab83b1d7", + "metadata": {}, + "outputs": [], + "source": [ + "# Parameter set with linear mixing in the reducing function\n", + "BIP_linear = [{\n", + " \"Name1\": \"CarbonDioxide\",\n", + " \"Name2\": \"Nitrogen\",\n", + " \"betaT\": 1.0,\n", + " \"gammaT\": 0.5*(Tc_K[0]+Tc_K[1])/(Tc_K[0]*Tc_K[1])**0.5,\n", + " \"betaV\": 1.0,\n", + " \"gammaV\": 4*(vc_m3mol[0]+vc_m3mol[1])/(vc_m3mol[0]**(1/3) + vc_m3mol[1]**(1/3))**3,\n", + " \"F\": 0.0\n", + "}]\n", + "\n", + "# Parameter set with Lorentz-Berthelot\n", + "BIP_LorentzBerthelot = [{\n", + " \"Name1\": \"CarbonDioxide\",\n", + " \"Name2\": \"Nitrogen\",\n", + " \"betaT\": 1.0,\n", + " \"gammaT\": 1.0,\n", + " \"betaV\": 1.0,\n", + " \"gammaV\": 1.0,\n", + " \"F\": 0.0\n", + "}]\n", + "\n", + "def get_model(BIPset=None):\n", + " if BIPset is not None and BIPset == 'linear':\n", + " return teqp.make_model({\n", + " 'kind': 'multifluid',\n", + " \"model\":{\n", + " \"components\": [\"CO2\", \"Nitrogen\"],\n", + " \"root\": teqp.get_datapath(),\n", + " \"BIP\": BIP_linear\n", + " }})\n", + "\n", + " # The Wilson activity coefficient model is from Lasala et al.: https://doi.org/10.1016/j.fluid.2016.05.015\n", + " # and then combined with the \n", + " j = {\n", + " \"kind\": \"multifluid-activity\",\n", + " \"model\": {\n", + " \"multifluid\":{\n", + " \"components\": [\"CO2\", \"Nitrogen\"],\n", + " \"root\": teqp.get_datapath(),\n", + " \"BIP\": BIPset\n", + " },\n", + " \"activity\": {\n", + " \"aresmodel\": {\n", + " \"type\": \"Wilson\", \n", + " \"m\": [[0.0, -3.4768], [3.5332, 0.0]], \n", + " \"n\": [[0.0, 825], [-585, 0.0]], \n", + " \"b\": get_bi()\n", + " },\n", + " \"options\": {\"u\": 1.17, \"b\": get_bi()}\n", + " }\n", + " }\n", + " }\n", + "\n", + " return teqp.make_model(j) if BIPset is not None else donor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "500aa1ae", + "metadata": {}, + "outputs": [], + "source": [ + "labels = ['linear+Wilson','Lorentz-Berthelot+Wilson','GERG-2008','linear']\n", + "for imod, model in enumerate([get_model(BIP_linear), get_model(BIP_LorentzBerthelot), get_model(), get_model('linear')]):\n", + " lw = 0.5+imod\n", + " for T in [223.15, 253.05, 273.08, 293.1]:\n", + " ipure = 0\n", + "\n", + " [rhoL0, rhoV0] = [CP.PropsSI('Dmolar','T',T,'Q',Q,'CO2') for Q in [0,1]]\n", + "\n", + " rhovecL0 = np.array([0.0, 0.0]); rhovecL0[ipure] = rhoL0\n", + " rhovecV0 = np.array([0.0, 0.0]); rhovecV0[ipure] = rhoV0\n", + "\n", + " J = model.trace_VLE_isotherm_binary(T, rhovecL0, rhovecV0)\n", + " df = pandas.DataFrame(J)\n", + " plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6,'k', lw=lw, label=labels[imod] if T == 223.15 else '')\n", + " plt.plot(df['xV_0 / mole frac.'], df['pV / Pa']/1e6,'k', lw=lw)\n", + " \n", + "plt.title('CO$_2$ (1) + N$_2$ (2)')\n", + "plt.plot(1-dat['x1'], dat['p/bar']/10, 'o')\n", + "plt.plot(1-dat['y1'], dat['p/bar']/10, '^')\n", + "plt.legend(loc='best')\n", + "plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / MPa', ylim=(0, 25))\n", + "plt.savefig(\"multifluid_Wilson_CO2_N2.pdf\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e1b7813a", + "metadata": {}, + "source": [ + "so we can see that the activity coefficient models are not too bad, clearly an improvement over pure linear mixing, but the model from GERG-2008 is unsurprisingly the most accurate" + ] + } + ], + "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/doc/source/models/index.rst b/doc/source/models/index.rst index 04a5b29c..213ae696 100644 --- a/doc/source/models/index.rst +++ b/doc/source/models/index.rst @@ -28,5 +28,6 @@ Models GERG ECS MultifluidPlusAssociation + MultifluidPlusAC IdealGas \ No newline at end of file diff --git a/include/teqp/models/activity/COSMOSAC.hpp b/include/teqp/models/activity/COSMOSAC.hpp new file mode 100644 index 00000000..cd05e504 --- /dev/null +++ b/include/teqp/models/activity/COSMOSAC.hpp @@ -0,0 +1,360 @@ +/** + This header is derived from the code in https://github.com/usnistgov/COSMOSAC, and the methodology is described in the paper: + + Ian H. Bell, Erik Mickoleit, Chieh-Ming Hsieh, Shiang-Tai Lin, Jadran Vrabec, Cornelia Breitkopf, and Andreas Jäger + Journal of Chemical Theory and Computation 2020 16 (4), 2635-2646 + DOI: 10.1021/acs.jctc.9b01016 + */ + +namespace teqp::activity::activity_models::COSMOSAC{ + +/** + A single sigma profile. In this data structure (for consistency with the 2005 Virginia Tech + database of COSMO-SAC parameters), the first column is the electron density in [e/A^2], and the second column + is the probability of finding a segment with this electron density, multiplied by the segment area, in A^2 + */ +class SigmaProfile { +public: + Eigen::ArrayXd m_sigma, m_psigmaA; + /// Default constructor + SigmaProfile() {}; + /// Copy-by-reference constructor + SigmaProfile(const Eigen::ArrayXd &sigma, const Eigen::ArrayXd &psigmaA) : m_sigma(sigma), m_psigmaA(psigmaA) {}; + /// Copy-by-reference constructor with std::vector + SigmaProfile(const std::vector &sigma, const std::vector &psigmaA) : m_sigma(Eigen::Map(&(sigma[0]), sigma.size())), m_psigmaA(Eigen::Map(&(psigmaA[0]), psigmaA.size())) {}; + /// Move constructor + SigmaProfile(Eigen::ArrayXd &&sigma, Eigen::ArrayXd &&psigmaA) : m_sigma(sigma), m_psigmaA(psigmaA) {}; + const Eigen::ArrayXd &psigmaA() const { return m_psigmaA; } + const Eigen::ArrayXd &sigma() const { return m_sigma; } + const Eigen::ArrayXd psigma(double A_i) const { return m_psigmaA/A_i; } +}; + +// A holder class for the sigma profiles for a given fluid +struct FluidSigmaProfiles { + SigmaProfile nhb, ///< The profile for the non-hydrogen-bonding segments + oh, ///< The profile for the OH-bonding segments + ot; ///< The profile for the "other" segments +}; + +struct COSMO3Constants { + double + AEFFPRIME = 7.25, // [A^2] + c_OH_OH = 4013.78, // [kcal A^4/(mol e^2)] + c_OT_OT = 932.31, // [kcal A^4 /(mol e^2)] + c_OH_OT = 3016.43, // [kcal A^4 /(mol e^2)] + A_ES = 6525.69, // [kcal A^4 /(mol e^2)] + B_ES = 1.4859e8, // [kcal A^4 K^2/(mol e^2)] + N_A = 6.022140758e23, // [mol^{-1}] + k_B = 1.38064903e-23, // [J K^{-1}] + R = k_B*N_A/4184, // [kcal/(mol*K)]; Universal gas constant of new redefinition of 2018, https://doi.org/10.1088/1681-7575/aa99bc + Gamma_rel_tol = 1e-8; // relative tolerance for Gamma in iterative loop + bool fast_Gamma = false; + std::string to_string() { + return "NOT IMPLEMENTED YET"; + //return "c_hb: " + std::to_string(c_hb) + " kcal A^4 /(mol*e^2) \nsigma_hb: " + std::to_string(sigma_hb) + " e/A^2\nalpha_prime: " + std::to_string(alpha_prime) + " kcal A^4 /(mol*e^2)\nAEFFPRIME: " + std::to_string(AEFFPRIME) + " A\nR: " + std::to_string(R) + " kcal/(mol*K)"; + } +}; + +enum class profile_type { NHB_PROFILE, OH_PROFILE, OT_PROFILE }; + +class COSMO3 { +private: + std::vector A_COSMO_A2; ///< The area per fluid, in \AA^2 + std::vector profiles; ///< The vector of profiles, one per fluid + COSMO3Constants m_consts; + Eigen::Index ileft, w; +public: + COSMO3(const std::vector& A_COSMO_A2, const std::vector &SigmaProfiles, const COSMO3Constants &constants = COSMO3Constants()) + : A_COSMO_A2(A_COSMO_A2), profiles(SigmaProfiles), m_consts(constants) { + Eigen::Index iL, iR; + std::tie(iL, iR) = get_nonzero_bounds(); + this->ileft = iL; this->w = iR - iL + 1; + }; + + /// Determine the left and right bounds for the psigma that is nonzero to accelerate the iterative calculations + std::tuple get_nonzero_bounds(){ + // Determine the range of entries in p(sigma) that are greater than zero, we + // will only calculate segment activity coefficients for these segments + Eigen::Index min_ileft = 51, max_iright = 0; + for (auto i = 0U; i < profiles.size(); ++i) { + const Eigen::ArrayXd psigma = profiles[i].nhb.psigma(A_COSMO_A2[i]); + Eigen::Index ileft = 0, iright = psigma.size(); + for (auto ii = 0; ii < psigma.size(); ++ii) { if (std::abs(psigma(ii)) > 1e-16) { ileft = ii; break; } } + for (auto ii = psigma.size() - 1; ii > ileft; --ii) { if (std::abs(psigma(ii)) > 1e-16) { iright = ii; break; } } + if (ileft < min_ileft) { min_ileft = ileft; } + if (iright > max_iright) { max_iright = iright; } + } + return std::make_tuple(min_ileft, max_iright); + } + + template + auto get_psigma_mix(const MoleFractions &z, profile_type type = profile_type::NHB_PROFILE) const { + Eigen::ArrayX> psigma_mix(51); psigma_mix.fill(0); + std::decay_t xA = 0.0; + for (auto i = 0; i < z.size(); ++i) { + switch (type) { + case profile_type::NHB_PROFILE: + psigma_mix += z[i] * profiles[i].nhb.psigmaA(); break; + case profile_type::OH_PROFILE: + psigma_mix += z[i] * profiles[i].oh.psigmaA(); break; + case profile_type::OT_PROFILE: + psigma_mix += z[i] * profiles[i].ot.psigmaA(); break; + } + xA += z[i] * A_COSMO_A2[i]; + } + psigma_mix /= xA; + return psigma_mix; + } + + /// Get access to the parameters that are in use + COSMO3Constants &get_mutable_COSMO_constants() { return m_consts; } + + std::tuple get_ileftw() const { return std::make_tuple(ileft, w);} ; + + double get_c_hb(profile_type type1, profile_type type2) const{ + + if (type1 == type2){ + if (type1 == profile_type::OH_PROFILE) { + return m_consts.c_OH_OH; + } + else if (type1 == profile_type::OT_PROFILE) { + return m_consts.c_OT_OT; + } + else { + return 0.0; + } + } + else if ((type1 == profile_type::OH_PROFILE && type2 == profile_type::OT_PROFILE) + || (type1 == profile_type::OT_PROFILE && type2 == profile_type::OH_PROFILE)) { + return m_consts.c_OH_OT; + } + else { + return 0.0; + } + } + template + auto get_DELTAW(const TType& T, profile_type type_t, profile_type type_s) const { + auto delta_sigma = 2*0.025/50; + Eigen::ArrayXX DELTAW(51, 51); + double cc = get_c_hb(type_t, type_s); + for (auto m = 0; m < 51; ++m) { + for (auto n = 0; n < 51; ++n) { + double sigma_m = -0.025 + delta_sigma*m, + sigma_n = -0.025 + delta_sigma*n, + c_hb = (sigma_m*sigma_n >= 0) ? 0 : cc; + auto c_ES = m_consts.A_ES + m_consts.B_ES/(T*T); + DELTAW(m, n) = c_ES*POW2(sigma_m + sigma_n) - c_hb*POW2(sigma_m-sigma_n); + } + } + return DELTAW; + } + template + auto get_DELTAW_fast(TType T, profile_type type_t, profile_type type_s) const { + auto delta_sigma = 2 * 0.025 / 50; + Eigen::ArrayXX DELTAW(51, 51); DELTAW.setZero(); + double cc = get_c_hb(type_t, type_s); + for (auto m = ileft; m < ileft+w+1; ++m) { + for (auto n = ileft; n < ileft+w+1; ++n) { + double sigma_m = -0.025 + delta_sigma*m, + sigma_n = -0.025 + delta_sigma*n, + c_hb = (sigma_m*sigma_n >= 0) ? 0 : cc; + auto c_ES = m_consts.A_ES + m_consts.B_ES / (T*T); + DELTAW(m, n) = c_ES*POW2(sigma_m + sigma_n) - c_hb*POW2(sigma_m - sigma_n); + } + } + return DELTAW; + } + + template + auto get_AA(const TType& T, PSigma psigmas){ + // Build the massive \Delta W matrix that is 153*153 in size + Eigen::ArrayXX DELTAW(51, 51); // Depends on temperature + double R = m_consts.R; + std::vector types = { profile_type::NHB_PROFILE, profile_type::OH_PROFILE, profile_type::OT_PROFILE }; + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + DELTAW.matrix().block(51 * i, 51 * j, 51, 51) = get_DELTAW(T, types[i], types[j]); + } + } + return Eigen::exp(-DELTAW/(R*T)).rowwise()*psigmas.transpose(); + } + /** + Obtain the segment activity coefficients \f$\Gamma\f$ for given a set of charge densities + + If fast_Gamma is enabled, then only segments that have some contribution are included in the iteration (can be a significant acceleration for nonpolar+nonpolar mixtures) + + \param T Temperature, in K + \param psigmas Charge densities, in the order of NHB, OH, OT. Length is 153. + */ + template + auto get_Gamma(const TType& T, PSigmaType psigmas) const { + auto startTime = std::chrono::high_resolution_clock::now(); + + using TXType = std::decay_t>; + + double R = m_consts.R; + Eigen::ArrayX Gamma(153), Gammanew(153); Gamma.setOnes(); Gammanew.setOnes(); + + // A convenience function to convert double values to string in scientific format + auto to_scientific = [](double val) { std::ostringstream out; out << std::scientific << val; return out.str(); }; + + auto max_iter = 2000; + if (!m_consts.fast_Gamma){ + // The slow and simple branch + + // Build the massive \Delta W matrix that is 153*153 in size + Eigen::ArrayXX DELTAW(153, 153); + std::vector types = { profile_type::NHB_PROFILE, profile_type::OH_PROFILE, profile_type::OT_PROFILE }; + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + DELTAW.matrix().block(51*i, 51*j, 51, 51) = get_DELTAW(T, types[i], types[j]); + } + } + + auto AA = (Eigen::exp(-DELTAW / (R*T)).template cast().rowwise()*psigmas.template cast().transpose()).eval(); + for (auto counter = 0; counter <= max_iter; ++counter) { + Gammanew = 1 / (AA.rowwise()*Gamma.transpose()).rowwise().sum(); + Gamma = (Gamma + Gammanew) / 2; + double maxdiff = getbaseval(((Gamma - Gammanew) / Gamma).cwiseAbs().real().maxCoeff()); + if (maxdiff < m_consts.Gamma_rel_tol) { + break; + } + if (counter == max_iter) { + throw std::invalid_argument("Could not obtain the desired tolerance of " + + to_scientific(m_consts.Gamma_rel_tol) + + " after " + + std::to_string(max_iter) + + " iterations in get_Gamma; current value is " + + to_scientific(maxdiff)); + } + } + return Gamma; + } + else { + // The fast branch! + // ---------------- + + // Build the massive AA matrix that is 153*153 in size + auto midTime = std::chrono::high_resolution_clock::now(); + std::vector types = { profile_type::NHB_PROFILE, profile_type::OH_PROFILE, profile_type::OT_PROFILE }; + std::vector offsets = {0*51, 1*51, 2*51}; + Eigen::ArrayXX AA(153, 153); + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + Eigen::Index rowoffset = offsets[i], coloffset = offsets[j]; + AA.matrix().block(rowoffset + ileft, coloffset + ileft, w, w) = Eigen::exp(-get_DELTAW_fast(T, types[i], types[j]).block(ileft, ileft, w, w).array() / (R*T)).template cast().rowwise()*psigmas.template cast().segment(coloffset+ileft,w).transpose(); + } + } + auto midTime2 = std::chrono::high_resolution_clock::now(); + + for (auto counter = 0; counter <= max_iter; ++counter) { + for (Eigen::Index offset : {51*0, 51*1, 51*2}){ + Gammanew.segment(offset + ileft, w) = 1 / ( + AA.matrix().block(offset+ileft,51*0+ileft,w,w).array().rowwise()*Gamma.segment(51*0+ileft, w).transpose() + + AA.matrix().block(offset+ileft,51*1+ileft,w,w).array().rowwise()*Gamma.segment(51*1+ileft, w).transpose() + + AA.matrix().block(offset+ileft,51*2+ileft,w,w).array().rowwise()*Gamma.segment(51*2+ileft, w).transpose() + ).rowwise().sum(); + } + for (Eigen::Index offset : {51 * 0, 51 * 1, 51 * 2}) { + Gamma.segment(offset + ileft, w) = (Gamma.segment(offset + ileft, w) + Gammanew.segment(offset + ileft, w)) / 2; + } + double maxdiff = getbaseval(((Gamma - Gammanew) / Gamma).cwiseAbs().real().maxCoeff()); + if (maxdiff < m_consts.Gamma_rel_tol) { + break; + } + if (counter == max_iter){ + throw std::invalid_argument("Could not obtain the desired tolerance of " + + to_scientific(m_consts.Gamma_rel_tol) + +" after " + +std::to_string(max_iter) + +" iterations in get_Gamma; current value is " + + to_scientific(maxdiff)); + } + } + auto endTime = std::chrono::high_resolution_clock::now(); + //std::cout << std::chrono::duration(midTime - startTime).count() << " s elapsed (DELTAW)\n"; + //std::cout << std::chrono::duration(midTime2 - midTime).count() << " s elapsed (AA)\n"; + //std::cout << std::chrono::duration(endTime - midTime2).count() << " s elapsed (comps)\n"; + //std::cout << std::chrono::duration(endTime - startTime).count() << " s elapsed (total)\n"; + return Gamma; + } + } + + /** + The residual part of ln(γ_i), for the i-th component + */ + template + auto get_lngamma_resid(std::size_t i, TType T, const Array &lnGamma_mix) const + { + double AEFFPRIME = m_consts.AEFFPRIME; + Eigen::ArrayX psigmas(3*51); // For a pure fluid, p(sigma) does not depend on temperature + double A_i = A_COSMO_A2[i]; + psigmas << profiles[i].nhb.psigma(A_i), profiles[i].oh.psigma(A_i), profiles[i].ot.psigma(A_i); + // double check_sum = psigmas.sum(); //// Should sum to 1.0 + auto lnGammai = get_Gamma(T, psigmas).log().eval(); + return A_i/AEFFPRIME*(psigmas*(lnGamma_mix - lnGammai)).sum(); + } + /** + This overload is a convenience overload, less computationally + efficient, but simpler to use and more in keeping with the other + contributions. It does not require the lnGamma_mix to be pre-calculated + and passed into this function + */ + + template + auto get_lngamma_resid(TType T, const MoleFracs& molefracs) const + { + using TXType = std::decay_t>; + + Eigen::Array psigmas; + psigmas << get_psigma_mix(molefracs, profile_type::NHB_PROFILE), get_psigma_mix(molefracs, profile_type::OH_PROFILE), get_psigma_mix(molefracs, profile_type::OT_PROFILE); + + Eigen::ArrayX lngamma(molefracs.size()); + // double check_sum = psigmas.sum(); //// Should sum to 1.0 + Eigen::ArrayX lnGamma_mix = get_Gamma(T, psigmas).log(); + for (Eigen::Index i = 0; i < molefracs.size(); ++i) { + lngamma(i) = get_lngamma_resid(i, T, lnGamma_mix); + } + return lngamma; + } + + // EigenArray get_lngamma_disp(const EigenArray &x) const { + // if (x.size() != 2){ throw std::invalid_argument("Multi-component mixtures not supported for dispersive contribution yet"); } + // double w = 0.27027; // default value + // auto cls0 = m_fluids[0].dispersion_flag, cls1 = m_fluids[1].dispersion_flag; + // using d = FluidProfiles::dispersion_classes; + // if ((cls0 == d::DISP_WATER && cls1 == d::DISP_ONLY_ACCEPTOR) + // | + // (cls1 == d::DISP_WATER && cls0 == d::DISP_ONLY_ACCEPTOR)){ + // w = -0.27027; + // } + // else if ((cls0 == d::DISP_WATER && cls1 == d::DISP_COOH) + // | + // (cls1 == d::DISP_WATER && cls0 == d::DISP_COOH)) { + // w = -0.27027; + // } + // else if ((cls0 == d::DISP_COOH && (cls1 == d::DISP_NHB || cls1 == d::DISP_DONOR_ACCEPTOR)) + // | + // (cls1 == d::DISP_COOH && (cls0 == d::DISP_NHB || cls0 == d::DISP_DONOR_ACCEPTOR))) { + // w = -0.27027; + // } + // + // double ekB0 = m_fluids[0].dispersion_eoverkB, ekB1 = m_fluids[1].dispersion_eoverkB; + // double A = w*(0.5*(ekB0+ekB1) - sqrt(ekB0*ekB1)); + // EigenArray lngamma_dsp(2); + // lngamma_dsp(0) = A*x[1]*x[1]; + // lngamma_dsp(1) = A*x[0]*x[0]; + // return lngamma_dsp; + // } + // EigenArray get_lngamma(double T, const EigenArray &x) const override { + // return get_lngamma_comb(T, x) + get_lngamma_resid(T, x) + get_lngamma_disp(x); + // } + + // Returns ares/RT + template + auto operator () (const TType& T, const MoleFractions& molefracs) const { + return contiguous_dotproduct(molefracs, get_lngamma_resid(T, molefracs)); + } +}; + +}; diff --git a/include/teqp/models/activity/activity_models.hpp b/include/teqp/models/activity/activity_models.hpp new file mode 100644 index 00000000..c5ea73d7 --- /dev/null +++ b/include/teqp/models/activity/activity_models.hpp @@ -0,0 +1,186 @@ +#pragma once + +#include "teqp/models/activity/COSMOSAC.hpp" + +namespace teqp::activity::activity_models{ + +/** + A residual Helmholtz term that returns nothing (the empty term) + */ +template +class NullResidualHelmholtzOverRT { +public: + template + auto operator () (const TType& /*T*/, const MoleFractions& molefracs) const { + std::common_type_t val = 0.0; + return val; + } +}; + +/** + + \f[ + \frac{a^{E,\gamma}_{total}}{RT} = -sum_iz_i\ln\left(\sum_jz_jOmega_{ij}(T)\right) + \f] + + \f[ + \frac{a^{E,\gamma}_{comb}}{RT} = -sum_iz_i\ln\left(\frac{\Omega_i}{z_i}\right) + \f] + + \f[ + \frac{a^{E,\gamma}_{res}}{RT} = \frac{a^{E,\gamma}_{total}}{RT} - \frac{a^{E,\gamma}_{comb}}{RT} + \f] + + Volume fraction of component \f$i\f$ +\f[ + \phi_i = \frac{z_iv_i}{\sum_j z_j v_j} + \f] + with \f$v_i = b_i\f$ + */ +template +class WilsonResidualHelmholtzOverRT { + +public: + const std::vector b; + const Eigen::ArrayXXd m, n; + WilsonResidualHelmholtzOverRT(const std::vector& b, const Eigen::ArrayXXd& m, const Eigen::ArrayXXd& n) : b(b), m(m), n(n) {}; + + template + auto combinatorial(const TType& /*T*/, const MoleFractions& molefracs) const { + if (b.size() != static_cast(molefracs.size())){ + throw teqp::InvalidArgument("Bad size of molefracs"); + } + + using TYPE = std::common_type_t; + // The denominator in Phi + TYPE Vtot = 0.0; + for (auto i = 0U; i < molefracs.size(); ++i){ + auto v_i = b[i]; + Vtot += molefracs[i]*v_i; + } + + TYPE summer = 0.0; + for (auto i = 0U; i < molefracs.size(); ++i){ + auto v_i = b[i]; + // The ratio phi_i/z_i is expressed like this to better handle + // the case of z_i = 0, which would otherwise be a divide by zero + // in the case that the composition of one component is zero + auto phi_i_over_z_i = v_i/Vtot; + summer += molefracs[i]*log(phi_i_over_z_i); + } + return summer; + } + + template + auto get_Aij(std::size_t i, std::size_t j, const TType& T) const{ + return forceeval(m(i,j)*T + n(i,j)); + } + + template + auto total(const TType& T, const MoleFractions& molefracs) const { + + using TYPE = std::common_type_t; + TYPE summer = 0.0; + for (auto i = 0U; i < molefracs.size(); ++i){ + auto v_i = b[i]; + TYPE summerj = 0.0; + for (auto j = 0U; j < molefracs.size(); ++j){ + auto v_j = b[j]; + auto Aij = get_Aij(i,j,T); + auto Omega_ji = v_j/v_i*exp(-Aij/T); + summerj += molefracs[j]*Omega_ji; + } + summer += molefracs[i]*log(summerj); + } + return forceeval(-summer); + } + + // Returns ares/RT + template + auto operator () (const TType& T, const MoleFractions& molefracs) const { + return forceeval(total(T, molefracs) - combinatorial(T, molefracs)); + } +}; + +///** +// \note This approach works well except that... the derivatives at the pure fluid endpoints don't. So this is more a record as a failed idea +// */ +//template +//class BinaryBetaResidualHelmholtzOverRT { +// +//public: +// const std::vector c; +// BinaryBetaResidualHelmholtzOverRT(const std::vector& c) : c(c) {}; +// +// // Returns ares/RT +// template +// auto operator () (const TType& /*T*/, const MoleFractions& molefracs) const { +// if (molefracs.size() != 2){ +// throw teqp::InvalidArgument("Must be size of 2"); +// } +// std::decay_t> out = c[0]*pow(molefracs[0], c[1])*pow(molefracs[1], c[2]); +// return out; +// } +//}; + +template +class BinaryInvariantResidualHelmholtzOverRT { + +public: + const std::vector c; + BinaryInvariantResidualHelmholtzOverRT(const std::vector& c) : c(c) {}; + + // Returns ares/RT + template + auto operator () (const TType& /*T*/, const MoleFractions& molefracs) const { + if (molefracs.size() != 2){ + throw teqp::InvalidArgument("Must be size of 2"); + } + std::decay_t> out = c[0]*molefracs[0]*molefracs[1]*(c[1] + c[2]*molefracs[1]); + return out; + } +}; + +using ResidualHelmholtzOverRTVariant = std::variant, WilsonResidualHelmholtzOverRT, BinaryInvariantResidualHelmholtzOverRT, COSMOSAC::COSMO3>; + +inline ResidualHelmholtzOverRTVariant ares_model_factory(const nlohmann::json& armodel) { + + std::string type = armodel.at("type"); + if (type == "Wilson"){ + std::vector b = armodel.at("b"); + auto mWilson = build_square_matrix(armodel.at("m")); + auto nWilson = build_square_matrix(armodel.at("n")); + return WilsonResidualHelmholtzOverRT(b, mWilson, nWilson); + } +// else if (type == "binaryBeta"){ +// std::vector c = armodel.at("c"); +// return BinaryBetaResidualHelmholtzOverRT(c); +// } + else if (type == "binaryInvariant"){ + std::vector c = armodel.at("c"); + return BinaryInvariantResidualHelmholtzOverRT(c); + } + else if (type == "COSMO-SAC-2010"){ + std::vector A_COSMOSAC_A2 = armodel.at("A_COSMOSAC / A^2"); + std::vector profiles; + for (auto& el : armodel.at("profiles")){ + COSMOSAC::FluidSigmaProfiles prof; + auto get_ = [](const auto& j){ + return COSMOSAC::SigmaProfile{ + j.at("sigma / e/A^2").template get>(), + j.at("p(sigma)*A / A^2").template get>() + }; + }; + prof.nhb = get_(el.at("nhb")); + prof.oh = get_(el.at("oh")); + prof.ot = get_(el.at("ot")); + profiles.push_back(prof); + } + return COSMOSAC::COSMO3(A_COSMOSAC_A2, profiles); + } + else{ + throw teqp::InvalidArgument("bad type of ares model: " + type); + } +}; + +} diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp index 10a985d4..eabdd6cb 100644 --- a/include/teqp/models/cubics.hpp +++ b/include/teqp/models/cubics.hpp @@ -17,6 +17,8 @@ Implementations of the canonical cubic equations of state #include "teqp/math/pow_templates.hpp" #include "nlohmann/json.hpp" +#include "teqp/models/activity/activity_models.hpp" +using namespace teqp::activity::activity_models; #include @@ -460,106 +462,6 @@ inline void from_json(const json& j, AdvancedPRaEOptions& o) { } } -/** - A residual Helmholtz term that returns nothing (the empty term) - */ -template -class NullResidualHelmholtzOverRT { -public: - template - auto operator () (const TType& /*T*/, const MoleFractions& molefracs) const { - std::common_type_t val = 0.0; - return val; - } -}; - -/** - - \f[ - \frac{a^{E,\gamma}_{total}}{RT} = -sum_iz_i\ln\left(\sum_jz_jOmega_{ij}(T)\right) - \f] - - \f[ - \frac{a^{E,\gamma}_{comb}}{RT} = -sum_iz_i\ln\left(\frac{\Omega_i}{z_i}\right) - \f] - - \f[ - \frac{a^{E,\gamma}_{res}}{RT} = \frac{a^{E,\gamma}_{total}}{RT} - \frac{a^{E,\gamma}_{comb}}{RT} - \f] - - Volume fraction of component \f$i\f$ -\f[ - \phi_i = \frac{z_iv_i}{\sum_j z_j v_j} - \f] - with \f$v_i = b_i\f$ - */ -template -class WilsonResidualHelmholtzOverRT { - -public: - const std::vector b; - const Eigen::ArrayXXd m, n; - WilsonResidualHelmholtzOverRT(const std::vector& b, const Eigen::ArrayXXd& m, const Eigen::ArrayXXd& n) : b(b), m(m), n(n) {}; - - template - auto combinatorial(const TType& /*T*/, const MoleFractions& molefracs) const { - if (b.size() != static_cast(molefracs.size())){ - throw teqp::InvalidArgument("Bad size of molefracs"); - } - - using TYPE = std::common_type_t; - // The denominator in Phi - TYPE Vtot = 0.0; - for (auto i = 0U; i < molefracs.size(); ++i){ - auto v_i = b[i]; - Vtot += molefracs[i]*v_i; - } - - TYPE summer = 0.0; - for (auto i = 0U; i < molefracs.size(); ++i){ - auto v_i = b[i]; - // The ratio phi_i/z_i is expressed like this to better handle - // the case of z_i = 0, which would otherwise be a divide by zero - // in the case that the composition of one component is zero - auto phi_i_over_z_i = v_i/Vtot; - summer += molefracs[i]*log(phi_i_over_z_i); - } - return summer; - } - - template - auto get_Aij(std::size_t i, std::size_t j, const TType& T) const{ - return forceeval(m(i,j)*T + n(i,j)); - } - - template - auto total(const TType& T, const MoleFractions& molefracs) const { - - using TYPE = std::common_type_t; - TYPE summer = 0.0; - for (auto i = 0U; i < molefracs.size(); ++i){ - auto v_i = b[i]; - TYPE summerj = 0.0; - for (auto j = 0U; j < molefracs.size(); ++j){ - auto v_j = b[j]; - auto Aij = get_Aij(i,j,T); - auto Omega_ji = v_j/v_i*exp(-Aij/T); - summerj += molefracs[j]*Omega_ji; - } - summer += molefracs[i]*log(summerj); - } - return forceeval(-summer); - } - - // Returns ares/RT - template - auto operator () (const TType& T, const MoleFractions& molefracs) const { - return forceeval(total(T, molefracs) - combinatorial(T, molefracs)); - } -}; - -using ResidualHelmholtzOverRTVariant = std::variant, WilsonResidualHelmholtzOverRT>; - /** Cubic EOS with advanced mixing rules, the EoS/aE method of Jaubert and co-workers diff --git a/include/teqp/models/multifluid/multifluid_activity.hpp b/include/teqp/models/multifluid/multifluid_activity.hpp new file mode 100644 index 00000000..7a287b35 --- /dev/null +++ b/include/teqp/models/multifluid/multifluid_activity.hpp @@ -0,0 +1,75 @@ +#pragma once + +#include "teqp/types.hpp" + +#include "teqp/models/multifluid.hpp" +#include "teqp/models/activity/activity_models.hpp" + +namespace teqp::multifluid::multifluid_activity { + +/** + Implementing the general approach of: + Jaeger et al., + A theoretically based departure function for multi-fluid mixture models, + https://doi.org/10.1016/j.fluid.2018.04.015 +*/ +class MultifluidPlusActivity{ +private: + using multifluid_t = decltype(multifluidfactory(nlohmann::json{})); + const multifluid_t m_multifluid; + const ResidualHelmholtzOverRTVariant m_activity; + const std::vector b; + const double u; +public: + MultifluidPlusActivity(const nlohmann::json &spec) : + m_multifluid(multifluidfactory(spec.at("multifluid"))), + m_activity(ares_model_factory(spec.at("activity").at("aresmodel"))), + b(spec.at("activity").at("options").at("b")), + u(spec.at("activity").at("options").at("u")){} + + /// Calculate the dimensionless value of \f$g_{\rm GE}^{\rm E,R}/RT\f$ from the AC model + auto calc_gER_over_RT(double T, const Eigen::ArrayXd& molefrac) const { + return std::visit([T, &molefrac](const auto& mod){return mod(T, molefrac); }, m_activity); + } + + template + auto R(const VecType& molefrac) const { + return get_R_gas(); + } + + template + auto alphar_activity(const TType& T, const RhoType& rho, const MoleFractions& molefrac) const { + auto gER_over_RT = std::visit([T, &molefrac](const auto& mod){return mod(T, molefrac); }, m_activity); // dimensionless + if (static_cast(b.size()) != molefrac.size()){ + throw teqp::InvalidArgument("Size of mole fractions is incorrect"); + } + + auto bm = contiguous_dotproduct(b, molefrac); + + const auto& Tcvec = m_multifluid.redfunc.Tc; + const auto& vcvec = m_multifluid.redfunc.vc; + + auto rhor = m_multifluid.redfunc.get_rhor(molefrac); + auto Tr = m_multifluid.redfunc.get_Tr(molefrac); + auto tau = forceeval(Tr/T); + auto delta_ref = forceeval(1.0/(u*bm*rhor)); + + std::decay_t> summer = 0.0; + for (auto i = 0; i < molefrac.size(); ++i){ + auto delta_i_ref = forceeval(1.0/(u*b[i]/vcvec(i))); + auto tau_i = forceeval(Tcvec(i)/T); + summer += molefrac(i)*(m_multifluid.alphar_taudeltai(tau, delta_ref, i) - m_multifluid.alphar_taudeltai(tau_i, delta_i_ref, i)); + } + return forceeval(log(1.0+rho*bm)/log(1.0+1.0/u)*(gER_over_RT - summer)); + } + + template + auto alphar(const TType& T, const RhoType& rho, const MoleFractions& molefrac) const { + return forceeval( + m_multifluid.alphar(T, rho, molefrac) + + alphar_activity(T, rho, molefrac) + ); + } +}; + +}; // namespace teqp diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp index a80e5f46..a7b55da9 100644 --- a/include/teqp/types.hpp +++ b/include/teqp/types.hpp @@ -26,6 +26,8 @@ namespace teqp { #include "boost/multiprecision/cpp_bin_float.hpp" #endif +#include "teqp/exceptions.hpp" + // autodiff include #include #include @@ -225,6 +227,41 @@ namespace teqp { } #endif + /** + \brief Take the dot-product of two vector-like objects that have contiguous memory and support the .size() method + + This allows to mix and match 1D Eigen arrays and vectors and STL vectors + */ + auto contiguous_dotproduct(const auto& x, const auto&y){ + using x_t = std::decay_t; + using y_t = std::decay_t; + using ret_t = std::common_type_t; + if (static_cast(x.size()) != y.size()){ + throw teqp::InvalidArgument("Arguments to contiguous_dotproduct are not the same size"); + } + ret_t summer = 0.0; + for (auto i = 0U; i < x.size(); ++i){ + summer += x[i]*y[i]; + } + return summer; +// auto get_ptr = [](auto& x){ +// using x_t = std::decay_t; +// const x_t* x_0_ptr; +// if constexpr (is_eigen_impl::value){ +// x_0_ptr = x.data(); +// } +// else{ +// x_0_ptr = &(x[0]); +// } +// return x_0_ptr; +// }; +// // element-wise product and then a sum +// return ( +// Eigen::Map>(get_ptr(x), x.size()).template cast() +// * Eigen::Map>(get_ptr(y), y.size()).template cast() +// ).sum(); + } + }; // namespace teqp #endif diff --git a/interface/CPP/model_factory_multifluid_activity.cpp b/interface/CPP/model_factory_multifluid_activity.cpp new file mode 100644 index 00000000..d4cb2421 --- /dev/null +++ b/interface/CPP/model_factory_multifluid_activity.cpp @@ -0,0 +1,11 @@ +#include "teqp/cpp/teqpcpp.hpp" +#include "teqp/cpp/deriv_adapter.hpp" +#include "teqp/models/multifluid/multifluid_activity.hpp" + +namespace teqp{ + namespace cppinterface{ + std::unique_ptr make_multifluid_activity(const nlohmann::json &j){ + return teqp::cppinterface::adapter::make_owned(multifluid::multifluid_activity::MultifluidPlusActivity(j)); + } + } +} diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp index 7dc365c4..22a20db1 100644 --- a/interface/CPP/teqp_impl_factory.cpp +++ b/interface/CPP/teqp_impl_factory.cpp @@ -24,6 +24,7 @@ namespace teqp { std::unique_ptr make_LKP(const nlohmann::json &); std::unique_ptr make_multifluid(const nlohmann::json &); std::unique_ptr make_multifluid_association(const nlohmann::json &); + std::unique_ptr make_multifluid_activity(const nlohmann::json &); std::unique_ptr make_multifluid_ECS_HuberEly1994(const nlohmann::json &); std::unique_ptr make_AmmoniaWaterTillnerRoth(); std::unique_ptr make_LJ126_TholJPCRD2016(); @@ -78,6 +79,7 @@ namespace teqp { {"multifluid", [](const nlohmann::json& spec){ return make_multifluid(spec);}}, {"multifluid-ECS-HuberEly1994", [](const nlohmann::json& spec){ return make_multifluid_ECS_HuberEly1994(spec);}}, {"multifluid-association", [](const nlohmann::json& spec){ return make_multifluid_association(spec);}}, + {"multifluid-activity", [](const nlohmann::json& spec){ return make_multifluid_activity(spec);}}, {"AmmoniaWaterTillnerRoth", [](const nlohmann::json& /*spec*/){ return make_AmmoniaWaterTillnerRoth();}}, {"LJ126_TholJPCRD2016", [](const nlohmann::json& /*spec*/){ return make_LJ126_TholJPCRD2016();}}, diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp index f938608b..8b28fc8f 100644 --- a/interface/pybind11_wrapper.cpp +++ b/interface/pybind11_wrapper.cpp @@ -17,6 +17,7 @@ #include "teqp/algorithms/ancillary_builder.hpp" #include "teqp/models/multifluid_ecs_mutant.hpp" #include "teqp/models/multifluid_association.hpp" +#include "teqp/models/multifluid/multifluid_activity.hpp" #include "teqp/models/saft/genericsaft.hpp" #include "teqp/algorithms/phase_equil.hpp" @@ -152,6 +153,7 @@ using CPA_t = decltype(teqp::CPA::CPAfactory("")); const std::type_index CPA_i{std::type_index(typeid(CPA_t))}; const std::type_index genericSAFT_i{std::type_index(typeid(teqp::saft::genericsaft::GenericSAFT))}; const std::type_index MultiFluidAssociation_i{std::type_index(typeid(MultifluidPlusAssociation))}; +const std::type_index MultiFluidActivity_i{std::type_index(typeid(teqp::multifluid::multifluid_activity::MultifluidPlusActivity))}; /** At runtime we can add additional model-specific methods that only apply for a particular model. We take in a Python-wrapped @@ -319,6 +321,11 @@ void attach_model_specific_methods(py::object& obj){ return get_typed(o).get_association().get_assoc_calcs(T, rhomolar, molefrac); }, "self"_a, "T"_a, "rhomolar"_a, "molefrac"_a), obj)); } + else if (index == MultiFluidActivity_i){ + setattr("calc_gER_over_RT", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ + return get_typed(o).calc_gER_over_RT(T, molefrac); + }, "self"_a, "T"_a, "molefrac"_a), obj)); + } }; /// Instantiate "instances" of models (really wrapped Python versions of the models), and then attach all derivative methods diff --git a/teqp/fluiddata/dev/fluids/Cyclopentane.json b/teqp/fluiddata/dev/fluids/Cyclopentane.json index 5a12d17d..cccae6ae 100644 --- a/teqp/fluiddata/dev/fluids/Cyclopentane.json +++ b/teqp/fluiddata/dev/fluids/Cyclopentane.json @@ -200,9 +200,9 @@ "T_units": "K", "hmolar": 38057.35612816187, "hmolar_units": "J/mol", - "p": 4571200, + "p": 4582800, "p_units": "Pa", - "rhomolar": 3820.000000000001, + "rhomolar": 3920.000000000001, "rhomolar_units": "mol/m^3", "smolar": 87.99746047649523, "smolar_units": "J/mol/K" @@ -240,26 +240,26 @@ "acentric_units": "-", "alpha0": [ { - "a1": 3.2489131288, - "a2": 2.6444166315, + "a1": -0.3946233253, + "a2": 2.4918910143, "type": "IdealGasHelmholtzLead" }, { - "a": 0.96, + "a": 3.0, "type": "IdealGasHelmholtzLogTau" }, { "n": [ - 3.34, - 18.6, - 13.9, - 4.86 + 1.34, + 13.4, + 17.4, + 6.65 ], "t": [ - 0.2345032439615415, - 2.540451809583366, - 5.276322989134683, - 10.35722660830141 + 0.4494645509262878, + 2.3059485656218244, + 4.299226139294927, + 10.161807238333463 ], "type": "IdealGasHelmholtzPlanckEinstein" } @@ -291,73 +291,79 @@ 1 ], "n": [ - 0.0536938, - 1.60394, - -2.41244, - -0.474009, - 0.203482, - -0.965616, - -0.344543, - 0.353975, - -0.231373, - -0.0379099 + 0.0630928, + 1.50365, + -2.37099, + -0.484886, + 0.191843, + -0.835582, + -0.435929, + 0.545607, + -0.209741, + -0.0387635 ], "t": [ - 1, + 1.0, 0.29, - 0.8, - 1.14, - 0.5, - 2, + 0.85, + 1.185, + 0.45, + 2.28, + 1.8, 1.5, - 1, - 3.36, - 0.95 - ], + 2.9, + 0.93], "type": "ResidualHelmholtzPower" }, { "beta": [ - 1.15, - 1.61, - 0.66, - 2.72 + 0.63, + 2.8, + 0.5, + 0.95, + 0.23 ], "d": [ 1, 1, 3, - 3 + 3, + 2 ], "epsilon": [ - 0.68, - 0.97, - 0.84, - 0.66 + 0.684, + 0.7, + 0.77, + 0.625, + 0.42 ], "eta": [ - 0.82, - 1.19, - 0.79, - 1.52 + 0.86, + 0.85, + 0.86, + 1.53, + 5.13 ], "gamma": [ - 1.08, - 0.36, - 0.09, - 1.48 + 1.22, + 0.32, + 0.22, + 1.94, + 1.21 ], "n": [ - 0.867586, - -0.381827, - -0.108741, - -0.0976984 + 0.677674, + -0.137043, + -0.0852862, + -0.128085, + -0.00389381 ], "t": [ - 1, - 2.5, - 2.5, - 1.5 + 1.05, + 4.0, + 2.33, + 1.5, + 1.0 ], "type": "ResidualHelmholtzGaussian" } @@ -426,7 +432,7 @@ "hmolar_units": "J/mol", "p": 4571200.0, "p_units": "Pa", - "rhomolar": 3820.000000000001, + "rhomolar": 3920.000000000001, "rhomolar_units": "mol/m^3", "smolar": 87.98268208465164, "smolar_units": "J/mol/K" @@ -557,4 +563,4 @@ "type": "Chung" } } -} \ No newline at end of file +}