Gravimetry in 2D - Part IΒΆ

Simple gravimetric field solution.

Calculate for the gravimetric potential \(u\)

\[\frac{\partial u}{\partial z}\]

along a profile for a cylindrical heterogeneity with different approaches.

import numpy as np
import pygimli as pg
from pygimli.meshtools import createCircle, createWorld, createMesh

from pygimli.physics.gravimetry import gradUCylinderHoriz, solveGravimetry

radius = 2.  # [m]
depth = 5.   # [m]
pos = [0., -depth]
dRho = 100

x = np.arange(-20, 20.1, .5)
pnts = np.array([x, np.zeros(len(x))]).T

Analytical solution first

gz_a = gradUCylinderHoriz(pnts, radius, dRho, pos)[:, 1]

Integration for a 2D polygon after [WB87]

circ = createCircle([0, -depth], radius=radius, marker=2, area=0.1,
                    segments=16)
gz_p = solveGravimetry(circ, dRho, pnts, complete=False)

Integration for complete 2D mesh after [WB87]

world = createWorld(start=[-20, 0], end=[20, -10], marker=1)
mesh = createMesh([world, circ])
dRhoC = pg.solver.parseMapToCellArray([[1, 0.0], [2, dRho]], mesh)
gc_m = solveGravimetry(mesh, dRhoC, pnts)

Finite Element solution for \(u\)

world = createWorld(start=[-200, 200], end=[200, -200], marker=1)

# Add some nodes to the measurement points to increase the accuracy a bit
[world.createNode(x_, 0.0,  1) for x_ in x]
mesh = createMesh([world, circ], quality=34)
mesh = mesh.createP2()

density = pg.solver.parseMapToCellArray([[1, 0.0], [2, dRho]], mesh)
u = pg.solver.solve(mesh, a=1, f=density, uB=[[-2, 0], [-1, 0]])

Calculate gradient of gravimetric potential \(\frac{\partial u}{\partial (x,z)}\)

dudz = np.zeros(len(pnts))

for i, p in enumerate(pnts):
    c = mesh.findCell(p)
    g = c.grad(p, u)
    dudz[i] = -g[1] * 4. * np.pi * pg.physics.constants.GmGal  # why 4 pi here?

Finishing the plots

ax1 = pg.plt.subplot(2, 1, 1)
ax1.plot(x, gz_a, '-b', marker='.', label='Analytical')
ax1.plot(x, gz_p, label='Integration: Polygon ')
ax1.plot(x, gc_m, label='Integration: Mesh')
ax1.plot(x, dudz, label=r'FEM: $\frac{\partial u}{\partial z}$')

ax2 = pg.plt.subplot(2, 1, 2)
pg.show([world,  circ], ax=ax2)
ax2.plot(x, x*0,  'bv')

ax1.set_ylabel(r'$\frac{\partial u}{\partial z}$ [mGal]')
ax1.set_xlabel('$x$-coordinate [m]')
ax1.grid()
ax1.legend()

ax2.set_aspect(1)
ax2.set_xlabel('$x$-coordinate [m]')
ax2.set_ylabel('$z$-coordinate [m]')
ax2.set_ylim((-9, 1))
ax2.set_xlim((-20, 20))

pg.wait()
../../_images/sphx_glr_plot_mod-gravimetry-2d_001.png

Total running time of the script: ( 0 minutes 2.968 seconds)

Generated by Sphinx-Gallery

News (Aug 2017):
pyGIMLi Paper out now! Download PDF here.


© 2017 - GIMLi Development Team
Created using Bootstrap and Sphinx. Last updated on Oct 15, 2017.