Skip to content

Commit

Permalink
Support SWU and SWLPC (#345) and allow cell-centre discrepancy
Browse files Browse the repository at this point in the history 
* Add support for SWU and SWLPC
* Add contact tolerance, allowing cell-centres to be slightly above or below contacts.
  • Loading branch information
@berland
berland authored Apr 27, 2021
1 parent de4594a commit 589f045
Show file tree
Hide file tree
Showing 4 changed files with 331 additions and 44 deletions.
81 changes: 56 additions & 25 deletions src/subscript/check_swatinit/check_swatinit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,10 +123,6 @@ def check_applicability(eclfiles):
"RPTRST not found in DATA-file, UNRST file is expected to be missing"
)

if "SWLPC" in deck:
logger.warning("SWLPC found in DATA file")
logger.warning("Not supported by check_swatinit. Do not trust output")

try:
eclfiles.get_rstfile()
except FileNotFoundError as exception:
Expand Down Expand Up @@ -240,8 +236,8 @@ def make_qc_gridframe(eclfiles):
"PCW",
"PPCW",
"SWL",
"SWU", # It is a feature request to process this
"SWLPC", # Extract in case it is there, but it is not supported yet.
"SWLPC",
"SWU",
],
rstdates="first",
)
Expand All @@ -257,12 +253,6 @@ def make_qc_gridframe(eclfiles):
logger.warning("Consider adding FILLEPS to the PROPS section")
grid_df["SWL"] = 0.0

if "SWU" in grid_df and (grid_df["SWU"] - 1).abs().sum() > 0:
logger.warning("SWU is less than 1, not yet supported by check_swatinit")

if "SWLPC" in grid_df and (grid_df["SWL"] - grid_df["SWLPC"]).abs().sum() > 0:
raise ValueError("SWLPC is in the data, but not supported by check_swatinit")

deck = eclfiles.get_ecldeck()
if "SWATINIT" in deck:
swatinit_deckdata = deck["SWATINIT"][0][0].get_raw_data_list()
Expand Down Expand Up @@ -337,9 +327,14 @@ def qc_flag(qc_frame):
else:
contact = "GWC"

# Eclipse and libecl does not calculate cell centres to the same decimals.
# Add some tolerance when testing towards fluid contacts.
contacttolerance = 1e-4

# SWATINIT ignored, water is lost if pc > 0, lost and/or gained if oil-wet pc-curve
qc_col[
(qc_frame["SWATINIT"] < 1) & (qc_frame["Z"] > qc_frame[contact])
(qc_frame["SWATINIT"] < 1)
& (qc_frame["Z"] > qc_frame[contact] - contacttolerance)
] = __HC_BELOW_FWL__

# SWATINIT accepted and PC is scaled.
Expand Down Expand Up @@ -375,14 +370,25 @@ def qc_flag(qc_frame):

# SWATINIT=1 above contact:
qc_col[
np.isclose(qc_frame["SWATINIT"], 1) & (qc_frame["Z"] < qc_frame[contact])
np.isclose(qc_frame["SWATINIT"], 1)
& (qc_frame["Z"] < qc_frame[contact] + contacttolerance)
] = __SWATINIT_1__

# If SWU is less than 1, SWATINIT is ignored whenever it is equal or larger
# than SWU. Behaviour is the same as SWATINIT=1; SWATINIT is ignored, and thus
# the same flag is reused:
if "SWU" in qc_frame:
qc_col[
(qc_frame["SWU"] < 1)
& ~np.isclose(qc_frame["SWATINIT"], qc_frame["SWAT"])
& (qc_frame["SWATINIT"] >= qc_frame["SWU"])
] = __SWATINIT_1__

# SWATINIT=1 below contact but with SWAT < 1, can happen with OIP_INIT:
if "OIP_INIT" in qc_frame:
qc_col[
(~np.isclose(qc_frame["OIP_INIT"], 0))
& (qc_frame["Z"] > qc_frame[contact])
& (qc_frame["Z"] > qc_frame[contact] - contacttolerance)
& (np.isclose(qc_frame["SWATINIT"], 1))
& (~np.isclose(qc_frame["SWATINIT"], qc_frame["SWAT"]))
] = __SWATINIT_1__
Expand All @@ -401,7 +407,18 @@ def qc_flag(qc_frame):
& (qc_frame["Z"] > qc_frame[contact])
] = __WATER__

