From afadb6f05cd0477976a04fad62122523c0c90420 Mon Sep 17 00:00:00 2001 From: Linus Torvalds Date: Thu, 3 Mar 2016 21:30:28 -0800 Subject: 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. Signed-off-by: Linus Torvalds Signed-off-by: Dirk Hohndel --- subsurface-core/gas-model.c | 98 +++++++++++++++++++-------------------------- 1 file changed, 42 insertions(+), 56 deletions(-) (limited to 'subsurface-core/gas-model.c') diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c index 420233612..81765e003 100644 --- a/subsurface-core/gas-model.c +++ b/subsurface-core/gas-model.c @@ -4,72 +4,58 @@ #include #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 least-square coefficients for O2/N2/He based on data from + * Cubic virial least-square 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 gas-specific 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; } -- cgit v1.2.3-70-g09d2