cross-hole forward question

[zmd142536]

Problem description

Regarding the vector division for cross-hole data, the lengths of the dividend and divisor are not equal.

  File ~/anaconda3/envs/pg/lib/python3.10/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File ~/Spyder/untitled25.py:88
    showSensitivity(data, mesh, nc=1)

  File ~/Spyder/untitled25.py:76 in showSensitivity
    normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)

  File ~/anaconda3/envs/pg/lib/python3.10/site-packages/pygimli/core/__init__.py:240 in __newRVectorTrueDiv__
    return __oldRVectorTrueDiv__(a, b)

RuntimeError: ./core/src/vector.h:689       GIMLI::Vector<ValueType>& GIMLI::Vector<ValueType>::operator/=(const GIMLI::Vector<ValueType>&) [with ValueType = double]  22 != 1525

Your environment

Ubuntu 20.04.6 pygimli 1.4.6

Steps to reproduce

Hello sir, I expected to find the optimal layout for cross-hole measurements in practice through forward modeling. However, I encountered an error when I added a circular geological body to the code of my case study. Let's take a look at my code first.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-v0_8")
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.meshtools import polytools as plc#新加用以创建园
from pygimli.physics import ert
xBH = [-2, 2]  # two boreholes, 4m spaced
zBH = -np.arange(1., 8.1, 1.0)
data = ert.DataContainer()
for x in xBH:
    for z in zBH:
        data.createSensor([x, z])

print(data)
world = mt.createWorld(start=[-4, -10], end=[4, 0], worldMarker=True, marker=1)
c1 = plc.createCircle(pos=[0.0, -5.0], radius=0.5, marker=2, area=.1)#创建园
geom = c1 + world
for pos in data.sensorPositions():
    world.createNode(pos)

mesh = mt.createMesh(geom, area=.1, quality=34, smooth=[1, 10])
ax, cb = pg.show(mesh, markers=True, showMesh=True)
ax.plot(pg.x(data), pg.y(data), "mo");

def plotABMN(ax, scheme, idx):
    """ Visualize four-point configuration on given axes. """
    def getABMN(scheme, idx):
        """ Get coordinates of four-point cfg with id `idx` from DataContainerERT
        `scheme`."""
        coords = {}
        for elec in "abmn":
            elec_id = int(scheme(elec)[idx])
            if elec_id >= 0:
                elec_pos = scheme.sensorPosition(elec_id)
                coords[elec] = elec_pos.x(), elec_pos.y()
        return coords

    coords = getABMN(scheme, idx)
    for elec in coords:
        x, y = coords[elec]
        if elec in "ab":
            color = "green"
        else:
            color = "magenta"
        ax.plot(x, y, marker=".", color=color, ms=10)
        ax.annotate(elec.upper(), xy=(x, y), size=12, ha="center", #fontsize=10, 
                    bbox=dict(boxstyle="round", fc=(0.8, 0.8, 0.8), ec=color), 
                    xytext=(20*np.sign(x), 0), textcoords='offset points', 
                    arrowprops=dict(arrowstyle="wedge, tail_width=.5", fc=color, ec=color,
                                    patchA=None, alpha=0.75))
        ax.plot(coords["a"][0],)

def showSensitivity(shm, mesh, idx=None, ax=None, nc=None):
    idx = idx or range(shm.size())
    nc = nc or len(idx)
    if ax is None:
        fig, ax = plt.subplots(len(idx)//nc, nc, sharex=True, sharey=True, figsize=(12, 8))

    fop = ert.ERTModelling(verbose=False)
    fop.setData(shm)
    fop.setMesh(mesh)
    model = np.ones(mesh.cellCount())
    fop.createJacobian(model)
    for i, ind in enumerate(idx):
        a = ax.flat[i]
        sens = fop.jacobian()[ind]
        normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)
        normsens /= np.max(normsens)
        pg.show(mesh, normsens, cMap="bwr", colorBar=False, ax=a,
                        label="sensitivity", nLevs=3, cMin=-1, cMax=1)
        plotABMN(a, shm, ind);
        # a.set_xlim(-2, 12)
        # a.set_ylim(-5, 1)
data.createFourPointData(data.size(), 1, 9, 2, 10)
data.createFourPointData(data.size(), 1, 14, 2, 15)
data.createFourPointData(data.size(), 1, 12, 7, 13)
data.createFourPointData(data.size(), 1, 9, 6, 14)
print(data)       
showSensitivity(data, mesh, nc=1)

