1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

16 
\section{Newtonian Potential} 
17 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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