Model computation is rapid, with a single block taking 2 us to compute on a Celeron 400. This allows models with hundreds of thousands or millions of elements to be computed in a few minutes on a normal desktop workstation. Some benchmarking results on various CPUs are available here.
The setup of the problem allows the computations to be easily parallelized, although not implemented here. As each point of computation and even the model elements are independant of each other, both the grid and computation points can be parallelized to speed the process in a nearly linear manner.
This package does little with model geometry generation; some example generators are included, but complex geometries will be difficult to create with these tools. The computation engine does no element bound or intersection checking; this should be handled in the geometry generator.
From the main README:
README FOR 3D GRAVITY BLOCK MODELING PACKAGE
--------------------------------------------
This package is a set of programs and C subroutines for computing the
gravity effect of a set of rectangular blocks. The blocks can be of
varying sizes, and need not be cubical. Virtually any density
structure can be approximated with this modelling scheme, although it
may be very ineffecient in some cases. It is, however, very general
and robust.
REQUIREMENTS
------------
* - ANSI C++ compiler, to compile the model computation, visualization,
and geometry generation code.
* - Python, for a model geometry and horizontal grid generator;
Python 1.5.2 has been tested, but anything after 1.5 should work.
* - VRML97/VRML2.0 viewer to use the visualization code
* - OpenGL library (libGL, libglut) for the OpenGL vis program
* - Some disk space; models take several hundred kB to several
tens of MB of space.
* - CPU cycles to burn.
LICENSE
-------
This code is released under the GPL. Have fun. Please send
improvements, bug fixes/reports, etc. back to the author
(gettings@mines.utah.edu).
There is no warranty, express or implied, of any kind.
UNITS
-----
For all measurements in the modelling code, units of meters and
kilograms are assumed. Hence, positions and dimensions should be in m,
and densities in kg/m^3. This means that results of the modelling are
returned as accelerations in m/s^2.
Careful use of other units can alter the result of the modelling to almost
arbitrary metric units. The only potential problem to the use of any
unit system is the definition of the universal gravitational constant
GAMMA; it is defined in model.h to be 6.671e-11, which is suitable for
m and kg. For other unit systems, this constant should be changed.
Alternatively, compute the effect in m/s/s and then apply a conversion
factor to the result.
Be careful with very small gravity effects to avoid underflow or loss
of precision. For models with gravity effects in the range of microGals
(1e-8 m/s/s), this is not a problem; doubles have a precision of roughly
1e-16 on a standard PC/workstation.
AXES
----
The model space defines a single set of right-handed Cartesian axes, with Z
positive being DOWN. The actual origin is arbitrary; the axes are not
referenced to any real-world systems.
All calculations, geometries, and points are referenced to this single
system.
DESCRIPTION OF GRAVITY ALGORITHM
--------------------------------
This algorithm depends on expanding the gravitational potential of
a block using Legendre polynomials to simplify the volume integral
evaluation into a series. It is then necessary to find the moments for
these polynomials for a block. Under the simplifying assumptions that
the coordinate origin is at the center of mass of the block, and that
only moments up to l=2 are important (see reference for justification), the
resulting series to compute has only 3 terms, with analytical moments
depending only on the block dimensions and density. By taking the
derivative with respect to z of the potential, the vertical component
of the gravitational attraction (Gz) is found. This can be done
analytically:
From Grant & West, Interpretation Theory in Applied Geophysics, 1965,
pp 222-227.
Start by defining a rectangular block with center of mass (and geometric
center) at the origin (0,0,0), uniform density p, and dimensions of (2a,
2b, 2c) in the (x, y, z) directions respectively.
Want to find the gravity effect at some position (x,y,z) that is
outside the block. Start by recalling that the gravity effect can be
found by differentiating the potential; so if the potential is found,
taking d/dz will yield the desired gravity effect.
Following Grant & West, expand the gravitational potential into a
series of Legendre polynomials with some set of moments. This
expansion is detailed in the text. The result is that the potential
is:
A B(3z^2 - r^2) 3C(x^2 - y^2)
-U(x,y,z) = - + ------------- + -------------
r 2*r^5 r^5
where r = (x^2 + y^2 + z^2)^0.5
and
0
A = 8Gpabc = B
0
A(2c^2 - a^2 - b^2) 0
B = ------------------- = B
6 2
A(a^2 - b^2) 2
C = ------------ = B
24 2
So need to find d/dz[U(x,y,z)]. After some algebra, the result is:
dU Az B 5z(3z^2 - r^2) 15Cz(x^2 - y^2)
-- = --- + ----[-------------- - 4z] + ---------------
dz r^3 2r^5 r^2 r^7