This example shows how the FMM implementation of pyGIMLi works. FMM is an alternative to the Shortest path method and utilized by the Refraction manager. We illustrate the effectiveness and compare the solutions.
The first-arrival travel-time \(t\) from a shot point at the origin is described by a Hamilton-Jacobi (HJ) type equation, known as Eikonal equation [PL91]
where \(s\) is the spatially distributed slowness (i.e. 1 / velocity). Under the assumption of piece-wise slowness (for individual cells) one can solve this equation by the so-called Fast-Marching method as first developed for 2D grids by Podvon & Lecompte (1991). Starting from the source location, a front is successively increased and filled by the travel times.
This example demonstrates the core features of an FMM Python implementation and compares the travel time for a two-layer model with the analytical solution (known for straight or sloped layers from the text books) and the shortest-path algorithm (Moser, 1991) implemented in the pyGIMLi core. [Mos91]
First we provide the analytical solution for a given offset vector x.
def analyticalSolution2Layer(x, zlay, v1, v2): """Analytical solution: minimum of direct and critically refracted wave.""" tdirect = np.abs(x) / v1 # direct wave alfa = asin(v1 / v2) # critically refracted wave angle xreflec = tan(alfa) * zlay * 2. # first critically refracted trefrac = (x - xreflec) / v2 + xreflec * v2 / v1**2 return np.minimum(tdirect, trefrac)
The model consists of two boxes, in the first is the source (1)
0--1--------2 | | 6-----------3 | | 5-----------4
We create a PLC (piece-wise linear complex) and insert the nodes.
xmin, xmax, zlay = -20., 150., 20. # model dimensions plc = pg.Mesh(2) nodes =  nodes.append(plc.createNode(xmin, 0., 0.)) # 0 nodes.append(plc.createNode(0.0, 0., 0.)) # 1 nodes.append(plc.createNode(xmax, 0., 0.)) # 2 nodes.append(plc.createNode(xmax, -zlay, 0.)) # 3 nodes.append(plc.createNode(xmax, -zlay * 2, 0.)) # 4 nodes.append(plc.createNode(xmin, -zlay * 2, 0.)) # 5 nodes.append(plc.createNode(xmin, -zlay, 0.)) # 6
The nodes are connected from from 0 to 6 and back to 0. An additional edge is drawn from 6 to 3. Node/edge markers do not matter.
<pygimli.core._pygimli_.Edge object at 0x7f40fc756e90>
We insert region markers (0 and 1) into the two layers and generate the mesh.
Mesh: Nodes: 1781 Cells: 3405 Boundaries: 5185
Next we generate a velocity model from the markers by using a map. The values are associated to the markers and stored as attributes.
We initialize the source position and the travel time vector initialize sets and tags and define the initial condition.
source = pg.RVector3(0., 0.) # does not have to be a node! times = pg.RVector(mesh.nodeCount(), 0.) upwind, downwind = set(), set() upTags, downTags = np.zeros(mesh.nodeCount()), np.zeros(mesh.nodeCount()) cell = mesh.findCell(source) for i, n in enumerate(cell.nodes()): times[n.id()] = cell.attribute() * n.pos().distance(source) upTags[n.id()] = 1 for i, n in enumerate(cell.nodes()): tmpNodes = pg.commonNodes(n.cellSet()) for nn in tmpNodes: if not upTags[nn.id()] and not downTags[nn.id()]: downwind.add(nn) downTags[nn.id()] = 1
Then we start marching until all fields are filled.
First, we plot the traveltime field and streamlines
<matplotlib.streamplot.StreamplotSet object at 0x7f40fc8e2208>
We compare the result with the analytical solution along the x axis
min(dt)=-0.0033639721106065723 ms max(dt)=1.6065909149150792 ms
In order to use the Dijkstra, we extract the surface positions >0
A data container with index arrays named s (shot) and g (geophones) is created and filled with the positions and shot/geophone indices.
data = pg.DataContainer() data.registerSensorIndex('s') data.registerSensorIndex('g') for pxi in px: data.createSensor(pg.RVector3(pxi, 0.0)) ndata = len(px) - 1 data.resize(ndata) data.set('s', pg.RVector(ndata, 0)) # only one shot at first sensor data.set('g', pg.utils.grange(1, ndata, 1)) # all others and geophones fop = pg.TravelTimeDijkstraModelling(mesh, data) tDijkstra = fop.response(mesh.cellAttributes())
We plot the calculated and measured travel times and relative differences
fig, ax = plt.subplots() ax.plot(x, tAna*1000, 'r-', label='analytical') ax.plot(x, tFMM*1000, 'b+-', label='FMM') ax.plot(px[1:], tDijkstra*1000, 'gx-', label='Dijkstra') ax.set_xlabel('x [m]') ax.set_ylabel('t [s]') ax.grid(True) ax.legend() ax2 = ax.twinx() dtFMM = (tFMM - tAna) * 1000 tAnaD = analyticalSolution2Layer(px[1:], zlay, v, v) dtDijkstra = (tDijkstra - tAnaD) * 1000 ax2.set_ylabel('dt [ms]') if 1: # relative differences in percent ax2.set_ylabel('dt [%]') dtFMM /= (tAna * 10) dtDijkstra /= (tAnaD * 10) ax2.plot(x, dtFMM, 'b.-') ax2.plot(px[1:], dtDijkstra, 'g.-') pg.wait()
Note that the Fast Marching Method is implemented in a modelling class (right now Python and very slow but to be replaced by fast C++ soon) FMModelling that can be more easily used with the Refraction manager by ra.useFMM(True)
Total running time of the script: ( 0 minutes 12.919 seconds)