From fa28bfcbc8c1d2313653384152fcc8ca8c127013 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 23 Jul 2024 16:54:16 -0400 Subject: [PATCH] Add a class for fitting of pure parameters, and example notebook --- doc/source/recipes/FitPureParameters.ipynb | 180 ++++++++++++++++ doc/source/recipes/index.rst | 3 +- .../algorithms/pure_param_optimization.hpp | 196 ++++++++++++++++++ interface/pybind11_wrapper.cpp | 46 ++++ src/boost_pool_test.cpp | 15 +- src/tests/catch_test_SAFTgeneric.cxx | 2 +- src/tests/catch_test_paramoptim.cxx | 126 +++++++++++ 7 files changed, 559 insertions(+), 9 deletions(-) create mode 100644 doc/source/recipes/FitPureParameters.ipynb create mode 100644 include/teqp/algorithms/pure_param_optimization.hpp create mode 100644 src/tests/catch_test_paramoptim.cxx diff --git a/doc/source/recipes/FitPureParameters.ipynb b/doc/source/recipes/FitPureParameters.ipynb new file mode 100644 index 00000000..c37434b8 --- /dev/null +++ b/doc/source/recipes/FitPureParameters.ipynb @@ -0,0 +1,180 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "0bace5e5", + "metadata": {}, + "outputs": [], + "source": [ + "import teqp\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import CoolProp.CoolProp as CP\n", + "import scipy.optimize\n", + "display(teqp.__version__)\n", + "\n", + "nonpolar = {\n", + " \"kind\": \"SAFT-VR-Mie\",\n", + " \"model\": {\n", + " \"coeffs\": [\n", + " {\n", + " \"name\": \"R32\",\n", + " \"BibTeXKey\": \"Bell\",\n", + " \"m\": 1.2476268271391935,\n", + " \"sigma_m\": 3.6080717234117107e-10,\n", + " \"epsilon_over_k\": 172.53065054286867,\n", + " \"lambda_r\": 14.634722358167384,\n", + " \"lambda_a\": 6\n", + " }\n", + " ]\n", + " }\n", + "}\n", + "template = { \n", + " 'kind': 'genericSAFT', \n", + " 'model': {\n", + " 'nonpolar': nonpolar\n", + " }\n", + "}\n", + "\n", + "# NOTE: '/' inside field names MUST be escaped as ~1; see https://datatracker.ietf.org/doc/html/rfc6901#section-3\n", + "pointers = [\n", + " '/model/nonpolar/model/coeffs/0/m',\n", + " '/model/nonpolar/model/coeffs/0/sigma_m',\n", + " '/model/nonpolar/model/coeffs/0/epsilon_over_k',\n", + " '/model/nonpolar/model/coeffs/0/lambda_r',\n", + " '/model/nonpolar/model/coeffs/0/lambda_a'\n", + "]\n", + "x0 = [1.5, 3e-10, 150, 19, 5.7]\n", + "bounds = [(1,5), (2e-10,5e-10), (100,400), (12,50), (5.1, 6.0)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5463428", + "metadata": {}, + "outputs": [], + "source": [ + "FLD = 'Methane'\n", + "ppo = teqp.PureParameterOptimizer(template, pointers)\n", + "\n", + "Ts = np.linspace(100, 170, 10)\n", + "\n", + "# Generate some artificial pseudo-experimental data from \n", + "# the reference EOS as implemented in CoolProp\n", + "for T in Ts:\n", + " pt = teqp.SatRhoLPWPoint()\n", + " rhoL, rhoV = [CP.PropsSI('Dmolar','Q',Q,'T',T,FLD) for Q in [0,1]]\n", + " p = CP.PropsSI('P','Q',0,'T',T,FLD)\n", + " w = CP.PropsSI('speed_of_sound','Q',0,'T',T,FLD)\n", + " pt.T = T\n", + " # Measurands (here, pseudo-experimental values)\n", + " pt.p_exp = p\n", + " pt.rhoL_exp = rhoL\n", + " pt.w_exp = w\n", + " \n", + " # Additional parameters\n", + " pt.rhoL_guess = rhoL\n", + " pt.rhoV_guess = rhoV\n", + " pt.R = 8.31446261815324\n", + " pt.M = CP.PropsSI('molemass',FLD)\n", + " AS = CP.AbstractState('HEOS',FLD)\n", + " AS.update(CP.DmolarT_INPUTS, 1e-10, T)\n", + " pt.Ao20 = AS.tau()**2*AS.d2alpha0_dTau2() # -cv0/R\n", + " \n", + " # Weights (multiplied by 100 to put on a percentage basis)\n", + " pt.weight_p = 1*100\n", + " pt.weight_rho = 1*100\n", + " pt.weight_w = 0.25*100\n", + " ppo.add_one_contribution(pt)\n", + "\n", + "def cost_function(x):\n", + "# return ppo.cost_function_threaded(x, 10) # This is an option if you have lots of threads\n", + " return ppo.cost_function(x)\n", + "\n", + "r = scipy.optimize.differential_evolution(cost_function, bounds=bounds, disp=True, maxiter=10000, popsize=8)\n", + "print(r)\n", + "x = r.x\n", + "model = teqp.make_model(ppo.build_JSON(x))\n", + "\n", + "Tc, rhoc = model.solve_pure_critical(400, 5000)\n", + "print(Tc, rhoc)\n", + "anc = teqp.build_ancillaries(model, Tc, rhoc, 0.5*Tc)\n", + "\n", + "def _get_SOS(model, T, rho, z, *, R, M, Ao20):\n", + " \"\"\" Helper function to calculate speed of sound \"\"\"\n", + " Ar0n = model.get_Ar02n(T, rho, z)\n", + " Ar01 = Ar0n[1]; Ar02 = Ar0n[2]\n", + " Ar11 = model.get_Ar11(T, rho, z)\n", + " Ar20 = model.get_Ar20(T, rho, z)\n", + "\n", + " # M*w^2/(R*T) where w is the speed of sound\n", + " # from the definition w = sqrt(dp/drho|s)\n", + " Mw2RT = 1 + 2*Ar01 + Ar02 - (1 + Ar01 - Ar11)**2/(Ao20 + Ar20)\n", + " if Mw2RT < 0:\n", + " return 1e6\n", + " w = (Mw2RT*R*T/M)**0.5\n", + " return w\n", + "\n", + "TcREF = CP.PropsSI('Tcrit', FLD)\n", + "Tsverify = np.linspace(100, TcREF*0.99999, 1000)\n", + "RHOL, RHOV, PPP, WWW = [],[],[],[]\n", + "z = np.array([1.0])\n", + "for T in Tsverify:\n", + " rhoL, rhoV = model.pure_VLE_T(T, anc.rhoL(T)*1.01, anc.rhoV(T)*0.9, 10)\n", + " pL = rhoL*model.get_R(z)*T*(1+model.get_Ar01(T, rhoL, z))\n", + " pV = rhoV*model.get_R(z)*T*(1+model.get_Ar01(T, rhoV, z))\n", + " RHOL.append(rhoL)\n", + " RHOV.append(rhoV)\n", + " PPP.append(pL)\n", + " AS = CP.AbstractState('HEOS',FLD)\n", + " AS.update(CP.DmolarT_INPUTS, 1e-10, T)\n", + " Ao20 = AS.d2alpha0_dTau2()*AS.tau()**2\n", + " WWW.append(_get_SOS(model, T, rhoL, z, R=8.314462618,M=CP.PropsSI('molemass',FLD),Ao20=Ao20))\n", + " \n", + "# Plot the T-rho for VLE\n", + "line, = plt.plot(np.array(RHOL), Tsverify)\n", + "plt.plot(np.array(RHOV), Tsverify, color=line.get_color())\n", + "for Q in [0, 1]:\n", + " D = CP.PropsSI('Dmolar','T',Tsverify,'Q',Q,FLD)\n", + " plt.plot(D, Tsverify, lw=2, color='k')\n", + "plt.gca().set(xlabel=r'$\\rho$ / mol/m$^3$', ylabel='$T$ / K')\n", + "\n", + "# And a deviation plot, much closer to the critical point\n", + "# than the fitting region\n", + "fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(10,5), sharex=True)\n", + "\n", + "ax1.plot(Tsverify, (np.array(PPP)/CP.PropsSI('P','T',Tsverify,'Q',0,FLD)-1)*100)\n", + "ax1.set(ylabel=r'$(p_{fit}/p_{\\rm pexp}-1)\\times 100$')\n", + "\n", + "ax2.plot(Tsverify, (np.array(RHOL)/CP.PropsSI('Dmolar','T',Tsverify,'Q',0,FLD)-1)*100)\n", + "ax2.set(ylabel=r'$(\\rho_{fit}/\\rho_{\\rm pexp}-1)\\times 100$')\n", + "\n", + "ax3.plot(Tsverify, (np.array(WWW)/CP.PropsSI('speed_of_sound','T',Tsverify,'Q',0,FLD)-1)*100)\n", + "ax3.set(ylabel=r'$(w_{fit}/w_{\\rm pexp}-1)\\times 100$', xlabel='$T$ / K')" + ] + } + ], + "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/recipes/index.rst b/doc/source/recipes/index.rst index 8806dbe1..35067a8f 100644 --- a/doc/source/recipes/index.rst +++ b/doc/source/recipes/index.rst @@ -8,4 +8,5 @@ The examples here demonstrate various capabilities of teqp to solve practical pr :caption: Contents: HeatPumpModel - \ No newline at end of file + IdealCurves + FitPureParameters \ No newline at end of file diff --git a/include/teqp/algorithms/pure_param_optimization.hpp b/include/teqp/algorithms/pure_param_optimization.hpp new file mode 100644 index 00000000..34c5e8d5 --- /dev/null +++ b/include/teqp/algorithms/pure_param_optimization.hpp @@ -0,0 +1,196 @@ +#pragma once + +#include "nlohmann/json.hpp" +#include "teqp/cpp/teqpcpp.hpp" +#include "teqp/exceptions.hpp" +#include "teqp/math/pow_templates.hpp" + +#include +#include + +#include + +namespace teqp::algorithms::pure_param_optimization { + +struct SatRhoLPoint{ + double T, rhoL_exp, rhoL_guess, rhoV_guess, weight=1.0; + auto check_fields() const{} + template + auto calculate_contribution(const Model& model) const{ + auto rhoLrhoV = model->pure_VLE_T(T, rhoL_guess, rhoV_guess, 10); + return std::abs(rhoLrhoV[0]-rhoL_exp)*weight; + } +}; + +#define stdstringify(s) std::string(#s) + +#define SatRhoLPPoint_optionalfields X(T) X(p_exp) X(rhoL_exp) X(rhoL_guess) X(rhoV_guess) +struct SatRhoLPPoint{ + #define X(field) std::optional field; + SatRhoLPPoint_optionalfields + #undef X + double weight_rho=1.0, weight_p=1.0, R=8.31446261815324; + Eigen::ArrayXd z = (Eigen::ArrayXd(1) << 1.0).finished(); + + auto check_fields() const{ + #define X(field) if (!field){ throw teqp::InvalidArgument("A field [" + stdstringify(field) + "] has not been initialized"); } + SatRhoLPPoint_optionalfields + #undef X + } + + template + auto calculate_contribution(const Model& model) const{ + auto rhoLrhoV = model->pure_VLE_T(T.value(), rhoL_guess.value(), rhoV_guess.value(), 10); + auto rhoL = rhoLrhoV[0]; + auto p = rhoL*R*T.value()*(1+model->get_Ar01(T.value(), rhoL, z)); +// std::cout << p << "," << p_exp << "," << (p_exp-p)/p_exp << std::endl; + return std::abs(rhoL-rhoL_exp.value())/rhoL_exp.value()*weight_rho + std::abs(p-p_exp.value())/p_exp.value()*weight_p; + } +}; + +#define SatRhoLPWPoint_optionalfields X(T) X(p_exp) X(rhoL_exp) X(w_exp) X(rhoL_guess) X(rhoV_guess) X(Ao20) X(M) X(R) +struct SatRhoLPWPoint{ + + #define X(field) std::optional field; + SatRhoLPWPoint_optionalfields + #undef X + double weight_rho=1.0, weight_p=1.0, weight_w=1.0; + Eigen::ArrayXd z = (Eigen::ArrayXd(1) << 1.0).finished(); + + auto check_fields() const{ + #define X(field) if (!field){ throw teqp::InvalidArgument("A field [" + stdstringify(field) + "] has not been initialized"); } + SatRhoLPWPoint_optionalfields + #undef X + } + + template + auto calculate_contribution(const Model& model) const{ + + auto rhoLrhoV = model->pure_VLE_T(T.value(), rhoL_guess.value(), rhoV_guess.value(), 10); + + // First part, density + auto rhoL = rhoLrhoV[0]; + + auto Ar0n = model->get_Ar02n(T.value(), rhoL, z); + double Ar01 = Ar0n[1], Ar02 = Ar0n[2]; + auto Ar11 = model->get_Ar11(T.value(), rhoL, z); + auto Ar20 = model->get_Ar20(T.value(), rhoL, z); + + // Second part, presure + auto p = rhoL*R.value()*T.value()*(1+Ar01); + + // Third part, speed of sound + // + // M*w^2/(R*T) where w is the speed of sound + // from the definition w = sqrt(dp/drho|s) + double Mw2RT = 1 + 2*Ar01 + Ar02 - POW2(1.0 + Ar01 - Ar11)/(Ao20.value() + Ar20); + double w = sqrt(Mw2RT*R.value()*T.value()/M.value()); + +// std::cout << p << "," << p_exp << "," << (p_exp-p)/p_exp << std::endl; + double cost_rhoL = std::abs(rhoL-rhoL_exp.value())/rhoL_exp.value()*weight_rho; + double cost_p = std::abs(p-p_exp.value())/p_exp.value()*weight_p; + double cost_w = std::abs(w-w_exp.value())/w_exp.value()*weight_w; + return cost_rhoL + cost_p + cost_w; + } +}; +using PureOptimizationContribution = std::variant; + +class PureParameterOptimizer{ +private: + auto make_pointers(const std::vector>>& pointerstrs){ + std::vector> pointers_; + for (auto & pointer: pointerstrs){ + if (std::holds_alternative(pointer)){ + const std::string& s= std::get(pointer); + pointers_.emplace_back(1, nlohmann::json::json_pointer(s)); + } + else{ + std::vector buffer; + for (const auto& s : std::get>(pointer)){ + buffer.emplace_back(s); + } + pointers_.push_back(buffer); + } + } + return pointers_; + } +public: + const nlohmann::json jbase; + std::vector> pointers; + std::vector contributions; + + PureParameterOptimizer(const nlohmann::json jbase, const std::vector>>& pointerstrs) : jbase(jbase), pointers(make_pointers(pointerstrs)){} + + void add_one_contribution(const PureOptimizationContribution& cont){ + std::visit([](const auto&c){c.check_fields();}, cont); + contributions.push_back(cont); + } + + auto prepare_helper_models(const auto& model) const { + return std::vector{1, 3}; + } + + template + nlohmann::json build_JSON(const T& x) const{ + + if (x.size() != pointers.size()){ + throw teqp::InvalidArgument("sizes don't match"); + }; + nlohmann::json j = jbase; + for (auto i = 0; i < x.size(); ++i){ + for (const auto& ptr : pointers[i]){ + j.at(ptr) = x[i]; + } + } + return j; + } + + template + auto prepare(const T& x) const { + auto j = build_JSON(x); + auto model = teqp::cppinterface::make_model(j); + auto helpers = prepare_helper_models(model); + return std::make_tuple(std::move(model), helpers); + } + + template + auto cost_function(const T& x) const{ + auto [model, helpers] = prepare(x); + double cost = 0.0; + for (const auto& contrib : contributions){ + cost += std::visit([&model](const auto& c){ return c.calculate_contribution(model); }, contrib); + } + if (!std::isfinite(cost)){ + return 1e30; + } + return cost; + } + + template + auto cost_function_threaded(const T& x, std::size_t Nthreads) { + boost::asio::thread_pool pool{Nthreads}; // Nthreads in the pool + const auto [model, helpers] = prepare(x); + std::valarray buffer(contributions.size()); + std::size_t i = 0; + for (const auto& contrib : contributions){ + auto& dest = buffer[i]; + auto payload = [&model, &dest, contrib] (){ + dest = std::visit([&model](const auto& c){ return c.calculate_contribution(model); }, contrib); + if (!std::isfinite(dest)){ dest = 1e30; } + }; + boost::asio::post(pool, payload); + i++; + } + pool.join(); + double summer = 0.0; + for (auto i = 0; i < contributions.size(); ++i){ +// std::cout << buffer[i] << std::endl; + summer += buffer[i]; + } +// std::cout << summer << std::endl; + return summer; + } + +}; + +} diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp index 69d1a799..13e3ad2e 100644 --- a/interface/pybind11_wrapper.cpp +++ b/interface/pybind11_wrapper.cpp @@ -18,6 +18,9 @@ #include "teqp/models/multifluid_ecs_mutant.hpp" #include "teqp/models/saft/genericsaft.hpp" +#include "teqp/algorithms/pure_param_optimization.hpp" +using namespace teqp::algorithms::pure_param_optimization; + namespace py = pybind11; using namespace py::literals; @@ -568,6 +571,49 @@ void init_teqp(py::module& m) { .def("get_rho", &NRIterator::get_rho) ; + py::class_(m, "SatRhoLPoint") + .def(py::init<>()) + .def_readwrite("T", &SatRhoLPoint::T) + .def_readwrite("rhoL_exp", &SatRhoLPoint::rhoL_exp) + .def_readwrite("rhoL_guess", &SatRhoLPoint::rhoL_guess) + .def_readwrite("rhoV_guess", &SatRhoLPoint::rhoV_guess) + .def_readwrite("weight", &SatRhoLPoint::weight) + ; + py::class_(m, "SatRhoLPPoint") + .def(py::init<>()) + .def_readwrite("T", &SatRhoLPPoint::T) + .def_readwrite("rhoL_exp", &SatRhoLPPoint::rhoL_exp) + .def_readwrite("p_exp", &SatRhoLPPoint::p_exp) + .def_readwrite("rhoL_guess", &SatRhoLPPoint::rhoL_guess) + .def_readwrite("rhoV_guess", &SatRhoLPPoint::rhoV_guess) + .def_readwrite("weight_rho", &SatRhoLPPoint::weight_rho) + .def_readwrite("weight_p", &SatRhoLPPoint::weight_p) + .def_readwrite("R", &SatRhoLPPoint::R) + ; + py::class_(m, "SatRhoLPWPoint") + .def(py::init<>()) + .def_readwrite("T", &SatRhoLPWPoint::T) + .def_readwrite("rhoL_exp", &SatRhoLPWPoint::rhoL_exp) + .def_readwrite("p_exp", &SatRhoLPWPoint::p_exp) + .def_readwrite("w_exp", &SatRhoLPWPoint::w_exp) + .def_readwrite("R", &SatRhoLPWPoint::R) + .def_readwrite("Ao20", &SatRhoLPWPoint::Ao20) + .def_readwrite("M", &SatRhoLPWPoint::M) + .def_readwrite("rhoL_guess", &SatRhoLPWPoint::rhoL_guess) + .def_readwrite("rhoV_guess", &SatRhoLPWPoint::rhoV_guess) + .def_readwrite("weight_rho", &SatRhoLPWPoint::weight_rho) + .def_readwrite("weight_p", &SatRhoLPWPoint::weight_p) + .def_readwrite("weight_w", &SatRhoLPWPoint::weight_w) + ; + + py::class_(m, "PureParameterOptimizer") + .def(py::init>>&>()) + .def("cost_function", &PureParameterOptimizer::cost_function) + .def("cost_function_threaded", &PureParameterOptimizer::cost_function_threaded) + .def("build_JSON", &PureParameterOptimizer::build_JSON) + .def("add_one_contribution", &PureParameterOptimizer::add_one_contribution) + ; + // // Some functions for timing overhead of interface // m.def("___mysummer", [](const double &c, const Eigen::ArrayXd &x) { return c*x.sum(); }); // using RAX = Eigen::Ref; diff --git a/src/boost_pool_test.cpp b/src/boost_pool_test.cpp index f13164e5..2f78113b 100644 --- a/src/boost_pool_test.cpp +++ b/src/boost_pool_test.cpp @@ -57,7 +57,7 @@ int evaluation() { auto rhovecV = (Eigen::ArrayXd(2) << 0.0, 2.20225704).finished(); using MultiFluid = decltype(teqp::build_multifluid_model(std::vector{"", ""}, "", "")); - std::vector models(60, teqp::build_multifluid_model({ "CarbonDioxide", "Water" }, "../mycp")); + std::vector models(40, teqp::build_multifluid_model({ "CarbonDioxide", "Water" }, "../teqp/fluiddata")); std::vector outputs(models.size()); auto serial = [&]() { @@ -70,15 +70,16 @@ int evaluation() { }; auto parallel = [&]() { // Launch the pool with four threads. - boost::asio::thread_pool pool(4); + boost::asio::thread_pool pool(40); std::size_t i = 0; for (auto& model : models) { + auto &model_ = models[0]; auto& o = outputs[i]; // Submit a lambda object to the pool. - boost::asio::post(pool, [&model, &o, &T, &rhovecL, &rhovecV]() { - using TDX = teqp::TDXDerivatives; - o = teqp::trace_VLE_isotherm_binary(model, T, rhovecL, rhovecV).dump(); + boost::asio::post(pool, [&model_, &o, &T, &rhovecL, &rhovecV]() { + using TDX = teqp::TDXDerivatives; + o = teqp::trace_VLE_isotherm_binary(model_, T, rhovecL, rhovecV).dump(); }); i++; } @@ -100,6 +101,6 @@ int evaluation() { int main() { - simple(); +// simple(); evaluation(); -} \ No newline at end of file +} diff --git a/src/tests/catch_test_SAFTgeneric.cxx b/src/tests/catch_test_SAFTgeneric.cxx index 8e829156..8b9265fa 100644 --- a/src/tests/catch_test_SAFTgeneric.cxx +++ b/src/tests/catch_test_SAFTgeneric.cxx @@ -96,7 +96,7 @@ TEST_CASE("Benchmark generic PC-SAFT+Association model", "[SAFTgeneric]"){ } -auto Dufal_contents = R"( +static auto Dufal_contents = R"( { "nonpolar": { "kind": "SAFT-VR-Mie", diff --git a/src/tests/catch_test_paramoptim.cxx b/src/tests/catch_test_paramoptim.cxx new file mode 100644 index 00000000..a492d0da --- /dev/null +++ b/src/tests/catch_test_paramoptim.cxx @@ -0,0 +1,126 @@ +#include +#include +using Catch::Approx; +#include + +#include "teqp/cpp/teqpcpp.hpp" +#include "teqp/cpp/deriv_adapter.hpp" +using namespace teqp; + +#include "teqp/json_tools.hpp" +#include "teqp/algorithms/pure_param_optimization.hpp" +#include "test_common.in" +#include "teqp/models/multifluid_ancillaries.hpp" + +using namespace teqp::algorithms::pure_param_optimization; + +static auto Dufal_contents = R"( +{ + "nonpolar": { + "kind": "SAFT-VR-Mie", + "model": { + "coeffs": [ + { + "name": "Water", + "BibTeXKey": "Dufal-2015", + "m": 1.0, + "sigma_Angstrom": 3.0555, + "epsilon_over_k": 418.00, + "lambda_r": 35.823, + "lambda_a": 6.0 + } + ] + } + } +} +)"_json; + +static auto Dufale_contents = R"( +{ + "nonpolar": { + "kind": "SAFT-VR-Mie", + "model": { + "coeffs": [ + { + "name": "Water", + "BibTeXKey": "Dufal-2015", + "m": 1.0, + "sigma_Angstrom": 3.0555, + "epsilon_over_k": 418.00, + "lambda_r": 35.823, + "lambda_a": 6.0 + } + ] + } + }, + "association": { + "kind": "canonical", + "model": { + "b / m^3/mol": [0.0000145], + "beta": [0.0692], + "Delta_rule": "CR1", + "epsilon / J/mol": [16655.0], + "molecule_sites": [["e","e","H","H"]], + "options": {"radial_dist": "CS"} + } + } +} +)"_json; + +auto gen_sat_data(){ + auto contents = R"({ + "kind": "multifluid", + "model": { + "components": ["R32"], + "root": "" + } + })"_json; + contents["model"]["root"] = FLUIDDATAPATH; + auto model = teqp::cppinterface::make_model(contents); + + auto jancillaries = load_a_JSON_file(FLUIDDATAPATH+"/dev/fluids/R32.json").at("ANCILLARIES"); + auto anc = teqp::MultiFluidVLEAncillaries(jancillaries); + + std::vector contribs; + for (auto T = 270.0; T < 310; T += 0.4){ + SatRhoLPoint pt; + pt.T = T; + auto rhoLrhoV = model->pure_VLE_T(T, anc.rhoL(T), anc.rhoV(T), 10); + pt.rhoL_exp = rhoLrhoV[0]; + pt.rhoL_guess = rhoLrhoV[0]; + pt.rhoV_guess = rhoLrhoV[1]; + contribs.push_back(pt); + } + return contribs; +} + +TEST_CASE("Benchmark param optimization", "[paramoptim]"){ + + nlohmann::json maincontents = { + {"kind", "genericSAFT"}, + {"model", Dufal_contents} + }; + + std::vector>> pointers = {"/model/nonpolar/model/coeffs/0/m"}; + PureParameterOptimizer ppo(maincontents, pointers); + for (auto pt : gen_sat_data()){ + ppo.add_one_contribution(pt); + } + + std::vector xx = {1.3}; + BENCHMARK("build JSON"){ + return ppo.build_JSON(xx); + }; + BENCHMARK("model building"){ + return ppo.prepare(xx); + }; + BENCHMARK("cost_function evaluation"){ + return ppo.cost_function(xx); + }; + BENCHMARK("cost_function evaluation threaded"){ + return ppo.cost_function_threaded(xx, 6); + }; + SECTION("check cost functions"){ + CHECK(ppo.cost_function_threaded(xx, 6) == ppo.cost_function(xx)); + } +}