Skip to content
Snippets Groups Projects
Commit 60330194 authored by Ian Bell's avatar Ian Bell
Browse files

Slight simplification of derivatives

parent 41ccad13
No related branches found
No related tags found
No related merge requests found
...@@ -17,8 +17,7 @@ caller(const FuncType& f, TType T, const ContainerType& rho) { ...@@ -17,8 +17,7 @@ caller(const FuncType& f, TType T, const ContainerType& rho) {
* respect to the first variable which here is temperature * respect to the first variable which here is temperature
*/ */
template <typename TType, typename ContainerType, typename FuncType> template <typename TType, typename ContainerType, typename FuncType>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type typename ContainerType::value_type derivT(const FuncType& f, TType T, const ContainerType& rho) {
derivT(const FuncType& f, TType T, const ContainerType& rho) {
double h = 1e-100; double h = 1e-100;
return f(std::complex<TType>(T, h), rho).imag() / h; return f(std::complex<TType>(T, h), rho).imag() / h;
} }
...@@ -28,8 +27,7 @@ derivT(const FuncType& f, TType T, const ContainerType& rho) { ...@@ -28,8 +27,7 @@ derivT(const FuncType& f, TType T, const ContainerType& rho) {
* respect to the first variable which here is temperature * respect to the first variable which here is temperature
*/ */
template <typename TType, typename ContainerType, typename FuncType> template <typename TType, typename ContainerType, typename FuncType>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type typename ContainerType::value_type derivTmcx(const FuncType& f, TType T, const ContainerType& rho) {
derivTmcx(const FuncType& f, TType T, const ContainerType& rho) {
using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>; using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
fcn_t wrapper = [&rho, &f](const MultiComplex<TType>& T_) {return f(T_, rho); }; fcn_t wrapper = [&rho, &f](const MultiComplex<TType>& T_) {return f(T_, rho); };
auto ders = diff_mcx1(wrapper, T, 1); auto ders = diff_mcx1(wrapper, T, 1);
...@@ -41,8 +39,7 @@ derivTmcx(const FuncType& f, TType T, const ContainerType& rho) { ...@@ -41,8 +39,7 @@ derivTmcx(const FuncType& f, TType T, const ContainerType& rho) {
* to the given composition variable * to the given composition variable
*/ */
template <typename TType, typename ContainerType, typename FuncType, typename Integer> template <typename TType, typename ContainerType, typename FuncType, typename Integer>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type typename ContainerType::value_type derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) {
derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) {
double h = 1e-100; double h = 1e-100;
using comtype = std::complex<ContainerType::value_type>; using comtype = std::complex<ContainerType::value_type>;
std::valarray<comtype> rhocom(rho.size()); std::valarray<comtype> rhocom(rho.size());
...@@ -57,8 +54,7 @@ derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) { ...@@ -57,8 +54,7 @@ derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) {
* \brief Calculate the Psir=ar*rho * \brief Calculate the Psir=ar*rho
*/ */
template <typename TType, typename ContainerType, typename Model> template <typename TType, typename ContainerType, typename Model>
typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type typename ContainerType::value_type get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
using container = decltype(rhovec); using container = decltype(rhovec);
auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0); auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
return model.alphar(T, rhovec)*model.R*T*rhotot_; return model.alphar(T, rhovec)*model.R*T*rhotot_;
...@@ -68,12 +64,8 @@ get_Psir(const Model& model, const TType T, const ContainerType& rhovec) { ...@@ -68,12 +64,8 @@ get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
* \brief Calculate the residual pressure from derivatives of alphar * \brief Calculate the residual pressure from derivatives of alphar
*/ */
template <typename Model, typename TType, typename ContainerType> template <typename Model, typename TType, typename ContainerType>
typename ContainerType::value_type get_pr( typename ContainerType::value_type get_pr(const Model& model, const TType T, const ContainerType& rhovec)
const Model& model,
const TType T,
const ContainerType& rhovec)
{ {
using container = decltype(rhovec);
auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0); auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
decltype(rhovec[0] * T) pr = 0.0; decltype(rhovec[0] * T) pr = 0.0;
for (auto i = 0; i < rhovec.size(); ++i) { for (auto i = 0; i < rhovec.size(); ++i) {
...@@ -82,16 +74,17 @@ typename ContainerType::value_type get_pr( ...@@ -82,16 +74,17 @@ typename ContainerType::value_type get_pr(
return pr*rhotot_*model.R*T; return pr*rhotot_*model.R*T;
} }
template <typename Model, typename TType, typename ContainerType>
typename ContainerType::value_type get_Ar10(const Model& model, const TType T, const ContainerType& rhovec){
return T*derivT([&model](const auto& T, const auto& rhovec) { return model.alphar(T, rhovec); }, T, rhovec);
}
/*** /***
* \brief Calculate the residual entropy (s^+ = -sr/R) from derivatives of alphar * \brief Calculate the residual entropy (s^+ = -sr/R) from derivatives of alphar
*/ */
template <typename Model, typename TType, typename ContainerType> template <typename Model, typename TType, typename ContainerType>
typename ContainerType::value_type get_splus( typename ContainerType::value_type get_splus(const Model& model, const TType T, const ContainerType& rhovec){
const Model& model, return model.alphar(T, rhovec) - get_Ar10(model, T, rhovec);
const TType T,
const ContainerType& rhovec)
{
return model.alphar(T, rhovec) + T*derivT([&model](const auto& T, const auto& rhovec) { return model.alphar(T, rhovec); }, T, rhovec);
} }
/*** /***
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment