Geoelectrics in 2.5DΒΆ

Aim is to perform geoelectrical (DC resistivity) forward modelling in 2.5 D, i.e. a 2D conductivity distribution and 3D point sources.

Let us start with the governing partial differential equation of Poisson type

\[\nabla\cdot(\sigma\nabla u)=-I\delta(\vec{r}-\vec{r}_{\text{s}}) \in R^3\]

The source term (point electrode) is three-dimensional, but the distribution of the electrical # conductivity \(\sigma(x,y)\) should be 2D so we apply a Fourier cosine transform from \(u(x,y,z) \mapsto u(x,k,z)\) with the wave number \(k\). (\(D^{(a)}(u(x,y,z)) \mapsto i^{|a|}k^a u(x,z)\))

\[\begin{split}\nabla\cdot( \sigma \nabla u ) - \sigma k^2 u &=-I\delta(\vec{r}-\vec{r}_{\text{s}}) \in R^2 \\ \frac{\partial }{\partial x}\left(\cdot \sigma \frac{\partial u}{\partial x}\right) + \frac{\partial }{\partial z}\left(\cdot\sigma \frac{\partial u}{\partial z}\right) - \sigma k^2 u & = -I\delta(x-x_{\text{s}})\delta(z-z_{\text{s}}) \in R^2 \\ \frac{\partial u}{\partial \vec{n}} & = 0 \quad\text{at the Surface}\quad (z=0) \\ \frac{\partial u}{\partial \vec{n}} & = a u \quad\text{in the Subsurface}\quad (z<0)\end{split}\]
import matplotlib
import numpy as np
import pygimli as pg

from pygimli.viewer.mpl import drawStreams

We know the exact solution by analytical formulas:

# .. math::
#     u = \frac{1}{2\pi\sigma} \cdot (K_0(\|r-r^+_s\| k)+K_0(\|r-r^-_s\| k))

