1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
7 
% 
8 
% Primary Business: Queensland, Australia 
9 
% Licensed under the Open Software License version 3.0 
10 
% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Newtonian Potential} 
15 

16 
In this chapter the gravitational potential field is developed for \esc. 
17 
Gravitational fields are present in many modelling scenarios, including 
18 
geophysical investigations, planetary motion and attraction and microparticle 
19 
interactions. Gravitational fields also present an opportunity to demonstrate 
20 
the saving and visualisation of vector data for Mayavi, and the construction of 
21 
variable sized meshes. 
22 

23 
The gravitational potential $U$ at a point $P$ due to a region with a mass 
24 
distribution of density $\rho(P)$, is given by Poisson's equation 
25 
\citep{Blakely1995} 
26 
\begin{equation} \label{eqn:poisson} 
27 
\nabla^2 U(P) = 4\pi\gamma\rho(P) 
28 
\end{equation} 
29 
where $\gamma$ is the gravitational constant. 
30 
Consider now the \esc general form, 
31 
\autoref{eqn:poisson} requires only two coefficients, 
32 
$A$ and $Y$, thus the relevant terms of the general form are reduced to; 
33 
\begin{equation} 
34 
\left(A_{jl} u_{,l} \right)_{,j} = Y 
35 
\end{equation} 
36 
One recognises that the LHS side is equivalent to 
37 
\begin{equation} \label{eqn:ex10a} 
38 
\nabla A \nabla u 
39 
\end{equation} 
40 
and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to 
41 
\begin{equation*} 
42 
\nabla^2 u 
43 
\end{equation*} 
44 
Thus Poisson's \autoref{eqn:poisson} satisfies the general form when 
45 
\begin{equation} 
46 
A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho 
47 
\end{equation} 
48 
The boundary condition is the last parameter that requires consideration. The 
49 
potential $U$ is related to the mass of a sphere by 
50 
\begin{equation} 
51 
U(P)=\gamma \frac{m}{r^2} 
52 
\end{equation} where $m$ is the mass of the sphere and $r$ is the distance from 
53 
the center of the mass to the observation point $P$. Plainly, the magnitude 
54 
of the potential is governed by an inversesquare distance law. Thus, in the 
55 
limit as $r$ increases; 
56 
\begin{equation} 
57 
\lim_{r\to\infty} U(P) = 0 
58 
\end{equation} 
59 
Provided that the domain being solved is large enough, and the source mass is 
60 
contained within a central region of the domain, the potential will decay to 
61 
zero. This is a dirichlet boundary condition where $U=0$. 
62 

63 
This boundary condition can be satisfied when there is some mass suspended in a 
64 
freespace. For geophysical models where the mass of interest may be an anomaly 
65 
inside a host rock, the anomaly can be isolated by subtracting the density of the 
66 
host rock from the model. This creates a fictitious free space model that will 
67 
satisfy the analytic boundary conditions. The result is that 
68 
$Y=4\pi\gamma\Delta\rho$, where $\Delta\rho=\rho\rho_0$ and $\rho_0$ is the 
69 
baseline or host density. This of course means that the true gravity response is 
70 
not being modelled, but the response due to an anomalous mass with a 
71 
density contrast $\Delta\rho$. 
72 

73 
For this example we have set all of the boundaries to zero but only one boundary 
74 
point needs to be set for the problem to be solvable. The normal flux condition 
75 
is also zero by default. Note, that for a more realistic and complicated models 
76 
it may be necessary to give careful consideration to the boundary conditions of the model, 
77 
which can have an influence upon the solution. 
78 

79 
Setting the boundary condition is relatively simple using the \verb!q! and 
80 
\verb!r! variables of the general form. First \verb!q! is defined as a masking 
81 
function on the boundary using 
82 
\begin{python} 
83 
q=whereZero(x[1]my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]mx) 
84 
mypde.setValue(q=q,r=0) 
85 
\end{python} 
86 
This identifies the points on the boundary and \verb!r! is simply 
87 
ser to \verb!r=0.0!. This is a Dirichlet boundary condition. 
88 

89 
\clearpage 
90 
\section{Gravity Pole} 
91 
\sslist{example10a.py} 
92 
A gravity pole is used in this example to demonstrate the vector characteristics 
93 
of gravity, and also to demonstrate how this information can be exported for 
94 
visualisation to Mayavi or an equivalent using the VTK data format. 
95 

