Skip to content

Commit

Permalink
Experiment with concepts to call alpha or alphar
Browse files Browse the repository at this point in the history
  • Loading branch information
@ianhbell
ianhbell committed Mar 2, 2024
1 parent 28db2bd commit 89b5967
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 69 deletions.
2 changes: 1 addition & 1 deletion include/teqp/cpp/deriv_adapter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ class DerivativeAdapter : public teqp::cppinterface::AbstractModel{
};

virtual EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z ) const override {
return DerivativeHolderSquare<2, AlphaWrapperOption::residual>(mp.get_cref(), T, rho, z).derivs;
return DerivativeHolderSquare<2>(mp.get_cref(), T, rho, z).derivs;
};
};

Expand Down
120 changes: 52 additions & 68 deletions include/teqp/derivs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <map>
#include <tuple>
#include <numeric>
#include <concepts>

#include "teqp/types.hpp"
#include "teqp/exceptions.hpp"
Expand Down Expand Up @@ -88,35 +89,6 @@ struct wrt_helper {
}
};

enum class AlphaWrapperOption {residual, idealgas};
/**
* \brief This class is used to wrap a model that exposes the generic
* functions alphar, alphaig, etc., and allow the desired member function to be
* called at runtime via perfect forwarding
*
* This class is needed because conventional binding methods (e.g., std::bind)
* require the argument types to be known, and they are not known in this case
* so we give the hard work of managing the argument types to the compiler
*/
template<AlphaWrapperOption o, class Model>
struct AlphaCallWrapper {
const Model& m_model;
AlphaCallWrapper(const Model& model) : m_model(model) {};

template <typename ... Args>
auto alpha(const Args& ... args) const {
if constexpr (o == AlphaWrapperOption::residual) {
// The alphar method is REQUIRED to be implemented by all
// models, so can just call it via perfect fowarding
return m_model.alphar(std::forward<const Args>(args)...);
}
else {
//throw teqp::InvalidArgument("Missing implementation for alphaig");
return m_model.alphaig(std::forward<const Args>(args)...);
}
}
};

