ERT shclumberger synthetic model, in order to forward simulate and invert a landslide

[AY131313]

Hi, I'm trying to make a 3 layered subsurface landslide geometry.
I provide the code I am running and would appreciate any insights and ideas to help me get through conducting a good inversion.
The array I'm trying to use is Schlumberger.

import matplotlib.pyplot as plt
import os
from pygimli.viewer.mpl import generateMatrix
from PIL import Image
from numpy import asarray
import numpy as np
from pygimli.meshtools import quality
import random 
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert

# Create geometry as provided
world = mt.createWorld(start=[0, 165], end=[117.5, 110])
layer1 = mt.createPolygon([[0.0, 137], [117.5, 164], [117.5, 162], [0.0, 135]],
                          isClosed=True, marker=1, area=1, boundaryMarker = 1)
layer2 = mt.createPolygon([[0.0, 126], [0.0, 135], [117.5, 162], [117.5, 153]],
                          isClosed=True, marker=2, boundaryMarker = 2)
layer3 = mt.createPolygon([[0.0, 110], [0.0, 126], [117.5, 153], [117.5, 110]],
                          isClosed=True, marker=3, boundaryMarker = 3)
freeAir = mt.createPolygon([[117.5,165] , [60,165] ,[0,165], [0,137]] , marker = 4, isClosed = True)
geom = layer1 + layer2 + layer3

# Define the number of electrodes and the desired positions
num_electrodes = 30
x_values = np.linspace(0.0, 117.5, num_electrodes)
y_values = np.linspace(137, 164, num_electrodes)
electrode_positions = np.column_stack((x_values, y_values))

# Create electrodes at the specified positions
for pos in electrode_positions:
    geom.createNode(pos[0], pos[1], 0.0)

ax, cb = pg.show(geom, showNodes=True)

#mesh creation
mesh = mt.createMesh(geom, quality = 32, area= 1)
pg.show(mesh);
plt.savefig("MODEL2.jpg",dpi=1200, bbox_inches="tight")


# Debugging: Check if any nodes have no cells
orphaned_nodes = []
for pos in electrode_positions:
    node = mesh.findNearestNode(pos)
    if node is None:
        orphaned_nodes.append(pos)

if orphaned_nodes:
    pg.warn("Orphaned nodes found:", len(orphaned_nodes))
    for i, pos in enumerate(orphaned_nodes):
        pg.warn(f"Orphaned Node {i+1}: Position {pos}")
        
        
# Create electrode data scheme
scheme = ert.createData(elecs=electrode_positions, schemeName="gr")
scheme['k'] = ert.geometricFactors(scheme) 
ert.show(scheme, 'k')
print(scheme)        


# Update electrode positions to match the mesh nodes
for i, pos in enumerate(scheme.sensorPositions()):
    nodeID = mesh.findNearestNode(pos)
    node = mesh.node(nodeID)
    electrode_positions[i] = [node.x(), node.y()]
pg.show(mesh)


# Debugging: Check if any nodes have no cells
orphaned_nodes = []
for pos in electrode_positions:
    node = mesh.findNearestNode(pos)
    if node is None:
        orphaned_nodes.append(pos)

if orphaned_nodes:
    pg.warn("Orphaned nodes found:", len(orphaned_nodes))
    for i, pos in enumerate(orphaned_nodes):
        pg.warn(f"Orphaned Node {i+1}: Position {pos}")
        
        
 #Define the resistivities for different regions using the region markers
resistivityMap = {
    1: 10.0,   # Secondary compacted regions with resistivity 10 ohm.m
    2: 75.0,   # Region with high moisture content with resistivity 75 ohm.m
    3: 150.0   # Hezar-Dareh alluvium with resistivity 150 ohm.m
}

# Assign resistivity values to the mesh using the region markers
resistivityValues = [resistivityMap[marker] for marker in mesh.cellMarkers()]

# Take a look at the mesh and the resistivity distribution
pg.show(mesh, data=resistivityValues, label=pg.unit('res'), showMesh=True)
plt.savefig("MODEL5.jpg",dpi=1200, bbox_inches="tight")


#forward simulation
data = ert.simulate(mesh, scheme=scheme, res=resistivityValues, noiseLevel = 0.03)
                    
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())