# with K0 being the Bessel function of first kind, and the normal and mirror
# sources r+ and r-. We define a function for it
def uAnalytical(p, sourcePos, k, sigma=1):
    """Calculates the analytical solution for the 2.5D geoelectrical problem.

    Solves the 2.5D geoelectrical problem for one wave number k.
    It calculates the normalized (for injection current 1 A and sigma=1 S/m)
    potential at position p for a current injection at position sourcePos.
    Injection at the subsurface is recognized via mirror sources along the
    surface at depth=0.

    p : pg.Pos
        Position for the sought potential
    sourcePos : pg.Pos
        Current injection position.
    k : float
        Wave number

    u : float
        Solution u(p)
    r1A = (p - sourcePos).abs()
    # Mirror on surface at depth=0
    r2A = (p - pg.Pos([1.0, -1.0])*sourcePos).abs()

    if r1A > 1e-12 and r2A > 1e-12:
        return  1 / (2.0 * np.pi) * 1/sigma * \
            (pg.math.besselK0(r1A * k) + pg.math.besselK0(r2A * k))
        return 0.

We assume the so-called mixed boundary conditions (Dey & Morrison, 1979). .. math:

\sigma k \frac{{\bf r}\cdot {\bf n}}{{|r|}} \frac{K_1(|r-r_s|k)}{K_0(|r-r_s|k)}
def mixedBC(boundary, userData):
    """Mixed boundary conditions.

    Define the derivative of the analytical solution regarding the outer normal
    direction :math:`\vec{n}`. So we can define the values for mixed boundary
    condition :math:`\frac{\partial u}{\partial \vec{n}} = -au`
    for the boundaries on the subsurface.

    ### ignore surface boundaries for wildcard boundary condition
    if boundary.norm()[1] == 1.0:
        return 0

    sourcePos = userData['sourcePos']
    k = userData['k']
    sigma = userData['s']
    r1 = - sourcePos

    # Mirror on surface at depth=0
    r2 = - pg.Pos(1.0, -1.0) * sourcePos
    r1A = r1.abs()
    r2A = r2.abs()

    n = boundary.norm()
    if r1A > 1e-12 and r2A > 1e-12:
        alpha = sigma * k * (( / r1A * pg.math.besselK1(r1A * k) +
                            ( / r2A * pg.math.besselK1(r2A * k)) / \
            (pg.math.besselK0(r1A * k) + pg.math.besselK0(r2A * k))

        return alpha

        # Note, the above is the same like:
        beta = 1.0
        return [alpha, beta, 0.0]

        return 0.0

We assemble the right-hand side (rhs) for the singular current term by hand since this cannot be done efficiently by pg.solve yet. We basically search for the cell containing the source and project the point using its shape functions N.

def rhsPointSource(mesh, source):
    """Define function for the current source term.

    :math:`\delta(x-pos), \int f(x) \delta(x-pos)=f(pos)=N(pos)`
    Right hand side entries will be shape functions(pos)
    rhs = pg.Vector(mesh.nodeCount())

    cell = mesh.findCell(source)
    rhs.setVal(cell.N(cell.shape().rst(source)), cell.ids())
    return rhs

Now we create a suitable mesh and solve the equation with pg.solve. Note that we use a mesh with quadratic shape functions by calling createP2.

mesh = pg.createGrid(x=np.linspace(-10.0, 10.0, 41),
                     y=np.linspace(-15.0,  0.0, 31))
mesh = mesh.createP2()

sourcePosA = [-5.25, -3.75]
sourcePosB = [+5.25, -3.75]

k = 1e-2
sigma = 1.0
bc={'Robin': {'1,2,4': mixedBC}}
u = pg.solve(mesh, a=sigma, b=-sigma * k*k,
             rhs=rhsPointSource(mesh, sourcePosA),
             bc=bc, userData={'sourcePos': sourcePosA, 'k': k, 's':sigma},

u -= pg.solve(mesh, a=sigma, b=-sigma * k*k,
              rhs=rhsPointSource(mesh, sourcePosB),
              bc=bc, userData={'sourcePos': sourcePosB, 'k': k, 's':sigma},

# The solution is shown by calling

ax =, data=u, cMap="RdBu_r", cMin=-1, cMax=1,
             orientation='horizontal', label='Potential $u$',
             nCols=16, nLevs=9, showMesh=True)[0]

# Additionally to the image of the potential we want to see the current flow.
# The current flows along the gradient of our solution and can be plotted as
# stream lines. By default, the drawStreams method draws one segment of a
# stream line per cell of the mesh. This can be a little confusing for dense
# meshes so we can give a second (coarse) mesh as a new cell base to draw the
# streams. If `drawStreams` gets scalar data, the gradients will be calculated.
gridCoarse = pg.createGrid(x=np.linspace(-10.0, 10.0, 20),
                           y=np.linspace(-15.0,   .0, 20))
drawStreams(ax, mesh, u, coarseMesh=gridCoarse, color='Black')
plot 04 ert 2 5d potential


Mesh:  Mesh: Nodes: 3741 Cells: 1200 Boundaries: 2470
Assembling time:  0.03
Solving time:  1.5030000000000001
Mesh:  Mesh: Nodes: 3741 Cells: 1200 Boundaries: 2470
Assembling time:  0.034
Solving time:  1.4060000000000001

We know the exact solution so we can compare it to the numerical results. Unfortunately, the point source singularity does not allow a good integration measure for the accuracy of the resulting field so we just look for the differences.

uAna = pg.Vector(list(map(lambda p__: uAnalytical(p__, sourcePosA, k, sigma),
uAna -= pg.Vector(list(map(lambda p__: uAnalytical(p__, sourcePosB, k, sigma),

ax =, data=pg.abs(uAna-u), cMap="Reds",
          orientation='horizontal', label='|$u_{exact}$ -$u$|',
          logScale=True, cMin=1e-7, cMax=1e-1,
          nCols=12, nLevs=7,

# print('l2:',
# print('L2:',, mesh)))
# print('H1:',, mesh)))
# np.testing.assert_approx_equal(pg.solver.normL2(uAna-u, mesh),
#                                0.02415, significant=3)
plot 04 ert 2 5d potential

Gallery generated by Sphinx-Gallery