1 

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

15 
\section{Newtonian Potential} 
16 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

247 
% \section{Gravity Surface over a fault model.} 
248 
% \sslist{example10c.py,example10m.py} 
249 
% This model demonstrates the gravity result for a more complicated domain which 
250 
% contains a fault. Additional information will be added when geophysical boundary 
251 
% conditions for a gravity scenario have been established. 
252 
% 
253 
% \begin{figure}[htp] 
254 
% \centering 
255 
% \subfigure[The geometry of the fault model in example10c.py.] 
256 
% {\label{fig:ex10cgeo} 
257 
% \includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\ 
258 
% \subfigure[The fault of interest from the fault model in 
259 
% example10c.py.] 
260 
% {\label{fig:ex10cmsh} 
261 
% \includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}} 
262 
% \end{figure} 
263 
% 
264 
% \begin{figure}[htp] 
265 
% \centering 
266 
% \includegraphics[width=0.8\textwidth]{figures/ex10cpot.png} 
267 
% \caption{Gravitational potential of the fault model with primary layers and 
268 
% faults identified as isosurfaces.} 
269 
% \label{fig:ex10cpot} 
270 
% \end{figure} 
271 
% \clearpage 
272 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

446 
There is a methodology which can help establish an appropriate zero mass region 
447 
to a domain. 
448 
\clearpage 