solve_1D2P_version1 crashes for some L, E, T values.

[felipeeler]

The solver crashes for some values of the parameters reported as feasible by the authors of the paper "SCA2005-32". The crash is at the .solve() call. It points to line 545 (_solver_result = solve_1D2Pversion)

The authors of LET mode say:

"Experience using the LET correlation indicates that the parameter L ≥ 1, E > 0 and T ≥ 0.5"

Using values like the ones listed below will cause this crash:

L1 = 2.00
E1 = 0.50
T1 = 0.85
L2 = 3.0
E2 = 0.60
T2 = 0.95

I'm pasting some sample code:

from scallib001.displacementmodel1D2P001 import DisplacementModel1D2P
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scallib001.relpermlib001 as rlplib

# Creating the sample data
core_name = 'sample'
core_diameter = 3.81 # cm
core_length = 5.5 # cm
core_perm = 100 # mD
core_porosity = 0.2 # fraction
core_pv = np.pi * np.power(core_diameter,2) * 0.25 * core_length * core_porosity  # cm3

# Creating Fluid Data
oil_mu = 5.0 # cP
oil_rho = 0.85 * 1000.0 # g/cm3 to kg/m3
water_mu = 0.5 # cP
water_rho = 1.05 * 1000.0 # g/cm3 to kg/m3

# Creating Ground Truth
swi_truth = 0.20
sor_truth = 0.15
krw_truth = 0.35
kro_truth = 0.80
p0_truth = 0.50  
beta_truth = 0.70
sd0_truth = 0.45 # Equivalent to 0.4925 without normalization
lw_truth = 2.0
ew_truth = 0.5
tw_truth = 0.85
lo_truth = 3.0
eo_truth = 0.60
to_truth = 0.95

# LET model
rlp_model_truth = rlplib.Rlp2PLET(
    swi_truth,
    sor_truth,
    lw_truth,
    ew_truth,
    tw_truth,
    lo_truth,
    eo_truth,
    to_truth,
    krw_truth,
    kro_truth,
)

# Log s^beta Pc model
class sbeta:

    """
    Log s^beta Pc model, as defined by Roland Lenormand at Cydar User Manual
    """

    def __init__(self, swi, sor, p0, beta, sd0):
        self.swi = swi
        self.sor = sor
        self.p0 = p0
        self.beta = beta
        self.sd0 = sd0
        self.a = 1 / (
            (beta * np.power(0.5, beta - 1)) / (1 - np.power(0.5, beta))
            + (beta * np.power(0.5, beta - 1)) / np.power(0.5, beta)
        )

    def calc(self, swv):
        swi = self.swi
        sor = self.sor
        swvd = (swv - swi) / (1 - swi - sor)

        # =A*Po*(LOG((1-[@Swd]^Beta)/[@Swd]^Beta)-LOG((1-S_Pc_0^Beta)/S_Pc_0^Beta))
        return (
            self.a
            * self.p0
            * (
                np.log((1 - np.power(swvd, self.beta)) / np.power(swvd, self.beta))
                - np.log(
                    (1 - np.power(self.sd0, self.beta)) / np.power(self.sd0, self.beta)
                )
            )
        )

# Pc model
sbeta_model_truth = sbeta(swi_truth, sor_truth, p0_truth, beta_truth, sd0_truth)
EPS = 1e-10
n = 101
swv_truth = np.linspace(sbeta_model_truth.swi + EPS, 1 - sbeta_model_truth.sor - EPS, n)
pcv_truth = sbeta_model_truth.calc(swv_truth)
krw_truth = rlp_model_truth.calc_kr1(swv_truth)
kro_truth = rlp_model_truth.calc_kr2(swv_truth)
cpr_model_truth = rlplib.CubicInterpolator(swv_truth, pcv_truth, lex=1, rex=1)

# Plotting the kr and Pc curves
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(swv_truth, krw_truth, "b-", label=f"$k_r$$_w$")
plt.plot(swv_truth, kro_truth, "r-", label=f"$k_r$$_o$")
plt.xlabel(f"$S_w$")
plt.ylabel(f"$k_r$")
plt.legend()
plt.subplot(1, 2, 2)
plt.plot(swv_truth, pcv_truth * 14.5038, "k-", label=f"$P_c$")
plt.xlabel(f"$S_w$")
plt.ylabel(f"$P_c$ (psi)")
plt.legend()
plt.tight_layout()
plt.show()

