gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
344 stars 131 forks source link

3D ERT cylinder meshing #702

Open HenryWHR opened 3 weeks ago

HenryWHR commented 3 weeks ago

Hello,

I am trying to simulate ERT measurements on a cylinder with two rings of electrodes. Based on the example '3D modelling in a closed geometry', I built the code below. However, it (the meshing) doesn't work unless I remove the nodes at electrode position.

If I remove the nodes at the electrode position, but keep nodes inserted 1 mm for refinement and consider these nodes as electrodes, the meshing and modelling can work.

I feel like the electrodes on the boundary are not in the domain of the cylinder, but I don't know why and how to solve this. Could you please have a look and help?

Many thanks!

Best regards Henry

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

# example configuration
pole = np.array((
    [0, 1, 2, 3],
    [0, 1, 3, 4], 
    [0, 1, 4, 5],
    [0, 1, 5, 6],
    [0, 1, 6, 7]))
shm = pg.DataContainerERT()

for i, elec in enumerate("abmn"):
    shm[elec] = pole[:,i]

r_container = 0.3 # radius of the container
h_container = 0.4 # height of the container

# electrode coordinates
# and also the points for refinement
N_elec_ring = 8 # number of electrode per ring
co_elec_ring1 = np.zeros((N_elec_ring,3)) # upper ring
co_elec_ring2 = np.zeros((N_elec_ring,3)) # lower ring
co_elec_ring1_rf = np.zeros((N_elec_ring,3)) # nodes for refinement
co_elec_ring2_rf = np.zeros((N_elec_ring,3)) # nodes for refinement

pi = np.pi
# additional nodes 1mm inside the cylinder to refine the grids around electrode
for n in range(N_elec_ring):
    co_elec_ring1[n,0] = (r_container)*np.sin(2*n*pi/N_elec_ring)
    co_elec_ring1[n,1] = (r_container)*np.cos(2*n*pi/N_elec_ring)
    co_elec_ring1[n,2] = 0.1 # position of the upper ring above 0
    co_elec_ring2[n,0] = (r_container)*np.sin(2*n*pi/N_elec_ring)
    co_elec_ring2[n,1] = (r_container)*np.cos(2*n*pi/N_elec_ring)
    co_elec_ring2[n,2] = -0.1 # position of the lower ring below 0
    co_elec_ring1_rf[n,0] = (r_container-1e-3)*np.sin(2*n*pi/N_elec_ring)
    co_elec_ring1_rf[n,1] = (r_container-1e-3)*np.cos(2*n*pi/N_elec_ring)
    co_elec_ring1_rf[n,2] = 0.1 # position of the upper ring above 0
    co_elec_ring2_rf[n,0] = (r_container-1e-3)*np.sin(2*n*pi/N_elec_ring)
    co_elec_ring2_rf[n,1] = (r_container-1e-3)*np.cos(2*n*pi/N_elec_ring)
    co_elec_ring2_rf[n,2] = -0.1 # position of the lower ring below 0

co_elec = np.concatenate((co_elec_ring1, co_elec_ring2), axis=0)
co_elec_rf = np.concatenate((co_elec_ring1_rf, co_elec_ring2_rf), axis=0)
# assign electrode coordinates in the scheme
shm.setSensorPositions(co_elec)

# create a Cylinder with center at [0,0,0]
plc = mt.createCylinder(radius=r_container, height=h_container,nSegments=32,boundaryMarker=1)

# creat nodes at electrode position and the additional refinement points
for s in co_elec:
    plc.createNode(s, marker=-99)
for s in co_elec_rf:
    plc.createNode(s)
# reference pole for dipole current injection
plc.createNode([0, 0, 0.15], marker=-999) # at the top
plc.createNode([0, 0, 0.15 - 1e-3])       # refinement
# calibration point for voltages
plc.createNode([0, 0, 0.2], marker=-1000) # at the top
# pg.show(plc, alpha=0.3)

mesh = mt.createMesh(plc)
# pg.show(mesh, markers=True, showMesh=True)
# pg.show(mesh, mesh.cellMarkers(), showMesh=True, filter={'clip':{'origin':(0, 0, 0)},})

hom = ert.simulate(mesh, res=1.0, scheme=shm, sr=False,
                    calcOnly=True, verbose=True)
HenryWHR commented 2 weeks ago

Update: Although I still don't know what is wrong with the previous code, the problem is partly solved now by removing the following lines from the code:

for s in co_elec:
    plc.createNode(s, marker=-99)
halbmy commented 2 weeks ago

The electrodes should really be added to have accurate results and they need to be inside of the modelling domain. Maybe this is a tolerance problem.

Often, the electrodes are not quite points but spheres or cylinders and it makes sense to put them a little bit (the dimension of the electrode) inside, i.e. using r_container-dx for the sin/cos terms. Maybe this helps.

HenryWHR commented 1 week ago

Thanks! The forward modelling is clear to me and it works well by inserting the electrodes into the mesh. r_container-dx helps.

When it comes to inversion, do you have any suggestions on how to set up the mesh? If I use the same idea, perhaps it's too fine for inversion grids around the electrodes.