output error:

  File ~/anaconda3/envs/pg/lib/python3.10/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File ~/Spyder/untitled25.py:88
    showSensitivity(data, mesh, nc=1)

  File ~/Spyder/untitled25.py:76 in showSensitivity
    normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)

  File ~/anaconda3/envs/pg/lib/python3.10/site-packages/pygimli/core/__init__.py:240 in __newRVectorTrueDiv__
    return __oldRVectorTrueDiv__(a, b)

RuntimeError: ./core/src/vector.h:689       GIMLI::Vector<ValueType>& GIMLI::Vector<ValueType>::operator/=(const GIMLI::Vector<ValueType>&) [with ValueType = double]  22 != 1525
...

> If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

In your case, you have two regions: a main region with marker=1 and a circle with marker=2 inside. In such cases, the region with the lower marker is automatically treated as background. This needs to be deactivated by

fop.setRegionProperties(1, background=False)

[zmd142536]

你好先生Halbmy，我遇到了一些问题超出了我的专业知识，当我完成建造前进的模型和编制计算。 能不能请你看一看在我的代码？

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-v0_8")
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.meshtools import polytools as plc#新加用以创建园
from pygimli.physics import ert
xBH = [-2, 2]  # two boreholes, 4m spaced
zBH = -np.arange(1., 8.1, 1.0)
data = ert.DataContainer()
for x in xBH:
    for z in zBH:
        data.createSensor([x, z])

print(data)
world = mt.createWorld(start=[-4, -10], end=[4, 0], worldMarker=True, marker=1)
c1 = plc.createCircle(pos=[0.0, -5.0], radius=0.5, marker=2, area=.1)#创建园
geom = c1 + world
for pos in data.sensorPositions():
    world.createNode(pos)

mesh = mt.createMesh(geom, area=.1, quality=34, smooth=[1, 10])
ax, cb = pg.show(mesh, markers=True, showMesh=True)
ax.plot(pg.x(data), pg.y(data), "mo");

def plotABMN(ax, scheme, idx):
    """ Visualize four-point configuration on given axes. """
    def getABMN(scheme, idx):
        """ Get coordinates of four-point cfg with id `idx` from DataContainerERT
        `scheme`."""
        coords = {}
        for elec in "abmn":
            elec_id = int(scheme(elec)[idx])
            if elec_id >= 0:
                elec_pos = scheme.sensorPosition(elec_id)
                coords[elec] = elec_pos.x(), elec_pos.y()
        return coords

    coords = getABMN(scheme, idx)
    for elec in coords:
        x, y = coords[elec]
        if elec in "ab":
            color = "green"
        else:
            color = "magenta"
        ax.plot(x, y, marker=".", color=color, ms=10)
        ax.annotate(elec.upper(), xy=(x, y), size=12, ha="center", #fontsize=10, 
                    bbox=dict(boxstyle="round", fc=(0.8, 0.8, 0.8), ec=color), 
                    xytext=(20*np.sign(x), 0), textcoords='offset points', 
                    arrowprops=dict(arrowstyle="wedge, tail_width=.5", fc=color, ec=color,
                                    patchA=None, alpha=0.75))
        ax.plot(coords["a"][0],)

def showSensitivity(shm, mesh, idx=None, ax=None, nc=None):

    idx = idx or range(shm.size())
    nc = nc or len(idx)
    if ax is None:
        fig, ax = plt.subplots(len(idx)//nc, nc, sharex=True, sharey=True, figsize=(12, 8)) 
    fop = ert.ERTModelling(verbose=False)
    fop.setData(shm)
    fop.setMesh(mesh)
    fop.setRegionProperties(1, background=False)  #确保低范围标记不被背景化
    model = np.ones(mesh.cellCount())
    fop.createJacobian(model)
    for i, ind in enumerate(idx):
        a = ax.flat[i]
        sens = fop.jacobian()[ind]
        normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)
        normsens /= np.max(normsens)
        pg.show(mesh, normsens, cMap="bwr", colorBar=False, ax=a,
                        label="sensitivity", nLevs=3, cMin=-1, cMax=1)
        plotABMN(a, shm, ind);
        # a.set_xlim(-2, 12)
        # a.set_ylim(-5, 1)
