The script subscript/interp_relperm use pyscal.scalrecommendation to interpolate between flow curves. Due to an issue in pyscal probably related to the add_fromtable function, Pc curves are set to zero as opposed to being set to whatever the table contain.

Even though the docstrings state that 'tag' is recognized. The tag should just be forwarded to the WaterOil and GasOil init's.

The crosspoint() function has been seen to crash with NotImplementedError when there are NaN's in the index. That is a bug that must be traced down and fixed, but further, we should probably catch some exceptions from crosspoint() because it only ends in a comment, which can be an error string, and never allow it to crash altogether.

This should probably not error:

> wo.SWFN()                                                                                                                                                                           
ERROR:pyscal.wateroil:krow data not found
> wo.table.columns
Index(['sw', 'swn', 'son', 'swnpc', 'krw', 'pc'], dtype='object')

Allow aliases for low-base-high in interpreted XLS/CSV files.

These should also be supported:
• Pessimistic
• Base
• Optimistic

• Pes
• Base
• Opt

It returns None, but modifies its first argument.

It modifies the second and third argument such that the caller is responsible for making deep copies.

PyscalFactory reuses the term kroend and kromax both for water-oil and for gas-oil.

First: does it make sense to have these different in a three-phase setting
Secondly: If so, we need to use different terms in pyscal factory, maybe krowend, krowmax, krogend, krogmax

The following code illustrates what happens when krgendanchor is ignored under interpolation. The segment from 1-sorg-swl to 1-swl becomes linear, and may break out of the envelope of the low and high curves.

    from pyscal import GasOil
    from pyscal.utils import interpolation
    go_low = GasOil(swl=0, sgcr=0.1, sorg=0.2, krgendanchor="")
    go_low.add_corey_gas(ng=5, krgend=0.5, krgmax=0.7)
    go_low.add_corey_oil(nog=2, kroend=0.6)
    go_high = GasOil(swl=0.02, sgcr=0.05, sorg=0.1, krgendanchor="")
    go_high.add_corey_gas(ng=4, krgend=0.5, krgmax=0.72)
    go_high.add_corey_oil(nog=2, kroend=0.6)
    # Interpolate to midpoint between the curves above
    go_ip = interpolation.interpolate_go(go_low, go_high, 0.5, h=0.1)

    from matplotlib import pyplot as plt
    _, mpl_ax = plt.subplots()
    go_low.plotkrgkrog(mpl_ax=mpl_ax, color="red")
    go_high.plotkrgkrog(mpl_ax=mpl_ax, color="blue")
    go_ip.plotkrgkrog(mpl_ax=mpl_ax, color="green")
    plt.show()

For debugging, it would have been great if all outputted tables from pyscal would include the pyscal version

This will allow mixing Corey and LET in a dataframe for different SATNUMs.

Right now, a Nan NW will give an assertionerror on epsilon < nw < MAX_EXPONENT

Only Python API is demonstrated

For consistency with the ecl2df package and the Eclipse world, it should be considered to make column names in dataframes uppercase-only.

Sample output from WaterOilGas.SOF3:

Header says SGOF, that is critical.

Then the tag is repeated - not critical.

 1 SGOF
  2 -- Garn
  3 -- Garn
  4 -- pyscal: v0.3.1+0.gfb21a5c.dirty
  5 -- swirr=0 swl=0.xx3098 swcr=0xx3098 sorw=0.x9
  6 -- swirr=0, sgcr=0.3, swl=0.xx098, sorg=0.xx, krgendanchor=sorg
  7 -- LET krow, l=2.x7, e=2.41, t=1.x3, kroend=1, kromax=1
  8 -- LET krog, l=2, e=x.65, t=1.x5, kroend=0.8, kromax=1
  9 -- SW     KROW      KROG      
 10 0.0000000 0.0000000 0.0000000
 11 0.0069025 0.0000000 0.0000000
 12 0.0169025 0.0000000 0.0000000

The factory module contains almost only static methods. It should probably not define a class, only provide functions.

When tables have been generated using add_fromtable(), knowledge of sorw and sorg is not present (they are assumed zero). However, these can be inferred from the data by separating data in linear and non-linear segments. Correct interpolation of relperm curves requires knowlegdge of this data.

 @--  ERROR  AT TIME        0.0   DAYS    (18-MAY-1999):
 @           ERROR IN SOF3 TABLE NUMBER   7                                  
 @           AT THE MAXIMUM OIL SATURATION ( 0.5811)                         
 @           KROW AND KROG SHOULD BOTH EQUAL THE OIL RELATIVE                
 @           PERMEABILITY FOR A SYSTEM WITH OIL AND CONNATE                  
 @           WATER ONLY - BUT IN THIS CASE THEY ARE DIFFERENT                
 @            ( KROW = 0.9000 AND KROG = 1.0000)  