96 
The solution script for this section is very simple. First the domain is 
97 
constructed, then the parameters of the model are set, and finally the steady 
98 
state solution is found. There are quite a few values that can now be derived 
99 
from the solution and saved to file for visualisation. 
100 

101 
The potential $U$ is related to the gravitational response $\vec{g}$ via 
102 
\begin{equation} 
103 
\vec{g} = \nabla U 
104 
\end{equation} 
105 
$\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where 
106 
\begin{equation} 
107 
g_{z}=\vec{g}\cdot\hat{z} 
108 
\end{equation} 
109 
Finally, there is the magnitude of the vertical component $g$ of 
110 
$g_{z}$ 
111 
\begin{equation} 
112 
g=g_{z} 
113 
\end{equation} 
114 
These values are derived from the \esc solution \verb!sol! to the potential $U$ 
115 
using the following commands 
116 
\begin{python} 
117 
g_field=grad(sol) #The gravitational acceleration g. 
118 
g_fieldz=g_field*[0,1] #The vertical component of the g field. 
119 
gz=length(g_fieldz) #The magnitude of the vertical component. 
120 
\end{python} 
121 
This data can now be simply exported to a VTK file via 
122 
\begin{python} 
123 
# Save the output to file. 
124 
saveVTK(os.path.join(save_path,"ex10a.vtu"),\ 
125 
grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz) 
126 
\end{python} 
127 

128 
It is quite simple to visualise the data from the gravity solution in Mayavi2. 
129 
With Mayavi2 open go to File, Load data, Open file \ldots as in 
130 
\autoref{fig:mayavi2openfile} and select the saved data file. The data will 
131 
have then been loaded and is ready for visualisation. Notice that under the data 
132 
object in the Mayavi2 navigation tree the 4 values saved to the VTK file are 
133 
available (\autoref{fig:mayavi2data}). There are two vector values, 
134 
\verbgfield and \verbgfieldz. Note that to plot both of these on the same 
135 
chart requires that the data object be imported twice. 
136 

137 
The point scalar data \verbgrav_pot is the gravitational potential and it is 
138 
easily plotted using a surface module. To visualise the cell data a filter is 
139 
required that converts to point data. This is done by right clicking on the data 
140 
object in the explorer tree and selecting the cell to point filter as in 
141 
\autoref{fig:mayavi2cell2point}. 
142 

143 
The settings can then be modified to suit personal taste. An example of the 
144 
potential and gravitational field vectors is illustrated in 
145 
\autoref{fig:ex10pot}. 
146 

147 
\begin{figure}[ht] 
148 
\centering 
149 
\includegraphics[width=0.75\textwidth]{figures/mayavi2_openfile.png} 
150 
\caption{Open a file in Mayavi2} 
151 
\label{fig:mayavi2openfile} 
152 
\end{figure} 
153 

154 
\begin{figure}[ht] 
155 
\centering 
156 
\includegraphics[width=0.75\textwidth]{figures/mayavi2_data.png} 
157 
\caption{The 4 types of data in the imported VTK file.} 
158 
\label{fig:mayavi2data} 
159 
\end{figure} 
160 

161 
\begin{figure}[ht] 
162 
\centering 
163 
\includegraphics[width=0.75\textwidth]{figures/mayavi2_cell2point.png} 
164 
\caption{Converting cell data to point data.} 
165 
\label{fig:mayavi2cell2point} 
166 
\end{figure} 
167 

168 
\begin{figure}[ht] 
169 
\centering 
170 
\includegraphics[width=0.75\textwidth]{figures/ex10apot.png} 
171 
\caption{Newtonian potential with $\vec{g}$ field directionality. The magnitude 
172 
of the field is reflected in the size of the vector arrows.} 
173 
\label{fig:ex10pot} 
174 
\end{figure} 
175 
\clearpage 
176 

177 
\section{Gravity Well} 
178 
\sslist{example10b.py} 
179 
Let us now investigate the effect of gravity in three dimensions. Consider a 
180 
volume which contains a spherical mass anomaly and a gravitational potential 
181 
which decays to zero at the base of the model. 
182 

