There will be a webinar on pyGIMLi hosted by SEG on March 19, 2024 at 4 pm CET. Register for free here.

Source code for pygimli.physics.petro.resistivity

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""For the manager look at BERT https://gitlab.com/resistivity-net/bert."""

import numpy as np
import matplotlib.pyplot as plt

import pygimli as pg
import pygimli.meshtools as mt

[docs] def resistivityArchie(rFluid, porosity, a=1.0, m=2.0, sat=1.0, n=2.0, mesh=None, meshI=None, fill=None, show=False): r"""Resistivity of rock for the petrophysical model from Archies law. Calculates resistivity of rock for the petrophysical model from Archie's law. :cite:`Archie1942` .. math:: \rho = a\rho_{\text{fl}}\phi^{-m} S^{-n} * :math:`\rho` - the electrical resistivity of the fluid saturated rock in :math:`\Omega\text{m}` * :math:`\rho_{\text{fl}}` - rFluid: electrical resistivity of the fluid in :math:`\Omega\text{m}` * :math:`\phi` - porosity 0.0 --1.0 * :math:`S` - fluid saturation 0.0 --1.0 [sat] * :math:`a` - Tortuosity factor. (common 1) * :math:`m` - Cementation exponent of the rock (usually in the range 1.3 -- 2.5 for sandstones) * :math:`n` - is the saturation exponent (usually close to 2) If mesh is not None the resulting values are calculated for each cell of the mesh. All parameter can be scalar, array of length mesh.cellCount() or callable(pg.cell). If rFluid is non-steady n-step distribution than rFluid can be a matrix of size(n, mesh.cellCount()) If meshI is not None the result is interpolated to meshI.cellCenters() and prolonged (if fill ==1). Notes ----- We experience some unstable nonlinear behavior. Until this is clarified all results are rounded to the precision 1e-6. Examples -------- >>> # WRITEME """ phi = porosity if isinstance(porosity, list): phi = np.array(porosity) if mesh is None: return rFluid * a * phi**(-m) * sat**(-n) rB = None if isinstance(rFluid, float): rB = pg.Matrix(1, mesh.cellCount()) rB[0] = pg.solver.parseArgToArray(rFluid, mesh.cellCount(), mesh) elif isinstance(rFluid, pg.Vector): rB = pg.Matrix(1, len(rFluid)) rB[0] = pg.solver.parseArgToArray(rFluid, mesh.cellCount(), mesh) elif hasattr(rFluid, 'ndim') and rFluid.ndim == 1: rB = pg.Matrix(1, len(rFluid)) rB[0] = pg.solver.parseArgToArray(rFluid, mesh.cellCount(), mesh) elif hasattr(rFluid, 'ndim') and rFluid.ndim == 2: rB = pg.Matrix(len(rFluid), len(rFluid[0])) for i, rFi in enumerate(rFluid): rB[i] = rFi phi = pg.solver.parseArgToArray(phi, mesh.cellCount(), mesh) a = pg.solver.parseArgToArray(a, mesh.cellCount(), mesh) m = pg.solver.parseArgToArray(m, mesh.cellCount(), mesh) S = pg.solver.parseArgToArray(sat, mesh.cellCount(), mesh) n = pg.solver.parseArgToArray(n, mesh.cellCount(), mesh) if show: pg.show(mesh, S, label='S') pg.show(mesh, phi, label='p') pg.wait() r = pg.Matrix(len(rB), len(rB[0])) for i, _ in enumerate(r): r[i] = rB[i] * a * phi**(-m) * S**(-n) r.round(1e-6) if meshI is None: if len(r) == 1: return r[0].copy() return r rI = pg.Matrix(len(r), meshI.cellCount()) if meshI: pg.interpolate(mesh, r, meshI.cellCenters(), rI) if fill: for i, ri_ in enumerate(rI): # slope == True produce unstable behavior .. check!!!!!! rI[i] = mt.fillEmptyToCellArray(meshI, ri_, slope=False) rI.round(1e-6) if len(rI) == 1: # copy here because of missing refcounter TODO return rI[0].array() return rI
[docs] def transFwdArchiePhi(rFluid=20, m=2): r"""Archies law transformation function for resistivity(porosity). .. math:: \rho & = a\rho_{\text{fl}}\phi^{-m}\S_w^{-n} \\ \rho & = \rho_{\text{fl}}\phi^(-m) = \left(\phi/\rho_{\text{fl}}^{-1/n}\right)^{-n} See also :py:mod:`pygimli.physics.petro.resistivityArchie` Returns ------- trans : :gimliapi:`GIMLI::RTransPower` Transformation function Examples -------- >>> from pygimli.physics.petro import * >>> phi = 0.3 >>> tFAPhi = transFwdArchiePhi(rFluid=20) >>> r1 = tFAPhi.trans(phi) >>> r2 = resistivityArchie(rFluid=20.0, porosity=phi, ... a=1.0, m=2.0, sat=1.0, n=2.0) >>> print(r1-r2 < 1e-12) True >>> phi = [0.3] >>> tFAPhi = transFwdArchiePhi(rFluid=20) >>> r1 = tFAPhi.trans(phi) >>> r2 = resistivityArchie(rFluid=20.0, porosity=phi, ... a=1.0, m=2.0, sat=1.0, n=2.0) >>> print((r1-r2 < 1e-12)[0]) True """ return pg.trans.TransPower(-m, rFluid**(1./m))
[docs] def transInvArchiePhi(rFluid=20, m=2): # phi(rho) """Inverse Archie transformation function porosity(resistivity). # rFluid/rho = phi^m ==> phi = (rFluid/rho)^(1/m) = (rho/rFluid)^(-1/m) See --- :py:mod:`pygimli.physics.petro.transFwdArchiePhi` """ return pg.trans.TransPower(-1/m, rFluid)
[docs] def transFwdArchieS(rFluid=20, phi=0.4, m=2, n=2): # rho(S) """Inverse Archie transformation function resistivity(saturation).""" # rho = rFluid * phi^(-m) S^(-n) return pg.trans.TransPower(-n, (rFluid*phi**(-m))**(1/n))
[docs] def transInvArchieS(rFluid=20, phi=0.4, m=2, n=2): # S(rho) """Inverse Archie transformation function saturation(resistivity).""" # rFluid/rho = phi^m S^n => S=(rFluid/rho/phi^m)^(1/n) # S = (rho/rFluid/phi^-m)^(-1/n) return pg.trans.TransPower(-1/n, rFluid*phi**(-m))
def test_Archie(): """Test Archie.""" import unittest dx = 0.01 phivec = np.arange(dx, 0.5, dx) swvec = np.arange(dx, 1, dx) phi0 = 0.4 # 40% porosity rhow = 20 # 20 Ohmm tap water tFAPhi = transFwdArchiePhi(rFluid=rhow) tFAS = transFwdArchieS(rFluid=rhow, phi=phi0) tIAPhi = transInvArchiePhi(rFluid=rhow) tIAS = transInvArchieS(rFluid=rhow, phi=phi0) ax = plt.subplots()[1] # direct function rA = resistivityArchie(rFluid=rhow, porosity=phivec) rS = resistivityArchie(rFluid=rhow, porosity=phi0, sat=swvec) ax.semilogy(phivec, rA, 'b-') ax.semilogy(swvec, rS, 'r-') # forward transformation fA = tFAPhi.trans(phivec) fS = tFAS.trans(swvec) np.testing.assert_allclose(rA, fA, rtol=1e-12) np.testing.assert_allclose(rS, fS, rtol=1e-12) ax.semilogy(phivec, fA, 'bx', markersize=10) ax.semilogy(swvec, fS, 'rx', markersize=10) # inverse transformation iA = tIAPhi.invTrans(phivec) iS = tIAS.invTrans(swvec) np.testing.assert_allclose(rA, iA, rtol=1e-12) np.testing.assert_allclose(rS, iS, rtol=1e-12) ax.semilogy(phivec, iA, 'bo') ax.semilogy(swvec, iS, 'ro') plt.show() if __name__ == "__main__": pass