From 91c4d09c870cb7489f7239dc28b1af7db6d8ce1f Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Thu, 11 Mar 2021 18:47:22 -0500
Subject: [PATCH] Switch Hessian evaluation to autodiff
Tracing is now crazy fast
---
CMakeLists.txt | 1 +
include/teqp/derivs.hpp | 27 ++++++++++++++++++++++++++-
src/test_autodiff.cpp | 1 -
src/trace_critical.cpp | 5 ++++-
4 files changed, 31 insertions(+), 3 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bec5188..72594c0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,6 +26,7 @@ add_executable(trace "${CMAKE_CURRENT_SOURCE_DIR}/src/trace_critical.cpp")
if (MSVC)
target_sources(trace PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis")
endif()
+target_link_libraries(trace autodiff)
add_executable(multifluid "${CMAKE_CURRENT_SOURCE_DIR}/src/multifluid.cpp")
if (MSVC)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 3e34171..8d106aa 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -6,6 +6,11 @@
#include "MultiComplex/MultiComplex.hpp"
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+#include <autodiff/forward/dual/eigen.hpp>
+using namespace autodiff;
+
template <typename TType, typename ContainerType, typename FuncType>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
caller(const FuncType& f, TType T, const ContainerType& rho) {
@@ -100,7 +105,7 @@ template <> void setval<std::valarray<std::valarray<double>>, std::size_t, doubl
/***
* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
-*
+*
* Requires the use of multicomplex derivatives to calculate second partial derivatives
*/
template<typename Model, typename TType, typename RhoType>
@@ -108,6 +113,26 @@ auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {
// Double derivatives in each component's concentration
// N^N matrix (symmetric)
+ dual2nd u; // the output scalar u = f(x), evaluated together with Hessian below
+ VectorXdual2nd g;
+ VectorXdual2nd rhovecc(2); rhovecc << rho[0], rho[1];
+ auto hfunc = [&model, &T](const VectorXdual2nd& rho_) {
+ auto rhotot_ = std::accumulate(std::begin(rho_), std::end(rho_), (decltype(rho_[0]))0.0);
+ return eval(model.alphar(T, rho_) * model.R * T * rhotot_);
+ };
+ return autodiff::hessian(hfunc, wrt(rhovecc), at(rhovecc), u, g); // evaluate the function value u, its gradient, and its Hessian matrix H
+}
+
+/***
+* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
+*
+* Requires the use of multicomplex derivatives to calculate second partial derivatives
+*/
+template<typename Model, typename TType, typename RhoType>
+auto build_Psir_Hessian_mcx(const Model& model, const TType T, const RhoType& rho) {
+ // Double derivatives in each component's concentration
+ // N^N matrix (symmetric)
+
// Lambda function for getting Psir with multicomplex concentrations
auto func = [&model, &T](const std::vector<MultiComplex<double>>& rhovec) {
auto N = rhovec.size();
diff --git a/src/test_autodiff.cpp b/src/test_autodiff.cpp
index de0718e..7bf7057 100644
--- a/src/test_autodiff.cpp
+++ b/src/test_autodiff.cpp
@@ -7,7 +7,6 @@
#include "MultiComplex/MultiComplex.hpp"
-
// autodiff include
#include <autodiff/forward/dual.hpp>
#include <autodiff/forward/dual/eigen.hpp>
diff --git a/src/trace_critical.cpp b/src/trace_critical.cpp
index 268c022..284cb1d 100644
--- a/src/trace_critical.cpp
+++ b/src/trace_critical.cpp
@@ -14,17 +14,20 @@ void trace() {
auto T = Tc_K[0];
const auto dT = 1;
std::valarray<double> rhovec = { rhoc0, 0.0 };
+ auto tic0 = std::chrono::steady_clock::now();
for (auto iter = 0; iter < 1000; ++iter){
auto drhovecdT = get_drhovec_dT_crit(vdW, T, rhovec);
rhovec += drhovecdT*dT;
T += dT;
int rr = 0;
auto z0 = rhovec[0] / rhovec.sum();
- std::cout << z0 << "," << rhovec[0] << "," << T << "," << get_splus(vdW, T, rhovec) << std::endl;
+ //std::cout << z0 << "," << rhovec[0] << "," << T << "," << get_splus(vdW, T, rhovec) << std::endl;
if (z0 < 0) {
break;
}
}
+ auto tic1 = std::chrono::steady_clock::now();
+ std::cout << std::chrono::duration<double>(tic1 - tic0).count() << " s to trace" << std::endl;
int rr = 0;
}
int main() {
--
GitLab