In [1]:

    

from fipy import *

# relperm parameters
swc = 0.1
sor = 0.1
krw0 = 0.3
kro0 = 1.0
nw = 2.0
no = 2.0

# domain and boundaries
k = 1e-12 # m^2
phi = 0.4
u = 1.e-5
p0 = 100e5 # Pa
Lx = 100.
Ly = 10.
nx = 100
ny = 10
dx = Lx/nx
dy = Ly/ny

# fluid properties
muo = 0.002
muw = 0.001

# define the fractional flow functions
def krw(sw):
    res = krw0*((sw-swc)/(1-swc-sor))**nw
    return res

def dkrw(sw):
    res = krw0*nw/(1-swc-sor)*((sw-swc)/(1-swc-sor))**(nw-1)
    return res


def kro(sw):
    res = kro0*((1-sw-sor)/(1-swc-sor))**no
    return res

def dkro(sw):
    res = -kro0*no/(1-swc-sor)*((1-sw-sor)/(1-swc-sor))**(no-1)
    return res

def fw(sw):
    res = krw(sw)/muw/(krw(sw)/muw+kro(sw)/muo)
    return res

def dfw(sw):
    res = (dkrw(sw)/muw*kro(sw)/muo-krw(sw)/muw*dkro(sw)/muo)/(krw(sw)/muw+kro(sw)/muo)**2
    return res

import matplotlib.pyplot as plt
import numpy as np

sw_plot = np.linspace(swc, 1-sor, 50)


    

In [2]:

    

krw_plot = [krw(sw) for sw in sw_plot]
kro_plot = [kro(sw) for sw in sw_plot]
fw_plot = [fw(sw) for sw in sw_plot]

plt.figure(1)
plt.plot(sw_plot, krw_plot, sw_plot, kro_plot)
plt.show()

plt.figure(2)
plt.plot(sw_plot, fw_plot)
plt.show()


    

In [3]:

    

# create the grid
mesh = Grid1D(dx = Lx/nx, nx = nx)
x = mesh.cellCenters

# create the cell variables and boundary conditions
sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = swc)
p = CellVariable(mesh=mesh, name="pressure", hasOld=True, value = p0)
# sw.setValue(1,where = x<=dx)
sw.constrain(1.0,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)
p.constrain(p0, mesh.facesRight)
p.constrain(1.02*p0, mesh.facesLeft)


    

In [4]:

    

eq_p = DiffusionTerm(var=p, coeff=-k*(krw(sw.faceValue)/muw+kro(sw.faceValue)/muo))- \
UpwindConvectionTerm(var=sw, coeff=-k*(dkrw(sw.faceValue)/muw+dkro(sw.faceValue)/muo)*p.faceGrad)- \
(k*(dkrw(sw.faceValue)/muw+dkro(sw.faceValue)/muo)*sw.faceValue*p.faceGrad).divergence == 0

eq_sw = TransientTerm(coeff=phi, var=sw) + \
DiffusionTerm(var=p, coeff=-k*krw(sw.faceValue)/muw)+ \
UpwindConvectionTerm(var=sw, coeff=-k*dkrw(sw.faceValue)/muw*p.faceGrad)- \
(-k*dkrw(sw.faceValue)/muw*p.faceGrad*sw.faceValue).divergence == 0

eq = eq_p & eq_sw
steps = 200
dt0 = 500.
dt = dt0
t_end = steps*dt0
t = 0.0
viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)
while t<t_end:
    swres = 1.0e6
    loop_count = 0
    while True:
        swres_new = eq.sweep(dt = dt)
        sw.value[sw.value>1-sor]=1-sor
        sw.value[sw.value<swc]=swc
        loop_count+=1
        if loop_count==1:
            sw_res = swres_new
        if swres_new>sw_res or loop_count>5:
            dt = dt/3
            sw.value = sw.old[:]
            p.value = p.old[:]
            continue
        swres=swres_new
        print(swres)
        if swres_new<1e-5:
            sw.updateOld()
            p.updateOld()
            t+=dt
            dt = dt0
            break
        
# Note: try to use the Appleyard method; the overflow is a result of wrong rel-perm values        
viewer.plot()


    

    




