Modelling an 2-D ERT on an 3-D unbounded domain to simulate a field acquisition

[LuFrag]

Greetings PyGimli team,

Following the discussion on issue #443, I am trying to model a 2-D ERT field acquisition using a 3-D model. In this case, an unbounded domain (?).

My goal is to simulate the effect of a 3-D object that is getting closer and closer to the ERT profile. I also setup a modelling with an homogenous 3-D world for comparison (see below).

I was wondering if my domain is correctly set or do I have to specify other type of boundary at the limits of my modelling to reach a better accuracy, specially for the homogenous case ? Or do I have to improve my mesh refinement ?

Thanks you all in advance,
Cheers,
Luis

Here's the result for an homogenous world with 100 Ohm.m as resistivity.

And the result for an homogenous world of 100 Ohm.m with an inhomogeneity of 1 Ohm.m at 1m meter of the ERT profile.

The same but with an inhomogeneity of 1 Ohm.m at 2m meter of the ERT profile.

And here is my code so far:

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
import matplotlib.pyplot as plt

#ERT electrode configuration
n_elec = 48
sep_elec = 2
sensors = np.zeros((n_elec,3))
sensors[:,0]= np.linspace(78,172,n_elec)
sensors[:,1]=[125]

#Creating the acquisition scheme in wenner alpha (wa)
shm =  ert.createData(elecs=sensors,schemeName='wa')

# Calculating geometric factors
data_k = ert.createGeometricFactors(shm)

#Creating the PLC world
world = mt.createWorld(start=[0,0,0], end=[250,250,-60],marker=1,worldMarker=True)

#Creating the boundary cells
domain = mt.createCube(start=[58,100,0], end=[192,150,-15],marker=2)

#Creating the 3D object (inhomogeneity)
inhom = mt.createCube(start=[125,120,-2], end=[130,124,-4],marker=3)

# Merging and visualizing
plc = world+domain+inhom
#fig,ax=plt.subplots()
pg.viewer.show(plc)
# pg.viewer.mpl.drawSensors(ax=ax,sensors=shm.sensors(), diam=1,
#                          facecolor='red', edgecolor='black')
#pg.viewer.show(plc,showNodes=True)

#Creating nodes in each electrode
for s in shm.sensors():
    plc.createNode(s,marker = -99)
# positioning reference electrode (should be labelled -999)
# plc.createNode([0,230,-11.25],marker=-999)

# positioning calibration electrode (should be labelled -1000)
# plc.createNode([250,250,-5], marker=-1000)

#Refining mesh around electrodes of the ERT (in the z direction)
for s in plc.positions(pg.find(plc.nodeMarkers() == -99)):
     plc.createNode(s-[0,0,0.1*sep_elec])


#Creating mesh
mesh = mt.createMesh(plc,smooth=True,area=25)

pg.show(mesh,mesh.cellMarkers())
# pg.show(mesh,mesh.cellMarkers(),filter={'clip':{'origin':(127, 0, 0.0)},})

# Populating cells at each marker domain with resistivity values
rhomap = [[1, 100.0], [2, 100.0],[3, 1.0]]
rhomap_hom = [[1, 100.0], [2, 100.0],[3, 100.0]]


#%% Simulating close 3D-block
inho = ert.simulate(mesh, res=rhomap, scheme=shm, sr=False,
                   calcOnly=True, verbose=True)

inho.save('inhomogeneous.ohm', 'a b m n u')
inho.set('rhoa', data_k * inho('u') / inho('i'))
inho.set('k',data_k)
pg.show(inho)
plt.title('Inhomogeneity')

#Simulating homogeneous
homo = ert.simulate(mesh, res=rhomap_hom, scheme=shm, sr=False,
                   calcOnly=True, verbose=True)

homo.save('homogeneous.ohm', 'a b m n u')
homo.set('rhoa', data_k * homo('u') / homo('i'))
homo.set('k',data_k)
pg.show(homo)
plt.title('homogeneous world')

PS.: When I try to use the keyword 'filter' using 'pg.show', I got the following error. I just would like to make a slice of the 3-D model like the one in the tutorial section to verify my geometry. Is it better to check on Paraview ?

  File "C:\Users\cavallui\Anaconda3\envs\pg\lib\site-packages\pyvista\plotting\plotting.py", line 3518, in __del__
    if not self._closed:
AttributeError: 'Plotter' object has no attribute '_closed'
03/11/22 - 15:56:59 - pyGIMLi - ERROR - None.showMesh3DVista(C:\Users\cavallui\Anaconda3\envs\pg\lib\site-packages\pygimli\viewer\vistaview.py:128)
fix pyvista bindings
__init__() got an unexpected keyword argument 'filter'

[halbmy]