pg.info('Simulated rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))
pg.info('Selected data noise %(min/max)', min(data['err'])*100, max(data['err'])*100)

after the last code snnipet related to forward simulation I always recieve an error like below:

RuntimeError Traceback (most recent call last)
Cell In[25], line 1
----> 1 data = ert.simulate(mesh, scheme=scheme, res=resistivityValues)
4 pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
5 pg.info('Simulated data', data)

File E:\Anaconda3\envs\pg-transform2022\lib\site-packages\pygimli\physics\ert\ert.py:121, in simulate(mesh, scheme, res, **kwargs)
118 # fop = self.createForwardOperator(useBert=True, # self.useBert,
119 # sr=sr, verbose=verbose)
120 fop.data = scheme
--> 121 fop.setMesh(mesh, ignoreRegionManager=True)
123 rhoa = None
124 phia = None

File E:\Anaconda3\envs\pg-transform2022\lib\site-packages\pygimli\frameworks\modelling.py:655, in MeshModelling.setMesh(self, mesh, ignoreRegionManager)
652 super(Modelling, self).setMesh(mesh, ignoreRegionManager=True)
653 pass
--> 655 self.setMeshPost(mesh)
656 return
657 self.clearRegionProperties()

File E:\Anaconda3\envs\pg-transform2022\lib\site-packages\pygimli\physics\ert\ertModelling.py:185, in ERTModelling.setMeshPost(self, mesh)
183 def setMeshPost(self, mesh):
184 """"""
--> 185 self._core.setMesh(mesh, ignoreRegionManager=True)

RuntimeError: C:/msys64/home/halbm/gimli/gimli/core/src/mesh.cpp:785 GIMLI::Cell* GIMLI::Mesh::findCell(const RVector3&, size_t&, bool) const no cells for this node. This is a corrupt mesh

P.S. I provide the pictures of the initial geometry, the Mesh, the data scheme, and the geometry merged with the mesh when injected resistivity values using the resistivity Map.

[halbmy]

I am getting an error already when I add the electrodes and want to show (i.e. in pg.show).

Plotting the mesh and the electrodes before (pg.show(geo)[0].plot(x_values, y_values))

it becomes clear that there are several issues

• the electrodes should not be on the side boundary
• some boundary around the electrodes is needed for boundary conditions of the forward problem
• creating a node where there is already one should be avoided by createNodeWithCheck
• there can be tolerance issues, i.e. the nodes are very close but below or above the topography line

I therefore suggest creating the mesh from the electrode positions so that they are included, using meshtools functions as done in several examples.

[halbmy]

Any progress?

[AY131313]

Hello Dr. Gunther. Actually yes. I get an answer. I will mail it .I really appreciate your help and support. An honor to be in touch with you. Regards. Amir Yazdanpanah

…

On Tue, 17 Oct 2023, 23:26 Thomas Günther, ***@***.***> wrote: Any progress? — Reply to this email directly, view it on GitHub <#578 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A7UQXQGPDHSSZEDXZ7FZHLDX73PFHAVCNFSM6AAAAAA4IW3S5OVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TGMBXG43TC> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[AY131313]

Now Im trying to add a topography to my final inversion results. Any ideas I can start with?? On Wed, 18 Oct 2023, 01:26 amir yazdanpanah, ***@***.***> wrote:

…

Hello Dr. Gunther. Actually yes. I get an answer. I will mail it .I really appreciate your help and support. An honor to be in touch with you. Regards. Amir Yazdanpanah On Tue, 17 Oct 2023, 23:26 Thomas Günther, ***@***.***> wrote: > Any progress? > > — > Reply to this email directly, view it on GitHub > <#578 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/A7UQXQGPDHSSZEDXZ7FZHLDX73PFHAVCNFSM6AAAAAA4IW3S5OVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TGMBXG43TC> > . > You are receiving this because you authored the thread.Message ID: > ***@***.***> >

[halbmy]

You should start right with the correct geometry, including the topography.

[AY131313]

You mean by creating a polygon?

…

On Wed, 18 Oct 2023, 10:20 Thomas Günther, ***@***.***> wrote: You should start right with the correct geometry, including the topography. — Reply to this email directly, view it on GitHub <#578 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A7UQXQFK5MJKMABEGNTWGZ3X7535DAVCNFSM6AAAAAA4IW3S5OVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TGMJRGQ2DK> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[halbmy]

By creating a mesh that already contains the correct topography, and also by using it in the data (from which the mesh can be created).

[AY131313]

can you please elaborate a little more? I mean how can I get provided with a mesh that already contains topography?

…

On Wed, Oct 18, 2023 at 10:32 PM Thomas Günther ***@***.***> wrote: By creating a mesh that already contains the correct topography, and also by using it in the data (from which the mesh can be created). — Reply to this email directly, view it on GitHub <#578 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A7UQXQDPZQSAE7WJRF6WRWTYAARWHAVCNFSM6AAAAAA4IW3S5OVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TGMJZGEZDM> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[AY131313]

Hello Dr.Gunther,
I have made my ERT and SRT models.
for inverting the ERT field data, I use pybert (pb.importData()).
in the inversion step my model is weird and not representable of my data (That belongs to a landslide and contains topography).

[halbmy]

1. Which format does your data have? Have you tried ert.load()?
2. It seems like the topography is not read properly (numerical are identical with analytical geometric factors, not topo in results).
3. Your pseudosection looks really weird? Why?

[AY131313]

hello. i used a .dat file in the format of res2dinv.
yes I did but ert.load() didn't read my data.
because it belongs to a landslide and doesn't represent it properly.

[halbmy]

ert.load() should be able to read res2dinv format and uses exactly the same function as pybert importData. Maybe we made an error, would be good to check if you can share the data.

[AY131313]

Hello Dr,
I tried ert.load() and I got an inversion result without any topography (the case is related to detecting the sliding surface of a landslide).
I have sent the field data of a profile in res2dinv format to your academic email.

[halbmy]

Something's wrong with the data import. I can't tell unless seing the data.

[AY131313]

the data has been sent to your email address Dr.

[halbmy]

I haven’t got anything. I haven’t even an idea who you are. Von: AY13 ***@***.***> Gesendet: Donnerstag, 16. November 2023 21:49 An: gimli-org/gimli ***@***.***> Cc: Guenther, Thomas ***@***.***>; Comment ***@***.***> Betreff: Re: [gimli-org/gimli] ERT shclumberger synthetic model, in order to forward simulate and invert a landslide (Discussion #578) the data has been sent to your email address Dr. — Reply to this email directly, view it on GitHub<#578 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABWOJEL3BE3WWX4OM74Y623YEZ33DAVCNFSM6AAAAAA4IW3S5OVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TKOJSGQ2DI>. You are receiving this because you commented.Message ID: ***@***.******@***.***>>

[AY131313]

[AY131313]

Hello dr,
any progress here? I mean, have u found the data I emailed to you?
If not, please send me an email address.
regards.