qc_col[qc_frame["SWL"] > qc_frame["SWATINIT"]] = __SWL_TRUNC__
qc_col[
np.isclose(qc_frame["SWAT"], qc_frame["SWL"])
& (qc_frame["SWL"] > qc_frame["SWATINIT"])
] = __SWL_TRUNC__

if "SWLPC" in qc_frame:
# SWLPC is not supported by OPM-flow, therefore we also check
# that SWAT == SWLPC before assigning this:
qc_col[
np.isclose(qc_frame["SWAT"], qc_frame["SWLPC"])
& (qc_frame["SWLPC"] > qc_frame["SWATINIT"])
] = __SWL_TRUNC__

# Tag the remainder with "unknown", when/if this happens, it is a bug or a
# feature request:
Expand Down Expand Up @@ -450,13 +467,15 @@ def qc_volumes(qc_frame):
return watergains


def _evaluate_pc(swats, scale_vert, swls, satfunc, sat_name="SW", pc_name="PCOW"):
def _evaluate_pc(swats, scale_vert, swls, swus, satfunc, sat_name="SW", pc_name="PCOW"):
"""Evaluate pc as a function of saturation on a scaled Pc-curve
Args:
swats (list): floats with water saturation values
scale_vert (list): floats with vertical scalers for pc
swls (list): List of SWL values for horizontal scaling
swls (list): List of SWL values for horizontal scaling.
If the model is using SWLPC, supply those values instead.
swus (list): List of SWU values for horizontal scaling.
satfunc (pd.DataFrame): Dataframe representing un-scaled
capillary pressure curve
sat_name (str): Column name for the column in the dataframe with the
Expand All @@ -465,21 +484,21 @@ def _evaluate_pc(swats, scale_vert, swls, satfunc, sat_name="SW", pc_name="PCOW"
values.
Returns:
pd.Series: computed capillary pressure values.
list: computed capillary pressure values.
"""
p_cap = []
sw_min = satfunc[sat_name].min()
sw_max = satfunc[sat_name].max()
if swls is None:
swls = [sw_min] * len(swats)
for swat, pc_scaling, swl in zip(swats, scale_vert, swls):
if swus is None:
swus = [sw_max] * len(swats)
for swat, pc_scaling, swl, swu in zip(swats, scale_vert, swls, swus):
p_cap.append(
np.interp(
swat,
swl
+ (satfunc[sat_name].values - sw_min)
/ (sw_max - sw_min)
* (sw_max - swl),
+ (satfunc[sat_name].values - sw_min) / (sw_max - sw_min) * (swu - swl),
satfunc[pc_name].values * pc_scaling,
)
)
Expand All @@ -496,13 +515,18 @@ def compute_pc(qc_frame, satfunc_df):
Likewise, we are not getting arbitrarily negative Pc values below contact
for the same reason, only slightly negative for an oil-wet curve.
Note: If SWLPC is present in the dataframe, it will be used instead
of SWL. As OPM-flow outputs SWLPC in the INIT files but otherwise
ignores SWLPC, this will result in an incorrect PC estimate for OPM-flow
when SWLPC is in use.
Args:
qc_frame (pd.DataFrame)
satfunc_df (pd.DataFrame)
Returns:
pd.Series, with capillary pressure values in bars (given Eclipse unit
is METRIC)
is METRIC)
"""
p_cap = pd.Series(index=qc_frame.index, dtype=np.float64)
p_cap[:] = np.nan
Expand All @@ -511,14 +535,21 @@ def compute_pc(qc_frame, satfunc_df):
return p_cap

for satnum, satnum_frame in qc_frame.groupby("SATNUM"):
if "SWL" in satnum_frame:
if "SWLPC" in satnum_frame:
swls = satnum_frame["SWLPC"].values
elif "SWL" in satnum_frame:
swls = satnum_frame["SWL"].values
else:
swls = None
if "SWU" in satnum_frame:
swus = satnum_frame["SWU"].values
else:
swus = None
p_cap[satnum_frame.index] = _evaluate_pc(
satnum_frame["SWAT"].values,
satnum_frame["PC_SCALING"].values,
swls,
swus,
satfunc_df[satfunc_df["SATNUM"] == satnum],
)
# Fix needed for OPM-flow above contact:
Expand Down
26 changes: 23 additions & 3 deletions src/subscript/check_swatinit/pillarmodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ def __init__(
swatinit=None,
satnum=None, # one value pr. cell
swl=None, # first saturation in swof, one pr. satnum
swlpc=None, # swl scaler for pc only, leaving kr unmodified
swu=None, # maximum saturation in swof, one pr. satnum
maxpc=None, # max pc in SWOF
minpc=None, # pc at sw=1 in SWOF
ppcwmax=None, # PPCWMAX keyword.
Expand Down Expand Up @@ -68,6 +70,20 @@ def __init__(
# SWL is pr. SATNUM, not pr. cell in this class
assert len(self.swl) == max(self.satnum)

if swlpc is None:
self.swlpc = None
else:
self.swlpc = swlpc
# SWLPC is pr. cell
assert len(self.swlpc) == self.cells

if swu is None:
self.swu = [1] * max(self.satnum)
else:
self.swu = swu
# SWU is pr. SATNUM, not pr. cell in this class
assert len(self.swu) == max(self.satnum)

if maxpc is None:
self.maxpc = [3] * max(self.satnum)
else:
Expand Down Expand Up @@ -220,7 +236,7 @@ def evaluate_sw(self, p_cap, scaling=1, satnum=1):
return np.interp(
p_cap,
[self.minpc[satnum - 1] * scaling, self.maxpc[satnum - 1] * scaling],
[1, self.swl[satnum - 1]],
[self.swu[satnum - 1], self.swl[satnum - 1]],
)

def props(self):
Expand All @@ -242,7 +258,7 @@ def props(self):
string += "SWOF\n"
string += "-- SW KRW KROW PC\n"
string += f"{self.swl[satnum]:g} 0 1 {self.maxpc[satnum]:g}\n"
string += f"1.0 1.0 0.0 {self.minpc[satnum]:g}\n/\n"
string += f"{self.swu[satnum]:g} 1.0 0.0 {self.minpc[satnum]:g}\n/\n"

if "GAS" in self.phases and "OIL" in self.phases:
for satnum in range(len(self.swl)):
Expand All @@ -256,14 +272,18 @@ def props(self):
string += "SWFN\n"
string += "-- SW KRW PC\n"
string += f"{self.swl[satnum]:g} 0 {self.maxpc[satnum]:g}\n"
string += f"1.0 1.0 {self.minpc[satnum]:g}\n/\n"
string += f"{self.swu[satnum]:g} 1.0 {self.minpc[satnum]:g}\n/\n"

for satnum in range(len(self.swl)):
string += "SGFN\n"
string += "-- SG KRG PC\n"
string += "0 0 0\n"
string += f"{1 - self.swl[satnum]:g} 1.0 0.0\n/\n"

if self.swlpc is not None:
string += "SWLPC\n"
string += _wrap(" ".join(map(str, self.swlpc)) + "/") + "\n"

string += """\n
DENSITY
800 1000 1.2 /
Expand Down
101 changes: 85 additions & 16 deletions tests/test_check_swatinit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,16 @@
@pytest.mark.parametrize(
"propslist, expected_flag",
[
([{"SWL": 0.3, "SWATINIT": 0.1}], __SWL_TRUNC__),
([{"SWL": 0.3, "SWATINIT": 0.1}], __UNKNOWN__),
([{"SWL": 0.3, "SWAT": 0.3, "SWATINIT": 0.1}], __SWL_TRUNC__),
(
[{"SWATINIT": 0.9, "SWAT": 1.0, "Z": 1000, "OWC": 900}],
__HC_BELOW_FWL__,
),
(
[{"SWATINIT": 0.9, "SWAT": 1.0, "Z": 899.99991, "OWC": 900}],
__HC_BELOW_FWL__,
),
([{"SWAT": 0.3, "SWATINIT": 0.3}], __PC_SCALED__),
(
[{"SWATINIT": 0.9, "SWAT": 0.8, "PPCW": 100, "PPCWMAX": 100}],
Expand All @@ -67,6 +72,10 @@
[{"SWATINIT": 1, "Z": 100, "OWC": 900}],
__SWATINIT_1__,
),
(
[{"SWATINIT": 1, "Z": 900.00001, "OWC": 900}],
__SWATINIT_1__,
),
(
[{"SWATINIT": 1, "SWAT": 0.9, "Z": 100, "OWC": 900}],
# In this case, E100 has ignored SWATINIT and found
Expand Down Expand Up @@ -224,6 +233,58 @@
],
__HC_BELOW_FWL__,
),
(
[
{
"SWATINIT": 0.9,
"SWAT": 0.9,
"SWU": 0.95,
"Z": 1100,
"OWC": 1200,
}
],
__PC_SCALED__,
),
(
[
{
"SWATINIT": 0.9,
"SWAT": 0.5,
"SWU": 0.9,
"Z": 1100,
"OWC": 1200,
}
],
__SWATINIT_1__,
),
(
[
{
"SWATINIT": 0.9,
"SWAT": 0.5,
"SWU": 0.7,
"Z": 1100,
"OWC": 1200,
}
],
__SWATINIT_1__,
),
(
[
{
"SWATINIT": 0.9,
"SWAT": 0.7,
"SWU": 0.7,
# What if below contact? Not sure what makes the most
# sense and if this is attainable from simulators:
"Z": 1300,
"OWC": 1200,
}
],
__SWATINIT_1__,
),
([{"SWL": 0.3, "SWAT": 0.1, "SWLPC": 0.05, "SWATINIT": 0.1}], __PC_SCALED__),
([{"SWL": 0.3, "SWAT": 0.4, "SWLPC": 0.4, "SWATINIT": 0.1}], __SWL_TRUNC__),
],
)
def test_qc_flag(propslist, expected_flag):
Expand Down Expand Up @@ -288,28 +349,36 @@ def test_qc_volumes(propslist, expected_dict):