183 
The script used to solve this model is very similar to that used for the gravity 
184 
pole in the previous section, but with an extra spatial dimension. As for all 
185 
the 3D problems examined in this cookbook, the extra dimension is easily 
186 
integrated into the \esc solution script. 
187 

188 
The domain is redefined from a rectangle to a Brick; 
189 
\begin{python} 
190 
domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz) 
191 
\end{python} 
192 
the source is relocated along $z$; 
193 
\begin{python} 
194 
x=x[mx/2,my/2,mz/2] 
195 
\end{python} 
196 
and, the boundary conditions are updated. 
197 
\begin{python} 
198 
q=whereZero(x[2]inf(x[2])) 
199 
\end{python} 
200 
No modifications to the PDE solution section are required. Thus the migration 
201 
from a 2D to a 3D problem is almost trivial. 
202 

203 
\autoref{fig:ex10bpot} illustrates the strength of a PDE solution. Three 
204 
different visualisation types help define and illuminate properties of the data. 
205 
The cut surfaces of the potential are similar to a 2D section for a given x or y 
206 
and z. The isosurfaces illuminate the 3D shape of the gravity field, as well as 
207 
its strength which is illustrated by the colour. Finally, the streamlines 
208 
highlight the directional flow of the gravity field in this example. 
209 

210 
The boundary conditions were discussed previously, but not thoroughly 
211 
investigated. It was stated, that in the limit as the boundary becomes more 
212 
remote from the source, the potential will reduce to zero. 
213 
\autoref{fig:ex10bpot2} is the solution to the same gravity problem 
214 
but with a slightly larger domain. It is obvious in this case that 
215 
the previous domain size was too small to accurately represent the 
216 
solution. The profiles in \autoref{fig:ex10p} further demonstrates the affect 
217 
the domain size has upon the solution. As the domain size increases, the 
218 
profiles begin to converge. This convergence is a good indicator of an 
219 
appropriately dimensioned model for the problem, and and sampling location. 
220 

221 
\begin{figure}[htp] 
222 
\centering 
223 
\includegraphics[width=0.75\textwidth]{figures/ex10bpot.png} 
224 
\caption{Gravity well with iso surfaces and streamlines of the vector 
225 
gravitational potential \textemdash small model.} 
226 
\label{fig:ex10bpot} 
227 
\end{figure} 
228 

229 
\begin{figure}[htp] 
230 
\centering 
231 
\includegraphics[width=0.75\textwidth]{figures/ex10bpot2.png} 
232 
\caption{Gravity well with iso surfaces and streamlines of the vector 
233 
gravitational potential \textemdash large model.} 
234 
\label{fig:ex10bpot2} 
235 
\end{figure} 
236 

237 
\begin{figure}[htp] 
238 
\centering 
239 
\includegraphics[width=0.85\textwidth]{figures/ex10p_boundeff.pdf} 
240 
\caption{Profile of the graviational provile along x where $y=0,z=250$ for 
241 
various sized domains.} 
242 
\label{fig:ex10p} 
243 
\end{figure} 
244 
\clearpage 
245 

246 
\section{Gravity Surface over a fault model.} 
247 
\sslist{example10c.py,example10m.py} 
248 
This model demonstrates the gravity result for a more complicated domain which 
249 
contains a fault. Additional information will be added when geophysical boundary 
250 
conditions for a gravity scenario have been established. 
251 

252 
\begin{figure}[htp] 
253 
\centering 
254 
\subfigure[The geometry of the fault model in example10c.py.] 
255 
{\label{fig:ex10cgeo} 
256 
\includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\ 
257 
\subfigure[The fault of interest from the fault model in 
258 
example10c.py.] 
259 
{\label{fig:ex10cmsh} 
260 
\includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}} 
261 
\end{figure} 
262 

263 
\begin{figure}[htp] 
264 
\centering 
265 
\includegraphics[width=0.8\textwidth]{figures/ex10cpot.png} 
266 
\caption{Gravitational potential of the fault model with primary layers and 
267 
faults identified as isosurfaces.} 
268 
\label{fig:ex10cpot} 
269 
\end{figure} 
270 
\clearpage 
271 

