Electrical modeling - Potential distribution

[Pbottelin]

Hi everyone, I am quite new to Pygimli, so forgive my post. I think this might help people like me to familiarize with the package and build the first simple models.

I am trying to compute the electrical potential distribution in the ground, as described in the 2.5D Geoelectrics example. In my example, I would like to model a pole-pole configuration (one injection electrode at the infinite, so it is not drawn in the model) for 'mise-à-la-masse' surveys.

I played around for quite long, but I have some questions, listed as comments in the code below (boundary conditions, solver, wavenumber, ...). Please feel free to comment my code! Best regards,

Pierre

# -*- coding: utf-8 -*-

"""
DESCRIPTION
This is a script which attempts to perform simple 2D modeling of electrical potential
distribution in a isotropic, layered, semi-infinite medium (i.e. layered ground).

It aims at estimating the opportunity to perform the "mise-a-la-masse" survey
to map the shape of a buried conductor. The method relies on the deformation
of the iso-potential lines along the buried target shape.
"""

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.solver import solve
from pygimli.viewer import show

""" **************************** PARAMETERS ****************************"""

sourcePosA = [0.0, 0.0]  #position of injection electrode (x, y) [in meters]

""" **************************** FUNCTIONS ****************************"""

def uAnalytical(p, sourcePos, k):
    """Analytical solution for one source location."""
    ## I don't really understand the meaning of this function
    ## describing the source term

    r1A = (p - sourcePos).abs()
    # Mirror on surface at depth=0
    r2A = (p - pg.RVector3(1.0, -1.0, 1.0) * sourcePos).abs() #what is RVector3

    # need rho here! Why this?
    if r1A > 1e-12 and r2A > 1e-12: #what is performed in this test?
        return (pg.besselK0(r1A * k) + pg.besselK0(r2A * k)) / (2.0 * np.pi)
    else:
        return 0.

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 value for the Neumann type
    Boundary conditions for the boundaries in the subsurface.
    """
    sourcePos = userData['sourcePos']
    k = userData['k']
    r1 = boundary.center() - sourcePos
    # Mirror on surface at depth=0
    r2 = boundary.center() - pg.RVector3(1.0, -1.0, 1.0) * sourcePos
    r1A = r1.abs()
    r2A = r2.abs()
    n = boundary.norm()

    # need rho here ! #why?
    if r1A > 1e-12 and r2A > 1e-12:
        return k * ((r1.dot(n)) / r1A * pg.besselK1(r1A * k) +
                    (r2.dot(n)) / r2A * pg.besselK1(r2A * k)) / \
            (pg.besselK0(r1A * k) + pg.besselK0(r2A * k))
    else:
        return 0.

def pointSource(cell, f, userData):
    """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)
    """
    sourcePos = userData['sourcePos']

    if cell.shape().isInside(sourcePos):
        f.setVal(cell.N(cell.shape().rst(sourcePos)), cell.ids())

""" **************************** MAIN ****************************"""

""" 1. mesh """
layer1 = pg.meshtools.createPolygon([[-50, 0], [+50, 0],
                           [+50, -1.0], [-50, -1.0]],
                          isClosed=True, marker=1, boundaryMarker=1)
layer2 = pg.meshtools.createPolygon([[-50, -1.0], [+50, -1.0],
                           [+50, -100], [-50, -100]],
                          isClosed=True, marker=2, boundaryMarker=2)

geom = pg.meshtools.mergePLC([layer1, layer2])

mesh = pg.meshtools.createMesh(geom, quality=30.0, area=0.1)

""" 2. Parametrization """
k = 1e-3 #wavenumber for the equation (see refs?)
#How to choose this value? what is its influence?

""" 3. model resistivities """
Rho1 = 1 #resistivity of first layer [Ohm.m]
Rho2 =  1000 #resistivity of second layer [Ohm.m]

a_list = [[1, 1/Rho1], [2, 1/Rho2]] #list of [region marker, resistivity]
b_list = [a * (k**2) for (i,a) in a_list] #is is correct?

""" 4. boundary conditions """

#condition based on du/dn (gradient of field along the boundary normal)
NeumannBC = [[mesh.findBoundaryByMarker(1), 0], #left
              [mesh.findBoundaryByMarker(2), 0], #right
              [mesh.findBoundaryByMarker(4), 0]] #bottom

#constant, scalar condition applied on boundary)
DirichletBC = [[mesh.findBoundaryByMarker(1), 20], # left
                [mesh.findBoundaryByMarker(2), 20], # right
             [mesh.findBoundaryByMarker(4), 20]]  # bottom

#mixed BC, with value derived from the analytical solution of potential diffusion in a
#semi infinite, homogeneous medium)
robBC = [[1, mixedBC],  # left boundary
         [2, mixedBC],  # right boundary
         [4, mixedBC]]  # bottom boundary

#top is not given because already taken into account in the source term?

""" 5. solver """
u = solve(mesh, a=a_list, b=b_list, f=pointSource,
          bc={'Robin': robBC},
          userData={'sourcePos': sourcePosA, 'k': k},
          verbose=True)

#how do I choose the potential at the injection electrode?
#the current intensity I is not asked for (only one electrode)
#the solution u is normalized?

# meaning of the following paragraph?
# uAna = pg.RVector(map(lambda p__: uAnalytical(p__, sourcePosA, k),
#                       grid.positions()))
# uAna -= pg.RVector(map(lambda p__: uAnalytical(p__, sourcePosB, k),
#                        grid.positions()))
# err = (1.0 -u/uAna) * 100.0
# print("error min max", min(err), max(err))

""" **************************** FIGURE ****************************"""

""" FIG1 """
#figure which draws the mesh
#pg.viewer.showMesh(geom, boundaryMarker=True, showMesh=True, marker=True, showBoundary=True)

""" FIG2 """
#drawing the results. Potential distribution in the 2-layered ground.
ax, cbar = show(mesh, data=u, fillContour=True, colorBar=True, cMap="RdBu_r",
          orientation='horizontal', label='U [V]', nLevs=51,
          logScale=False, hold=True, showMesh=True, marker=True)
ax.set_xlabel('Distance [m]')
ax.set_ylabel('Depth [m]')
ax.set_aspect('equal')
#ax.set_xlim(-20, 20)
#ax.set_ylim(-20, 0)

# Additional to the image of the potential we want to see the current flow too.
# The current flows along the gradient of our solution and can be plotted as
# stream lines. On 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 basis to draw the
# streams. If the drawStreams get 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, grid, u, coarseMesh=gridCoarse, color='Black')
plt.show()
pg.wait()

[carsten-forty2]

Hello Pierre,

I just try to answer your questions without modifying your script.

The example in examples/dc_and_ip/plot_mod-dc-2d.py is supposed to show a simple way on how the pg finite element solver can be used to solve a little part of the 2.5D geoelectrics problem. Its not supposed to being a complete and flexible solution to the ERT forward problem. If I find the time I will improve the example, so thank you for pointing out important spots that needs clarification. A more advanced reference implementation of the ERT problem you can find in pygimli/physics/ert/ert.py`. The most advanced solution to solve any ERT problems would be the BERT package: https://gitlab.com/resistivity-net/bert

To your question as far I can see them:

1. ## I don't really understand the meaning of this function ## describing the source term

   uAnalytical(p, sourcePos, k)

   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.

2. #what is performed in this test? This test handles the singular 1/r case for a distance r=0 at the source position.

3. k = 1e-3 #wave number for the equation (see refs?) #How to choose this value? what is its influence?

   Its an arbitrary number to fulfill the equation to be solved. The 2.5D wave number approach is to handle the 3D characteristics of the current density due to a point source in 2 dimensions and solves for u(x,y,k). If you want to have your potential distribution u(x,y, z=0) you need so solve for multiple wave numbers and integrate all solutions. But the latter is beyond the focus of this simple example. You can find a possible solution in the reference implementation.

4. b_list = [a * (k**2) for (i,a) in a_list] #is is correct? Parameter b for pg.solver should be only k*k regarding the equation to be solved and the api reference for pg.solve. https://www.pygimli.org/pygimliapi/_generated/pygimli.solver.html#pygimli.solver.solveFiniteElements

5. #top is not given because already taken into account in the source term?

   Top is probably ground/air surface and hence be treated a homogeneous Neumann (no flow) condition which is default for all FEM boundaries.

6. #how do I choose the potential at the injection electrode? #the current intensity I is not asked for (only one electrode) #the solution u is normalized?

   u is normalized to 1 A injection current. The potential at the injection electrode is 1/(r=0) and not defined and hence set to zero. You can maybe estimated it to some higher value if you like nice images for post processing.

7. # meaning of the following paragraph? you can compare the analytical exact solution with the numerical solution (for constant sigma=1) here, I need to think about to make some tests in the script with it.

Hope I could clarify things a bit.

Cheers, Carsten

[florian-wagner]

Hello @Pbottelin,

could you please give a quick feedback if Carstens answer helped and if this issue can be closed?

Thanks.

[Pbottelin]

Hi Florian, Carsten, thanks for the help and the reminder. Yes I think I'll manage to go on based on Carsten answers. This issue can be closed! Best regards,

Pierre

[rvb3]

Dear Carsten and Florian,

could you please specify the answer to question 6? You say, the potential at the injection electrode can be changed. But I do not understand how.

Thanks and best regards, Rhea

[carsten-forty2]

Assuming the Point Electrode Model (PEM), the potential u at the injection electrode is 1/r with r=0, so it is undefined and cannot be determined. If you use total field calculation it is set through the solution of the finite element linear solver to a meaningless high value. For the solution with the singularity removal technique, it is determined as 1/r(half distance to the next nearest node)

However, if you want nice images of your potential distribution, you can set it to any arbitrary high value.