diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 7f0ce9d38bbc334ce6961f32d900715ce1c4f021..154960f1ff95cde6ceadba361f8ebce382cfa9ba 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -81,6 +81,7 @@ struct wrt_helper {
}
};
+enum class AlphaWrapperOption {residual, idealgas};
/**
* \brief This class is used to wrap a model that exposes the generic
* functions alphar, alphaig, etc., and allow the desired member function to be
@@ -90,14 +91,14 @@ struct wrt_helper {
* require the argument types to be known, and they are not known in this case
* so we give the hard work of managing the argument types to the compiler
*/
-template<int i, class Model>
+template<AlphaWrapperOption o, class Model>
struct AlphaCallWrapper {
const Model& m_model;
AlphaCallWrapper(const Model& model) : m_model(model) {};
template <typename ... Args>
auto alpha(const Args& ... args) const {
- if constexpr (i == 0) {
+ if constexpr (o == AlphaWrapperOption::residual) {
// The alphar method is REQUIRED to be implemented by all
// models, so can just call it via perfect fowarding
return m_model.alphar(std::forward<const Args>(args)...);
@@ -212,7 +213,7 @@ struct TDXDerivatives {
*/
template<int iT, int iD, ADBackends be>
static auto get_Arxy(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
- auto wrapper = AlphaCallWrapper<0, decltype(model)>(model);
+ auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
if constexpr (iT == 0 && iD == 0) {
return wrapper.alpha(T, rho, molefrac);
}
@@ -231,7 +232,7 @@ struct TDXDerivatives {
*/
template<int iT, int iD, ADBackends be>
static auto get_Aigxy(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
- auto wrapper = AlphaCallWrapper<1, decltype(model)>(model);
+ auto wrapper = AlphaCallWrapper<AlphaWrapperOption::idealgas, decltype(model)>(model);
if constexpr (iT == 0 && iD == 0) {
return wrapper.alpha(T, rho, molefrac);
}
@@ -338,7 +339,7 @@ struct TDXDerivatives {
*/
template<int iT, ADBackends be = ADBackends::autodiff>
static auto get_Arn0(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
- auto wrapper = AlphaCallWrapper<0, decltype(model)>(model);
+ auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
return get_Agenn0<iT, be>(wrapper, T, rho, molefrac);
}
@@ -352,7 +353,7 @@ struct TDXDerivatives {
*/
template<int iD, ADBackends be = ADBackends::autodiff>
static auto get_Ar0n(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
- auto wrapper = AlphaCallWrapper<0, decltype(model)>(model);
+ auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
return get_Agen0n<iD, be>(wrapper, T, rho, molefrac);
}
@@ -904,4 +905,30 @@ struct IsochoricDerivatives{
}
};
+template<int Nderivsmax, AlphaWrapperOption opt>
+class DerivativeHolderSquare{
+
+public:
+ Eigen::Array<double, Nderivsmax+1, Nderivsmax+1> derivs;
+
+ template<typename Model, typename Scalar, typename VecType>
+ DerivativeHolderSquare(const Model& model, const Scalar& T, const Scalar& rho, const VecType& z) {
+ using tdx = TDXDerivatives<decltype(model), Scalar, VecType>;
+ static_assert(Nderivsmax == 2, "It's gotta be 2 for now");
+ AlphaCallWrapper<opt, Model> wrapper(model);
+
+ auto AX02 = tdx::template get_Agen0n<2>(wrapper, T, rho, z);
+ derivs(0, 0) = AX02[0];
+ derivs(0, 1) = AX02[1];
+ derivs(0, 2) = AX02[2];
+
+ auto AX20 = tdx::template get_Agenn0<2>(wrapper, T, rho, z);
+ derivs(0, 0) = AX20[0];
+ derivs(1, 0) = AX20[1];
+ derivs(2, 0) = AX20[2];
+
+ derivs(1, 1) = tdx::template get_Agenxy<1,1>(wrapper, T, rho, z);
+ }
+};
+
}; // namespace teqp
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 61e0d08d3bc013a8eb1c32196e2074232b42ed04..dcd96be9df907fdd75ae212bfabf7c3647832b29 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -164,6 +164,7 @@ void init_teqp(py::module& m) {
auto alphaig = py::class_<IdealHelmholtz>(m, "IdealHelmholtz").def(py::init<const nlohmann::json&>());
alphaig.def_static("convert_CoolProp_format", [](const std::string &path, int index){return convert_CoolProp_idealgas(path, index);});
add_ig_derivatives<IdealHelmholtz>(m, alphaig);
+ alphaig.def("get_deriv_mat2", [](const IdealHelmholtz &ig, double T, double rho, const Eigen::ArrayXd& z){return DerivativeHolderSquare<2, AlphaWrapperOption::idealgas>(ig, T, rho, z).derivs;});
add_vdW(m);
add_PCSAFT(m);
diff --git a/interface/pybind11_wrapper.hpp b/interface/pybind11_wrapper.hpp
index e223e712795ed47c0c81c0d109bee5aad9629793..a9579fc05c3ec695c6efe3f9a15587eb917462e2 100644
--- a/interface/pybind11_wrapper.hpp
+++ b/interface/pybind11_wrapper.hpp
@@ -106,5 +106,7 @@ void add_derivatives(py::module &m, Wrapper &cls) {
cls.def("get_Ar05n", &(tdx::template get_Ar0n<5>), py::arg("T"), py::arg("rho"), py::arg("molefrac").noconvert());
cls.def("get_Ar06n", &(tdx::template get_Ar0n<6>), py::arg("T"), py::arg("rho"), py::arg("molefrac").noconvert());
cls.def("get_neff", &(tdx::template get_neff<ADBackends::autodiff>), py::arg("T"), py::arg("rho"), py::arg("molefrac").noconvert());
+
+ cls.def("get_deriv_mat2", [](const Model &model, double T, double rho, const Eigen::ArrayXd& z){return DerivativeHolderSquare<2, AlphaWrapperOption::residual>(model, T, rho, z).derivs;});
}
diff --git a/src/bench_derivbox.cpp b/src/bench_derivbox.cpp
index 81ef2fff4ef451613a22d91d0dc828643e203b4b..233b18d1c4abc22af0be4f6576af82e9f69a6c52 100644
--- a/src/bench_derivbox.cpp
+++ b/src/bench_derivbox.cpp
@@ -8,32 +8,6 @@
using namespace teqp;
-template<int Nderivsmax, int k>
-class DerivativeHolderSquare{
-
-public:
- Eigen::Array<double, Nderivsmax+1, Nderivsmax+1> derivs;
-
- template<typename Model, typename Scalar, typename VecType>
- DerivativeHolderSquare(const Model& model, const Scalar& T, const Scalar& rho, const VecType& z) {
- using tdx = TDXDerivatives<decltype(model), Scalar, VecType>;
- static_assert(Nderivsmax == 2, "It's gotta be 2");
- AlphaCallWrapper<k, Model> wrapper(model);
-
- auto AX02 = tdx::template get_Agen0n<2>(wrapper, T, rho, z);
- derivs(0, 0) = AX02[0];
- derivs(0, 1) = AX02[1];
- derivs(0, 2) = AX02[2];
-
- auto AX20 = tdx::template get_Agenn0<2>(wrapper, T, rho, z);
- derivs(0, 0) = AX20[0];
- derivs(1, 0) = AX20[1];
- derivs(2, 0) = AX20[2];
-
- derivs(1, 1) = tdx::template get_Agenxy<1,1>(wrapper, T, rho, z);
- }
-};
-
TEST_CASE("multifluid derivatives", "[mf]")
{
std::vector<std::string> names = { "Propane" };
@@ -52,10 +26,10 @@ TEST_CASE("multifluid derivatives", "[mf]")
};
BENCHMARK("All residual derivatives needed for first derivatives of h,s,u,p w.r.t. T&rho") {
- return DerivativeHolderSquare<2,0>(model, T, rho, z).derivs;
+ return DerivativeHolderSquare<2,AlphaWrapperOption::residual>(model, T, rho, z).derivs;
};
BENCHMARK("All ideal-gas derivatives needed for first derivatives of h,s,u,p w.r.t. T&rho") {
- return DerivativeHolderSquare<2,1>(aig, T, rho, z).derivs;
+ return DerivativeHolderSquare<2,AlphaWrapperOption::idealgas>(aig, T, rho, z).derivs;
};
}
diff --git a/src/tests/catch_test_alphaig.cxx b/src/tests/catch_test_alphaig.cxx
index a5230d242cce982b52cabe5f99bc45ebf94af092..8473c3978916a981075563db97a901a31fdaf522 100644
--- a/src/tests/catch_test_alphaig.cxx
+++ b/src/tests/catch_test_alphaig.cxx
@@ -29,7 +29,7 @@ TEST_CASE("alphaig derivative", "[alphaig]") {
j.push_back(demo_pure_term(a_1, a_2));
IdealHelmholtz ih(j);
auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
- auto wih = AlphaCallWrapper<1, decltype(ih)>(ih);
+ auto wih = AlphaCallWrapper<AlphaWrapperOption::idealgas, decltype(ih)>(ih);
wih.alpha(T, rho, molefrac);
using tdx = TDXDerivatives<decltype(ih), double, Eigen::ArrayXd>;
SECTION("All cross derivatives should be zero") {
@@ -59,8 +59,10 @@ TEST_CASE("Ammonia derivative", "[alphaig][NH3]") {
nlohmann::json j = {{ {"R", 8.31446261815324}, {"terms", j0terms} }};
IdealHelmholtz ih(j);
auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
- auto wih = AlphaCallWrapper<1, decltype(ih)>(ih);
+ auto wih = AlphaCallWrapper<AlphaWrapperOption::idealgas, decltype(ih)>(ih);
auto calc = wih.alpha(T, rho, molefrac);
auto expected = -5.3492909452728545;
REQUIRE(calc == Approx(expected));
+
+ DerivativeHolderSquare<2, AlphaWrapperOption::idealgas> dhs(ih, T, rho, molefrac);
}