272 
\section{Variable meshelement sizes} 
273 
\sslist{example10m.py} 
274 
We saw in a previous section that the domain needed to be sufficiently 
275 
large for the boundary conditions to be satisfied. This can be troublesome when 
276 
trying to solve problems that require a dense mesh, either for solution 
277 
resolution of stability reasons. The computational cost of solving a large 
278 
number of elements can be prohibitive.S 
279 

280 
With the help of Gmsh, it is possible to create a mesh for \esc, which has a 
281 
variable element size. Such an approach can significantly reduce the number of 
282 
elements that need to be solved, and a large domain can be created that has 
283 
sufficient resolution close to the source and extends to distances large enough 
284 
that the infinity clause is satisfied. 
285 

286 
To create a variable size mesh, multiple meshing domains are required. The 
287 
domains must share points, boundaries and surfaces so that they are joined; and 
288 
no subdomains are allowed to overlap. Whilst this initially seems complicated, 
289 
it is quite simple to implement. 
290 

291 
This example creates a mesh which contains a high resolution subdomain at its 
292 
center. We begin by defining two curve loops which describe the large or 
293 
big subdomain and the smaller subdomain which is to contain the high 
294 
resolution portion of the mesh. 
295 
\begin{python} 
296 
################################################BIG DOMAIN 
297 
#ESTABLISHING PARAMETERS 
298 
width=10000. #width of model 
299 
depth=10000. #depth of model 
300 
bele_size=500. #big element size 
301 
#DOMAIN CONSTRUCTION 
302 
p0=Point(0.0, 0.0) 
303 
p1=Point(width, 0.0) 
304 
p2=Point(width, depth) 
305 
p3=Point(0.0, depth) 
306 
# Join corners in anticlockwise manner. 
307 
l01=Line(p0, p1) 
308 
l12=Line(p1, p2) 
309 
l23=Line(p2, p3) 
310 
l30=Line(p3, p0) 
311 

312 
cbig=CurveLoop(l01,l12,l23,l30) 
313 

314 
################################################SMALL DOMAIN 
315 
#ESTABLISHING PARAMETERS 
316 
xwidth=2000.0 #x width of model 
317 
zdepth=2000.0 #y width of model 
318 
sele_size=10. #small element size 
319 
#TRANSFORM 
320 
xshift=width/2xwidth/2 
321 
zshift=depth/2zdepth/2 
322 
#DOMAIN CONSTRUCTION 
323 
p4=Point(xshift, zshift) 
324 
p5=Point(xwidth+xshift, zshift) 
325 
p6=Point(xwidth+xshift, zdepth+zshift) 
326 
p7=Point(xshift, zdepth+zshift) 
327 
# Join corners in anticlockwise manner. 
328 
l45=Line(p4, p5) 
329 
l56=Line(p5, p6) 
330 
l67=Line(p6, p7) 
331 
l74=Line(p7, p4) 
332 

333 
csmall=CurveLoop(l45,l56,l67,l74) 
334 
\end{python} 
335 
The small subdomain curve can then be used to create a surface. 
336 
\begin{python} 
337 
ssmall=PlaneSurface(csmall) 
338 
\end{python} 
339 
However, so that the two domains do not overlap, when the big subdomain 
340 
curveloop is used to create a surface it must contain a hole. The hole is 
341 
defined by the small subdomain curveloop. 
342 
\begin{python} 
343 
sbig=PlaneSurface(cbig,holes=[csmall]) 
344 
\end{python} 
345 
The two subdomains now have a common geometry and no overlaping features as 
346 
per \autoref{fig:ex10mgeo}. Notice, that both domains have a normal in the 
347 
same direction. 
348 

349 
The next step, is exporting each subdomain individually, with an appropriate 
350 
element size. This is carried out using the \pycad Design command. 
351 
\begin{python} 
352 
# Design the geometry for the big mesh. 
353 
d1=Design(dim=2, element_size=bele_size, order=1) 
354 
d1.addItems(sbig) 
355 
d1.addItems(PropertySet(l01,l12,l23,l30)) 
356 
d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo")) 
357 
MakeDomain(d1) 
358 

