[PATCH] gas-model: add R compressibility script

[Linus Torvalds]

From: Linus Torvalds <torvalds at linux-foundation.org>
Date: Wed, 9 Mar 2016 13:14:05 -0800
Subject: [PATCH] gas-model: add R compressibility script

It annoyed me that we hand-waved a bit about how the virial factors were
actually computed in the gas-model.c file, so here's an actual R script
that computes them and plots the results.

You can run it with (for example):

    R --vanilla < compressibility.r

and it will generate a Rplots.pdf of the plots, and the coefficients
will be shown on stdout.

The result actually differs in insignificant ways from the values that
Lubomir computed, which is likely just due to tools.  I used R, Lubomir
seems to have used

    http://polynomialregression.drque.net/online.php

but the actual curve is pretty much the same.

NOTE! R is not entirely happy about the non-linear fit of the Helium
curve: the fit is *so* precise that it failes the R relative-offset
convergence criterion.  That is apparently generally a sign of
artificial data.

That is probably because Lubomir generated them from the linear mix of
two polynomial fits, rather than a linear mix of the original data.  But
maybe the original data was artificial?

Signed-off-by: Linus Torvalds <torvalds at linux-foundation.org>
---

This generates no code, nor did I even update the existing coefficients 
with the slightly different coefficients I get using 'R'. However, it 
really did bother me that I had no way of reproducing the magic that 
Lubomir had done, and that our source code just had those magic numbers in 
them.

At least this _explains_ the magic numbers, and if we ever use more or 
different pressure data, this has the scripts to generate the virial 
factors.

NOTE! I'm not an R person. I've never used R in my life before this 
exercise. A real R person would probably do this much better. This is a 
lot of trial and error and reading examples.

 subsurface-core/compressibility.r | 78 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)
 create mode 100644 subsurface-core/compressibility.r

diff --git a/subsurface-core/compressibility.r b/subsurface-core/compressibility.r
new file mode 100644
index 000000000000..15aca6ea27c9
--- /dev/null
+++ b/subsurface-core/compressibility.r
@@ -0,0 +1,78 @@
+# Compressibility data gathered by Lubomir I Ivanov:
+#
+#  "Data obtained by finding two books online:
+#
+#   [1]
+#    PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION
+#    pretty serious book, from which the wiki AIR values come from!
+#
+#     http://www.unhas.ac.id/rhiza/arsip/kuliah/Sistem-dan-Tekn-Kendali-Proses/PDF_Collections/REFERENSI/Perrys_Chemical_Engineering_Handbook.pdf
+#     page 2-165
+#
+#     [*](Computed from pressure-volume-temperature tables in Vasserman monographs)
+#     ^ i have no idea idea what this means, but the values might not be exactly
+#     experimental?!
+#
+#   the only thing this book is missing is helium, thus [2]!
+#
+#   [2]
+#    VOLUMETRIC BEHAVIOR OF HELIUM-ARGON MIXTURES AT HIGH PRESSURE AND MODERATE TEMPERATURE.
+#
+#     https://shareok.org/bitstream/handle/11244/2062/6614196.PDF?sequence=1
+#     page 108
+#
+#
+#     the book has some tables with pressure values in atmosphere units. i'm
+#     converting them bars. one of the relevant tables is for 323K and one for 273K
+#     (both almost equal distance from 300K).
+#
+#     this again is a linear mix operation between isotherms, which is probably not
+#     the most accurate solution but it works.
+#
+# all data sets contain Z values at 300k, while the pressures are in bars in
+# the 1 to 500 range
+#
+#
+
+x  = c(1, 5, 10, 20, 40, 60, 80, 100, 200, 300, 400, 500)
+o2 = c(0.9994, 0.9968, 0.9941, 0.9884, 0.9771, 0.9676, 0.9597, 0.9542, 0.9560, 0.9972, 1.0689, 1.1572)
+n2 = c(0.9998, 0.9990, 0.9983, 0.9971, 0.9964, 0.9973, 1.0000, 1.0052, 1.0559, 1.1422, 1.2480, 1.3629)
+he = c(1.0005, 1.0024, 1.0048, 1.0096, 1.0191, 1.0286, 1.0381, 1.0476, 1.0943, 1.1402, 1.1854, 1.2297)
+
+
+options(digits=15)
+
+#
+# Get the O2 virial coefficients
+#
+plot(x,o2)
+o2fit = nls(o2 ~ 1.0 + p1*x + p2 *x^2 + p3*x^3, start=list(p1=0,p2=0,p3=0))
+summary(o2fit)
+
+new = data.frame(x = seq(min(x),max(x),len=200))
+lines(new$x,predict(o2fit,newdata=new))
+
+#
+# Get the N2 virial coefficients
+#
+plot(x,n2)
+n2fit = nls(n2 ~ 1.0 + p1*x + p2 *x^2 + p3*x^3, start=list(p1=0,p2=0,p3=0))
+summary(n2fit)
+
+new = data.frame(x = seq(min(x),max(x),len=200))
+lines(new$x,predict(n2fit,newdata=new))
+
+#
+# Get the He virial coefficients
+#
+# NOTE! This will not confirm convergence, thus the warnOnly.
+# That may be a sign that the data is possibly artificial.
+#
+plot(x,he)
+hefit = nls(he ~ 1.0 + p1*x + p2 *x^2 + p3*x^3,
+	start=list(p1=0,p2=0,p3=0),
+	control=nls.control(warnOnly=TRUE))
+summary(hefit)
+
+new = data.frame(x = seq(min(x),max(x),len=200))
+lines(new$x,predict(hefit,newdata=new))
-- 
2.7.2.334.g7c0da37