0.0005369273839403227
0.0005369273839403227






    




/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/cellTerm.py:118: RuntimeWarning: overflow encountered in true_divide
  L.addAt(coeffVectors['new value'].ravel() / dt, ids.ravel(), ids.swapaxes(0, 1).ravel())
/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/cellTerm.py:116: RuntimeWarning: overflow encountered in true_divide
  b += (oldArray.value[numerix.newaxis] * coeffVectors['old value']).sum(-2).ravel() / dt
/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/cellTerm.py:116: RuntimeWarning: divide by zero encountered in true_divide
  b += (oldArray.value[numerix.newaxis] * coeffVectors['old value']).sum(-2).ravel() / dt
/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/cellTerm.py:118: RuntimeWarning: divide by zero encountered in true_divide
  L.addAt(coeffVectors['new value'].ravel() / dt, ids.ravel(), ids.swapaxes(0, 1).ravel())






    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-4-986a9bfede1b> in <module>
     19     loop_count = 0
     20     while True:
---> 21         swres_new = eq.sweep(dt = dt)
     22         sw.value[sw.value>1-sor]=1-sor
     23         sw.value[sw.value<swc]=swc

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/term.py in sweep(self, var, solver, boundaryConditions, dt, underRelaxation, residualFn, cacheResidual, cacheError)
    212             `Term`
    213         """
--> 214         solver = self._prepareLinearSystem(var=var, solver=solver, boundaryConditions=boundaryConditions, dt=dt)
    215         solver._applyUnderRelaxation(underRelaxation=underRelaxation)
    216         residual = solver._calcResidual(residualFn=residualFn)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/term.py in _prepareLinearSystem(self, var, solver, boundaryConditions, dt)
    128                                                            transientGeomCoeff=self._getTransientGeomCoeff(var),
    129                                                            diffusionGeomCoeff=self._getDiffusionGeomCoeff(var),
--> 130                                                            buildExplicitIfOther=self._buildExplcitIfOther)
    131 
    132         self._buildCache(matrix, RHSvector)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/coupledBinaryTerm.py in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     86                                                                                      transientGeomCoeff=uncoupledTerm._getTransientGeomCoeff(tmpVar),
     87                                                                                      diffusionGeomCoeff=uncoupledTerm._getDiffusionGeomCoeff(tmpVar),
---> 88                                                                                      buildExplicitIfOther=buildExplicitIfOther)
     89 
     90                 termMatrix += tmpMatrix

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/binaryTerm.py in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     32                                                                         transientGeomCoeff=transientGeomCoeff,
     33                                                                         diffusionGeomCoeff=diffusionGeomCoeff,
---> 34                                                                         buildExplicitIfOther=buildExplicitIfOther)
     35 
     36             matrix += tmpMatrix

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/binaryTerm.py in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     32                                                                         transientGeomCoeff=transientGeomCoeff,
     33                                                                         diffusionGeomCoeff=diffusionGeomCoeff,
---> 34                                                                         buildExplicitIfOther=buildExplicitIfOther)
     35 
     36             matrix += tmpMatrix

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/unaryTerm.py in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     65                                                        dt=dt,
     66                                                        transientGeomCoeff=transientGeomCoeff,
---> 67                                                        diffusionGeomCoeff=diffusionGeomCoeff)
     68         elif buildExplicitIfOther:
     69             _, matrix, RHSvector = self._buildMatrix(self.var,

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/abstractDiffusionTerm.py in _buildMatrix(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff)
    321             ids = self._reshapeIDs(var, numerix.arange(mesh.numberOfCells))
    322             L.addAt(self.constraintL.ravel(), ids.ravel(), ids.swapaxes(0, 1).ravel())
--> 323             b += numerix.reshape(self.constraintB.ravel(), ids.shape).sum(-2).ravel()
    324 
    325         return (var, L, b)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in ravel(self)
   1394 
   1395     def ravel(self):
-> 1396         return self.value.ravel()
   1397 
   1398     def _axisClass(self, axis):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _getValue(self)
    491 
    492         if self.stale or not self._isCached() or self._value is None:
--> 493             value = self._calcValue()
    494             if self._isCached():
    495                 self._setValueInternal(value=value)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/operatorVariable.py in _calcValue(self)
     53                     return self._execInline(comment=self.comment)
     54                 else:
---> 55                     return self._calcValue_()
     56 
     57         def _calcValue_(self):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/binaryOperatorVariable.py in _calcValue_(self)
     46                 val1 = self.var[1]
     47 
---> 48             return self.op(self.var[0].value, val1)
     49 
     50         @property

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _getValue(self)
    491 
    492         if self.stale or not self._isCached() or self._value is None:
--> 493             value = self._calcValue()
    494             if self._isCached():
    495                 self._setValueInternal(value=value)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/operatorVariable.py in _calcValue(self)
     53                     return self._execInline(comment=self.comment)
     54                 else:
---> 55                     return self._calcValue_()
     56 
     57         def _calcValue_(self):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/binaryOperatorVariable.py in _calcValue_(self)
     46                 val1 = self.var[1]
     47 
---> 48             return self.op(self.var[0].value, val1)
     49 
     50         @property

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _getValue(self)
    491 
    492         if self.stale or not self._isCached() or self._value is None:
--> 493             value = self._calcValue()
    494             if self._isCached():
    495                 self._setValueInternal(value=value)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/operatorVariable.py in _calcValue(self)
     53                     return self._execInline(comment=self.comment)
     54                 else:
---> 55                     return self._calcValue_()
     56 
     57         def _calcValue_(self):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/unaryOperatorVariable.py in _calcValue_(self)
     34     class unOp(operatorClass):
     35         def _calcValue_(self):
---> 36             return self.op(self.var[0].value)
     37 
     38         @property

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _getValue(self)
    491 
    492         if self.stale or not self._isCached() or self._value is None:
--> 493             value = self._calcValue()
    494             if self._isCached():
    495                 self._setValueInternal(value=value)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/addOverFacesVariable.py in _calcValue(self)
     29             return self._calcValueInline()
     30         else:
---> 31             return self._calcValueNoInline()
     32 
     33     def _calcValueInline(self):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/addOverFacesVariable.py in _calcValueNoInline(self)
     70         ids = self.mesh.cellFaceIDs
     71 
---> 72         contributions = numerix.take(self.faceVariable, ids, axis=-1)
     73 
     74         # FIXME: numerix.MA.filled casts away dimensions

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/tools/numerix.py in take(a, indices, axis, fill_value)
    600 
    601     if _isPhysical(a):
--> 602         taken = a.take(indices, axis=axis)
    603     elif isinstance(indices, type(MA.array((0)))):
    604         ## Replaces `MA.take`. `MA.take` does not always work when

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in take(self, ids, axis)
   1463 
   1464     def take(self, ids, axis=0):
-> 1465         return numerix.take(self.value, ids, axis)
   1466 
   1467     def allclose(self, other, rtol=1.e-5, atol=1.e-8):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _getValue(self)
    495                 self._setValueInternal(value=value)
    496             else:
--> 497                 self._setValueInternal(value=None)
    498             self._markFresh()
    499         else:

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _setValueInternal(self, value, unit, array)
    622 
    623     def _setValueInternal(self, value, unit=None, array=None):
--> 624         self._value = self._makeValue(value=value, unit=unit, array=array)
    625 
    626     def _makeValue(self, value, unit=None, array=None):

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/variables/variable.py in _makeValue(self, value, unit, array)
    659                 array[:] = value
    660                 value = array
--> 661             elif type(value) not in (type(None), type(numerix.array(1)), type(numerix.MA.array(1))):
    662                 value = numerix.array(value)
    663 ##                 # numerix does strange things with really large integers.

KeyboardInterrupt: 





    

In [9]:

    

import fractional_flow as ff
xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil = ff.frac_flow_wf(muw=muw, muo=muo, ut=u, phi=1.0, \
  k=1e-12, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0, \
  nw=nw, sw0=swc, sw_inj=1.0, L=Lx, pv_inj=5.0)


    

In [10]:

    

plt.figure()
plt.plot(xt_prf, sw_prf)
plt.plot(x.value.squeeze()/(steps*dt), sw.value)
plt.show()


    

    




/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in true_divide
  This is separate from the ipykernel package so we can avoid doing imports until