enum class ADBackends { autodiff
#if defined(TEQP_MULTICOMPLEX_ENABLED)
,multicomplex
Expand All @@ -126,13 +98,37 @@ enum class ADBackends { autodiff
#endif
};

template<typename T, typename U, typename V, typename W>
concept CallableAlpha = requires(T t, U u, V v, W w) {
{ t.alpha(u,v,w) };
};

template<typename T, typename U, typename V, typename W>
concept CallableAlphar = requires(T t, U u, V v, W w) {
{ t.alphar(u,v,w) };
};

template<typename T, typename U, typename V, typename W>
concept CallableAlpharTauDelta = requires(T t, U u, V v, W w) {
{ t.alphar_taudelta(u,v,w) };
};

template<typename Model, typename Scalar = double, typename VectorType = Eigen::ArrayXd>
struct TDXDerivatives {

static auto get_Ar00(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
return model.alphar(T, rho, molefrac);
}

template<typename AlphaWrapper, typename S1, typename S2, typename Vec>
static auto AlphaCaller(const AlphaWrapper& w, const S1& T, const S2& rho, const Vec& molefrac) requires CallableAlpha<AlphaWrapper, S1, S2, Vec>{
return w.alpha(T, rho, molefrac);
}
template<typename AlphaWrapper, typename S1, typename S2, typename Vec>
static auto AlphaCaller(const AlphaWrapper& w, const S1& T, const S2& rho, const Vec& molefrac) requires CallableAlphar<AlphaWrapper, S1, S2, Vec>{
return w.alphar(T, rho, molefrac);
}

/**
* Calculate the derivative \f$\Lambda_{xy}\f$, where
* \f[
Expand All @@ -144,25 +140,27 @@ struct TDXDerivatives {
template<int iT, int iD, ADBackends be = ADBackends::autodiff, class AlphaWrapper>
static auto get_Agenxy(const AlphaWrapper& w, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {

static_assert(iT > 0 || iD > 0);
if constexpr (iT == 0 && iD > 0) {
if constexpr (iT == 0 && iD == 0){
return AlphaCaller(w, T, rho, molefrac);
}
else if constexpr (iT == 0 && iD > 0) {
if constexpr (be == ADBackends::autodiff) {
// If a pure derivative, then we can use autodiff::Real for that variable and Scalar for other variable
autodiff::Real<iD, Scalar> rho_ = rho;
auto f = [&w, &T, &molefrac](const auto& rho__) { return w.alpha(T, rho__, molefrac); };
auto f = [&w, &T, &molefrac](const auto& rho__) { return AlphaCaller(w, T, rho__, molefrac); };
return powi(rho, iD)*derivatives(f, along(1), at(rho_))[iD];
}
#if defined(TEQP_COMPLEXSTEP_ENABLED)
else if constexpr (iD == 1 && be == ADBackends::complex_step) {
double h = 1e-100;
auto rho_ = std::complex<Scalar>(rho, h);
return powi(rho, iD) * w.alpha(T, rho_, molefrac).imag() / h;
return powi(rho, iD) * AlphaCaller(w, T, rho_, molefrac).imag() / h;
}
#endif
#if defined(TEQP_MULTICOMPLEX_ENABLED)
else if constexpr (be == ADBackends::multicomplex) {
using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
fcn_t f = [&](const auto& rhomcx) { return w.alpha(T, rhomcx, molefrac); };
fcn_t f = [&](const auto& rhomcx) { return AlphaCaller(w, T, rhomcx, molefrac); };
auto ders = diff_mcx1(f, rho, iD, true /* and_val */);
return powi(rho, iD)*ders[iD];
}
Expand All @@ -176,7 +174,7 @@ struct TDXDerivatives {
if constexpr (be == ADBackends::autodiff) {
// If a pure derivative, then we can use autodiff::Real for that variable and Scalar for other variable
autodiff::Real<iT, Scalar> Trecipad = Trecip;
auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(forceeval(1.0/Trecip__), rho, molefrac); };
auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return AlphaCaller(w, forceeval(1.0/Trecip__), rho, molefrac); };
return powi(Trecip, iT)*derivatives(f, along(1), at(Trecipad))[iT];
}
#if defined(TEQP_COMPLEXSTEP_ENABLED)
Expand All @@ -189,7 +187,7 @@ struct TDXDerivatives {
#if defined(TEQP_MULTICOMPLEX_ENABLED)
else if constexpr (be == ADBackends::multicomplex) {
using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
fcn_t f = [&](const auto& Trecipmcx) { return w.alpha(1.0/Trecipmcx, rho, molefrac); };
fcn_t f = [&](const auto& Trecipmcx) { return AlphaCaller(w, 1.0/Trecipmcx, rho, molefrac); };
auto ders = diff_mcx1(f, Trecip, iT, true /* and_val */);
return powi(Trecip, iT)*ders[iT];
}
Expand All @@ -204,7 +202,7 @@ struct TDXDerivatives {
adtype Trecipad = 1.0 / T, rhoad = rho;
auto f = [&w, &molefrac](const adtype& Trecip, const adtype& rho_) {
adtype T_ = 1.0/Trecip;
return eval(w.alpha(T_, rho_, molefrac)); };
return eval(AlphaCaller(w, T_, rho_, molefrac)); };
auto wrts = std::tuple_cat(build_duplicated_tuple<iT>(std::ref(Trecipad)), build_duplicated_tuple<iD>(std::ref(rhoad)));
auto der = derivatives(f, std::apply(wrt_helper(), wrts), at(Trecipad, rhoad));
return powi(forceeval(1.0 / T), iT) * powi(rho, iD) * der[der.size() - 1];
Expand All @@ -214,7 +212,7 @@ struct TDXDerivatives {
using fcn_t = std::function< mcx::MultiComplex<double>(const std::valarray<mcx::MultiComplex<double>>&)>;
const fcn_t func = [&w, &molefrac](const auto& zs) {
auto Trecip = zs[0], rhomolar = zs[1];
return w.alpha(1.0 / Trecip, rhomolar, molefrac);
return AlphaCaller(w, 1.0 / Trecip, rhomolar, molefrac);
};
std::vector<double> xs = { 1.0 / T, rho};
std::vector<int> order = { iT, iD };
Expand Down Expand Up @@ -244,7 +242,7 @@ struct TDXDerivatives {
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
return eval(AlphaCaller(w, 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)));
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];
Expand Down Expand Up @@ -304,7 +302,7 @@ struct TDXDerivatives {
Eigen::ArrayX<adtype> molefracdual = molefrac.template cast<adtype>();
molefracdual[i] = xi_;
molefracdual[j] = xj_;
return eval(w.alphar(T_, rho_, molefracdual)); };
return eval(AlphaCaller(w, 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];
Expand All @@ -327,7 +325,7 @@ struct TDXDerivatives {
molefracdual[i] = xi_;
molefracdual[j] = xj_;
molefracdual[k] = xk_;
return eval(w.alphar(T_, rho_, molefracdual)); };
return eval(AlphaCaller(w, 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 All @@ -343,13 +341,7 @@ struct TDXDerivatives {
*/
template<int iT, int iD, ADBackends be = ADBackends::autodiff>
static auto get_Arxy(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
if constexpr (iT == 0 && iD == 0) {
return wrapper.alpha(T, rho, molefrac);
}
else {
return get_Agenxy<iT, iD, be>(wrapper, T, rho, molefrac);
}
return get_Agenxy<iT, iD, be>(model, T, rho, molefrac);
}

/**
Expand All @@ -362,13 +354,7 @@ struct TDXDerivatives {
*/
template<int iT, int iD, ADBackends be = ADBackends::autodiff>
static auto get_Aigxy(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
auto wrapper = AlphaCallWrapper<AlphaWrapperOption::idealgas, decltype(model)>(model);
if constexpr (iT == 0 && iD == 0) {
return wrapper.alpha(T, rho, molefrac);
}
else {
return get_Agenxy<iT, iD, be>(wrapper, T, rho, molefrac);
}
return get_Agenxy<iT, iD, be>(model, T, rho, molefrac);
}

template<ADBackends be = ADBackends::autodiff>
Expand Down Expand Up @@ -427,7 +413,7 @@ struct TDXDerivatives {
if constexpr (be == ADBackends::autodiff) {
// If a pure derivative, then we can use autodiff::Real for that variable and Scalar for other variable
autodiff::Real<Nderiv, Scalar> rho_ = rho;
auto f = [&w, &T, &molefrac](const auto& rho__) { return w.alpha(T, rho__, molefrac); };
auto f = [&w, &T, &molefrac](const auto& rho__) { return AlphaCaller(w, T, rho__, molefrac); };
auto ders = derivatives(f, along(1), at(rho_));
for (auto n = 0; n <= Nderiv; ++n) {
o[n] = forceeval(powi(rho, n) * ders[n]);
Expand All @@ -438,7 +424,7 @@ struct TDXDerivatives {
else {
using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
bool and_val = true;
fcn_t f = [&w, &T, &molefrac](const auto& rhomcx) { return w.alpha(T, rhomcx, molefrac); };
fcn_t f = [&w, &T, &molefrac](const auto& rhomcx) { return AlphaCaller(w, T, rhomcx, molefrac); };
auto ders = diff_mcx1(f, rho, Nderiv, and_val);
for (auto n = 0; n <= Nderiv; ++n) {
o[n] = powi(rho, n) * ders[n];
Expand All @@ -456,7 +442,7 @@ struct TDXDerivatives {
if constexpr (be == ADBackends::autodiff) {
// If a pure derivative, then we can use autodiff::Real for that variable and Scalar for other variable
autodiff::Real<Nderiv, Scalar> Trecipad = Trecip;
auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(forceeval(1.0/Trecip__), rho, molefrac); };
auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return AlphaCaller(w, forceeval(1.0/Trecip__), rho, molefrac); };
auto ders = derivatives(f, along(1), at(Trecipad));
for (auto n = 0; n <= Nderiv; ++n) {
o[n] = powi(Trecip, n) * ders[n];
Expand All @@ -465,7 +451,7 @@ struct TDXDerivatives {
#if defined(TEQP_MULTICOMPLEX_ENABLED)
else if constexpr (be == ADBackends::multicomplex) {
using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
fcn_t f = [&](const auto& Trecipmcx) { return w.alpha(1.0/Trecipmcx, rho, molefrac); };
fcn_t f = [&](const auto& Trecipmcx) { return AlphaCaller(w, 1.0/Trecipmcx, rho, molefrac); };
auto ders = diff_mcx1(f, Trecip, Nderiv+1, true /* and_val */);
for (auto n = 0; n <= Nderiv; ++n) {
o[n] = powi(Trecip, n) * ders[n];
Expand All @@ -488,8 +474,7 @@ struct TDXDerivatives {
*/
template<int iT, ADBackends be = ADBackends::autodiff>
static auto get_Arn0(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
return get_Agenn0<iT, be>(wrapper, T, rho, molefrac);
return get_Agenn0<iT, be>(model, T, rho, molefrac);
}

/**
Expand All @@ -502,8 +487,7 @@ struct TDXDerivatives {
*/
template<int iD, ADBackends be = ADBackends::autodiff>
static auto get_Ar0n(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
return get_Agen0n<iD, be>(wrapper, T, rho, molefrac);
return get_Agen0n<iD, be>(model, T, rho, molefrac);
}


Expand Down Expand Up @@ -1356,7 +1340,7 @@ struct IsochoricDerivatives{
}
};

template<int Nderivsmax, AlphaWrapperOption opt>
template<int Nderivsmax>
class DerivativeHolderSquare{

public:
Expand All @@ -1366,19 +1350,19 @@ class DerivativeHolderSquare{
DerivativeHolderSquare(const Model& model, const Scalar& T, const Scalar& rho, const VecType& z) {
using tdx = TDXDerivatives<decltype(model), Scalar, VecType>;
static_assert(Nderivsmax == 2, "It's gotta be 2 for now");
AlphaCallWrapper<opt, Model> wrapper(model);

auto AX02 = tdx::template get_Agen0n<2>(wrapper, T, rho, z);

auto AX02 = tdx::template get_Agen0n<2>(model, T, rho, z);
derivs(0, 0) = AX02[0];
derivs(0, 1) = AX02[1];
derivs(0, 2) = AX02[2];

auto AX20 = tdx::template get_Agenn0<2>(wrapper, T, rho, z);
auto AX20 = tdx::template get_Agenn0<2>(model, T, rho, z);
derivs(0, 0) = AX20[0];
derivs(1, 0) = AX20[1];
derivs(2, 0) = AX20[2];

derivs(1, 1) = tdx::template get_Agenxy<1,1>(wrapper, T, rho, z);
derivs(1, 1) = tdx::template get_Agenxy<1,1>(model, T, rho, z);
}
};

Expand Down

0 comments on commit 89b5967

Please sign in to comment.