data.createFourPointData(data.size(), 1, 9, 2, 10)
data.createFourPointData(data.size(), 1, 9, 4, 12)
data.createFourPointData(data.size(), 1, 9, 6, 14)
data.createFourPointData(data.size(), 2, 10, 3, 11)
data.createFourPointData(data.size(), 2, 10, 5, 13)
data.createFourPointData(data.size(), 2, 10, 7, 15)
data.createFourPointData(data.size(), 3, 11, 4, 12)
data.createFourPointData(data.size(), 3, 11, 6, 14)
data.createFourPointData(data.size(), 4, 12, 5, 13)
data.createFourPointData(data.size(), 4, 12, 7, 15)
data.createFourPointData(data.size(), 5, 13, 6, 14)
data.createFourPointData(data.size(), 2, 10, 3, 11)
data.createFourPointData(data.size(), 3, 11, 4, 12)
data.createFourPointData(data.size(), 4, 12, 5, 13)
data.createFourPointData(data.size(), 5, 13, 6, 14)
data.createFourPointData(data.size(), 6, 14, 7, 15)
print(data)     
rhomap = [[1, 100.],
          [2, 500.]]

mgr = ert.ERTManager()
data = ert.simulate(mesh, scheme=plotABMN, res=rhomap, noiseLevel=0.01, noiseAbs=1e-6, seed=1337)
pg.viewer.mpl.drawSensors(ax, data, diam=1.0, facecolor='white', edgecolor='black')
raise SystemExit

Error output： `Traceback (most recent call last):

File ~/anaconda3/envs/pg/lib/python3.10/site-packages/spyder_kernels/py3compat.py:356 in compat_exec exec(code, globals, locals)

File ~/Spyder/kuakong-youdizhiti(baocuo).py:105 data = ert.simulate(mesh, scheme=plotABMN, res=rhomap, noiseLevel=0.01, noiseAbs=1e-6, seed=1337)

File ~/anaconda3/envs/pg/lib/python3.10/site-packages/pygimli/physics/ert/ert.py:163 in simulate if not scheme.allNonZero('k') and not calcOnly:

AttributeError: 'function' object has no attribute 'allNonZero'`

[halbmy]

You cannot just paste your mess of code and expect other people to start thinking. Just have a look at the error message: Obviously, plotABMN is of the wrong type, as it is a function. ert.simulate expects a data container as the docstring says, in your case data.

[zmd142536]

I apologize for this, Mr. Halbmy. Please forgive my foolishness.But I do need professional help indeed. I've made the code modifications and it seems to have run successfully. However, I encountered difficulties during the data plotting step. I'm unsure if it's due to my electrode configuration being unreasonable or if the data plotting scheme across the boreholes isn't suitable, resulting in less than ideal data plotting.

Steps to reproduce

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-v0_8")
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.meshtools import polytools as plc#Newly added for creating circles.
from pygimli.physics import ert
from pygimli.viewer.mpl import showDataContainerAsMatrix
#The function to read data from a file is defined as follows
def read_txt_file(filename):
    """Read data from a TXT file and return it."""
    with open(filename, 'r') as file:
        lines = file.readlines()

    # Parse the data and store it as a list.
    data_list = []
    for line in lines:
        # Assuming the data is comprised of four integers separated by spaces
        values = line.strip().split()
        data_list.append([int(value) for value in values])

    return data_list
#Create sensors and data containers
xBH = [-2, 2]  # two boreholes, 4m spaced
zBH = -np.arange(1., 16.1, 1.0)
data = ert.DataContainer()
for x in xBH:
    for z in zBH:
        data.createSensor([x, z])

print(data)
world = mt.createWorld(start=[-4, -20], end=[4, 0], worldMarker=True, marker=1)
c1 = plc.createCircle(pos=[0.0, -8.0], radius=0.5, marker=2, area=.1)#Create a circle.
geom = c1 + world
for pos in data.sensorPositions():
    world.createNode(pos)

mesh = mt.createMesh(geom, area=.1, quality=34, smooth=[1, 10])
#Read the data file and pass the data to the createFourPointData function
data_from_file = read_txt_file('/home/zmd/Spyder/paojifangshi1/paojifangshi2.txt')
for data_point in data_from_file:
    data.createFourPointData(data.size(), *data_point)
ax, cb = pg.show(mesh, markers=True, showMesh=True)
ax.plot(pg.x(data), pg.y(data), "mo");

