From 4d3f8fd2c8f1239eea1331f3269e21145e3f70f8 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Tue, 8 Nov 2022 15:35:20 -0700
Subject: [PATCH] Add Newton-Raphson iterator for homogeneous iterative
calculations
---
include/teqp/algorithms/iteration.hpp | 72 +++++++++++++++++++++++++++
interface/pybind11_wrapper.cpp | 13 +++++
src/bench_derivbox.cpp | 36 ++++++++++++--
3 files changed, 117 insertions(+), 4 deletions(-)
create mode 100644 include/teqp/algorithms/iteration.hpp
diff --git a/include/teqp/algorithms/iteration.hpp b/include/teqp/algorithms/iteration.hpp
new file mode 100644
index 0000000..9d05044
--- /dev/null
+++ b/include/teqp/algorithms/iteration.hpp
@@ -0,0 +1,72 @@
+#pragma once
+
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/cpp/derivs.hpp"
+
+namespace teqp {
+namespace iteration {
+
+using namespace cppinterface;
+/**
+ A class for doing Newton-Raphson steps to solve for two unknown thermodynamic variables
+
+ */
+class NRIterator{
+private:
+ const std::shared_ptr<AbstractModel> ar, aig;
+ const std::vector<char> vars;
+ const Eigen::Ref<const Eigen::ArrayXd> vals;
+ double T, rho;
+ const Eigen::Ref<const Eigen::ArrayXd> z;
+
+public:
+ NRIterator(const std::shared_ptr<AbstractModel> &ar, const std::shared_ptr<AbstractModel> &aig, const std::vector<char>& vars, const Eigen::Ref<const Eigen::ArrayXd>& vals, double T, double rho, const Eigen::Ref<const Eigen::ArrayXd>& z) : ar(ar), aig(aig), vars(vars), vals(vals), T(T), rho(rho), z(z){}
+
+ /// Return the variables that are being used in the iteration
+ std::vector<char> get_vars() const { return vars; }
+ /// Return the target values to be obtained
+ Eigen::ArrayXd get_vals() const { return vals; }
+ /// Return the current temperature
+ auto get_T() const { return T; }
+ /// Return the current molar density
+ auto get_rho() const { return rho; }
+ /// Return the current mole fractions
+ Eigen::ArrayXd get_molefrac() const { return z; }
+
+ /** Do the calculations needed for the step and return the step and the other data
+ * In C++, the copy will be elided (the return value will be moved)
+ * \param T Temperature
+ * \param rhomolar Molar density
+ */
+ auto calc_step(double T, double rho){
+ auto Ar = ar->get_deriv_mat2(T, rho, z);
+ auto Aig = aig->get_deriv_mat2(T, rho, z);
+ auto R = ar->get_R(z);
+ auto im = build_iteration_Jv(vars, Ar, Aig, R, T, rho, z);
+ return std::make_tuple((im.J.matrix().colPivHouseholderQr().solve((-(im.v-vals)).matrix())).eval(), im);
+ }
+
+ /// Take one step, return the residuals
+ auto take_step(){
+ auto [dx, im] = calc_step(T, rho);
+ T += dx(0);
+ rho += dx(1);
+ return (im.v-vals).eval();
+ }
+
+ /** Take a given number of steps
+ * \param N The number of steps to take
+ */
+ auto take_steps(int N){
+ if (N <= 0){
+ throw teqp::InvalidArgument("N must be greater than 0");
+ }
+ for (auto i = 0; i < N; ++i){
+ take_step();
+ }
+ }
+};
+
+
+}
+}
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 06588c6..6e6d571 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -11,6 +11,7 @@
#include "teqp/derivs.hpp"
#include "teqp/cpp/teqpcpp.hpp"
#include "teqp/models/multifluid_ancillaries.hpp"
+#include "teqp/algorithms/iteration.hpp"
namespace py = pybind11;
using namespace py::literals;
@@ -353,6 +354,18 @@ void init_teqp(py::module& m) {
m.def("_make_model", &teqp::cppinterface::make_model);
m.def("attach_model_specific_methods", &attach_model_specific_methods);
+ using namespace teqp::iteration;
+ py::class_<NRIterator>(m, "NRIterator")
+ .def(py::init<const std::shared_ptr<AbstractModel> &, const std::shared_ptr<AbstractModel> &, const std::vector<char>&, const Eigen::Ref<const Eigen::ArrayXd>&, double, double, const Eigen::Ref<const Eigen::ArrayXd>&>())
+ .def("take_step", &NRIterator::take_step)
+ .def("take_steps", &NRIterator::take_steps)
+ .def("get_vars", &NRIterator::get_vars)
+ .def("get_vals", &NRIterator::get_vals)
+ .def("get_molefrac", &NRIterator::get_molefrac)
+ .def("get_T", &NRIterator::get_T)
+ .def("get_rho", &NRIterator::get_rho)
+ ;
+
// // 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<const Eigen::ArrayXd>;
diff --git a/src/bench_derivbox.cpp b/src/bench_derivbox.cpp
index e79ffc7..7223d7e 100644
--- a/src/bench_derivbox.cpp
+++ b/src/bench_derivbox.cpp
@@ -6,11 +6,14 @@
#include "teqp/derivs.hpp"
#include "teqp/ideal_eosterms.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/algorithms/iteration.hpp"
+
using namespace teqp;
TEST_CASE("multifluid derivatives", "[mf]")
{
- std::vector<std::string> names = { "Propane" };
+ std::vector<std::string> names = { "Ethane" };
auto model = build_multifluid_model(names, "../mycp");
double T = 300, rho = 2;
@@ -21,6 +24,15 @@ TEST_CASE("multifluid derivatives", "[mf]")
nlohmann::json jigs = nlohmann::json::array(); jigs.push_back(jig);
auto aig = teqp::IdealHelmholtz(jigs);
+ auto amm = teqp::cppinterface::make_multifluid_model(names, "../mycp");
+ auto aigg = teqp::cppinterface::make_model({{"kind","IdealHelmholtz"}, {"model", jigs}});
+
+ std::vector<char> vars = {'P', 'S'};
+ const auto vals = (Eigen::ArrayXd(2) << 300.0, 400.0).finished();
+
+ Eigen::Ref<const Eigen::ArrayXd> rvals = vals, rz = z;
+ teqp::iteration::NRIterator NR(amm, aigg, vars, rvals, T, rho, rz);
+
BENCHMARK("alphar") {
return model.alphar(T, rho, z);
};
@@ -33,8 +45,24 @@ TEST_CASE("multifluid derivatives", "[mf]")
return DerivativeHolderSquare<2,AlphaWrapperOption::idealgas>(aig, T, rho, z).derivs;
};
- BENCHMARK("IterationStepper matrix") {
- std::vector<char> vars = {'P','S','T','D'};
- return build_iteration_Jv(vars, model, aig, T, rho, z);
+ BENCHMARK("All residual derivatives (via AbstractModel) needed for first derivatives of h,s,u,p w.r.t. T&rho") {
+ return amm->get_deriv_mat2(T, rho, z);
+ };
+ BENCHMARK("All residual derivatives (via AbstractModel) needed for first derivatives of h,s,u,p w.r.t. T&rho") {
+ return aigg->get_deriv_mat2(T, rho, z);
};
+
+ BENCHMARK("Newton-Raphson construction") {
+ return teqp::iteration::NRIterator(amm, aigg, vars, rvals, T, rho, rz);
+ };
+ BENCHMARK("Newton-Raphson calc_step") {
+ return NR.calc_step(T, rho);
+ };
+ BENCHMARK("Newton-Raphson take_step") {
+ return NR.take_step();
+ };
+ BENCHMARK("Newton-Raphson take_steps") {
+ return NR.take_steps(5);
+ };
+
}
--
GitLab