1. Well, if you're working in an unbounded domain, there is no reason for total field computations, so you should use ert.simulate(..., sr=True) (which is the default). This makes the computations way more accurate. If insisting on total field, you would have to use quadratic shape functions (actually it would be good to have a keyword argument to do that on the fly).
2. You can directly put the geometric factors into the scheme by shm["k"] = ert.createGeometricFactors(shm) and simplify the call.
3. Computing homogeneous resistivity can also be done by passing a scalar, e.g. ert.simulate(mesh, res=1.0, scheme=shm).
4. Obviously, mt.createWorld decides on start and end arguments (and not analysing the z values) which face obtains the Neumann boundary condition marker, so better use world = mt.createWorld(start=[0, 0, -60], end=[250, 250, 0]) instead. This is also influencing the boundary conditions and could cause error messages.

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
import matplotlib.pyplot as plt

# ERT electrode configuration
sep_elec = 2
xe = np.arange(78, 172.1, sep_elec)
n_elec = len(xe)
sensors = np.zeros((n_elec, 3))
sensors[:, 0] = xe
sensors[:, 1] = 125
# Creating the acquisition scheme in wenner alpha (wa)
shm = ert.createData(elecs=sensors, schemeName='wa')
# Calculating geometric factors
shm["k"] = ert.createGeometricFactors(shm)
# Creating the PLC world
world = mt.createWorld(start=[0, 0, -60], end=[250, 250, 0],
                       marker=1, worldMarker=True)
# Creating the boundary cells
domain = mt.createCube(start=[58, 100, -15], end=[192, 150, 0], marker=2)
# Creating the 3D object (inhomogeneity)
inhom = mt.createCube(start=[125, 120, -4], end=[130, 124, -2], marker=3)
# Merging and visualizing
plc = world + domain + inhom
# Creating nodes in each electrode
for s in shm.sensors():
    plc.createNode(s, marker=-99)
# Refining mesh around electrodes of the ERT (in the z direction)
for s in plc.positions(pg.find(plc.nodeMarkers() == -99)):
    plc.createNode(s - [0, 0, 0.1*sep_elec])
# Creating mesh
mesh = mt.createMesh(plc, smooth=True, area=25)
# Simulating homogeneous
homo = ert.simulate(mesh, res=100.0, scheme=shm, verbose=True)
pg.show(homo)
#
rhomap = [[1, 100.0], [2, 100.0], [3, 1.0]]
inho = ert.simulate(mesh, res=rhomap, scheme=shm, verbose=True)
pg.show(inho)

The homogeneous case of course gains constant 100 Ohmm and the inhomogeneous is as expected:

[LuFrag]

Thanks for your detailed answer/recommandations :) Works pretty fine.

I confirm the error messages concerning the boundary conditions of my first try. They were badly defined by my createWorld start and end arguments.

1. To avoid this problem in the future, is there a way to verify the boundary conditions (specially the Neumann one) of a 3-D model using PyGimli and visualize slices from it? I managed to export the mesh for a visualisation in Paraview, but I can only see the 'cellMarkers' values.

2. To consider the effect of the 3-D anomaly on an 2-D ERT inversion scheme, I added the following lines to the code. I reorganized the electrode positions, zeroying the previous y coordinates. I will be able to evaluate the effects of a close 3-D anomaly like this (hopefully).

Thanks again for your help.

#%% Inverting considering the Earth as a 2-D model instead of the original 3-D

ERT_mgr = pg.physics.ert.ERTManager(homo)

homo['err'] = ert.estimateError(homo, relativeError=0.05)
inho['err'] = ert.estimateError(inho, relativeError=0.05)

# Retrieving the 3-D position of the electrodes
pos =np.array([homo.sensors()])

# Switching y and z on the sensors positions for a 2-D inversion
for i in range(homo.sensorCount()):
    homo.setSensorPosition(i,[pos[0,i,0], 0.0 ,0.0])
    inho.setSensorPosition(i,[pos[0,i,0], 0.0 ,0.0])

# Verifying the switched electrodes
pos =np.array([homo.sensors()])
print(pos)

#%% Inversion
mod_ho = ERT_mgr.invert(data=homo,lam=5,verbose=True,maxIter=5)

ERT_mgr.showResultAndFit()

mod_in = ERT_mgr.invert(data=inho,lam=5,verbose=True,maxIter=5)

ERT_mgr.showResultAndFit()

Here is the result for the homogeneous case (I could improve the visualisation):

And for the inhomogenous 3-D anomaly close to the ERT:

[halbmy]

1. for exporting boundary conditions you can use mesh.exportBoundaryVTU("boundary.vtu")
2. inverting a homogeneous model does not really make sense. For the inhomogeneous model, I strongly suggest adding noise.

[LuFrag]

Thanks for your answer and recommandations !

I visualized the boundary conditions of the 3-D mesh and yes, adding noise is the next step in my study.