Description
Problem description
I installed the latest version of pygimly in a fresh environment (named pg15+ on my machine) using the command given on the website . My issue is that when I try to do a topographic inversion of ERT field measure there is a Fatal Python error: Segmentation fault.
I've run the same code on windows and the inversion just worked perfectly fine, i've also created a new environment with pygimly 1.2.1 and the inversion was able to be processed easily on my linux machine.
Your environment
Date: mar. déc. 17 21:00:40 2024 CET
OS : Linux (Ubuntu 24.04)
CPU(s) : 12
Machine : x86_64
Architecture : 64bit
RAM : 30.7 GiB
Environment : IPython
File system : ext4
Python 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:36:13) [GCC
12.3.0]
pygimli : 1.5.3
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.9.1
scipy : 1.14.1
tqdm : 4.67.1
IPython : 8.30.0
pyvista : 0.44.2
Steps to reproduce
So here is the code with the latest set of data I've tried (I've also done it with data formatted for pygimly and some data using the command pg.getExampleData)
r"""
ERT field data with topography
==============================
Simple example of data measured over a slagdump demonstrating:
- 2D inversion with topography
- geometric factor generation
- topography effect
"""
# sphinx_gallery_thumbnail_number = 5
import pygimli as pg
from pygimli.physics import ert
import pygimli.meshtools as mt
from pygimli.frameworks.resolution import resolutionMatrix
###############################################################################
# Get some example data with, typically by a call like
data = ert.load('/home/lena/Documents/Documents/Universite/S9/04_Near_surface_geophy/04_Field_trip_data/G1-22112024/20241119_ddp64.dat')
# that supports various file formats
print(data)
sensors=data.sensorPositions()
# %%
###############################################################################
# The data file does not contain geometric factors (token field 'k'),
# so we create them based on the given topography.
data['k'] = ert.createGeometricFactors(data, numerical=True)
data['rhoa'] = data['k'] * data['r']
###############################################################################
# We initialize the ERTManager for further steps and eventually inversion.
mgr = ert.ERTManager(sr=False)
###############################################################################
# It might be interesting to see the topography effect, i.e the ratio between
# the numerically computed geometry factor and the analytical formula
k0 = ert.createGeometricFactors(data)
ert.showData(data, vals=k0/data['k'], label='Topography effect')
###############################################################################
# The data container has no apparent resistivities (token field 'rhoa') yet.
# We can let the Manager fix this later for us (as we now have the 'k' field),
# or we do it manually.
mgr.checkData(data)
print(data)
# %%
###############################################################################
# The data container does not necessarily contain data errors data errors
# (token field 'err'), requiring us to enter data errors. We can let the
# manager guess some defaults for us automaticly or set them manually
#data['err'] = ert.estimateError(data, relativeError=0.03, absoluteUError=5e-5)
# or manually:
# data['err'] = data_errors # somehow
ert.show(data, data['rhoa'])
# %%
plc = mt.createParaMeshPLC(data,
paraDX=0.25, #Résolution entre capteur (nb triangles entre capteur)
paraMaxCellSize=15, #Taille max (surface) du triangle
surfaceMeshQuality=33.6, #qualité du maillage
#surfaceMeshArea=2000,
paraBoundary=5, #de combien on déborde aux extrémités du profil (à gauche et à droite, en %)
#boudary=0.1,
#addTopo=topo,
paraDepth=30, #prof du profil
#addNodes=2, #ça raffine. Pas trop certain de l'utilité pour ces données-ci
)
for pos in sensors:
plc.createNode(pos)
mesh = mt.createMesh(plc)
pg.show(mesh)
###############################################################################
# Now the data have all necessary fields ('rhoa', 'err' and 'k') so we can run
# the inversion. The inversion mesh will be created with some optional values
# for the parametric mesh generation.
#
mod = mgr.invert(data,
#lam=10, # Facteur de régularisation
verbose=True,
mesh = mesh,
robust = True, #Norm L1 (False = L2) avec un fort contraste entre les couches L1 marche mieux car lisse moins
zWeight= 0.3, #Anisotropie (milieu stratifié < 1 ; variations latérales > 1)
)
mgr.showResult()
###############################################################################
# We can view the resulting model in the usual way.
mgr.showResultAndFit()
# np.testing.assert_approx_equal(ert.inv.chi2(), 1.10883, significant=3)
# %%
###############################################################################
# Or just plot the model only using your own options.
mgr.showResult(mod, cMin=10, cMax=1000, cMap="Spectral_r", logScale=True)
Expected behavior
Get all the result of the inversion ? Eg the figures
Actual behavior
The kernel just die :
%runfile /home/lena/Documents/Documents/Universite/S9/04_Near_surface_geophy/04_Field_trip_data/ERT_inversion.py --wdir
17/12/24 - 21:04:24 - pyGIMLi - INFO - could not read unified data format for ERT ... try res2dinv
17/12/24 - 21:04:24 - pyGIMLi - INFO - Create default mesh for geometric factor calculation.
ModellingBase::setMesh() copying new mesh ... Found datafile: 64 electrodes
Found: 64 node-electrodes
rMin = 1, rMax = 252
NGauLeg + NGauLag for inverse Fouriertransformation: 14 + 4
Found non-Neumann domain
0.0105294 s
FOP updating mesh dependencies ... 1.956e-06 s
/tmp/spyder-lena/tmpot2fkmkm: line 3: 42782 Segmentation fault (core dumped) /home/lena/anaconda3/envs/pg15+/bin/python -Xfrozen_modules=off -m spyder_kernels.console -f /home/lena/.local/share/jupyter/runtime/kernel-89c49f7d7589.json
The kernel died, restarting...
Fatal Python error: Segmentation fault
Main thread:
Thread 0x00007960c4337740 (most recent call first):
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/pygimli/physics/ert/ert.py", line 192 in simulate
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/pygimli/physics/ert/ert.py", line 385 in createGeometricFactors
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/pygimli/utils/cache.py", line 277 in wrapper
File "/home/lena/Documents/Documents/Universite/S9/04_Near_surface_geophy/04_Field_trip_data/ERT_inversion.py", line 34 in <module>
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/spyder_kernels/customize/utils.py", line 209 in exec_encapsulate_locals
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/spyder_kernels/customize/code_runner.py", line 506 in _exec_code
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/spyder_kernels/customize/code_runner.py", line 336 in _exec_file
File "/home/lena/anaconda3/envs/pg15+/lib/python3.11/site-packages/spyder_kernels/customize/code_runner.py", line 165 in runfile
I thank you in advance for you help ! I'd really like to be able to use this latest version ! :)
Best,
Lena