From bbb33836faa5d9587fb512ae4cffb99737634e24 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Wed, 24 Feb 2021 09:45:25 -0500
Subject: [PATCH] More generic derivative functions
---
include/teqp/derivs.hpp | 24 ++++++++----------------
src/main.cpp | 12 ++++++------
2 files changed, 14 insertions(+), 22 deletions(-)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index c2bd36d..0cd4358 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -6,37 +6,29 @@ caller(const FuncType& f, TType T, const ContainerType& rho) {
return f(T, rho);
}
-/// Given a
+/// Given a function, use complex step derivatives to calculate the derivative with respect to the first variable
+/// which here is temperature
template <typename TType, typename ContainerType, typename FuncType>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-deriv1(const FuncType& f, TType T, const ContainerType& rho) {
+derivT(const FuncType& f, TType T, const ContainerType& rho) {
double h = 1e-100;
return f(std::complex<TType>(T, h), rho).imag() / h;
}
-template <typename TType, typename ContainerType, typename FuncType>
+/// Given a function, use complex step derivatives to calculate the derivative with respect to the first composition variable
+template <typename TType, typename ContainerType, typename FuncType, typename Integer>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-deriv2(const FuncType& f, TType T, const ContainerType& rho) {
+derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer j) {
double h = 1e-100;
using comtype = std::complex<ContainerType::value_type>;
std::valarray<comtype> rhocom(rho.size());
- rhocom[0] = comtype(rho[0], h);
- for (auto i = 1; i < rho.size(); ++i) {
+ for (auto i = 0; i < rho.size(); ++i) {
rhocom[i] = comtype(rho[i], 0.0);
}
+ rhocom[j] = comtype(rho[j], h);
return f(T, rhocom).imag() / h;
}
-template <typename TType, typename ContainerType, typename FuncType>
-typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-deriv3(const FuncType& f, TType T, const ContainerType& rho) {
- double h = 1e-100;
- using comtype = std::complex<ContainerType::value_type>;
- std::valarray<comtype> rhocom(rho.size());
- rhocom[0] = comtype(rho[0], 0.0);
- rhocom[1] = comtype(rho[1], h);
- return f(T, rhocom).imag() / h;
-}
template<typename Model, typename TType, typename RhoType>
auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {
diff --git a/src/main.cpp b/src/main.cpp
index 9a42ddc..5f58eea 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -33,8 +33,8 @@ void test_vdW() {
auto t21 = std::chrono::steady_clock::now();
auto Psir = vdW.Psir(T, rhovec);
- auto dPsirdrho0 = rhovec[0] * deriv2([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
- auto dPsirdrho1 = rhovec[1] * deriv3([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
+ auto dPsirdrho0 = rhovec[0] * derivrhoi([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 0);
+ auto dPsirdrho1 = rhovec[1] * derivrhoi([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 1);
auto pfromderiv = rho * R * T - Psir + dPsirdrho0 + dPsirdrho1;
auto t31 = std::chrono::steady_clock::now();
@@ -58,12 +58,12 @@ void test_vdwMix() {
auto t2 = std::chrono::steady_clock::now();
auto Psir = vdW.Psir(T, rhovec);
- auto dPsirdrho0 = rhovec[0] * deriv2([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
- auto dPsirdrho1 = rhovec[1] * deriv3([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
+ auto dPsirdrho0 = rhovec[0] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 0);
+ auto dPsirdrho1 = rhovec[1] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 1);
auto pfromderiv = rho*R*T - Psir + dPsirdrho0 + dPsirdrho1;
{
- auto term0 = rhovec[0] * deriv2([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec);
- auto term1 = rhovec[1] * deriv3([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec);
+ auto term0 = rhovec[0] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec, 0);
+ auto term1 = rhovec[1] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec, 1);
auto pr = (term0 + term1)*rhotot*R*T;
auto pfromderiv2 = rho*R*T + pr;
auto err2 = pfromderiv / pfromderiv2 - 1;
--
GitLab