From 8a3ef27de35670bf30b22c0ad26e097cfe3aaf61 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Tue, 2 Aug 2022 17:51:15 -0400
Subject: [PATCH] Mixtures work now too for alphaig
Just need to fill in the remaining terms and test, test, test
---
include/teqp/ideal_eosterms.hpp | 49 ++++++++++++++++++++++++++++----
src/tests/catch_test_alphaig.cxx | 8 ++++--
2 files changed, 50 insertions(+), 7 deletions(-)
diff --git a/include/teqp/ideal_eosterms.hpp b/include/teqp/ideal_eosterms.hpp
index 6d22c14..537ee66 100644
--- a/include/teqp/ideal_eosterms.hpp
+++ b/include/teqp/ideal_eosterms.hpp
@@ -30,11 +30,11 @@ namespace teqp {
using IdealHelmholtzTerms = std::variant<IdealHelmholtzLead>;
- class IdealHelmholtz {
+ class PureIdealHelmholtz {
public:
std::vector<IdealHelmholtzTerms> contributions;
- IdealHelmholtz(const nlohmann::json& j) {
- for (auto& term : j) {
+ PureIdealHelmholtz(const nlohmann::json& jpure) {
+ for (auto& term : jpure) {
if (term.at("type") == "Lead") {
contributions.emplace_back(IdealHelmholtzLead(term.at("a_1"), term.at("a_2")));
}
@@ -43,8 +43,8 @@ namespace teqp {
}
}
}
- template<typename TType, typename RhoType, typename MoleFrac>
- auto alphaig(const TType& T, const RhoType &rho, const MoleFrac& /**/) const{
+ template<typename TType, typename RhoType>
+ auto alphaig(const TType& T, const RhoType &rho) const{
std::common_type_t <TType, RhoType> ig = 0.0;
for (const auto& term : contributions) {
auto contrib = std::visit([&](auto& t) { return t.alphaig(T, rho); }, term);
@@ -54,4 +54,43 @@ namespace teqp {
}
};
+ /**
+ * @brief Ideal-gas Helmholtz energy container
+ *
+ * \f[ \alpha^{\rm ig} = \sum_i x_i[\alpha^{\rm ig}_{oi}(T,\rho) + x_i] \f]
+ *
+ * where x_i are mole fractions
+ *
+ */
+ class IdealHelmholtz {
+
+ public:
+
+ std::vector<PureIdealHelmholtz> pures;
+
+ IdealHelmholtz(const nlohmann::json &jpures){
+ for (auto &jpure : jpures){
+ pures.emplace_back(jpure);
+ }
+ }
+
+ template<typename TType, typename RhoType, typename MoleFrac>
+ auto alphaig(const TType& T, const RhoType &rho, const MoleFrac &molefrac) const {
+ std::common_type_t <TType, RhoType, decltype(molefrac[0])> ig = 0.0;
+ if (molefrac.size() != pures.size()){
+ throw teqp::InvalidArgument("molefrac and pures are not the same length");
+ }
+ std::size_t i = 0;
+ for (auto &pure : pures){
+ if (molefrac[i] != 0){
+ ig += molefrac[i]*(pure.alphaig(T, rho) + log(molefrac[i]));
+ }
+ else{
+ // lim_{x\to 0} x*ln(x) => 0
+ }
+ i++;
+ }
+ return ig;
+ }
+ };
}
\ No newline at end of file
diff --git a/src/tests/catch_test_alphaig.cxx b/src/tests/catch_test_alphaig.cxx
index ab1a40e..7a00f5e 100644
--- a/src/tests/catch_test_alphaig.cxx
+++ b/src/tests/catch_test_alphaig.cxx
@@ -10,8 +10,10 @@ using namespace teqp;
TEST_CASE("Simplest case","[alphaig]") {
double a_1 = 1, a_2 = 2, T = 300;
+ nlohmann::json j0 = nlohmann::json::array();
+ j0.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } });
nlohmann::json j = nlohmann::json::array();
- j.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } });
+ j.push_back(j0);
IdealHelmholtz ih(j);
std::valarray<double> molefrac{1.0};
REQUIRE(ih.alphaig(T, 1, molefrac) == log(1) + a_1 + a_2 / T);
@@ -19,8 +21,10 @@ TEST_CASE("Simplest case","[alphaig]") {
TEST_CASE("alphaig derivative", "[alphaig]") {
double a_1 = 1, a_2 = 2, T = 300, rho = 1;
+ nlohmann::json j0 = nlohmann::json::array();
+ j0.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } }); // For the first component
nlohmann::json j = nlohmann::json::array();
- j.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } });
+ j.push_back(j0);
IdealHelmholtz ih(j);
auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
auto wih = AlphaCallWrapper<1, decltype(ih)>(ih);
--
GitLab