Description
Problem description
Hi everyone,
I am running ERT-inversions with custom constraint matrices to account for the specific geometry in my setup. Running these inversions on the my local HPC, I run into segmentation fault errors when the constraint matrix is being build. I have attached a small minimum working example, where an identity matrix as constraint matrix crashes the jupyter kernel before it can properly throw an exception. Running as a script I encounter the same error in a mwe on my local machine.
Your environment
Date: Wed Feb 28 11:09:28 2024 CET
OS : Linux
CPU(s) : 64
Machine : x86_64
Architecture : 64bit
RAM : 251.5 GiB
Environment : Jupyter
File system : ext4
Python 3.10.13 | packaged by conda-forge | (main, Oct 26 2023, 18:07:37)
[GCC 12.3.0]
pygimli : 1.4.4
pgcore : 1.4.0
numpy : 1.24.4
matplotlib : 3.7.2
scipy : 1.11.3
IPython : 8.17.2
pyvista : 0.42.3
Steps to reproduce
Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:
import pygimli as pg
import pygimli.meshtools as mt
import numpy as np
import pygimli.physics.ert as ert
mesh_inner_plc = pg.meshtools.createWorld(np.array([10,4]), np.array((20, -4)), worldMarker=True, marker=1)
mesh_inner_plc.setBoundaryMarkers([-2]*4)
mesh_inner = pg.meshtools.createMesh(mesh_inner_plc, quality=32, area=.1, smooth=True)
final_mesh=mesh_inner
pg.show(final_mesh, markers=True, showMesh=True)
### Create data
sensorPositions_x = np.linspace(12,18,11)
sensorPositions_top = np.vstack((sensorPositions_x, np.ones_like(sensorPositions_x)*2)).T
sensorPositions_bot = np.vstack((sensorPositions_x, np.ones_like(sensorPositions_x)*(-2))).T
scheme = pg.DataContainerERT()
for pos in sensorPositions_top:
scheme.createSensor(pos)
for pos in sensorPositions_bot:
scheme.createSensor(pos)
for sen_top in range(len(sensorPositions_top)):
for sen_bot in range(len(sensorPositions_top),len(sensorPositions_top)+len(sensorPositions_bot)):
scheme.addFourPointData(sen_top,-1,sen_bot,-1)
scheme["k"] +=1
res = np.array([h.center().y() for h in final_mesh.cells()])**2+100
pg.show(final_mesh, res)
syn_data = ert.simulate(mesh=final_mesh, scheme=scheme, res=res, noiseLevel=1)
pg.show(final_mesh, markers=True, showMesh=True)
### Invert with ERT manager
sM=100
manager = ert.Manager()
manager.setData(syn_data)
manager.setMesh(final_mesh)
C=pg.Matrix(np.eye(len(final_mesh.cells())))
weights = pg.Vector(np.ones(len(final_mesh.cells())))
weights = np.ones(len(final_mesh.cells()))
# weights = 1
manager.inv.setRegularization(C=C)
manager.inv.setConstraintWeights(weights)
result = manager.invert(
data=syn_data,
mesh=final_mesh,
lam =1e2,
startModel=sM,
verbose=True
)
Expected behavior
Run properly.
Actual behavior
...
Building constraints matrix
Segmentation fault (core dumped)
ModellingBase::setMesh() copying new mesh ... Found datafile: 22 electrodes
Found: 22 free-electrodes
rMin = 0.3, rMax = 14.4222
NGauLeg + NGauLag for inverse Fouriertransformation: 10 + 4
Found non-Neumann domain
0.00672143 s
FOP updating mesh dependencies ... 3.85e-06 s
Calculating response for model: min = 100 max = 115.381
Allocating memory for primary potential...... 0.000443408
No primary potential for secondary field calculation. Calculating analytically...
Forward: time: 0.0617532s
Response: min = 1.22867 max = 2.19441 mean = 1.89735
Reciprocity rms(modelReciprocity) 0.0936967%, max: 0.240693%
relativeError set to a value > 0.5 .. assuming this is a percentage Error level dividing them by 100
Data error estimate (min:max) 0.010045556265465905 : 0.010081375519050945
28/02/24 - 11:11:12 - pyGIMLi - INFO - Found 1 regions.
28/02/24 - 11:11:12 - pyGIMLi - INFO - Found 1 regions.
28/02/24 - 11:11:12 - pyGIMLi - INFO - Creating forward mesh from region infos.
28/02/24 - 11:11:12 - pyGIMLi - INFO - Creating refined mesh (H2) to solve forward task.
28/02/24 - 11:11:12 - pyGIMLi - INFO - Mesh for forward task: Mesh: Nodes: 2740 Cells: 5300 Boundaries: 4064
min/max(dweight) = 99.1928/99.5465
28/02/24 - 11:11:12 - pyGIMLi - INFO - Starting inversion.
fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x7fa2ee027560>
Data transformation: <pgcore._pygimli_.RTransLogLU object at 0x7fa2e549f100>
Model transformation: <pgcore._pygimli_.RTransLog object at 0x7fa2e549e5c0>
min/max (data): 1.23/2.24
min/max (error): 1%/1.01%
min/max (start model): 100/100
--------------------------------------------------------------------------------
Calculating response for model: min = 100 max = 100
Allocating memory for primary potential...... 0.00057357
No primary potential for secondary field calculation. Calculating analytically...
Forward: time: 0.251483s
Response: min = 1.10435 max = 1.99086 mean = 1.71883
Reciprocity rms(modelReciprocity) 0.0226251%, max: 0.0574602%
min/max(dweight) = 99.1928/99.5465
Building constraints matrix
Segmentation fault (core dumped)
bug