TravelTimeManager crashes with strong topography

GregBievre commented 3 years ago

Hi pyGIMLi team,

Problem description

Seismic refraction inversion using the TravelTimeManager seems to crash when marked topography is present. When removing topography it works. When, increasing progressively the slope/topography, it works until reaching topographic slope values of around 35%.

I tried to reproduce the problem using the plot_04_koenigsee.py script provided on the website (instead of using my own script) associated with a linear slope along the profile and I eventually have the same issues: it crashes when reaching slopes of around 35% and higher.

Other profiles with lower slopes work very well (on this site or other sites). To specify that I think it might not be related to my pickings only (well, this was my first guess), I also emphasize that we conducted both P and S-wave measurements along the same profile and that they both crash with the same slope threshold. Finally, I also mention that all the process (with the same data) worked when I used the refraction manager with an older version of pyGIMLI.

I attach the data file sp01p2. sp01p2.zip

Your environment

Operating system: Windows 10 64 Python version: 3.7.8 pyGIMLi version: 1.1.0+31.g47c8d810 Way of installation: Conda

Steps to reproduce

import pygimli as pg from pygimli.physics import TravelTimeManager

data = pg.physics.traveltime.load('sp01p2.dat') print(data)

fig, ax = pg.plt.subplots() pg.physics.traveltime.drawFirstPicks(ax, data)

mgr = TravelTimeManager() model=mgr.invert(data, secNodes=3, verbose=1)

Expected behavior

It works with those data when the slope is lower than 35%

Actual behavior

The kernel crashes. I attach the data file with topography values producing the crash.

11/11/20 - 16:02:03 - pyGIMLi - INFO - Found 1 regions.
11/11/20 - 16:02:04 - pyGIMLi - ERROR - <class 'pygimli.physics.traveltime.TravelTimeManager.TravelTimeManager'>._ensureError(C:\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py:339)
DataContainer has no "err" values. Fallback set to 3%
11/11/20 - 16:02:04 - pyGIMLi - INFO - Creating forward mesh from region infos.
11/11/20 - 16:02:04 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
11/11/20 - 16:02:04 - pyGIMLi - INFO - Starting inversion.
11/11/20 - 16:02:04 - pyGIMLi - INFO - Create gradient starting model. 500: 5000

Kernel died, restarting

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

halbmy commented 3 years ago

To create a mesh, internally the function

plc = pg.meshtools.createParaMesh(data, paraDepth=depth, boundary=0)

is called. If you use depth=50, you end up in Apparently, the lower model boundary is a horizontal line at the depth below the first sensor and not (like in BERT) a sloped line from 50m below the first sensor to 50m below the last sensor. We should change this accordingly.

By default, the depth is determined by the maximum shot-geophone offset divided by 3. If the slow is too high, it cuts the topography and produce an invalid mesh.

halbmy commented 3 years ago

Actually, createParaMesh has a parameter balanceDepth=True. If set to False, it produces a nicer mesh: Maybe this should be False by default. In my opinion its interpretation is a bit confusing (does balance mean to introduce the slope or to avoid the slope?).

Anyway, you can pass this argument to the invert call:

mgr.invert(secNodes=3, verbose=1, balanceDepth=False)
mgr.showResult()

GregBievre commented 3 years ago

Wow, It works fine. Many thanks Thomas for your efficience and your quick answer! It actually did not come to me to investigate the createParaMesh parameters since I did not explicitely write it in the script (and despite it is actually used...). Best regards, Greg

gimli-org / gimli