359 
# Design the geometry for the small mesh. 
360 
d2=Design(dim=2, element_size=sele_size, order=1) 
361 
d2.addItems(ssmall) 
362 
d2.setScriptFileName(os.path.join(save_path,"example10m_small.geo")) 
363 
MakeDomain(d2) 
364 
\end{python} 
365 
Finally, a system call to Gmsh is required to merge and then appropriately 
366 
mesh the two subdomains together. 
367 
\begin{python} 
368 
# Join the two meshes using Gmsh and then apply a 2D meshing algorithm. 
369 
# The small mesh must come before the big mesh in the merging call!!@!!@! 
370 
sp.call("gmsh 2 "+ 
371 
os.path.join(save_path,"example10m_small.geo")+" "+ 
372 
os.path.join(save_path,"example10m_big.geo")+" o "+ 
373 
os.path.join(save_path,"example10m.msh"),shell=True) 
374 
\end{python} 
375 
The ``2'' option is responsible for the 2D meshing, and the ``o'' option 
376 
provides the output path. The resulting mesh is depicted in 
377 
\autoref{fig:ex10mmsh} 
378 

379 
To use the Gmsh ``*.msh'' file in the solution script, the mesh reading function 
380 
``ReadGmsh'' is required. It can be imported via; 
381 
\begin{python} 
382 
from esys.finley import ReadGmsh 
383 
\end{python} 
384 
To read in the file the function is called 
385 
\begin{python} 
386 
domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain 
387 
\end{python} 
388 
where the integer argument is the number of domain dimensions. 
389 
% 
390 
\begin{figure}[ht] 
391 
\centering 
392 
\includegraphics[width=0.8\textwidth]{figures/ex10m_geo.png} 
393 
\caption{Geometry of two surfaces for a single domain.} 
394 
\label{fig:ex10mgeo} 
395 
\end{figure} 
396 

397 
\begin{figure}[ht] 
398 
\centering 
399 
\includegraphics[width=0.8\textwidth]{figures/ex10m_msh.png} 
400 
\caption{Mesh of merged surfaces, showing variable element size. Elements 
401 
range from 10m in the centroid to 500m at the boundary.} 
402 
\label{fig:ex10mmsh} 
403 
\end{figure} 
404 
\clearpage 
405 

406 
\section{Unbounded problems} 
407 
With a variable elementsize, it is now possible to solve the potential problem 
408 
over a very large mesh. To test the accuracy of the solution, we will compare 
409 
the \esc result with the analytic solution for the vertical gravitational 
410 
acceleration $g_z$ of an infinite horizontal cylinder. 
411 

412 
For a horizontal cylinder with a circular crosssection with infinite strike, 
413 
the analytic solution is give by 
414 
\begin{equation} 
415 
g_z = 2\gamma\pi R^2 \Delta\rho \frac{z}{(x^2+z^2)} 
416 
\end{equation} 
417 
where $\gamma$ is the gravitational constant (as defined previously), $R$ is the 
418 
radius of the cylinder, $\Delta\rho$ is the density contrast and $x$ and $z$ are 
419 
geometric factors, relating the observation point to the center of the source 
420 
via the horizontal and vertical displacements respectively. 
421 

422 
The accuracy of the solution was tested using a square domain. For each test the 
423 
dimensions of the domain were modified, being set to 5, 10, 20 and 40 Km. The 
424 
results are compared with the analytic solution and are depicted in 
425 
\autoref{fig:ex10q boundeff} and \autoref{fig:ex10q boundeff zoom}. Clearly, 
426 
as the domain size increases, the results are valid at greater 
427 
distances from the source. The same is true at the anomaly peak, where the 
428 
variation around the source diminishes with an increasing domain size. 
429 

430 
\begin{figure}[ht] 
431 
\centering 
432 
\includegraphics[width=0.8\textwidth]{figures/ex10q_boundeff.pdf} 
433 
\caption{Solution profile 1000.0 meters from the source as the domain size 
434 
increases.} 
435 
\label{fig:ex10q boundeff} 
436 
\end{figure} 
437 

438 
\begin{figure}[ht] 
439 
\centering 
440 
\includegraphics[width=0.8\textwidth]{figures/ex10q_boundeff_zoom.pdf} 
441 
\caption{Magnification of \autoref{fig:ex10q boundeff}.} 
442 
\label{fig:ex10q boundeff zoom} 
443 
\end{figure} 
444 
\clearpage 
445 

446 
\subsection{Inversion using scipy} 