From input data:

SATNUM	TAG	SWL	SORW	KRWEND	KROWEND	SORG	KRGEND	SGCR	KROGEND
7	CUSTOM	0.418948342	0.16	0.35	0.9	0.01	0.72	0	1

There should be proper unit handling for any function dealing with data with units.

Find a suitable Python module first for the task.

This gets emitted from add_fromtable()

/home/berland/projects/pyscal/pyscal/wateroil.py:299: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  dframe[nonlinearpart][swcolname].astype(float),
/home/berland/projects/pyscal/pyscal/wateroil.py:300: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  dframe[nonlinearpart][krowcolname].astype(float),

SOF2 and SGWFN are currently not supported. SOF2 is relevant for two-phase oil-water decks and SGWFN is relevant for two-phase gas-water. It is unclear if there is a need for this support.

The following code illustrates a bug with WaterOil.add_fromtable():

from pyscal import WaterOil
import matplotlib.pyplot as plt

(fig, ax) = plt.subplots()

wo = WaterOil(sorw=0.3, h=0.01)
wo.add_corey_water(nw=2, krwend=0.8)
wo.add_corey_oil(now=2)
wo.plotkrwkrow(ax=ax)

wo_fine = WaterOil(h=0.0001)
wo_fine.add_fromtable(wo.table, swcolname="sw")
wo_fine.plotkrwkrow(ax=ax, color="green")
plt.show()

The green curve should follow the blue in this plot.

The reason is that add_fromtable() is using a monotone cubic interpolator, which makes a lot of sense for sw < 1 - sorw, but not above.

The PyscalList container class should have a __getitem__() function so that it can be accessed using the []construct.

The question is then, should it be zero-indexed or 1-indexed?

After having initialized a PyscalList with SATNUM 1 to SATNUM 7, should pyscallist[1] refer to SATNUM 1 or 2?

How should the API for add_corey_oil() be? To simplify usage in interactive sessions, it is tempting to add **kwargs to allow for sending in a dict where only a subset of the args are used, then one can do:

p = dict(swl=.., sorg=.., nw=.., now=.., ng=.., etc)
wog = WaterOilGas(p)
wog.wateroil.add_corey_oil(p)
wog.wateroil.add_corey_water(p)
wog.gasoil.add_corey_gas(p)
wog.gasoil.add_corey_oil(p)

If kwargs is not there to swallow superfluous arguments, we need to be explicit:

wog.wateroil.add_corey_water(nw=p['nw'], krwend=p['krwend'], krwmax=p['krwmax'])

(and we cannot use the default values provided inadd_corey_water function signature)

This is neat, but allowing deliberately superfluous arguments to a function might not be good practice.

The docstrings for SCAL recommendation objects state that it should work for two-phase problems, but the code does not really support that.

Pretty-print the formula.

Show equivalence to other formulas, relate to RMS and to standard petrophysics formula/parameters.

The interpolation code handles saturation endpoints naively.

The current code illustrates how bad it will look in an extreme setting:

from pyscal import *
import matplotlib.pyplot as plt 
import copy
wo0 = WaterOil(swl=0.00)
wo1 = WaterOil(swl=0.5)
wo0.add_corey_water(nw=1)
wo1.add_corey_water(nw=1)
wo0.add_corey_oil(now=3)
wo1.add_corey_oil(now=3)
ip = WaterOil(swl=0.05)
utils.interpolator(ip, copy.deepcopy(wo0), copy.deepcopy(wo1), 0.5)
(fig, ax) = plt.subplots()
wo0.plotkrwkrow(ax=ax, color='blue')
wo1.plotkrwkrow(ax=ax, color='red')
ip.plotkrwkrow(ax=ax, color='green')
plt.show()

The green curves are here the interpolants between red and blue, with an unwanted discontinuity at sw=0.5.

The test codes do not test how kromax and krgmax (+ krgendanchor) behaves through interpolation. No reason yet to believe there are errors, but cannot be trusted before test code is present.

The normalized oil saturation in WaterOil is not using swcr, leading to wrong results, and no effect of kroend/kromax for add_corey_oil. Check also add_LET_oil, has potentially the same issue.

credit: cott

It takes about 0.07-0.10 seconds from initializing a WaterOilGas object to writing out both SWOF and SGOF. In most use cases this is completely fine. 👍

