ERROR - There are cells with marker -1. Did you define a boundary region? (This is not needed).

[makeabhishek]

Traveltime tomography

I'm trying to invert a model with circular geometry. I have created a mesh with two regions and assign the markers to each. However, while running the inversion it is giving the error Core - ERROR - There are cells with marker -1. Did you define a boundary region? (This is not needed).

Operating system: Linux Python version: Python 3.9.16 pyGIMLi version: 1.3.1

Steps to reproduce

import matplotlib.pyplot as plt
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import TravelTimeManager
import pygimli.physics.traveltime as tt
import math
mr = TravelTimeManager()

c0 = mt.createCircle(pos=(0,0), radius=10, nSegments=64, marker=1)
c1 = mt.createCircle(pos=(0, 0), radius=5, nSegments=64, marker=2)
geom = c0 + c1
mesh = mt.createMesh(geom, quality=34, area=2, smooth=[1, 10])
ax, cb = pg.show(geom, markers=True)

fig = ax.get_figure()
for nr, marker in enumerate(geom.regionMarkers()):
    print('Position marker number {}:'.format(nr + 1), marker.x(), marker.y(),
          marker.z())
    ax.scatter(marker.x(), marker.y(), s=(4 - nr) * 20, color='k')
ax.set_title('marker positions - working example')
fig.show()

pg.show(mesh, markers=True, showMesh=True);

numberSources = 16
sensors = np.linspace(0., 360, numberSources) # x position
scheme = pg.physics.traveltime.createRAData(sensors) 

# Sensor positions
pos = np.zeros((numberSources, 3))
pos = np.array(scheme.sensors())
for x in (np.arange(numberSources)):
    angle = 360/numberSources*x
    pos[x,0] = math.sin(angle*math.pi/180)*9.99
    pos[x,1] = math.cos(angle*math.pi/180)*9.99
scheme.setSensors(pos)

mgr = TravelTimeManager()
vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 1000
vp[vp == 2] = 3000

ax, _ = pg.show(mesh, vp, colorBar=True, logScale=False, label='v [m/s]')
pg.viewer.mpl.drawSensors(ax, scheme.sensors(), diam=0.5,
                         facecolor='red', edgecolor='black')

ttData = mgr.simulate(slowness=1.0 / vp, scheme=scheme, mesh=mesh,
                    noiseLevel=0.001, noiseAbs=0.001, seed=1337,
                    verbose=True)

# Check at the traveltime data matrix: show first arrival traveltime curves
mgr.showData(ttData)

TT = tt.TravelTimeManager(verbose=False)
velInv = TT.invert(ttData, mesh=mesh, lam=100, useGradient=0, zWeight=1.0)
TT.inv.echoStatus()

Error

*** /home/carsten/miniconda3/envs/build/conda-bld/pgcore_1661190390662/work/gimli/core/src/ttdijkstramodelling.cpp:500  

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], line 8
      1 #velInv = mgr.invert(ttdata, useGradient=True, secNodes=2, paraMaxCellSize=15.0,
      2 #                  maxIter=10, verbose=0)
      3 #np.testing.assert_array_less(mgr.inv.inv.chi2(), 1.1)
      7 TT = tt.TravelTimeManager(verbose=False)
----> 8 velInv = TT.invert(ttData, mesh=mesh, lam=100, useGradient=0, zWeight=1.0)
      9 TT.inv.echoStatus()

File ~/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/physics/traveltime/TravelTimeManager.py:247, in TravelTimeManager.invert(self, data, useGradient, vTop, vBottom, secNodes, **kwargs)
    244 else:
    245     self.fop._useGradient = None
--> 247 slowness = super().invert(data, mesh, **kwargs)
    248 velocity = 1.0 / slowness
    249 self.fw.model = velocity

File ~/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/frameworks/methodManager.py:793, in MeshMethodManager.invert(self, data, mesh, startModel, **kwargs)
    790     self.fw.blockyModel = kwargs["blockyModel"]
    792 self.preRun(**kwargs)
--> 793 self.fw.run(dataVals, errorVals, **kwargs)
    794 self.postRun(**kwargs)
    795 return self.paraModel(self.fw.model)

File ~/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/frameworks/inversion.py:587, in Inversion.run(self, dataVals, errorVals, **kwargs)
    584 if self._preStep and callable(self._preStep):
    585     self._preStep(0, self)
--> 587 self.inv.start()
    588 self.maxIter = maxIterTmp
    590 if self._postStep and callable(self._postStep):

File ~/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/physics/traveltime/modelling.py:84, in TravelTimeDijkstraModelling.createJacobian(self, par)
     81 if not self.mesh():
     82     pg.critical("no mesh")
---> 84 return self._core.createJacobian(par)

RuntimeError: /home/carsten/miniconda3/envs/build/conda-bld/pgcore_1661190390662/work/gimli/core/src/sparsematrix.h:342     GIMLI::SparseMapMatrix<ValueType, IndexType>::MatElement GIMLI::SparseMapMatrix<ValueType, IndexType>::Aux::operator[](IndexType) [with ValueType = double; IndexType = long unsigned int; GIMLI::SparseMapMatrix<ValueType, IndexType>::MatElement = GIMLI::MatrixElement<double, long unsigned int, std::map<std::pair<long unsigned int, long unsigned int>, double, std::less<std::pair<long unsigned int, long unsigned int> >, std::allocator<std::pair<const std::pair<long unsigned int, long unsigned int>, double> > > >]  idx = 18446744073709551615, 0 maxcol = 412 stype: 0

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".

