gimli-org / gimli

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

Problem With ERT Sensor in Slope Model #698

Closed tonapawilliam closed 6 months ago

tonapawilliam commented 6 months ago

Excuse me Sir, can you fix my code? I try to make model with slope but there is problem with simulate data ERT.

import numpy as np
import pygimli as pg
from pygimli.physics import ert
import pygimli.meshtools as mt
from pygimli.physics.ert.importData import importRes2dInv
from pygimli.physics.ert.importData import importAsciiColumns
layer1 = mt.createPolygon([[0, -20], [345, 0], [345, -10], [0, -30]],
                          isClosed=True, marker=1, area=1)
layer2 = mt.createPolygon([[0.0, -30], [345, -10], [345, -30], [300, -40], [250, -80], [200, -80], [50, -50], [0, -50]],
                          isClosed=True, marker=2)
layer3 = mt.createPolygon([[0.0, -50], [50, -50], [200, -80], [250, -80], [300, -40], [345, -30], [345, -90], [0, -90]],
                          isClosed=True, marker=3)
anomaly1 = mt.createCircle(pos=[150, -40], radius=[25,9], area=5, marker=4)
anomaly2 = mt.createCircle(pos=[225, -50], radius=[25, 15], marker=5)
anomaly3 = mt.createCircle(pos=[100,-45], radius=10, marker=6)
geom = layer1 + layer2 + layer3 + anomaly1 + anomaly2 + anomaly3
pg.show(geom)
scheme = ert.createData(elecs=np.linspace(start=0, stop=345, num=71), schemeName ='gr')
mesh = mt.createMesh(geom, quality = 34, smooth =[0,10], area = 0.5, verbose=True)
rhomap = [[1, 100], #limonate
          [2, 50],#limonite ore
          [3, 500], #Bedrock
          [4, 1000], #Boulder
          [5, 1000], #Boulder
          [6, 1000]] #Boulder
pg.show(mesh, data=rhomap, showMesh=False, cMap = "Spectral_r")
data = ert.simulate(mesh, scheme=scheme, res=rhomap)
data ['err'] = ert.estimateError(data)
ert.show(data)]

The Error:

Traceback (most recent call last):

  File ~\miniconda3\envs\pg\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File d:\skripsi\code\hasil\model syntetik vale\model_syntetic.py:57
    data = ert.simulate(mesh, scheme=scheme, res=rhomap)

  File ~\miniconda3\envs\pg\lib\site-packages\pygimli\physics\ert\ert.py:124 in simulate
    fop.setMesh(mesh, ignoreRegionManager=True)

  File ~\miniconda3\envs\pg\lib\site-packages\pygimli\frameworks\modelling.py:731 in setMesh
    self.setMeshPost(mesh)

  File ~\miniconda3\envs\pg\lib\site-packages\pygimli\physics\ert\ertModelling.py:211 in setMeshPost
    self._core.setMesh(mesh, ignoreRegionManager=True)

RuntimeError: There is a requested electrode that does not match the given mesh. 

./core/src/bert/dcfemmodelling.cpp:931      virtual void GIMLI::DCMultiElectrodeModelling::searchElectrodes_()  0   0   0 mesh BoundingBox [0   -90 0, 345  0   0]

Constructing Delaunay triangulation by divide-and-conquer method.
Delaunay milliseconds:  17
Recovering segments in Delaunay triangulation.
Segment milliseconds:  0
Removing unwanted triangles.
Spreading regional attributes and area constraints.
Hole milliseconds:  0
Adding Steiner points to enforce quality.
Quality milliseconds:  141

Writing vertices.
Writing triangles.
Writing segments.
Writing edges.

Output milliseconds:  8
Total running milliseconds:  169

Statistics:

  Input vertices: 48
  Input segments: 50
  Input holes: 0

  Mesh vertices: 55961
  Mesh triangles: 110834
  Mesh edges: 166794
  Mesh exterior boundary edges: 1086
  Mesh interior boundary edges: 1281
  Mesh subsegments (constrained edges): 2367

ModellingBase::setMesh() copying new mesh ... Found Neumann domain. Setting topography=1.
Found Neumann domain. but 2.5D -> neumann: false
andieie commented 6 months ago

It looks like your sensors are located at y=0. You can check your sensor positions by scheme.sensorPositions().array(). Since you are modelling a slope you should put the sensors along the slope. Also, the sensors are at the start and end, right at 0 and 345 which are the beginning and end of your mesh. I would put them more centrally in your mesh just to make sure it is not that. You can also add the sensors to your geom and then make your mesh. You can do this by :

for i in scheme.sensorPositions():
        geom.createNode(i)

Hope this helps.

tonapawilliam commented 6 months ago

How to put the value slope in my code? Can You Help me?

andieie commented 6 months ago

You can use this as an example: https://www.pygimli.org/_examples_auto/2_seismics/plot_01_refraction_manager.html

tonapawilliam commented 6 months ago

Thank You, It's work GRD12