gimli-org / gimli

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

Setting up regions of fixed slowness for travel-time inversion #487

Open Whitt1985 opened 1 year ago

Whitt1985 commented 1 year ago

Problem description

Aiming to use travel time tomography in a small scale environment. Problem is that we have two borehole liners that act as a conduit for the signal and fixing these known velocity regions will allow the problem to converge to a correct model.

Your environment

Operating system: Windows Python version: 3.8 pyGIMLi version: 1.2.5 Way of installation: e.g. Conda package, manual compilation from source, etc.

Steps to reproduce

_```
_# -*- coding: utf-8 -*-
"""
Created on Tue Oct  5 13:44:27 2021

@author: james.whittaker
"""

import matplotlib.pyplot as plt
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt

data = tt.load("survey_tom.TOM")

from pygimli.viewer.mpl import showVecMatrix

def showCrossholeData(data, vals, **kwargs):
    d = -pg.y(data)  # sensor depth
    ds = d[data["s"]]  # shot depth
    dg = d[data["g"]]  # geophone depth
    ax, cb = showVecMatrix(ds, dg, vals, label="apparent velocity (m/s)", **kwargs);
    ax.set_xlabel("shot depth");
    ax.set_ylabel("geophone depth");

showCrossholeData(data, data["s"])

left = 0
right = 3100
depth = -1000

rect1 = mt.createRectangle(
    start=[left, 0.0],
    end=[right, depth],
    isClosed=True,
    marker=1,
    markerPosition=[500,-400],
)

rect3 = mt.createRectangle(
    start=[0.0, 0.0],
    end=[right, depth-100],
    isClosed=True,
    marker=1,
    markerPosition=[500,-1100],
)

rect2 = mt.createRectangle(
    start=[0.0, 0.0],
    end=[right, -100],
    isClosed=True,
    marker=2,
)

plc = rect1 + rect2 + rect3

ax, cb = pg.show(plc, markers=True)

fig = ax.get_figure()
for nr, marker in enumerate(plc.regionMarkers()):
    print('Position marker number {}:'.format(nr + 1), marker.x(), marker.y(),
          marker.z())
    ax.scatter(marker.x(), marker.y(), s=(4 - nr) * 20, color='k')
ax.set_title('marker positions - working example')
fig.show()

mesh = mt.createMesh(plc, quality=34.4)
for b in mesh.boundaries():
    if b.marker() == -1 and not b.outside():
        b.setMarker(2)

print(mesh)
pg.show(mesh, markers=True, showMesh=True);

mgr = tt.TravelTimeManager(data, verbose=True)
mgr.setMesh(mesh)

mgr.inv.setRegularization(0,fix=1/2300)
mgr.inv.setRegularization(2,fix=1/2300)
mgr.invert(data,useGradient=False,startModel=0.005, verbose=True, isReference=True)
cov = pg.Vector(mgr.model.size(), 1.0)
kw = dict(cMin=200, cMax=3000, logScale=False, cMap="Spectral_r", coverage=cov)
ax, cb = mgr.showResult(**kw)__

Output errors in console

Expected behavior

Aim is to have a mesh with the top and bottom regions fixed at a velocity of 2300 m/s leaving the data in between more able to reproduce a result of varying velocity

Actual behavior

Inversion fails and doesn't run.

Reset region parameter
RegionManager copying mesh ...0.003 s
create NeighborInfos ... 0.001 s
analysing mesh ... 3 regions.
creating para domain ... 0.002 s
Reset region parameter
RegionManager copying mesh ...0.003 s
create NeighborInfos ... 0 s
analysing mesh ... 3 regions.
creating para domain ... 0.002 s
Reset region parameter
RegionManager copying mesh ...0.002 s
create NeighborInfos ... 0.001 s
analysing mesh ... 3 regions.
creating para domain ... 0.002 s
creating para domain ... 0.002 s
Warning: Region Nr: 0  is background and should not get a model transformation.
Warning: Region Nr: 0  is background and should not get a model control.
ModellingBase::setMesh() copying new mesh ... 0.002 s
FOP updating mesh dependencies ... 0.001 s
min/max(dweight) = 20.3445/46.2534
*** 0 C:/msys64/home/halbm/gimli/gimli/core/src/modellingbase.cpp:465
Warning: ** TMP HACk ** fixing region:  0  to:  17.2414
min/max(dweight) = 20.3445/46.2534
Building constraints matrix
constraint matrix of size(nBounds x nModel) 434 x 324
check Jacobian: wrong dimensions: (0x0) should be (1181x324)  force: 1
jacobian size invalid, forced recalc: 1
calculating jacobian matrix (forced=1)...*** C:/msys64/home/halbm/gimli/gimli/core/src/ttdijkstramodelling.cpp:500
Error!: There are cells with marker -1. Did you define a boundary region? (This is not needed).

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question". Figure 2023-02-07 153536

halbmy commented 1 year ago

To reproduce the error, I would need the data file survey_tom.TOM

halbmy commented 1 year ago

The function drawPLC (called by pg.show) has a keyword argument regionMarker (True by default) that leads to coloring the regions. From the docstring Show region marker one would expect to see the position of the marker as well: image

At any rate, there are only markers 0, 1, 2 so I don't understand Error!: There are cells with marker -1.

The messages Region Nr: 0 is background are clear as it is a fixed region, so that'a ok. I will continue once I have the data file.

halbmy commented 1 year ago

Actually, I noticed that the lowest region marker has also a marker of 1: image However it is on the boundary and does therefore not affect the cells.

Whitt1985 commented 1 year ago

Hi Thomas, I've forwarded the file onto your email. Github wouldn't let me upload the file. As I said regions 2 and 0 should be set at a velocity of 2300 m/s in the inversion. Thanks for looking into this. survey_tom.TOM.gz