@pytest.mark.parametrize(
"swats, scale_vert, swls, expected_pc",
"swats, scale_vert, swls, swus, expected_pc",
[
([0], [1], None, [3]),
([0.1], [1], None, [3]),
([1], [1], None, [0]),
([0, 1], [1, 1], None, [3, 0]),
([2], [1], None, [0]), # constant extrapolation
([0.55], [1], None, [1.5]),
([0.55], [2], None, [3]),
([0.55], [2], [0], [3 - 3 / 10.0]),
([0], [1], [0], [3]),
([0.1], [1], [0], [3 - 3 / 10.0]),
([1], [1], [0], [0]),
([0.5], [1], [0.5], [3]),
# First without any SWL/SWU scaling:
([0], [1], None, None, [3]),
([0.1], [1], None, None, [3]),
([1], [1], None, None, [0]),
([0, 1], [1, 1], None, None, [3, 0]),
([2], [1], None, None, [0]), # constant extrapolation
([0.55], [1], None, None, [1.5]),
([0.55], [2], None, None, [3]),
# Test SWL scaling:
([0.55], [2], [0], None, [3 - 3 / 10.0]),
([0], [1], [0], None, [3]),
([0.1], [1], [0], None, [3 - 3 / 10.0]),
([1], [1], [0], None, [0]),
([0.5], [1], [0.5], None, [3]),
# Test SWU scaling:
([0.5], [1], None, [0.5], [0]),
([0.5], [10], None, [0.5], [0]),
([0.5], [1], [0.3], [0.5], [0]),
([0.1], [1], None, [0.5], [3]),
([0.1], [1], [0], [0.9], [3 - 1 / 3]),
],
)
def test_evaluate_pc(swats, scale_vert, swls, expected_pc):
def test_evaluate_pc(swats, scale_vert, swls, swus, expected_pc):
"""Test that we can evaluate capillary pressure from a saturation table
when we also allow the pc in the table to be scaled as Eclipse is doing it"""
satfunc_df = pd.DataFrame([{"SW": 0.1, "PCOW": 3}, {"SW": 1, "PCOW": 0}])
assert np.isclose(
_evaluate_pc(swats, scale_vert, swls, satfunc_df), expected_pc
_evaluate_pc(swats, scale_vert, swls, swus, satfunc_df), expected_pc
).all()


Expand Down
Loading

0 comments on commit 589f045

Please sign in to comment.