Skip to content

Commit

Permalink
Add docs for comp derivs; add a few more derivatives
Browse files Browse the repository at this point in the history 
Closes #101
  • Loading branch information
@ianhbell
ianhbell committed Mar 3, 2024
1 parent 2699b9c commit dc8031d
Show file tree
Hide file tree
Showing 3 changed files with 294 additions and 2 deletions.
278 changes: 278 additions & 0 deletions doc/source/derivs/compderivs.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,278 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "f751a5e0",
"metadata": {},
"source": [
"# Composition derivatives\n",
"\n",
"For other mixture calculations composition derivatives of the form\n",
"\n",
"$$ \\Lambda^r_{xyz_i} = (1/T)^x(\\rho)^y\\left(\\frac{\\partial^{x+y+z_i}(\\alpha^r)}{\\partial (1/T)^x\\partial \\rho^y \\partial \\mathbf{Z}_i^{z_i}} \\right) $$\n",
"\n",
"are needed. This function is exposed in teqp (as of version 0.19) as the function ``get_ATrhoXi``. In order to limit the binary size and compilation time, x has a max of 2 and y does as well. ``z_i`` can be up to 3, and must be at least 1, otherwise you can use the other derivative functions that do not require any composition derivatives.\n",
"\n",
"The mixed composition derivative of the form\n",
"\n",
"$$\\Lambda^r_{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)$$\n",
"\n",
"supports $x$ and $y$ either 0 or 1, with at most two composition derivatives. The triple composition derivative of the form\n",
"\n",
"$$\\Lambda^r_{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)$$\n",
"\n",
"supports $x$ and $y$ either 0 or 1, with up to first derivatives in each composition variable. If this is not enough derivatives, open a feature request here : https://github.com/usnistgov/teqp/issues\n",
"\n",
"## $\\tau$ and $\\delta$ derivatives\n",
"\n",
"In the multi-fluid modeling approach used in NIST REFPROP and the GERG-2004 GERG-2008 models, the derivatives are in the form\n",
"\n",
"$$ \\Lambda^r_{xyz_i} = \\tau^x\\delta^y\\left(\\frac{\\partial^{x+y+z_i}(\\alpha^r)}{\\partial \\tau^x\\partial \\delta^y \\partial \\mathbf{Z}_i^{z_i}} \\right) $$\n",
"\n",
"with $\\tau=T_{\\rm red}(\\mathbf{Z})/T$ and $\\delta=\\rho/\\rho_{\\rm red}(\\mathbf{Z})$. The higher derivatives are similarly equal to\n",
"\n",
"$$\\Lambda^r_{xyz_i z_j } = \\tau^x\\delta^y\\left(\\frac{\\partial^{x+y+z_i+z_j}(\\alpha^r)}{\\partial \\tau^x\\partial \\delta^y \\partial \\mathbf{Z}_i^{z_i} \\partial \\mathbf{Z}_j^{z_j} }\\right)$$\n",
"\n",
"$$\\Lambda^r_{xyz_i z_j z_k} = \\tau^x\\delta^y\\left(\\frac{\\partial^{x+y+z_i+z_j+z_k}(\\alpha^r)}{\\partial \\tau^x\\partial \\delta^y \\partial \\mathbf{Z}_i^{z_i} \\partial \\mathbf{Z}_j^{z_j} \\partial \\mathbf{Z}_k^{z_k} }\\right)$$\n",
"\n",
"The same limitations on the numbers of derivatives are used for the derivatives with $(1/T)$ and $\\rho$ as independent variables."
]
},
{
"cell_type": "raw",
"id": "59ede80a",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
"The Python methods are documented here: \n",
"\n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_ATrhoXi` \n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_ATrhoXiXj` \n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_ATrhoXiXjXk` \n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_AtaudeltaXi` \n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_AtaudeltaXiXj` \n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_AtaudeltaXiXjXk`"
]
},
{
"cell_type": "markdown",
"id": "d46932aa",
"metadata": {},
"source": [
"## xN (in)dependent\n",
"\n",
"Let's suppose that some quantity $\\Upsilon$ depends on mole fractions. If all the mole fractions are considered to be independent, the total differential is obtained from\n",
"\n",
"$\\newcommand{\\deriv}[3]{\\left(\\frac{\\partial #1}{\\partial #2}\\right)_{#3}}$\n",
"\n",
"$$\n",
"{\\rm d} \\Upsilon = \\sum_j\\deriv{\\Upsilon}{x_j}{x_{i\\neq j}} {\\rm d}x_j\n",
"$$\n",
"\n",
"If instead the last mole fraction is defined to be dependent on the others via\n",
"$$\n",
"x_N = 1-\\sum_{i=1}^{N-1}x_i\n",
"$$\n",
"then the total differential is obtained from\n",
"\n",
"$$\n",
" {\\rm d} \\Upsilon = \\sum_{j=1}^{N-1}\\deriv{\\Upsilon}{x_j}{x_{i\\neq j}} {\\rm d}x_j + \\deriv{\\Upsilon}{x_N}{x_{i\\neq j}} {\\rm d}{x_N}\n",
"$$\n",
"where $x_N$ is considered to be an independent variable in the derivative $\\deriv{\\Upsilon}{x_N}{x_{i\\neq j}}$. Thus derivatives with respect to one of the dependent mole fractions ($x_k$ with $k < N$) would be equal to \n",
"$$\n",
"\\deriv{\\Upsilon}{x_k}{x_{j\\neq k}} =\\deriv{\\Upsilon}{x_k}{x_{i\\neq k}} - \\deriv{\\Upsilon}{x_N}{x_{i\\neq N}} \n",
"$$\n",
"because \n",
"$$\n",
"\\deriv{x_N}{x_i}{} = -1\n",
"$$\n",
"\n",
"So if the library (e.g., CoolProp and TREND) allows for the fractions to be independent, you can use the molar composition derivatives with the mole fractions treated as being indepenent to obtain derivatives with one of the mole fractions dependent on the other $N-1$ fractions."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "2ec44ccf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.18.1.dev0'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import teqp, numpy as np\n",
"teqp.__version__"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "53ba18fe",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-0.04358738425351122\n",
"-0.21188579988125827\n",
"-0.03650566667904926\n",
"-0.07488856488580686\n",
"-0.20693890096529247\n",
"0.014468933385218725\n",
"-0.005978809921279951\n",
"-0.0027918555000108207\n",
"0.0\n"
]
}
],
"source": [
"j = {\n",
" 'components': [\"Methane\", \"Nitrogen\", \"Oxygen\"],\n",
" 'root': teqp.get_datapath(),\n",
" 'BIP': '',\n",
" 'departure': ''\n",
"}\n",
"model = teqp.make_model({'kind':'multifluid', 'model': j})\n",
"\n",
"\n",
"T = 300 # K\n",
"rhomolar = 3000 # mol/m^3\n",
"z = np.array([0.3, 0.5, 0.2]) # mole fractions\n",
"\n",
"Tr = model.get_Tr(z)\n",
"rhor = model.get_rhor(z)\n",
"tau = Tr/T\n",
"delta = rhomolar/rhor\n",
"\n",
"Ntau = 0\n",
"Ndelta = 0\n",
"Nxi = 1\n",
"print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))\n",
"\n",
"Ntau = 1\n",
"Ndelta = 0\n",
"Nxi = 1\n",
"print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))\n",
"\n",
"Ntau = 0\n",
"Ndelta = 1\n",
"Nxi = 1\n",
"print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))\n",
"\n",
"Ntau = 2\n",
"Ndelta = 0\n",
"Nxi = 1\n",
"print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))\n",
"\n",
"Ntau = 1\n",
"Ndelta = 1\n",
"Nxi = 1\n",
"print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))\n",
"\n",
"Ntau = 0\n",
"Ndelta = 2\n",
"Nxi = 1\n",
"print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))\n",
"\n",
"Ntau = 1\n",
"Ndelta = 0\n",
"Nxi = 1\n",
"Nxj = 1\n",
"print(model.get_AtaudeltaXiXj(tau, Ntau, delta, Ndelta, z, 0, Nxi, 1, Nxj))\n",
"\n",
"Ntau = 0\n",
"Ndelta = 1\n",
"Nxi = 1\n",
"Nxj = 1\n",
"print(model.get_AtaudeltaXiXj(tau, Ntau, delta, Ndelta, z, 0, Nxi, 1, Nxj))\n",
"\n",
"Ntau = 0\n",
"Ndelta = 0\n",
"Nxi = 1\n",
"Nxj = 1\n",
"Nxk = 1\n",
"print(model.get_AtaudeltaXiXjXk(tau, Ntau, delta, Ndelta, z, 0, Nxi, 1, Nxj, 2, Nxk))"
]
},
{
"cell_type": "markdown",
"id": "594e48be",
"metadata": {},
"source": [
"With CoolProp, version 6.6.0, the following script in C++:\n",
"```C++\n",
"#include \"AbstractState.h\"\n",
"#include \"Backends/Helmholtz/MixtureDerivatives.h\"\n",
"\n",
"int main(){\n",
" std::shared_ptr<CoolProp::AbstractState> AS(\n",
" CoolProp::AbstractState::factory(\"HEOS\",\"Methane&Nitrogen&Oxygen\")\n",
" );\n",
" AS->set_mole_fractions({0.3, 0.5, 0.2});\n",
" AS->specify_phase(CoolProp::iphase_gas);\n",
" AS->update(CoolProp::DmolarT_INPUTS, 3000, 300);\n",
" auto& HEOS = *dynamic_cast<CoolProp::HelmholtzEOSMixtureBackend*>(AS.get());\n",
" auto xN = CoolProp::x_N_dependency_flag::XN_INDEPENDENT;\n",
" using md = CoolProp::MixtureDerivatives;\n",
" std::cout << md::dalphar_dxi(HEOS, 0, xN) << std::endl;\n",
" \n",
" std::cout << md::d2alphar_dxi_dTau(HEOS, 0, xN)*AS->tau() << std::endl;\n",
" std::cout << md::d2alphar_dxi_dDelta(HEOS, 0, xN)*AS->delta() << std::endl;\n",
" std::cout << md::d3alphar_dxi_dTau2(HEOS, 0, xN)*pow(AS->tau(), 2) << std::endl;\n",
" std::cout << md::d3alphar_dxi_dDelta_dTau(HEOS, 0, xN)*AS->tau()*AS->delta() << std::endl;\n",
" std::cout << md::d3alphar_dxi_dDelta2(HEOS, 0, xN)*pow(AS->delta(), 2) << std::endl;\n",
" \n",
" std::cout << md::d3alphar_dxi_dxj_dTau(HEOS, 0, 1, xN)*AS->tau() << std::endl;\n",
" std::cout << md::d3alphar_dxi_dxj_dDelta(HEOS, 0, 1, xN)*AS->delta() << std::endl;\n",
" std::cout << md::d3alphardxidxjdxk(HEOS, 0, 1, 2, xN) << std::endl;\n",
"}\n",
"```\n",
"yields the output:\n",
"```\n",
"-0.0435874\n",
"-0.211886\n",
"-0.0365057\n",
"-0.0748886\n",
"-0.206939\n",
"0.0144689\n",
"-0.00597881\n",
"-0.00279186\n",
"0\n",
"```\n",
"which is the same as the above"
]
}
],
"metadata": {
"celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
1 change: 1 addition & 0 deletions doc/source/derivs/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,4 +6,5 @@ Derivatives
:caption: Contents:

derivs
compderivs
ideal_gas
17 changes: 15 additions & 2 deletions include/teqp/derivs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -258,15 +258,25 @@ struct TDXDerivatives {
return powi(forceeval(1.0 / T), iT) * powi(rho, iD) * der[der.size() - 1];
}

#define get_ATrhoXi_runtime_combinations X(0,0,1) \
#define get_ATrhoXi_runtime_combinations \
X(0,0,1) \
X(0,0,2) \
X(0,0,3) \
X(1,0,1) \
X(1,0,2) \
X(1,0,3) \
X(0,1,1) \
X(0,1,2) \
X(0,1,3)
X(0,1,3) \
X(2,0,1) \
X(2,0,2) \
X(2,0,3) \
X(1,1,1) \
X(1,1,2) \
X(1,1,3) \
X(0,2,1) \
X(0,2,2) \
X(0,2,3)

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){
Expand Down Expand Up @@ -305,12 +315,15 @@ struct TDXDerivatives {
X(0,0,0,1,1) \
X(0,0,1,0,1) \
X(0,0,1,1,0) \
X(0,0,1,1,1) \
X(1,0,0,1,1) \
X(1,0,1,0,1) \
X(1,0,1,1,0) \
X(1,0,1,1,1) \
X(0,1,0,1,1) \
X(0,1,1,0,1) \
X(0,1,1,1,0) \
X(0,1,1,1,1)

template<typename AlphaWrapper>
static auto get_ATrhoXiXjXk_runtime(const AlphaWrapper& w, const Scalar& T, int iT, const Scalar& rho, int iD, const VectorType& molefrac, int i, int iXi, int j, int iXj, int k, int iXk){
Expand Down

0 comments on commit dc8031d

Please sign in to comment.