derivs.hpp

Go to the documentation of this file.

#pragma once
 
namespace teqp{
namespace cppinterface{
 
struct IterationMatrices{
    std::vector<char> vars; 
    Eigen::ArrayXXd J; 
    Eigen::ArrayXd v; 
};
 
template<typename Array>
auto build_iteration_Jv(const std::vector<char>& vars, const Eigen::Array<double, 3, 3>& Ar, const Eigen::Array<double, 3, 3>& Aig, const double R, const double T, const double rho, const Array &z){
    IterationMatrices im; im.J.resize(vars.size(), 2); im.v.resize(vars.size()); im.vars = vars;
    
    auto A = Ar + Aig;
    Eigen::ArrayXd& v = im.v;
    Eigen::ArrayXXd& J = im.J;
    auto Trecip = 1.0/T;
    auto dTrecipdT = -Trecip*Trecip;
    
    // Define some lambda functions for things that *might* be needed
    // The lambdas are used to get a sort of lazy evaluation. The lambdas are
    // only evaluated on an as-needed basis.  If the lambda is not called,
    // no cost is incurred.
    //
    // Probably the compiler will inline these functions anyhow.
    //
    // Derivatives of total alpha; alpha = a/(R*T) = a*R/Trecip
    auto alpha = [&](){ return A(0,0); };
    auto dalphadTrecip = [&](){ return A(1,0)/Trecip; };
    auto dalphadrho = [&](){ return A(0,1)/rho; };
    auto d2alphadTrecip2 = [&](){ return A(2,0)/(Trecip*Trecip); };
    auto d2alphadTrecipdrho = [&](){ return A(1,1)/(Trecip*rho); };
    auto d2alphadrho2 = [&](){ return A(0,2)/(rho*rho); };
    //
    // Derivatives of total Helmholtz energy a in terms of derivatives of alpha
    auto dadTrecip = [&](){ return R/(Trecip*Trecip)*(Trecip*dalphadTrecip()-alpha());};
    auto d2adTrecip2 = [&](){ return R/(Trecip*Trecip*Trecip)*(Trecip*Trecip*d2alphadTrecip2()-2*Trecip*dalphadTrecip()+2*alpha());};
    auto dadrho = [&](){return R/Trecip*(dalphadrho());};
    auto d2adrho2 = [&](){return R/Trecip*(d2alphadrho2());};
    auto d2adTrecipdrho = [&](){ return R/(Trecip*Trecip)*(Trecip*d2alphadTrecipdrho()-dalphadrho());};
    
    for (auto i = 0; i < vars.size(); ++i){
        switch(vars[i]){
            case 'T':
                v(i) = T;
                J(i, 0) = 1.0;
                J(i, 1) = 0.0;
                break;
            case 'D':
                v(i) = rho;
                J(i, 0) = 0.0;
                J(i, 1) = 1.0;
                break;
            case 'P':
                v(i) = rho*R*T*(1 + Ar(0,1));
                J(i, 0) = rho*R*(1 + Ar(0,1) - Ar(1,1));
                J(i, 1) = R*T*(1 + 2*Ar(0,1) + Ar(0,2));
                break;
            case 'S':
                v(i) = Trecip*Trecip*dadTrecip();
                J(i, 0) = (Trecip*Trecip*d2adTrecip2() + 2*Trecip*dadTrecip())*dTrecipdT;
                J(i, 1) = Trecip*Trecip*d2adTrecipdrho();
                break;
            default:
                throw std::invalid_argument("bad var: " + std::to_string(vars[i]));
        }
    }
    return im;
}
 
}
};