Newer
Older
template<ADBackends be = ADBackends::autodiff>
static auto get_ln_fugacity_coefficients_Trhomolefracs(const Model& model, const Scalar& T, const Scalar& rhotot, const VectorType& molefrac) {
auto R = model.R(molefrac);
using tdx = TDXDerivatives<Model, Scalar, VectorType>;
auto Z = 1.0 + tdx::template get_Ar01<be>(model, T, rhotot, molefrac);
auto rhovec = (rhotot*molefrac).eval();
auto grad = build_Psir_gradient<be>(model, T, rhovec).eval();
auto RT = R * T;
auto lnphi = ((grad / RT).array() - log(Z)).eval();
return forceeval(lnphi.eval());
}
template<ADBackends be = ADBackends::autodiff>
static auto get_ln_fugacity_coefficients1(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
auto R = model.R(molefrac);
auto grad = build_Psir_gradient<be>(model, T, rhovec).eval();
auto RT = R * T;
auto lnphi = ((grad / RT).array()).eval();
return forceeval(lnphi);
}
template<ADBackends be = ADBackends::autodiff>
static auto get_ln_fugacity_coefficients2(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
using tdx = TDXDerivatives<Model, Scalar, VectorType>;
auto Z = 1.0 + tdx::template get_Ar01<be>(model, T, rhotot, molefrac);
return forceeval(-log(Z));
}
static Eigen::ArrayXd build_d2PsirdTdrhoi_autodiff(const Model& model, const Scalar& T, const VectorType& rho) {
Eigen::ArrayXd deriv(rho.size());
// d^2psir/dTdrho_i
for (auto i = 0; i < rho.size(); ++i) {
auto psirfunc = [&model, &rho, i](const auto& T, const auto& rhoi) {
ArrayXdual2nd rhovecc(rho.size()); for (auto j = 0; j < rho.size(); ++j) { rhovecc[j] = rho[j]; }
rhovecc[i] = rhoi;
auto rhotot_ = rhovecc.sum();
auto molefrac = (rhovecc / rhotot_).eval();
return eval(model.alphar(T, rhotot_, molefrac) * model.R(molefrac) * T * rhotot_);
};
dual2nd Tdual = T, rhoidual = rho[i];
auto [u00, u10, u11] = derivatives(psirfunc, wrt(Tdual, rhoidual), at(Tdual, rhoidual));
deriv[i] = u11;
}
return deriv;
}
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
static Eigen::ArrayXd build_d2alphardrhodxi_constT(const Model& model, const Scalar& T, const Scalar& rhomolar, const VectorType& molefrac) {
Eigen::ArrayXd deriv(molefrac.size());
// d^2alphar/drhodx_i|T
for (auto i = 0; i < molefrac.size(); ++i) {
auto alpharfunc = [&model, &T, &molefrac, i](const auto& rho, const auto& xi) {
ArrayXdual2nd molefracdual = molefrac.template cast<autodiff::dual2nd>();
molefracdual[i] = xi;
return forceeval(model.alphar(T, rho, molefracdual));
};
dual2nd rhodual = rhomolar, xidual = molefrac[i];
auto [u00, u10, u11] = derivatives(alpharfunc, wrt(rhodual, xidual), at(rhodual, xidual));
deriv[i] = u11;
}
return deriv;
}
/***
* \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. temperature at constant mole concentrations (implying constant mole fractions and density)
*
* Uses autodiff to calculate derivatives by default
* \f[
* \deriv{ \ln\vec\phi}{T}{\vec \rho}
* \f]
*/
template<ADBackends be = ADBackends::autodiff>
static auto get_d_ln_fugacity_coefficients_dT_constrhovec(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
auto R = model.R(molefrac);
using tdx = TDXDerivatives<Model, Scalar, VectorType>;
auto Z = 1.0 + tdx::template get_Ar01<be>(model, T, rhotot, molefrac);
auto dZdT_Z = tdx::template get_Ar11<be>(model, T, rhotot, molefrac)/(-T)/Z; // Note: (1/T)dX/d(1/T) = -TdX/dT, the deriv in RHS is what we want, the left is what we get, so divide by -T
VectorType grad = build_Psir_gradient<be>(model, T, rhovec).eval();
VectorType Tgrad = build_d2PsirdTdrhoi_autodiff(model, T, rhovec);
return forceeval((1/(R*T)*(Tgrad - 1.0/T*grad)-dZdT_Z).eval());
}
/***
* \brief Calculate ln(Z), Z, and dZ/drho at constant temperature and mole fractions
*
* Uses autodiff to calculate derivatives by default
*/
template<ADBackends be = ADBackends::autodiff>
static auto get_lnZ_Z_dZdrho(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
using tdx = TDXDerivatives<Model, Scalar, VectorType>;
auto Ar01 = tdx::template get_Ar01<be>(model, T, rhotot, molefrac);
auto Ar02 = tdx::template get_Ar02<be>(model, T, rhotot, molefrac);
auto Z = 1.0 + Ar01;
auto dZdrho = (Ar01 + Ar02)/rhotot; // (dZ/rhotot)_{T,x}
return std::make_tuple(log(Z), Z, dZdrho);
}
/***
* \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. molar density at constant temperature and mole fractions
*
* \f[
* \deriv{ \ln\vec\phi}{\rho}{T,\vec x} = \frac{1}{RT}H(\Psi_r)\vec x - \frac{1}{Z}\deriv{Z}{\rho}{T,\vec x}
* \f]
*/
template<ADBackends be = ADBackends::autodiff>
static auto get_d_ln_fugacity_coefficients_drho_constTmolefracs(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
auto R = model.R(molefrac);
auto [lnZ, Z, dZdrho] = get_lnZ_Z_dZdrho(model, T, rhovec);
auto hessian = build_Psir_Hessian_autodiff(model, T, rhovec);
return forceeval((1/(R*T)*(hessian*molefrac.matrix()).array() - dZdrho/Z).eval());
}
/***
* \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. molar volume at constant temperature and mole fractions
*
* \f[
* \deriv{ \ln\vec\phi}{v}{T,\vec x} = \deriv{ \ln\vec\phi}{\rho}{T,\vec x}\deriv{ \rho}{v}{}
* \f]
*/
template<ADBackends be = ADBackends::autodiff>
static auto get_d_ln_fugacity_coefficients_dv_constTmolefracs(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto drhodv = -rhotot*rhotot; // rho = 1/v; drho/dv = -1/v^2 = -rho^2
return get_d_ln_fugacity_coefficients_drho_constTmolefracs(model, T, rhovec)*drhodv;
}
/***
* \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. mole fraction of each component, at constant temperature and molar density
*
* \f[
* \deriv{ \ln\vec\phi}{\vec x}{T, \rho}
* \f]
*/
template<ADBackends be = ADBackends::autodiff>
static auto get_d_ln_fugacity_coefficients_dmolefracs_constTrho(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
auto R = model.R(molefrac);
using tdx = TDXDerivatives<Model, Scalar, VectorType>;
auto Ar01 = tdx::template get_Ar01<be>(model, T, rhotot, molefrac);
auto Z = 1.0 + Ar01;
VectorType dZdx = rhotot*build_d2alphardrhodxi_constT(model, T, rhotot, molefrac);
Eigen::RowVector<decltype(rhotot), Eigen::Dynamic> dZdx_Z = dZdx/Z;
// Starting matrix is from the first term
auto hessian = build_Psir_Hessian_autodiff(model, T, rhovec);
Eigen::ArrayXXd out = rhotot/(R*T)*hessian;
// Then each row gets the second part
out.rowwise() -= dZdx_Z.array();
return out;
}
template<ADBackends be = ADBackends::autodiff>
static auto get_d_ln_fugacity_coefficients_dmolefracs_constTrho1(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = forceeval(rhovec.sum());
auto molefrac = (rhovec / rhotot).eval();
auto R = model.R(molefrac);
auto hessian = build_Psir_Hessian_autodiff(model, T, rhovec);
// Starting matrix is from the first term
Eigen::ArrayXXd out = 1/(R*T)*rhotot*hessian;
return out;
}
/***
* \brief Calculate the temperature derivative of the chemical potential of each component
* \note: Some contributions to the ideal gas part are missing (reference state and cp0), but are not relevant to phase equilibria
*/
static auto get_dchempotdT_autodiff(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = rhovec.sum();
auto molefrac = (rhovec / rhotot).eval();
auto rhorefideal = 1.0;
return (build_d2PsirdTdrhoi_autodiff(model, T, rhovec) + model.R(molefrac)*(rhorefideal + log(rhovec/rhorefideal))).eval();
}
/***
* \brief Calculate the temperature derivative of the pressure at constant molar concentrations
*/
static auto get_dpdT_constrhovec(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = rhovec.sum();
auto molefrac = (rhovec / rhotot).eval();
auto dPsirdT = get_dPsirdT_constrhovec(model, T, rhovec);
return rhotot*model.R(molefrac) - dPsirdT + rhovec.matrix().dot(build_d2PsirdTdrhoi_autodiff(model, T, rhovec).matrix());
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
/***
* \brief Calculate the molar concentration derivatives of the pressure at constant temperature
*/
static auto get_dpdrhovec_constT(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = rhovec.sum();
auto molefrac = (rhovec / rhotot).eval();
auto RT = model.R(molefrac)*T;
auto [func, grad, hessian] = build_Psir_fgradHessian_autodiff(model, T, rhovec); // The hessian matrix
return (RT + (hessian*rhovec.matrix()).array()).eval(); // at constant temperature
}
/***
* \brief Calculate the partial molar volumes of each component
*
* \f[
* \hat v_i = \left(\frac{\partial V}{\partial n_i}\right)_{T,V,n_{j \neq i}}
* \f]
*/
static auto get_partial_molar_volumes(const Model& model, const Scalar& T, const VectorType& rhovec) {
auto rhotot = rhovec.sum();
auto molefrac = (rhovec / rhotot).eval();
using tdx = TDXDerivatives<Model, Scalar, VectorType>;
auto RT = model.R(molefrac)*T;
auto dpdrhovec = get_dpdrhovec_constT(model, T, rhovec);
auto numerator = -rhotot*dpdrhovec;
auto ders = tdx::template get_Ar0n<2>(model, T, rhotot, molefrac);
auto denominator = -pow2(rhotot)*RT*(1 + 2*ders[1] + ders[2]);
return (numerator/denominator).eval();
}
static VectorType get_Psir_sigma_derivs(const Model& model, const Scalar& T, const VectorType& rhovec, const VectorType& v) {
autodiff::Real<4, double> sigma = 0.0;
auto rhovecad = rhovec.template cast<decltype(sigma)>(), vad = v.template cast<decltype(sigma)>();
auto wrapper = [&rhovecad, &vad, &T, &model](const auto& sigma_1) {
auto rhovecused = (rhovecad + sigma_1 * vad).eval();
auto rhotot = rhovecused.sum();
auto molefrac = (rhovecused / rhotot).eval();
return forceeval(model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot);
};
auto der = derivatives(wrapper, along(1), at(sigma));
for (auto i = 0; i < ret.size(); ++i){ ret[i] = der[i];}
return ret;
}
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
template<int Nderivsmax, AlphaWrapperOption opt>
class DerivativeHolderSquare{
public:
Eigen::Array<double, Nderivsmax+1, Nderivsmax+1> derivs;
template<typename Model, typename Scalar, typename VecType>
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);
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);
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);
}
};