3D magnetics inversion can't run "inversion.run()" function

[liutingli]

Problem description

Jacobian matrix size was wrong. There are  information :
min/max(dweight) = 2.14119e+09/8.33918e+13
Building constraints matrix
constraint matrix of size(nBounds x nModel) 156 x 100
check Jacobian: wrong dimensions: (8000x2670) should be (8000x100)  force: 1
jacobian size invalid, forced recalc: 1
Calculating Jacobian matrix (forced=1)...... 4.4e-06 s
Traceback (most recent call last):
  File "\pygimli\Two_magnetic_model.py", line 58, in <module>
    invmodel = inv.run(data, absoluteError=absError)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "\miniconda3\envs\pg\Lib\site-packages\pygimli\frameworks\inversion.py", line 692, in run
    self.inv.start()
RuntimeError: ./core/src/inversion.cpp:95		double GIMLI::RInversion::getPhiD(const Vec&) const  getPhiD == nan

Your environment

            OS : Windows (10 10.0.22631 SP0 Multiprocessor Free)
        CPU(s) : 16
       Machine : AMD64
  Architecture : 64bit
   Environment : Python

Python 3.11.10 | packaged by conda-forge | (main, Sep 22 2024, 14:00:36)
[MSC v.1941 64 bit (AMD64)]
pygimli : 1.5.2
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.9.2
scipy : 1.14.1
tqdm : 4.66.5
pyvista : 0.44.1

Steps to reproduce

import numpy as np
import pygimli as pg
from pygimli.meshtools import readGmsh
from pygimli.viewer import pv
from pygimli.physics.gravimetry import MagneticsModelling
import matplotlib


susceptibility1 = 1
susceptibility2 = 2.5
matplotlib.use("Qt5Agg")
# In processing of save mesh.msh file（version 2-ASCII）,don't choose "save all elements" button.
mesh = readGmsh(r"t4", verbose=True)
magconductivity = [[99, susceptibility1]]
model = pg.solver.parseMapToCellArray(magconductivity, mesh)
mesh["magnetic conductivity(SI)"] = model
# pl, _ = pg.show(mesh, label="magnetic conductivity(SI)", showMesh=True, cMap='jet', alpha=0.4)
# _ = pl.show()

F, I, D = 50000e-9, 90, 0
H = F * np.cos(np.deg2rad(I))
X = H * np.cos(np.deg2rad(D))
Y = H * np.sin(np.deg2rad(D))
Z = F * np.sin(np.deg2rad(I))
igrf = [D, I, H, X, Y, Z, F]  

xx = np.arange(-200, 200, 10)
yy = np.arange(-200, 200, 10)
py, px = np.meshgrid(xx, yy)

px = px.ravel()
py = py.ravel()
points = np.column_stack((px, py, np.ones_like(px)*30))
nxx = xx.shape
nyy = yy.shape
# The forward operator
cmp = ['Bxy', 'Bxz', 'Byy', 'Byz', 'Bzz']
# 'Bxy', 'Bxz', 'Byy', 'Byz', 'Bzz'
fop = MagneticsModelling(mesh, points, cmp, igrf)
data = fop.response(mesh["magnetic conductivity(SI)"])

# depth weighting
bz = np.array([abs(b.center().z()) for b in mesh.boundaries() if not b.outside()])
z0 = 25
wz = 100 / (bz+z0)**1.5
print(min(wz), max(wz))

inv = pg.Inversion(fop=fop, verbose=True, stopAtChi1=True)

inv.setRegularization(limits=[0, 1.2])  # to limit values and add relative weight for vertical boundaries
# inv.setConstraintWeights(wz)
# relaError = 0.001 * np.random.randn(len(data))
absError = 0.002 * np.abs(data)
data += np.random.randn(len(data)) * absError
# relaError = absError / data
# data += absError
invmodel = inv.run(data, absoluteError=absError)
mesh["inv"] = invmodel

pl, _ = pg.show(mesh, label="magnetic conductivity(SI)", style="wireframe", hold=True,
                filter={"threshold": dict(value=0.25, scalars="magnetic conductivity(SI)")})
pv.drawMesh(pl, mesh, label="inv", style="surface", cMap="Spectral_r",
            filter={"threshold": dict(value=0.25, scalars="inv")})
pv.drawMesh(pl, mesh, label="inv", style="surface", cMap="Spectral_r",
            filter={"slice": dict(normal=[-1, 0, 0],origin=[0, 0, 80])})
pl.camera_position = "yz"
pl.camera.roll = 90
pl.camera.azimuth = 180 - 15
pl.camera.elevation = 10
pl.camera.zoom(1.2)
_ = pl.show()

magnetic_inversion.zip

[halbmy]

Some of your data are near-zero. If you consider a relative error transformed into an absolute one

absError = 0.002 * np.abs(data)

the minimum of this vector becomes practically zero (1e-27) and this makes the inversion crash. I strongly suggest adding a fixed absolute error (e.g. magnetometer resolution) and expect this to resolve the issue.

[liutingli]

Thank you very much, sir! I have tried to add a fixed absolute error before, but it doesn't work. I want to know if there is a problem with my initial grid settings in the inversion function: inv = pg.Inversion(fop=fop, verbose=True, stopAtChi1=True).

[halbmy]

No, that's ok.

I am pretty sure the error message is different if adding an absolute value to the error.