From bbac7dc5e1c2d10cc5f79054c8154e70f17d03f2 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Tue, 1 Nov 2022 16:14:07 -0400
Subject: [PATCH] Transition more things into the AbstractModel framework
---
include/teqp/algorithms/critical_tracing.hpp | 25 ++-----------
.../algorithms/critical_tracing_types.hpp | 27 ++++++++++++++
include/teqp/models/fwd.hpp | 6 ++--
interface/CPP/teqpcpp.cpp | 36 ++++++++++++++++---
interface/CPP/teqpcpp.hpp | 22 ++++++++++--
interface/CPP/test_teqpcpp.cpp | 12 ++++++-
6 files changed, 95 insertions(+), 33 deletions(-)
create mode 100644 include/teqp/algorithms/critical_tracing_types.hpp
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 3cdf044..72c4e33 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -8,6 +8,7 @@
#include <Eigen/Dense>
#include "teqp/algorithms/rootfinding.hpp"
#include "teqp/algorithms/critical_pure.hpp"
+#include "teqp/algorithms/critical_tracing_types.hpp"
#include "teqp/exceptions.hpp"
// Imports from boost
@@ -16,28 +17,6 @@
#include <boost/numeric/odeint/stepper/euler.hpp>
namespace teqp {
-// This has to be outside the critical tracing struct so that the pybind11 wrapper doesn't fight with the types
-struct TCABOptions {
- double abs_err = 1.0e-6,
- rel_err = 1.0e-6,
- init_dt = 10, ///< The initial step size
- max_dt = 10000000000,
- T_tol = 1e-6, ///< The tolerance on temperature to indicate that it is converged
- init_c = 1.0; ///< The c parameter which controls the initial search direction for the first step. Choices are 1 or -1
- int small_T_count = 5; ///< How many small temperature steps indicates convergence
- int integration_order = 5; ///< The order of integration, either 1 for simple Euler or 5 for adaptive RK45
- int max_step_count = 1000; ///< Maximum number of steps allowed
- int skip_dircheck_count = 1; ///< Only start checking the direction dot product after this many steps
- bool polish = false; ///< If true, polish the solution at every step
- double polish_reltol_T = 0.01; ///< The maximum allowed change in temperature when polishing
- double polish_reltol_rho = 0.05; ///< The maximum allowed change in any molar concentration when polishing
- bool terminate_negative_density = true; ///< Stop the tracing if the density is negative
- bool calc_stability = false; ///< Calculate the local stability with the method of Deiters and Bell
- double stability_rel_drho = 0.001; ///< The relative size of the step (relative to the sum of the molar concentration vector) to be used when taking the step in the direction of \f$\sigma_1\f$ when assessing local stability
- int verbosity = 0; ///< The greater the verbosity, the more output you will get, especially about polishing failures
- bool polish_exception_on_fail = false; ///< If true, when polishing fails, throw an exception, otherwise, terminate tracing
- bool pure_endpoint_polish = false; ///< If true, if the last step crossed into negative concentrations, try to interpolate to find the pure fluid endpoint hiding in the data
-};
template<typename Model, typename Scalar = double, typename VecType = Eigen::ArrayXd>
struct CriticalTracing {
@@ -776,4 +755,4 @@ struct CriticalTracing {
}; // namespace VecType
-}; // namespace teqp
\ No newline at end of file
+}; // namespace teqp
diff --git a/include/teqp/algorithms/critical_tracing_types.hpp b/include/teqp/algorithms/critical_tracing_types.hpp
new file mode 100644
index 0000000..48703d7
--- /dev/null
+++ b/include/teqp/algorithms/critical_tracing_types.hpp
@@ -0,0 +1,27 @@
+# pragma once
+
+namespace teqp {
+
+struct TCABOptions {
+ double abs_err = 1.0e-6,
+ rel_err = 1.0e-6,
+ init_dt = 10, ///< The initial step size
+ max_dt = 10000000000,
+ T_tol = 1e-6, ///< The tolerance on temperature to indicate that it is converged
+ init_c = 1.0; ///< The c parameter which controls the initial search direction for the first step. Choices are 1 or -1
+ int small_T_count = 5; ///< How many small temperature steps indicates convergence
+ int integration_order = 5; ///< The order of integration, either 1 for simple Euler or 5 for adaptive RK45
+ int max_step_count = 1000; ///< Maximum number of steps allowed
+ int skip_dircheck_count = 1; ///< Only start checking the direction dot product after this many steps
+ bool polish = false; ///< If true, polish the solution at every step
+ double polish_reltol_T = 0.01; ///< The maximum allowed change in temperature when polishing
+ double polish_reltol_rho = 0.05; ///< The maximum allowed change in any molar concentration when polishing
+ bool terminate_negative_density = true; ///< Stop the tracing if the density is negative
+ bool calc_stability = false; ///< Calculate the local stability with the method of Deiters and Bell
+ double stability_rel_drho = 0.001; ///< The relative size of the step (relative to the sum of the molar concentration vector) to be used when taking the step in the direction of \f$\sigma_1\f$ when assessing local stability
+ int verbosity = 0; ///< The greater the verbosity, the more output you will get, especially about polishing failures
+ bool polish_exception_on_fail = false; ///< If true, when polishing fails, throw an exception, otherwise, terminate tracing
+ bool pure_endpoint_polish = false; ///< If true, if the last step crossed into negative concentrations, try to interpolate to find the pure fluid endpoint hiding in the data
+};
+
+}
diff --git a/include/teqp/models/fwd.hpp b/include/teqp/models/fwd.hpp
index 9a65966..1894738 100644
--- a/include/teqp/models/fwd.hpp
+++ b/include/teqp/models/fwd.hpp
@@ -21,12 +21,14 @@ namespace teqp {
using vad = std::valarray<double>;
+ using PCSAFTType = decltype(PCSAFT::PCSAFTfactory(nlohmann::json{}));
+
// Define the EOS types by interrogating the types returned by the respective factory function
using AllowedModels = std::variant<
vdWEOS1,
decltype(canonical_PR(vad{}, vad{}, vad{})),
decltype(CPA::CPAfactory(nlohmann::json{})),
- decltype(PCSAFT::PCSAFTfactory(nlohmann::json{})),
+ PCSAFTType,
decltype(multifluidfactory(nlohmann::json{}))
>;
-}
\ No newline at end of file
+}
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index 2e309a9..ecac905 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -2,6 +2,7 @@
#include "teqp/derivs.hpp"
#include "teqp/json_builder.hpp"
#include "teqp/algorithms/critical_tracing.hpp"
+#include "teqp/algorithms/VLE.hpp"
using namespace teqp;
using namespace teqp::cppinterface;
@@ -12,21 +13,48 @@ namespace teqp {
class ModelImplementer : public AbstractModel {
protected:
const AllowedModels m_model;
+
public:
ModelImplementer(AllowedModels&& model) : m_model(model) {};
- double get_Arxy(const int NT, const int ND, const double T, const double rho, const Eigen::ArrayXd& molefracs) const override {
+ double get_Arxy(const int NT, const int ND, const double T, const double rho, const EArrayd& molefracs) const override {
return std::visit([&](const auto& model) {
- using tdx = teqp::TDXDerivatives<decltype(model), double, Eigen::ArrayXd>;
+ using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
return tdx::template get_Ar(NT, ND, model, T, rho, molefracs);
}, m_model);
}
- nlohmann::json trace_critical_arclength_binary(const double T0, const Eigen::ArrayXd& rhovec0) const override {
+ nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename_, const std::optional<TCABOptions> &options_) const override {
return std::visit([&](const auto& model) {
using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec0)>>;
return crit::trace_critical_arclength_binary(model, T0, rhovec0, "");
}, m_model);
}
+ EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const override {
+ return std::visit([&](const auto& model) {
+ return teqp::pure_VLE_T(model, T, rhoL, rhoV, maxiter);
+ }, m_model);
+ }
+ EArrayd get_fugacity_coefficients(const double T, const EArrayd& rhovec) const override {
+ return std::visit([&](const auto& model) {
+ using id = IsochoricDerivatives<decltype(model), double, EArrayd>;
+ return id::get_fugacity_coefficients(model, T, rhovec);
+ }, m_model);
+ }
+
+
+ // Methods only available for PC-SAFT
+ EArrayd get_m() const override {
+ return std::get<PCSAFTType>(m_model).get_m();
+ }
+ EArrayd get_sigma_Angstrom() const override {
+ return std::get<PCSAFTType>(m_model).get_sigma_Angstrom();
+ }
+ EArrayd get_epsilon_over_k_K() const override {
+ return std::get<PCSAFTType>(m_model).get_m();
+ }
+ double max_rhoN(const double T, const EArrayd& z) const override {
+ return std::get<PCSAFTType>(m_model).max_rhoN(T, z);
+ }
};
std::unique_ptr<AbstractModel> make_model(const nlohmann::json& j) {
@@ -37,4 +65,4 @@ namespace teqp {
return std::make_unique<ModelImplementer>(build_multifluid_model(components, coolprop_root, BIPcollectionpath, flags, departurepath));
}
}
-}
\ No newline at end of file
+}
diff --git a/interface/CPP/teqpcpp.hpp b/interface/CPP/teqpcpp.hpp
index cd8b3e7..b86440f 100644
--- a/interface/CPP/teqpcpp.hpp
+++ b/interface/CPP/teqpcpp.hpp
@@ -4,13 +4,29 @@
#include <Eigen/Dense>
#include "nlohmann/json.hpp"
+// The only headers that can be included here are
+// ones that define and use POD (plain ole' data) types
+#include "teqp/algorithms/critical_tracing_types.hpp"
+
+using EArray2 = Eigen::Array<double, 2, 1>;
+using EArrayd = Eigen::ArrayX<double>;
+
namespace teqp {
namespace cppinterface {
class AbstractModel {
public:
- virtual double get_Arxy(const int, const int, const double, const double, const Eigen::ArrayXd&) const = 0;
- virtual nlohmann::json trace_critical_arclength_binary(const double T0, const Eigen::ArrayXd& rhovec0) const = 0;
+ virtual double get_Arxy(const int, const int, const double, const double, const EArrayd&) const = 0;
+ virtual nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>&, const std::optional<TCABOptions> &) const = 0;
+ virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
+ virtual EArrayd get_fugacity_coefficients(const double T, const EArrayd& rhovec) const = 0;
+
+ // Methods only available for PCSAFT
+ virtual EArrayd get_m() const = 0;
+ virtual EArrayd get_sigma_Angstrom() const = 0;
+ virtual EArrayd get_epsilon_over_k_K() const = 0;
+ virtual double max_rhoN(const double, const EArrayd&) const = 0;
+
virtual ~AbstractModel() = default;
};
@@ -27,4 +43,4 @@ namespace teqp {
const std::string& departurepath = {}
);
}
-}
\ No newline at end of file
+}
diff --git a/interface/CPP/test_teqpcpp.cpp b/interface/CPP/test_teqpcpp.cpp
index 801c4e1..7ddb9ba 100644
--- a/interface/CPP/test_teqpcpp.cpp
+++ b/interface/CPP/test_teqpcpp.cpp
@@ -26,4 +26,14 @@ int main() {
auto z = (Eigen::ArrayXd(2) << 0.5, 0.5).finished();
double Ar01 = am->get_Arxy(0, 1, 300, 3, z);
std::cout << Ar01 << std::endl;
-}
\ No newline at end of file
+
+ auto fugcoeff = am->get_fugacity_coefficients(300.0, z*300.0);
+ std::cout << fugcoeff << std::endl;
+
+ try{
+ std::cout << am->get_m() << std::endl;
+ }
+ catch(...){
+ std::cout << "This fails, as it should because the model is not a PCSAFT one" << std::endl;
+ }
+}
--
GitLab