ERT 4-point sensitivities with heterogenous subspace

[RobertDormilon]

Dear pygimli Team,

in order to visualize the sensitivity distributions of a heterogenous subspace I very slightly modified the example code under https://www.pygimli.org/_examples_auto/3_ert/plot_05_ert_4_point_sensitivities.html . In the example code a homogenous halfspace is assumed. To make it heterogenous I added some layers in line 43 where I included the kwarg
layers = [-0.1,-1.1]
to the mt.createWorld() command. However this results in an error message at the end of the code, where a sensitivity array is attempted to be displayed on a mesh of a different cell number than the sensitivity array length.
Also, differently to the example code, I get the message, that the region with the lowest marker (0) is set to background. That seems to cause the forward operation to be carried out on a different mesh than the one that has been previously created.
I have tried a number of things to solve the issue, but couldn't find a solution.

Kind regards,
Robert

PS: Also, I am new to Github and couldn't find the form that should be filled out in order to post to your repository. I apologize for that.

Sensitivities_try4.txt

[halbmy]

There is this (sometimes helpful, sometimes nasty) default behaviour to set the lowest region to background.
You can avoid this to set all regions to background by

fop.setRegionProperties('*', background=False)

Then the length of the Jacobian rows match the mesh, however, the forward operator renumbers the model cells according to increasing region markers. You can retrieve the renumberd (inversion) mesh by fop.paraDomain that should now match the Jacobian matrix, so that fop.paraDomain.show(fop.jacobian().row(0)) should work.

[RobertDormilon]

Thank you very much for your help!
I included

fop.setRegionProperties('*', background=False)

which did make the number of cells in the mesh match the jacobian matrix. However, the renumbering of the mesh cells still seems to be an issue, because the result doesn't match the output from the example code. I retrieved the renumbered mesh with

fop.paraDomain

and displayed the cell ids, but they seem to be the same as the ones from the initially created mesh.

fop.paraDomain.show(fop.jacobian().row(0))

didn't work because fop.paraDomain has now show() method, so I tried

pg.show(fop.paraDomain,fop.jacobian().row(0))

instead, but the in the resulting plot it was apparent again, that the entries of the jacobian matrix were displayed over the wrong cell ids. Is there maybe a way to force the forward operator to not renumber the cells of the mesh?

I use pygimli version 1.5.2 on python 3.11

Sensitivities_try4.txt

[halbmy]

Yes, the mesh.show(x) way was added after v1.5.3 (and does the same like pg.show(mesh, x).

Indeed, there is something wrong with the numbering we still have to dig in why.

If the forward operator is initialized with different markers, region-wise renumbering will be done. However, you can create the resistivity vector from the markers and then reset the markers before initializing the forward operator:

res = mesh.populate("res", {1:100, 2:1000, 3:300})
fop = ert.ERTModelling()
fop.setData(data)
fop.setMesh(mesh)
fop.createJacobian(res)
mesh.show(fop.jacobian().row(1))