gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
344 stars 131 forks source link

Inversion manager ignoring set background. #651

Closed soehag closed 1 month ago

soehag commented 3 months ago

Problem description

Hello everyone! I have been trying to run some "simple" ERT inversions with a specific region of interest. To only invert for the specific region, I created a mesh (region of interest) with a triangle boundary around it with associated region markers. When using the ERTmanager and setting the regularisation via manager.inv.setRegularization, the set backgrounds etc. get ignored.

Your environment


Date: Fri Feb 16 14:43:27 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

## create meshes
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 = pg.meshtools.appendTriangleBoundary(mesh_inner, xbound=10, ybound = 2 , quality=32, 
                                                    area=1, isSubSurface=False, smooth=True, addNodes=50, 
                                                    marker = 2)
pg.show(final_mesh, markers=True, showMesh=True)

## create sensors and data
sensorPositions_x = np.linspace(5,15,11)
sensorPositions_top = np.vstack((sensorPositions_x, np.ones_like(sensorPositions_x)*4)).T
sensorPositions_bot = np.vstack((sensorPositions_x, np.ones_like(sensorPositions_x)*(-4))).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 background set to marker and additional regularisation constraints

manager = ert.Manager()
manager.setData(np.random.random)
manager.setMesh(final_mesh)

manager.inv.setRegularization(2, background=True)
manager.inv.setRegularization(1, cType=1)

result = manager.invert(
    data=syn_data,
    mesh=final_mesh,
    lam =1e4
)
pg.show(manager.paraDomain, result)
pg.show(final_mesh, res)

Expected behavior

The background should be with marker 2, i.e. the inner box.

Actual behavior

The background is overwritten by the manager to the lowest marker value.

In addition it would be a nice feature, to have flexibility with the startModel input (e.g. startModel can be define on the full grid even if parts is omitted as background, fixed to a certain value, or treated as homogenous) and model output (retrieve the model on the full mesh, even if background is omitted during inversion). Cheers and thank you!

halbmy commented 2 months ago

background is a property that every region has (there can be several or no background regions), so setting it for one region does not change it for the other.

By default, region 1 (the lowest in case of multiple regions) background is set to True, so you will have to deactivate this explicitly:

manager.inv.setRegularization(1, cType=1, background=False)
soehag commented 2 months ago

Hi @halbmy, thanky you for your answer - that perfectly explains the behaviour.

However I still think, that the intuitive behaviour is to not assign an additional background region, if there is a background region already set in the manager. Especially, if a specific regularisation is set for the region with the lowest marker, overwriting the background property is unintuitive in my opinion.

I have a further question, regarding the background marker. In my understanding, setting background=True should completely remove the cells with the specific marker from the inversion, i.e. that paraDomain. But manually setting background=True doesnt seem to achieve that besides for the marker with the lowest number.

Cheers and have a nice weekend, Hagen

halbmy commented 2 months ago

I agree that this default behaviour is not very intuitive, it just matches the default mesh generation done by the method managers. So we have to place it into other abstraction levels and have to better document this behaviour.

To your question: setting background=True for a specific region excludes it from inversion. In your case, the problem was passing the mesh after setting it

manager.setMesh(final_mesh)
manager.inv.setRegularization(2, background=True)
manager.inv.setRegularization(1, cType=1, background=False)
result = manager.invert(data=syn_data, mesh=final_mesh, lam=1e4)
manager.showResult()

which will override the regularization set before. If you remove the mesh argument from invert() all will be well. Except that the inversion fails, but this is due to a wrong geometry since your electrodes are outside of the inversion domain. image

soehag commented 2 months ago

Hi,

thank you for your answer! That makes sense - starting the inversion with an adjusted geometry (such that the electrodes are well into the domain) and not passing the mesh to manager.inv() still throws a message, that the region with marker 1 is set to the background. I assume, that this is a "ghost" message?

To add onto the documentation - it would also be great to add at this place, how proper (non homogeneous) starting models can be supplied, when the flag "single=True" is set.

Cheers and thank you Hagen