Skip to content

Commit

Permalink
Add the T, rho, x combined derivative and expose
Browse files Browse the repository at this point in the history
  • Loading branch information
@ianhbell
ianhbell committed Mar 2, 2024
1 parent af95aad commit 8169bf1
Show file tree
Hide file tree
Showing 3 changed files with 111 additions and 1 deletion.
3 changes: 3 additions & 0 deletions include/teqp/cpp/deriv_adapter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,9 @@ class DerivativeAdapter : public teqp::cppinterface::AbstractModel{
virtual double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const override {
return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_dmBnvirdTm_runtime(Nderiv, NTderiv, mp.get_cref(), T, molefrac);
};
virtual double get_ATrhoXi(const double T, int NT, const double rhomolar, int ND, const EArrayd& molefrac, int i, int NX) const override {
return TDXDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_ATrhoXi_runtime(mp.get_cref(), T, NT, rhomolar, ND, molefrac, i, NX);
};

// Derivatives from isochoric thermodynamics (all have the same signature within each block), and they differ by their output argument
#define X(f) virtual double f(const double T, const EArrayd& rhovec) const override { return IsochoricDerivatives<decltype(mp.get_cref()), double, EArrayd>::f(mp.get_cref(), T, rhovec); };
Expand Down
3 changes: 3 additions & 0 deletions include/teqp/cpp/teqpcpp.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,9 @@ namespace teqp {
virtual double get_B12vir(const double T, const EArrayd& z) const = 0;
virtual double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& z) const = 0;

// Composition derivatives
virtual double get_ATrhoXi(const double T, int NT, const double rhomolar, int ND, const EArrayd& molefrac, int i, int NX) const = 0;

// Derivatives from isochoric thermodynamics (all have the same signature whithin each block)
#define X(f) virtual double f(const double T, const EArrayd& rhovec) const = 0;
ISOCHORIC_double_args
Expand Down
106 changes: 105 additions & 1 deletion include/teqp/derivs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,111 @@ struct TDXDerivatives {
throw std::invalid_argument("algorithmic differentiation backend is invalid in get_Agenxy for iD > 0 and iT > 0");
}
}
// return static_cast<Scalar>(-999999999*T); // This will never hit, only to make compiler happy because it doesn't know the return type
// return static_cast<Scalar>(-999999999*T); // This will never hit, only to make compiler happy because it doesn't know the return type
}

/**
Calculate the derivative
\f[
\Lambda_{xyz_i} = (1/T)^x(\rho)^y\deriv{^{x+y+z_i}(\alpha^r)}{(1/T)^x\partial \rho^y \partial \mathbf{Z}_i^{z_i}}}{}
\f]
in which all the compositions are treated as being independent
*/
template<int iT, int iD, int iXi, typename AlphaWrapper>
static auto get_ATrhoXi(const AlphaWrapper& w, const Scalar& T, const Scalar& rho, const VectorType& molefrac, int i){
using adtype = autodiff::HigherOrderDual<iT + iD + iXi, double>;
adtype Trecipad = 1.0 / T, rhoad = rho, xi = molefrac[i];
auto f = [&w, &molefrac, &i](const adtype& Trecip, const adtype& rho_, const adtype& xi_) {
adtype T_ = 1.0/Trecip;
Eigen::ArrayX<adtype> molefracdual = molefrac.template cast<adtype>();
molefracdual[i] = xi_;
return eval(w.alphar(T_, rho_, molefracdual)); }; // TODO: generalize to alphar, alpha, or alphaig
auto wrts = std::tuple_cat(build_duplicated_tuple<iT>(std::ref(Trecipad)), build_duplicated_tuple<iD>(std::ref(rhoad)), build_duplicated_tuple<iXi>(std::ref(xi)));
auto der = derivatives(f, std::apply(wrt_helper(), wrts), at(Trecipad, rhoad, xi));
return powi(forceeval(1.0 / T), iT) * powi(rho, iD) * der[der.size() - 1];
}