# schedules
T_END = 192  # hours
Movie_times = np.linspace(0, T_END, 2304 + 1)  # 1 frame per 5 minutes

Schedule = pd.DataFrame(
    [
        [0.0, 0.25, 1.0],
        [24.001, 0.50, 1.0],
        [48.001, 0.75, 1.0],
        [72.001, 1.00, 1.0],
        [96.001, 2.00, 1.0],
        [144.001, 4.00, 1.0],
    ],
    columns=["StartTime", "InjRate", "FracFlow"],
)

# Define 1D2P displacement model
model_truth = DisplacementModel1D2P(
    NX=50,
    gravity_multiplier=0.0,  # For horizontal displacements, the gravity does not affect dP
    core_length=core_length,
    core_area=np.pi * np.power(core_diameter, 2) * 0.25,
    permeability=core_perm,
    porosity=core_porosity,
    sw_initial=swi_truth,
    viscosity_w=water_mu,
    viscosity_n=oil_mu,
    density_w=water_rho,
    density_n=oil_rho,
    rlp_model=rlp_model_truth,
    cpr_model=cpr_model_truth,
    time_end=T_END,
    rate_schedule=Schedule,
    movie_schedule=Movie_times,
)

# Solve the model with the initial gues relperm and capcurve models
# note that first time code needs to be called by numba, takes some time
model_truth.solve()

# Keep the results
result_truth = model_truth.results

# Plotting the results
plt.figure(figsize=(15, 4))
tss_truth = result_truth.tss_table
oil_balance_truth = tss_truth.CumOIL - tss_truth.CumOILInj
plt.subplot(1, 2, 1)
plt.plot(tss_truth.TIME, oil_balance_truth, "k", label="Ground Truth")
plt.legend()
plt.ylabel("Oil Production[cm3]")
plt.xlabel("Time [hour]")
plt.subplot(1, 2, 2)
plt.plot(tss_truth.TIME, tss_truth.delta_P * 14.5038, "k", label="Ground Truth")
plt.legend()
plt.ylabel("Pressure drop [psi]")
plt.xlabel("Time [hour]")
plt.show()

The main difference is that I'm using another Pc function, but I think it does not influence the behavior I'm reporting. This code works fine with L, E, T values ≥ 1.

[steffenbergshell]

The LET relperm can have infinite saturation derivatives at endpoints (Sw = Swc or Sw = 1 - Sorw) when one of its Corey exponents (L or T) is smaller than 1. Infinite valued derivatives cause the Newton solver to fail. We discussed the issue with the LET model authors. They acknowledged the problem and suggested methods to regularise the derivative. We are currently developing a more robust implementation of the LET relperm and expect to provide an update soon.

[felipeeler]

Thank you very much for the reply!

[harmdijk]

We have updated the implementation of the LET relperm such that the function and its derivative are always finite. For the implementation we followed the suggestion of the LET authors.

In case you are interested in implementation details, the key update is the introduction of the 'power_eps' function in the relperm library. When the exponent (L or T in the LET relperm model) is less than 1, power_eps returns a regularized version of the power function by using a small parameter 'eps',

power_eps(s, N) = ((s + eps)**N - eps**N)/((1+eps)**N - eps**N) for N < 1.

When the exponent (L or T) is greater or equal 1, power_eps is equivalent to the standard power function,

power_eps(s, N) = s**N for N >= 1.

The eps parameter can be specified as the last input parameter for the Rlp2PLET and Rlp2PCorey models; eps = 0.0001 by default.

In addition, to further increase robustness, we introduced checks on input parameters for Corey and LET relperm models. Parameters will be accepted only, if they are within a physically acceptable range (0 <= Sr1 < 1, 0 < N1, Sr1 + Sr2 < 1, etc.).

With the updated implementation, your example problem now runs without any problem. Thanks for sharing the issue!

[felipeeler]

I really appreciate it! Thank you!