However, in some special use cases, with many relative permeability curves being created, this time quickly adds up. (In my case I'm creating 1000 curves in 100 realization with 4 iterations = 11 hrs compute time).

It would be great to see if there are any (obvious) performance gains to be made.

Given the functionality to infer the endpoints from tabulated data, it is low-hanging fruit to be able to infer (via regression, in scipy) Corey parameters for a tabulated curve.

LET is probably also possible, but wil incur some more numerical issues as the formulation is not as unique.

SGOF
-- SATNUM 1 interpolant at 0.7
-- Sg Krg Krog Pc
-- swirr=0, sgcr=0.01, swl=0.3, sorg=0.01
, krgendanchor=""-- krg = krog @ sw=0.39839
0.0000000 0.0000000 1.0000000 0

could be around here:

hypothesis managed to find a corner case where at least SLGOF output is invalid:

sgcr = 0.04995000000000013
wog = WaterOilGas(swl=0, sorg=0, sgcr=sgcr, h=0.05)
wog.gasoil.add_corey_gas(krgmax=1)
wog.gasoil.add_corey_oil()

First data point in SLGOF should be zero, but it is not. Probably related to 0.00005 and 0.00000 being regarded as equal with the default setting of SWINTEGERS.

This function needs to be added. It should get proper documentation, test-code and a name - which?

def calcPc(SwJ, PERM, PORO, a, b, sigma_costau):
    """
        Calculate capillary pressure based on the water saturation, permeability, porosity, a, b and sigma*cos(tau)

        <-- Input
        SwJ             = Water saturation as a fraction
        PERM            = Permeability in mD (Typical value = formation average)
        PORO            = Porosity as a fraction (Typical value = formation average)
        a & b           = petrophysical J-function fitting parameters (Typical value, see below)
        sigma_costau    = Interfacial tension in mN/m (Typical value = 30 mN/m)
        Result -->l
        Pc              = Capillary pressure in bars

      """
    Pc = (
        (
            ((SwJ / a) ** (1 / b))
            / (math.sqrt((float(0.0000000000009869233) * PERM / float(1000)) / PORO))
        )
        * (sigma_costau / float(1000))
    ) / float(101325)

    return Pc

I noticed that there might be a bug for LET imbibition capillary pressure.

Here is the screenshot of the output from pyscal. As you see, I have set PcMin to -2 bar but the SWOF table generated pressure down to 2 bar.

Propose fix:

    def add_LET_pc_imb(self, Ls, Es, Ts, Lf, Ef, Tf, Pcmax, Pcmin, Pct):
        """Add an imbition LET capillary pressure curve.

        Docs: https://wiki.equinor.com/wiki/index.php/Res:The_LET_correlation_for_capillary_pressure
        """

        # Normalized water saturation including sorw
        self.table["swnpco"] = (self.table.sw - self.swirr) / (
                1 - self.sorw - self.swirr
        )

        # The "forced part"
        # ---- Original code
        # self.table["Fficow"] = (1 - self.table.swnpco) ** Ls / (
        #         (1 - self.table.swnpco) ** Ls + Es * self.table.swnpco ** Ts
        # )
        # ------------------
        self.table["Fficow"] = self.table.swnpco ** Lf / (
                self.table.swnpco ** Lf + Ef * (1 - self.table.swnpco) ** Tf
        )

        # The spontaneous part:
        # Original code:
        # self.table["Fticow"] = self.table.swnpco ** Lf / (
        #         self.table.swnpco ** Lf + Ef * (1 - self.table.swnpco) ** Tf
        # )
        # ------------------
        self.table["Fsicow"] = (1 - self.table.swnpco) ** Ls / (
                (1 - self.table.swnpco) ** Ls + Es * self.table.swnpco ** Ts
        )

        # Putting it together:
        # ---- Original code
        # self.table["pc"] = (
        #         (Pcmax + Pct) * self.table.Fficow - (Pcmin - Pct) * self.table.Fticow - Pct
        # )
        # ------------------
        self.table["pc"] = (
                (Pcmax - Pct) * self.table.Fsicow + (Pcmin - Pct) * self.table.Fficow + Pct
        )

        # Special handling of the interval [0,swirr]
        self.table.loc[self.table.swnpco < epsilon, "pc"] = Pcmax
        # and [1-sorw,1]
        # ---- Original code
        # self.table.loc[self.table.swnpco > 1 - epsilon, "pc"] = self.table.pc.min()
        # -------------------
        self.table.loc[self.table.swnpco > 1 - epsilon, "pc"] = Pcmin
        self.pccomment = (
                "-- LET correlation for imbibition Pc;\n"
                + "-- Ls=%g, Es=%g, Ts=%g, Lf=%g, Ef=%g, Tf=%g, Pcmax=%g, Pcmin=%g, Pct=%g\n"
                % (Ls, Es, Ts, Lf, Ef, Tf, Pcmax, Pcmin, Pct)
        )

The following code crashes pyscal:

data = {'sgcr': 0.273770080949912, 'Lg': 1, 'Eg': 1, 'Tg': 1,'Log': 1, 'Eog': 1, 'Tog': 1.1, 'swl': 0.051} 
from pyscal import PyscalFactory 
factory = PyscalFactory() 
go = factory.create_gas_oil(data) 
go.table.head()                                                        
WARNING:pyscal.gasoil:sorg was close to zero, set to zero
WARNING:pyscal.gasoil:krgmax ignored when not anchoring to sorg
Out[28]: 
        sg       sl       sgn       son       krg      krog
0  0.00000  1.00000 -0.405447  1.405447  0.000000  1.000000
1  0.27377  0.72623  0.000000  1.000000  0.000000       NaN
2  0.28377  0.71623  0.014810  0.985190  0.014810  0.990232
3  0.29377  0.70623  0.029620  0.970380  0.029620  0.978983
4  0.30377  0.69623  0.044429  0.955571  0.044429  0.967067

Capillary pressure is not interpolated from SCALrecommendation objects, zero pc is returned.

This is a regression in v0.3.0. But not present before that.

Avoid filling the log with WARNING:root:message, instead by WARNING:pyscal:message

Capillary pressure is currently interpolated linearly, while it might make more sense to interpolate this logarithmically.

According to E100 doc, SGCR

is the largest gas saturation for which the gas relative permeability is zero

The estimate_sgcr() is probably not fully in sync with this. Probably from practial reasons it estimates based on krog by default assuming there is a linear section in krog to get it from, but it really should be taken from krg-data. But if krg is linear zero, it will return the right endpoint according to docs, which is not according to the definition of sgcr.

It returns objects with empty strings as tag. Improve to something more sensible, conserve the comment if it identical in input curves, or show both if they are different. Don't repeat SATNUM index.

Ensure PyscalList are able to include the satnum index.

There should be a container class that holds a list of WaterOilGas objects.

This is relevant for a user that then can call .SWOF(filename or filehandle) to dump
SWOF data for all SATNUMS to disk, instead of doing the looping themselves.

The delegator pattern can be used in WaterOilGas.

f.ex as in https://www.michaelcho.me/article/method-delegation-in-python

This is probably the square symbol in the pc-comment. Triggered when adding simple J to a pc curve. Why not caught in tests??

Traceback (most recent call last):
  File "create_relperm_smb.py", line 13, in <module>
    swof_fh.write(wog.SWOF(header=False))
  File "/prog/res/komodo/2020.02.rc3-py27/root/lib/python2.7/site-packages/pyscal/wateroilgas.py", line 91, in SWOF
    return self.wateroil.SWOF(header, dataincommentrow)
  File "/prog/res/komodo/2020.02.rc3-py27/root/lib/python2.7/site-packages/pyscal/wateroil.py", line 985, in SWOF
    string += self.pccomment
UnicodeDecodeError: 'ascii' codec can't decode byte 0xc2 in position 103: ordinal not in range(128)

The interpolate() functions in scalrecommendation.py belongs better in their own module than bundled together with the SCAL recommendation object.

A script meant for end-users should be made:

• Input should contain LET or Corey parameters from a dataframe (csv/excel/more?)
• Multiple SATNUMs should be recognized
• Script should make an effort to explain inconsistencies in user input.
• SCAL-recommendations should be supported, if so it must be possible to supply interpolation parameters (through command line and/or parameter files)
• Eclipse keyword family 1 at least. Make room in command line options for family 2.
• Simple J for capillary pressure.

from pyscal import WaterOil
wo = WaterOil(swirr=0, swl=0.025, swcr=0.07, sorw=0.15)
wo.add_corey_oil(now=3, kroend=1, kromax=1)
wo.add_corey_water(nw=4, krwend=0.55, krwmax=1)
print(wo.SWOF())
SWOF
-- 
-- Sw Krw Krow Pc
-- swirr=0 swl=0.025 swcr=0.07 sorw=0.15
-- Corey krw, nw=4, krwend=0.55, krwmax=1
-- Corey krow, now=3, kroend=1, kromax=1
-- krw = krow @ sw=0.53206
-- Zero capillary pressure
0.0250000 0.0000000 1.0000000 0
0.0350000 0.0000000 1.0000000 0
0.0450000 0.0000000 1.0000000 0
0.0550000 0.0000000 1.0000000 0
0.0650000 0.0000000 1.0000000 0
0.0700000 0.0000000 1.0000000 0
0.0750000 0.0000000 0.9808922 0
0.0850000 0.0000001 0.9434100 0
...

where the kro values are not monotonely decreasing, which is required by Eclipse 100. This regression happened in release 0.1.4