+
Skip to content

Segmentation Fault with custom Constraint Matrix #659

Closed
@soehag

Description

@soehag

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

Metadata

Metadata

Labels

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions

    点击 这是indexloc提供的php浏览器服务,不要输入任何密码和下载