Add a 'fast'-mode, surpassing robustness code (#59)

In normal mode: WaterOilGas construction and SWOF+SGOF: 66 ms
In fast mode: 36 ms

pyscal/gasoil.py

@@ -53,10 +53,21 @@ class GasOil(object):
             it does not matter.
         h (float): Saturation step-length in the outputted table.
         tag (str): Optional string identifier, only used in comments.
+        fast (bool): Set to true to prioritize speed over robustness. Not guaranteed
+            to give you include files that Eclipse does not crash on - you are essentially
+            on your own. Default false.
     """
 
     def __init__(
-        self, swirr=0, sgcr=0.0, h=0.01, swl=0.0, sorg=0.0, tag="", krgendanchor="sorg"
+        self,
+        swirr=0,
+        sgcr=0.0,
+        h=0.01,
+        swl=0.0,
+        sorg=0.0,
+        tag="",
+        krgendanchor="sorg",
+        fast=False,
     ):
         assert epsilon < h <= 1
         assert -epsilon < swirr < 1.0 + epsilon
@@ -88,6 +99,9 @@ def __init__(
         else:
             logging.warning("Unknown krgendanchor %s, ignored", str(krgendanchor))
             self.krgendanchor = ""
+
+        self.fast = fast
+
         if np.isclose(sorg, 0.0) and self.krgendanchor == "sorg":
             # This is too much info..
             self.krgendanchor = ""  # This is critical to avoid bugs due to numerics.
@@ -104,34 +118,39 @@ def __init__(
             map(int, list(map(round, self.table["sg"] * SWINTEGERS)))
         )
         self.table.drop_duplicates("sgint", inplace=True)
-        # Now sg=1-sorg-swl might be accidentally dropped, so make sure we
-        # have it by replacing the closest value by 1 - sorg exactly
-        sorgindex = (
-            (self.table["sg"] - (1 - self.sorg - self.swl)).abs().sort_values().index[0]
-        )
-        self.table.loc[sorgindex, "sg"] = 1 - self.sorg - self.swl
-
-        # Same for sg=sgcr
-        sgcrindex = (self.table["sg"] - (self.sgcr)).abs().sort_values().index[0]
-        self.table.loc[sgcrindex, "sg"] = self.sgcr
-        if sgcrindex == 0 and sgcr > 0.0:
-            # Need to conserve sg=0
-            zero_row = pd.DataFrame({"sg": 0}, index=[0])
-            self.table = pd.concat([zero_row, self.table], sort=False).reset_index(
-                drop=True
+
+        if not self.fast:
+            # Now sg=1-sorg-swl might be accidentally dropped, so make sure we
+            # have it by replacing the closest value by 1 - sorg exactly
+            sorgindex = (
+                (self.table["sg"] - (1 - self.sorg - self.swl))
+                .abs()
+                .sort_values()
+                .index[0]
             )
+            self.table.loc[sorgindex, "sg"] = 1 - self.sorg - self.swl
+
+            # Same for sg=sgcr
+            sgcrindex = (self.table["sg"] - (self.sgcr)).abs().sort_values().index[0]
+            self.table.loc[sgcrindex, "sg"] = self.sgcr
+            if sgcrindex == 0 and sgcr > 0.0:
+                # Need to conserve sg=0
+                zero_row = pd.DataFrame({"sg": 0}, index=[0])
+                self.table = pd.concat([zero_row, self.table], sort=False).reset_index(
+                    drop=True
+                )
 
-        # If sg=1-swl was dropped, then sorg was close to zero:
-        if not np.isclose(self.table["sg"].max(), 1 - self.swl):
-            # Add it as an extra row:
-            self.table.loc[len(self.table) + 1, "sg"] = 1 - self.swl
-        # Ensure the value closest to 1-swl is actually 1-swl:
-        swl_right_index = (
-            (self.table["sg"] - (1 - self.swl)).abs().sort_values().index[0]
-        )
-        self.table.loc[swl_right_index, "sg"] = 1 - self.swl
+            # If sg=1-swl was dropped, then sorg was close to zero:
+            if not np.isclose(self.table["sg"].max(), 1 - self.swl):
+                # Add it as an extra row:
+                self.table.loc[len(self.table) + 1, "sg"] = 1 - self.swl
+            # Ensure the value closest to 1-swl is actually 1-swl:
+            swl_right_index = (
+                (self.table["sg"] - (1 - self.swl)).abs().sort_values().index[0]
+            )
+            self.table.loc[swl_right_index, "sg"] = 1 - self.swl
 
-        self.table.sort_values(by="sg", inplace=True)
+            self.table.sort_values(by="sg", inplace=True)
         self.table.reset_index(inplace=True)
         self.table = self.table[["sg"]]
         self.table["sl"] = 1 - self.table["sg"]
@@ -609,14 +628,14 @@ def SGOF(self, header=True, dataincommentrow=True):
         as Eclipse comments.
 
         Args:
-            header: boolean for whether the SGOF string should be emitted.
+            header (bool): Whether the SGOF string should be emitted.
                 If you have multiple satnums, you should have True only
                 for the first (or False for all, and emit the SGOF yourself).
                 Defaults to True.
-            dataincommentrow: boolean for wheter metadata should be printed,
+            dataincommentrow (bool): Whether metadata should be printed,
                 defaults to True.
         """
-        if not self.selfcheck():
+        if not self.fast and not self.selfcheck():
             return
         string = ""
         if "pc" not in self.table:
@@ -630,7 +649,8 @@ def SGOF(self, header=True, dataincommentrow=True):
             string += self.sgcomment
             string += self.krgcomment
             string += self.krogcomment
-            string += "-- krg = krog @ sw=%1.5f\n" % self.crosspoint()
+            if not self.fast:
+                string += "-- krg = krog @ sw=%1.5f\n" % self.crosspoint()
             string += self.pccomment
         string += self.table[["sg", "krg", "krog", "pc"]].to_csv(
             sep=" ", float_format="%1.7f", header=None, index=False

pyscal/wateroil.py

@@ -34,12 +34,18 @@ class WaterOil(object):
     parametrizations, or from a dataframe (will incur interpolation).
 
     Can be dumped as include files for Eclipse/OPM and Nexus simulators.
+
+    Args:
+        fast (bool): Set to True if in order to skip some integrity checks
+            and nice-to-have features. Not needed to set for normal pyscal
+            runs, as speed is seldom crucial. Default False
     """
 
-    def __init__(self, swirr=0.0, swl=0.0, swcr=0.0, sorw=0.0, h=0.01, tag=""):
+    def __init__(
+        self, swirr=0.0, swl=0.0, swcr=0.0, sorw=0.0, h=0.01, tag="", fast=False
+    ):
         """Sets up the saturation range. Swirr is only relevant
         for the capillary pressure, not for relperm data."""
-
         assert -epsilon < swirr < 1.0 + epsilon
         assert -epsilon < swl < 1.0 + epsilon
         assert -epsilon < swcr < 1.0 + epsilon
@@ -62,29 +68,33 @@ def __init__(self, swirr=0.0, swl=0.0, swcr=0.0, sorw=0.0, h=0.01, tag=""):
         self.sorw = sorw
         self.h = h
         self.tag = tag
+        self.fast = fast
         sw = list(np.arange(self.swl, 1 - sorw, h)) + [self.swcr] + [1 - sorw] + [1]
         self.table = pd.DataFrame(sw, columns=["sw"])
+
         # Ensure that we do not have sw values that are too close
         # to each other, determined rougly by the distance 1/10000
         self.table["swint"] = list(
             map(int, list(map(round, self.table["sw"] * SWINTEGERS)))
         )
         self.table.drop_duplicates("swint", inplace=True)
-        # Now, sw=1-sorw might be accidentaly dropped, so make sure we
-        # have it by replacing the closest value by 1-sorw exactly
-        sorwindex = (self.table["sw"] - (1 - self.sorw)).abs().sort_values().index[0]
-        self.table.loc[sorwindex, "sw"] = 1 - self.sorw
 
-        # Same for sw=swcr:
-        swcrindex = (self.table["sw"] - (self.swcr)).abs().sort_values().index[0]
-        self.table.loc[swcrindex, "sw"] = self.swcr
+        if not self.fast:
+            # Now, sw=1-sorw might be accidentaly dropped, so make sure we
+            # have it by replacing the closest value by 1-sorw exactly
+            sorwindex = (self.table["sw"] - (1 - self.sorw)).abs().sort_values().index[0]
+            self.table.loc[sorwindex, "sw"] = 1 - self.sorw
 
-        # If sw=1 was dropped, then sorw was close to zero:
-        if not np.isclose(self.table["sw"].max(), 1.0):
-            # Add it as an extra row:
-            self.table.loc[len(self.table) + 1, "sw"] = 1.0
+            # Same for sw=swcr:
+            swcrindex = (self.table["sw"] - (self.swcr)).abs().sort_values().index[0]
+            self.table.loc[swcrindex, "sw"] = self.swcr
 
-        self.table.sort_values(by="sw", inplace=True)
+            # If sw=1 was dropped, then sorw was close to zero:
+            if not np.isclose(self.table["sw"].max(), 1.0):
+                # Add it as an extra row:
+                self.table.loc[len(self.table) + 1, "sw"] = 1.0
+
+            self.table.sort_values(by="sw", inplace=True)
         self.table.reset_index(inplace=True)
         self.table = self.table[["sw"]]  # Drop the swint column
 
@@ -885,7 +895,24 @@ def selfcheck(self):
             return True
 
     def SWOF(self, header=True, dataincommentrow=True):
-        if not self.selfcheck():
+        """
+        Produce SWOF input for Eclipse reservoir simulator.
+
+        The columns sw, krw, krow and pc are outputted and
+        formatted accordingly.
+
+        Meta-information for the tabulated data are printed
+        as Eclipse comments.
+
+        Args:
+            header (bool): Indicate whether the SWOF string should
+                be emitted. If you have multiple SATNUMs, you should
+                set this to True only for the first (or False for all,
+                and emit the SWOF yourself). Default True
+            dataincommentrow (bool): Wheter metadata should
+                be printed. Defualt True
+        """
+        if not self.fast and not self.selfcheck():
             return
         string = ""
         if "pc" not in self.table.columns:
@@ -899,7 +926,8 @@ def SWOF(self, header=True, dataincommentrow=True):
             string += self.swcomment
             string += self.krwcomment
             string += self.krowcomment
-            string += "-- krw = krow @ sw=%1.5f\n" % self.crosspoint()
+            if not self.fast:
+                string += "-- krw = krow @ sw=%1.5f\n" % self.crosspoint()
             string += self.pccomment
         string += self.table[["sw", "krw", "krow", "pc"]].to_csv(
             sep=" ", float_format="%1.7f", header=None, index=False
@@ -921,7 +949,7 @@ def SWFN(self, header=True, dataincommentrow=True):
         if dataincommentrow:
             string += self.swcomment
             string += self.krwcomment
-            if "krow" in self.table.columns:
+            if "krow" in self.table.columns and not self.fast:
                 string += "-- krw = krow @ sw=%1.5f\n" % self.crosspoint()
             string += self.pccomment
         string += self.table[["sw", "krw", "pc"]].to_csv(
@@ -943,7 +971,8 @@ def WOTABLE(self, header=True, dataincommentrow=True):
             string += self.swcomment.replace("--", "!")
             string += self.krwcomment.replace("--", "!")
             string += self.krowcomment.replace("--", "!")
-            string += "! krw = krow @ sw=%1.5f\n" % self.crosspoint()
+            if not self.fast:
+                string += "! krw = krow @ sw=%1.5f\n" % self.crosspoint()
             string += self.pccomment.replace("--", "!")
         string += self.table[["sw", "krw", "krow", "pc"]].to_csv(
             sep=" ", float_format="%1.7f", header=None, index=False

pyscal/wateroilgas.py

@@ -32,16 +32,30 @@ class WaterOilGas(object):
         sgcr (float): Critical gas saturation, gas is immobile below this
         h (float): Saturation intervals in generated tables.
         tag (str): Optional text that will be included as comments.
+        fast (bool): Set to True if you prefer speed over robustness. Not recommended,
+            pyscal will not guarantee valid output in this mode.
    """
 
     def __init__(
-        self, swirr=0, swl=0.0, swcr=0.0, sorw=0.00, sorg=0, sgcr=0, h=0.01, tag=""
+        self,
+        swirr=0,
+        swl=0.0,
+        swcr=0.0,
+        sorw=0.00,
+        sorg=0,
+        sgcr=0,
+        h=0.01,
+        tag="",
+        fast=False,
     ):
         """Sets up the saturation range for three phases"""
+        self.fast = fast
         self.wateroil = WaterOil(
-            swirr=swirr, swl=swl, swcr=swcr, sorw=sorw, h=h, tag=tag
+            swirr=swirr, swl=swl, swcr=swcr, sorw=sorw, h=h, tag=tag, fast=fast
+        )
+        self.gasoil = GasOil(
+            swirr=swirr, sgcr=sgcr, sorg=sorg, swl=swl, h=h, tag=tag, fast=fast
         )
-        self.gasoil = GasOil(swirr=swirr, sgcr=sgcr, sorg=sorg, swl=swl, h=h, tag=tag)
 
     def selfcheck(self):
         """Run selfcheck on both wateroil and gasoil.

tests/benchme.py

@@ -0,0 +1,43 @@
+"""Example code for benchmarking of the "fast" feature"""
+import timeit
+
+from pyscal import WaterOilGas
+
+
+def benchme(fast=False, doprint=False):
+    """Test function for benchmarking the "fast"-feature
+
+    Pipe the following code snip into "ipython" on the shell:
+    > echo "import benchme
+      %timeit benchme.benchme()
+      %timeit benchme.benchme(fast=True)
+      " | ipython
+
+    or run this module directly
+    """
+    wog = WaterOilGas(swl=0.1, h=0.1, fast=fast)
+    wog.wateroil.add_corey_oil(now=3, kroend=0.9)
+    wog.wateroil.add_corey_water(nw=2, krwend=0.2)
+    wog.gasoil.add_corey_oil()
+    wog.gasoil.add_corey_gas()
+    if not doprint:
+        len(wog.SWOF())
+        len(wog.SGOF())
+    else:
+        print(wog.SWOF())
+        print(wog.SGOF())
+
+
+if __name__ == "__main__":
+    print("Running in robust and slow mode:")
+    print(
+        timeit.timeit(
+            stmt="benchme(fast=False)", setup="from benchme import benchme", number=100
+        )
+    )
+    print("Running in fast mode:")
+    print(
+        timeit.timeit(
+            stmt="benchme(fast=True)", setup="from benchme import benchme", number=100
+        )
+    )

tests/test_gasoil.py

@@ -191,6 +191,7 @@ def check_endpoints(go, krgend, krgmax, kroend, kromax):
 
 
 def test_gasoil_kromax():
+    """Manual test of kromax behaviour"""
     go = GasOil(h=0.1, sgcr=0.1)
     go.add_corey_oil(2, 0.5)  # Default for kromax
     assert float_df_checker(go.table, "sg", 0.0, "krog", 1.0)
@@ -206,6 +207,23 @@ def test_gasoil_kromax():
     assert float_df_checker(go.table, "sg", 0.0, "krog", 0.5)
 
 
+def test_gasoil_kromax_fast():
+    """Test some code in fast mode"""
+    go = GasOil(h=0.1, sgcr=0.1, fast=True)
+    go.add_corey_oil(2, 0.5)  # Default for kromax
+    assert float_df_checker(go.table, "sg", 0.0, "krog", 1.0)
+    assert float_df_checker(go.table, "sg", 0.1, "krog", 0.5)
+    go.add_corey_oil(2, 0.5, 0.7)
+    assert float_df_checker(go.table, "sg", 0.0, "krog", 0.7)
+    assert float_df_checker(go.table, "sg", 0.1, "krog", 0.5)
+
+    go = GasOil(h=0.1, sgcr=0.0)
+    go.add_corey_oil(2, 0.5)
+    assert float_df_checker(go.table, "sg", 0.0, "krog", 0.5)
+    go.add_corey_oil(2, 0.5, 1)  # A warning will be given
+    assert float_df_checker(go.table, "sg", 0.0, "krog", 0.5)
+
+
 def test_gasoil_krgendanchor():
     """Test behaviour of the krgendanchor"""
     gasoil = GasOil(krgendanchor="sorg", sorg=0.2, h=0.1)