Add Newton-Raphson iterator for homogeneous iterative calculations

include/teqp/algorithms/iteration.hpp

+#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();
+        }
+    }
+};
+
+
+}
+}

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>;

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);
+    };
+    
 }