[makeabhishek]

Regarding the above issue, if I create a parametric mesh, it eliminate the error. Below is the set of lines I added after generating the mesh. But I wondering what is the function parametric mesh?

# Create the parametric mesh that only reflect the domain geometry
paraMesh = mt.createMesh(c0, q=34.0, smooth=[1, 10])
paraMesh.scale(10.0)

Furthermore, the inverted results are not as close if I follow the petrophysical inversion example.

May I seek your advise?

Regards AB

[halbmy]

Using the same mesh for forward calculation for inversion is considered an "inverse crime" in general. In your special case, you are already providing the information about the velocity distribution in the two given regions, so an inversion does not make much sense. You should rather create a new mesh for inversion as you already do in your second post.

For some methods (like ERT) you need a background region to account for boundary conditions. If more than one region is used for inversion, the one with the lowest marker is used as background, which is indeed a very bad default for traveltime tomography that should be avoided by default, at least for traveltime. The paraMesh returns the mesh, where inversion is done.

[makeabhishek]

Yeah, that make sense. Thanks for the explanation.
However, when I'm trying to use field data with 16 # shot/geophone points. It is giving me same error. Although, as I followed the field data example and read that the mesh is generated by coordinates. Here is the error I'm getting,

- pyGIMLi - ERROR - <class 'pygimli.physics.traveltime.TravelTimeManager.TravelTimeManager'>.checkError(TravelTimeManager.py:93)
DataContainer has no "err" values. Fallback to 3%
18/04/23 - 11:59:22 - pyGIMLi - INFO - Create constant starting model: 347.45950000000005
18/04/23 - 11:59:22 - pyGIMLi - INFO - Created startmodel from forward operator: 35 [347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005, 347.45950000000005]
18/04/23 - 11:59:22 - pyGIMLi - INFO - Starting inversion.
18/04/23 - 11:59:22 - Core - ERROR - There are cells with marker -1. Did you define a boundary region? (This is not needed).

As, I'm considering circular geometry, do i need to define the geometry separately? Furthermore, How can I create a initial velocity model, as I found that pygimli is considering arbitrary initial model.

[halbmy]

What does pg.show(mesh, markers=True) show? Are there cells with marker -1?

[makeabhishek]

/I have not created any mesh for field data. It was mentioned in the petrophysics example that "Finally, we call the invert method and plot the result.The mesh is created based on the sensor positions on-the-fly."

However, when i create a mesh its generating a non-close circular geometry. I'm expecting a circular geometry.

import pygimli as pg
from pygimli.physics import TravelTimeManager
import pygimli.physics.traveltime as tt
#data = tt.load("srtfieldline1.dat")
data = pg.load("traveltimeTestData.sgt")
mgr = TravelTimeManager(data)
mesh=mgr.createMesh(data,quality=34)
sm = tt.createGradientModel2D(data, mesh, vTop=200,vBot=1000)
pg.show(mesh, data=1/sm, label=pg.unit('vel'), showMesh=True, markers=True)

Here are the co-ordinates in .sgt file ``

[makeabhishek]

I also realized that for circular geometry, the coordinate sequence is changing the mesh.

[halbmy]

mt.createMesh is intended to work for surface sensors but not circular geometries. Please have a look at the examples with circular geometry (e.g. petrophysical joint inversion) and use mt.createPolygon(data.sensors(), isClosed=True) instead.

[makeabhishek]

Thanks for your advise. May I also know, how to create an initial model with a given velocity for inversion of field data? Is this the correct way to initialise the starting model?

TT.invert(data, mesh=paraMesh, lam=100, useGradient=0, zWeight=1.0, secNodes=2,
                  startModel=1500)

[halbmy]

For traveltime tomography in a circular geometry, you can start with a homogeneous velocity derived from the mean value of apparent velocities. For seismic refraction on the Earths surface, one usually starts with a gradient model, for which upper and lower velocity can be specified in the traveltime manager.

[makeabhishek]

So to start with a homogeneous velocity, do I need to define the apparent velocity explicitly or it’s automatically obtain from data? If need to define explicitly, is this the correct code line

TT.invert(data, mesh=paraMesh, lam=100, useGradient=0, zWeight=1.0, secNodes=2,
                  startModel=1500)

[halbmy]

Yes, the starting model can be specified by TT.invert(..., startModel= or by TT.startModel=, however, you will have to specify slowness instead of velocity, i.e. 1/1500.

I don't understand why you deleted a significant part of your comment, the one with the image enclosed.

A main purpose of these issues is to help other users solving their problems by finding solutions of previous problems, thus avoiding to repeat ourselves. Please keep this in mind.

[makeabhishek]

Actually the image was not showing in my computer, I thought there is a glitch in uploading. So I removed that.

[halbmy]

There was a problem with backquotes, but I had edited your comment to make it visible. Klick on "Edited" button to see the versions.