Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
Whitt1985 opened this issue Feb 7, 2023 · 4 comments
Open

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

Whitt1985 opened this issue Feb 7, 2023 · 4 comments
Assignees
Labels

Comments

@Whitt1985
Copy link

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 halbmy self-assigned this Feb 8, 2023
@halbmy
Copy link
Contributor

halbmy commented Feb 8, 2023

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

@halbmy
Copy link
Contributor

halbmy commented Feb 8, 2023

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
Copy link
Contributor

halbmy commented Feb 8, 2023

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
Copy link
Author

Whitt1985 commented Feb 9, 2023

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

@Whitt1985 Whitt1985 reopened this Feb 9, 2023
@halbmy halbmy added the bug label Oct 24, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants