[RFC PATCH] gas model: use virial cubic polynomial form
Linus Torvalds
torvalds at linuxfoundation.org
Thu Mar 3 21:46:19 PST 2016
From: Linus Torvalds <torvalds at linuxfoundation.org>
Date: Thu, 3 Mar 2016 21:30:28 0800
Subject: [PATCH] gas model: use virial cubic polynomial form
The "virial" form of the Z compression factor is of the form
Z = 1.0 + A*p + B*p^2 + C*p^3 + ..
and it's considered the "right" polynomial form to use. It happens to
also make for one constant less per gas (since the 1.0 can be added
later), and can be used to simplify the expression and avoid a few
floating point operations.
However, in order for that kind of expression simplification to make
sense, we need to make sure that we don't calculate the powers of the
pressure multiple times either, and that means we have to inline all the
actual calculations.
Our compiler options still mean that the generated code isn't optimal,
but that's a separate issue. And it is a lot better than it used to be.
Being clever about this does potentially make the code a tiny bit less
legible, but maybe that's not too bad for something that we'd expect to
not ever touch once we get it right.
Signedoffby: Linus Torvalds <torvalds at linuxfoundation.org>

NOTE THE "RFC" part. I'm not sure this is worth it. It came from the G+
discussion: the virial form is apparently the preferred form for the Z
calculation and the difference to the pure cubic leastsquare is minimal.
And it does generate better code, just because there are fewer constants
and fewer FP operations and we can delay adding the 1.0 term until after
the gas mix, rather than doing it for each gas.
At the same time, I'm not exactly claiming that this really helps make the
code obvious. Maybe it's premature optimization. And really, considering
that we don't even inline the "get_o2()" trivial function etc, the code
generation really isn't very optional.
So for this to make sense, we'd need to feel that performance matters. We
*do* end up doing the gas use for each plot point, so it's not like this
function might not get called a lot.
But compiler options will change the code generation even more than this
ever will, and aggressive compiler optimizations would have made the old
code better too (ie just inlining and CSE on the pressure power series
would do good things to the thing without this patch too).
In other words: I think this patch is correct, and it *does* make the code
faster. But do we even care at this point?
subsurfacecore/gasmodel.c  98 +++++++++++++++++++
1 file changed, 42 insertions(+), 56 deletions()
diff git a/subsurfacecore/gasmodel.c b/subsurfacecore/gasmodel.c
index 420233612cb9..81765e003002 100644
 a/subsurfacecore/gasmodel.c
+++ b/subsurfacecore/gasmodel.c
@@ 4,72 +4,58 @@
#include <stdlib.h>
#include "dive.h"
/*
 * Generic cubic polynomial
 */
static double cubic(double bar, const double coefficient[])
{
 double x0 = 1.0,
 x1 = bar,
 x2 = x1*x1,
 x3 = x2*x1;

 return x0 * coefficient[0] +
 x1 * coefficient[1] +
 x2 * coefficient[2] +
 x3 * coefficient[3];
}
+/* "Virial minus one"  the virial cubic form without the initial 1.0 */
+#define virial_m1(C, x1, x2, x3) (C[0]*x1+C[1]*x2+C[2]*x3)
/*
 * Cubic leastsquare coefficients for O2/N2/He based on data from
+ * Cubic virial leastsquare coefficients for O2/N2/He based on data from
*
* PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION
*
* with the lookup and curve fitting by Lubomir.
*
+ * The "virial" form of the compression factor polynomial is
+ *
+ * Z = 1.0 + C[0]*P + C[1]*P^2 + C[2]*P^3 ...
+ *
+ * and these tables do not contain the initial 1.0 term.
+ *
* NOTE! Helium coefficients are a linear mix operation between the
* 323K and one for 273K isotherms, to make everything be at 300K.
*/
static const double o2_coefficients[4] = {
 +1.00117935180448264158,
 0.00074149079841471884,
 +0.00000291901111247509,
 0.00000000162092217461
};

static const double n2_coefficients[4] = {
 +1.00030344355797817778,
 0.00022528077251905598,
 +0.00000295430303276288,
 0.00000000210649996337
};

static const double he_coefficients[4] = {
 +1.00000137322788319261,
 +0.000488393024887620665,
 0.000000054210166760015,
 +0.000000000010908069275
};

static double o2_compressibility_factor(double bar) { return cubic(bar, o2_coefficients); }
static double n2_compressibility_factor(double bar) { return cubic(bar, n2_coefficients); }
static double he_compressibility_factor(double bar) { return cubic(bar, he_coefficients); }

/*
 * We end up using a simple linear mix of the gasspecific functions.
 */
double gas_compressibility_factor(struct gasmix *gas, double bar)
{
 double o2, n2, he; // Z factors
 double of, nf, hf; // gas fractions ("partial pressures")

 o2 = o2_compressibility_factor(bar);
 n2 = n2_compressibility_factor(bar);
 he = he_compressibility_factor(bar);

 of = get_o2(gas) / 1000.0;
 hf = get_he(gas) / 1000.0;
 nf = 1.0  of  nf;

 return o2*of + n2*nf + he*hf;
+ static const double o2_coefficients[3] = {
+ 0.00071809207370164567,
+ +0.00000281852572807643,
+ 0.00000000150290620491
+ };
+ static const double n2_coefficients[3] = {
+ 0.00021926035329221337,
+ +0.00000292844845531647,
+ 0.00000000207613482075
+ };
+ static const double he_coefficients[3] = {
+ +0.00047961098687979363,
+ 0.00000004077670019935,
+ +0.00000000000077707035
+ };
+ double o2, he;
+ double x1, x2, x3;
+ double Z;
+
+ o2 = get_o2(gas) / 1000.0;
+ he = get_he(gas) / 1000.0;
+
+ x1 = bar; x2 = x1*x1; x3 = x2*x1;
+
+ Z = virial_m1(o2_coefficients, x1, x2, x3) * o2 +
+ virial_m1(he_coefficients, x1, x2, x3) * he +
+ virial_m1(n2_coefficients, x1, x2, x3) * (1.0  o2  he);
+
+ /*
+ * We add the 1.0 at the very end  the linear mixing of the
+ * three 1.0 terms is still 1.0 regardless of the gas mix.
+ */
+ return Z + 1.0;
}

2.7.2.334.g7c0da37
More information about the subsurface
mailing list