def plotABMN(ax, scheme, idx):
    """ Visualize four-point configuration on given axes. """
    def getABMN(scheme, idx):
        """ Get coordinates of four-point cfg with id `idx` from DataContainerERT
        `scheme`."""
        coords = {}
        for elec in "abmn":
            elec_id = int(scheme(elec)[idx])
            if elec_id >= 0:
                elec_pos = scheme.sensorPosition(elec_id)
                coords[elec] = elec_pos.x(), elec_pos.y()
        return coords

    coords = getABMN(scheme, idx)
    for elec in coords:
        x, y = coords[elec]
        if elec in "ab":
            color = "green"
        else:
            color = "magenta"
        ax.plot(x, y, marker=".", color=color, ms=10)
        ax.annotate(elec.upper(), xy=(x, y), size=12, ha="center", #fontsize=10, 
                    bbox=dict(boxstyle="round", fc=(0.8, 0.8, 0.8), ec=color), 
                    xytext=(20*np.sign(x), 0), textcoords='offset points', 
                    arrowprops=dict(arrowstyle="wedge, tail_width=.5", fc=color, ec=color,
                                    patchA=None, alpha=0.75))
        ax.plot(coords["a"][0],)

def showSensitivity(shm, mesh, idx=None, ax=None, nc=None):

    idx = idx or range(shm.size())
    nc = nc or len(idx)
    if ax is None:
        fig, ax = plt.subplots(len(idx)//nc, nc, sharex=True, sharey=True, figsize=(12, 8)) 
    fop = ert.ERTModelling(verbose=False)
    fop.setData(shm)
    fop.setMesh(mesh)
    fop.setRegionProperties(1, background=False)  #Ensure that the low-range markers are not backgrounded.
    model = np.ones(mesh.cellCount())
    fop.createJacobian(model)
    for i, ind in enumerate(idx):
        a = ax.flat[i]
        sens = fop.jacobian()[ind]
        normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)
        normsens /= np.max(normsens)
        pg.show(mesh, normsens, cMap="bwr", colorBar=False, ax=a,
                        label="sensitivity", nLevs=3, cMin=-1, cMax=1)
        plotABMN(a, shm, ind);
        # a.set_xlim(-2, 12)
        # a.set_ylim(-5, 1)

rhomap = [[1, 100.],
          [2, 500.]]
mgr = ert.ERTManager()
data = ert.simulate(mesh, data, res=rhomap, noiseLevel=0.01, noiseAbs=1e-6, seed=1337)
pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
pg.info('Simulated data', data)
pg.info('The data contains:', data.dataMap().keys())
ert.show(data)

### error output ： 20/03/24 - 16:30:52 - pyGIMLi - WARNING - Something gone wrong while drawing data. Try fallback with equidistant electrodes.

There are 32 sensors in the data, but the number of nonzero entries is 0. Does this have anything to do with it?

[zmd142536]

[Uploading paojifangshi2.txt…]()

[halbmy]

The print data command comes directly after generation of the sensors

data = ert.DataContainer()
for x in xBH:
    for z in zBH:
        data.createSensor([x, z])

print(data)

So of course there are no data at this point but only sensors. The arrays are added after it

data_from_file = read_txt_file('paojifangshi2.txt')
for data_point in data_from_file:
    data.createFourPointData(data.size(), *data_point)

I cannot download your file this as something went wrong when you uploaded them.

[zmd142536]

paojifangshi2.txt This is my way of extreme running, thank you.

[halbmy]

So what exactly is the problem? The pseudosection image?

There is no unique way to plot pseudosections for crosshole cases. The default ert.show(data) (data.show()) is only good for surface data. You can try data.show(style="A-M").

[halbmy]

(I just used "AB-MN")

[zmd142536]

Thank you, sir, as you said, this is related to the running mode. It seems I've got the answer I wanted, now I just need to optimize the running mode of the electrode to get good results.

here's one last question: regarding the mesh after inversion, I used unstructured mesh for modeling, but after inversion, I got a square grid. What controls this?

# %% inversion outmesh 
inner = mt.createGrid(np.arange(-4, 4, 0.5), np.arange(-20.0, .1, 0.5),
                      marker=2)
mesh = mt.appendTriangleBoundary(inner, xbound=30, ybound=82)
ax, _ = pg.show(mesh, markers=True)
# %% inversion calculation 
mgr = ert.ERTManager('AMBN1.dat')
inv = mgr.invert(data, mesh, secNodes=4, lam=10,verbose=True)
mgr.showResult()

[halbmy]

You control the inversion mesh by yourself, in this case a grid (surrounded by a trianglular outer mesh for forward calculation), but of course you can create a triangle mesh as done in other examples.

[halbmy]

As you "seem to got the answer wanted", I am closing the issue.