template<typename AlphaWrapper>
static auto get_ATrhoXi_runtime(const AlphaWrapper& w, const Scalar& T, int iT, const Scalar& rho, int iD, const VectorType& molefrac, int i, int iXi){
if (iT == 0 and iD == 0){
if (iXi == 1){
return get_ATrhoXi<0, 0, 1>(w, T, rho, molefrac, i);
}
else if (iXi == 2){
return get_ATrhoXi<0, 0, 2>(w, T, rho, molefrac, i);
}
else if (iXi == 3){
return get_ATrhoXi<0, 0, 3>(w, T, rho, molefrac, i);
}
}
else if (iT == 1 and iD == 0){
if (iXi == 1){
return get_ATrhoXi<1, 0, 1>(w, T, rho, molefrac, i);
}
else if (iXi == 2){
return get_ATrhoXi<1, 0, 2>(w, T, rho, molefrac, i);
}
else if (iXi == 3){
return get_ATrhoXi<1, 0, 3>(w, T, rho, molefrac, i);
}
}
else if (iT == 0 and iD == 1){
if (iXi == 1){
return get_ATrhoXi<0, 1, 1>(w, T, rho, molefrac, i);
}
else if (iXi == 2){
return get_ATrhoXi<0, 1, 2>(w, T, rho, molefrac, i);
}
else if (iXi == 3){
return get_ATrhoXi<0, 1, 3>(w, T, rho, molefrac, i);
}
}
throw teqp::InvalidArgument("Can't match these derivative counts");
}

/**
Calculate the derivative
\f[
\Lambda_{xyz_i z_j } = (1/T)^x(\rho)^y\left(\frac{\partial^{x+y+z_i+z_j}(\alpha^r)}{\partial (1/T)^x\partial \rho^y \partial \mathbf{Z}_i^{z_i} \partial \mathbf{Z}_j^{z_j} \}\right)
\f]
in which all the compositions are treated as being independent
*/
template<int iT, int iD, int iXi, int iXj, typename AlphaWrapper>
static auto get_ATrhoXiXj(const AlphaWrapper& w, const Scalar& T, const Scalar& rho, const VectorType& molefrac, int i, int j){
using adtype = autodiff::HigherOrderDual<iT + iD + iXi + iXj, double>;
adtype Trecipad = 1.0 / T, rhoad = rho, xi = molefrac[i], xj = molefrac[j];
auto f = [&w, &molefrac, i, j](const adtype& Trecip, const adtype& rho_, const adtype& xi_, const adtype& xj_) {
adtype T_ = 1.0/Trecip;
Eigen::ArrayX<adtype> molefracdual = molefrac.template cast<adtype>();
molefracdual[i] = xi_;
molefracdual[j] = xj_;
return eval(w.alphar(T_, rho_, molefracdual)); };
auto wrts = std::tuple_cat(build_duplicated_tuple<iT>(std::ref(Trecipad)), build_duplicated_tuple<iD>(std::ref(rhoad)), build_duplicated_tuple<iXi>(std::ref(xi)), build_duplicated_tuple<iXj>(std::ref(xj)));
auto der = derivatives(f, std::apply(wrt_helper(), wrts), at(Trecipad, rhoad, xi, xj));
return powi(forceeval(1.0 / T), iT) * powi(rho, iD) * der[der.size() - 1];
}

/**
Calculate the derivative
\f[
\Lambda_{xyz_i z_j z_k} = (1/T)^x(\rho)^y\left(\frac{\partial^{x+y+z_i+z_j+z_k}(\alpha^r)}{\partial (1/T)^x\partial \rho^y \partial \mathbf{Z}_i^{z_i} \partial \mathbf{Z}_j^{z_j} \partial \mathbf{Z}_k^{z_k} \}\right
\f]
in which all the compositions are treated as being independent
*/
template<int iT, int iD, int iXi, int iXj, int iXk, typename AlphaWrapper>
static auto get_ATrhoXiXjXk(const AlphaWrapper& w, const Scalar& T, const Scalar& rho, const VectorType& molefrac, int i, int j, int k){
using adtype = autodiff::HigherOrderDual<iT + iD + iXi + iXj + iXk, double>;
adtype Trecipad = 1.0 / T, rhoad = rho, xi = molefrac[i], xj = molefrac[j], xk = molefrac[k];
auto f = [&w, &molefrac, i, j, k](const adtype& Trecip, const adtype& rho_, const adtype& xi_, const adtype& xj_, const adtype& xk_) {
adtype T_ = 1.0/Trecip;
Eigen::ArrayX<adtype> molefracdual = molefrac.template cast<adtype>();
molefracdual[i] = xi_;
molefracdual[j] = xj_;
molefracdual[k] = xk_;
return eval(w.alphar(T_, rho_, molefracdual)); };
auto wrts = std::tuple_cat(build_duplicated_tuple<iT>(std::ref(Trecipad)), build_duplicated_tuple<iD>(std::ref(rhoad)), build_duplicated_tuple<iXi>(std::ref(xi)), build_duplicated_tuple<iXj>(std::ref(xj)), build_duplicated_tuple<iXk>(std::ref(xk)));
auto der = derivatives(f, std::apply(wrt_helper(), wrts), at(Trecipad, rhoad, xi, xj, xk));
return powi(forceeval(1.0 / T), iT) * powi(rho, iD) * der[der.size() - 1];
}

/**
Expand Down

0 comments on commit 8169bf1

Please sign in to comment.