1 
ahallam 
3153 

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 
ahallam 
3232 
In this chapter the gravitational potential field is developed for \esc. 
17 


Gravitational fields are present in many modelling scenarios, including 
18 
ahallam 
3410 
geophysical investigations, planetary motion and attraction and microparticle 
19 
ahallam 
3373 
interactions. Gravitational fields also present an opportunity to demonstrate 
20 
ahallam 
3410 
the saving and visualisation of vector data for Mayavi, and the construction of 
21 


variable sized meshes. 
22 
ahallam 
3153 

23 


The gravitational potential $U$ at a point $P$ due to a region with a mass 
24 
ahallam 
3232 
distribution of density $\rho(P)$, is given by Poisson's equation 
25 


\citep{Blakely1995} 
26 
ahallam 
3153 
\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 
ahallam 
3410 
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 
ahallam 
3153 
\begin{equation} 
34 
jfenwick 
3308 
\left(A_{jl} u_{,l} \right)_{,j} = Y 
35 
ahallam 
3153 
\end{equation} 
36 
ahallam 
3373 
One recognises that the LHS side is equivalent to 
37 
ahallam 
3153 
\begin{equation} \label{eqn:ex10a} 
38 


\nabla A \nabla u 
39 


\end{equation} 
40 
ahallam 
3392 
and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to 
41 
ahallam 
3153 
\begin{equation*} 
42 


\nabla^2 u 
43 


\end{equation*} 
44 
ahallam 
3373 
Thus Poisson's \autoref{eqn:poisson} satisfies the general form when 
45 
ahallam 
3153 
\begin{equation} 
46 
jfenwick 
3308 
A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho 
47 
ahallam 
3153 
\end{equation} 
48 
ahallam 
3373 
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 
ahallam 
3410 
of the potential is governed by an inversesquare distance law. Thus, in the 
55 
ahallam 
3373 
limit as $r$ increases; 
56 


\begin{equation} 
57 


\lim_{r\to\infty} U(P) = 0 
58 


\end{equation} 
59 
ahallam 
3410 
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 
ahallam 
3373 

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 
ahallam 
3410 
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 
ahallam 
3373 

73 
ahallam 
3410 
For this example we have set all of the boundaries to zero but only one boundary 
74 
ahallam 
3373 
point needs to be set for the problem to be solvable. The normal flux condition 
75 
ahallam 
3410 
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 
ahallam 
3232 
which can have an influence upon the solution. 
78 
ahallam 
3153 

79 


Setting the boundary condition is relatively simple using the \verb!q! and 
80 
ahallam 
3232 
\verb!r! variables of the general form. First \verb!q! is defined as a masking 
81 


function on the boundary using 
82 
ahallam 
3153 
\begin{python} 
83 


q=whereZero(x[1]my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]mx) 
84 
ahallam 
3410 
mypde.setValue(q=q,r=0) 
85 
ahallam 
3153 
\end{python} 
86 
ahallam 
3232 
This identifies the points on the boundary and \verb!r! is simply 
87 
ahallam 
3410 
ser to \verb!r=0.0!. This is a Dirichlet boundary condition. 
88 
ahallam 
3153 

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

96 
ahallam 
3232 
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 
ahallam 
3373 
$\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where 
106 
ahallam 
3232 
\begin{equation} 
107 
jfenwick 
3308 
g_{z}=\vec{g}\cdot\hat{z} 
108 
ahallam 
3232 
\end{equation} 
109 


Finally, there is the magnitude of the vertical component $g$ of 
110 
jfenwick 
3308 
$g_{z}$ 
111 
ahallam 
3232 
\begin{equation} 
112 
jfenwick 
3308 
g=g_{z} 
113 
ahallam 
3232 
\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 
ahallam 
3410 
g_field=grad(sol) #The gravitational acceleration g. 
118 
ahallam 
3232 
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 
ahallam 
3373 
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 
ahallam 
3410 
object in the explorer tree and selecting the cell to point filter as in 
141 
ahallam 
3373 
\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 
ahallam 
3153 
\begin{figure}[ht] 
148 


\centering 
149 
ahallam 
3373 
\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 
ahallam 
3153 
\includegraphics[width=0.75\textwidth]{figures/ex10apot.png} 
171 
ahallam 
3392 
\caption{Newtonian potential with $\vec{g}$ field directionality. The magnitude 
172 


of the field is reflected in the size of the vector arrows.} 
173 
ahallam 
3153 
\label{fig:ex10pot} 
174 


\end{figure} 
175 
ahallam 
3373 
\clearpage 
176 
ahallam 
3153 

177 


\section{Gravity Well} 
178 


\sslist{example10b.py} 
179 
ahallam 
3232 
Let us now investigate the effect of gravity in three dimensions. Consider a 
180 
ahallam 
3410 
volume which contains a spherical mass anomaly and a gravitational potential 
181 
ahallam 
3373 
which decays to zero at the base of the model. 
182 
ahallam 
3153 

183 
ahallam 
3232 
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 
ahallam 
3410 
integrated into the \esc solution script. 
187 
ahallam 
3232 

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 
ahallam 
3410 
No modifications to the PDE solution section are required. Thus the migration 
201 


from a 2D to a 3D problem is almost trivial. 
202 
ahallam 
3232 

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 
ahallam 
3410 
its strength which is illustrated by the colour. Finally, the streamlines 
208 


highlight the directional flow of the gravity field in this example. 
209 
ahallam 
3232 

210 
ahallam 
3406 
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 
ahallam 
3410 
solution. The profiles in \autoref{fig:ex10p} further demonstrates the affect 
217 
ahallam 
3406 
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 
ahallam 
3410 
appropriately dimensioned model for the problem, and and sampling location. 
220 
ahallam 
3392 

221 
ahallam 
3153 
\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 
ahallam 
3392 
gravitational potential \textemdash small model.} 
226 
ahallam 
3153 
\label{fig:ex10bpot} 
227 


\end{figure} 
228 



229 
ahallam 
3392 
\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 
ahallam 
3153 
\section{Gravity Surface over a fault model.} 
247 


\sslist{example10c.py,example10m.py} 
248 
ahallam 
3232 
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 
ahallam 
3153 
\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 
ahallam 
3392 
\clearpage 
271 
ahallam 
3153 

272 
ahallam 
3392 
\section{Variable meshelement sizes} 
273 
ahallam 
3406 
\sslist{example10m.py} 
274 


We saw in a previous section that the domain needed to be sufficiently 
275 
ahallam 
3410 
large for the boundary conditions to be satisfied. This can be troublesome when 
276 
ahallam 
3406 
trying to solve problems that require a dense mesh, either for solution 
277 
ahallam 
3410 
resolution of stability reasons. The computational cost of solving a large 
278 


number of elements can be prohibitive.S 
279 
ahallam 
3406 

280 
ahallam 
3392 
With the help of Gmsh, it is possible to create a mesh for \esc, which has a 
281 
ahallam 
3406 
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 
ahallam 
3392 

286 


To create a variable size mesh, multiple meshing domains are required. The 
287 
ahallam 
3410 
domains must share points, boundaries and surfaces so that they are joined; and 
288 
ahallam 
3392 
no subdomains are allowed to overlap. Whilst this initially seems complicated, 
289 


it is quite simple to implement. 
290 



291 
ahallam 
3410 
This example creates a mesh which contains a high resolution subdomain at its 
292 
ahallam 
3406 
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 
ahallam 
3392 
\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 
ahallam 
3410 
The two subdomains now have a common geometry and no overlaping features as 
346 
ahallam 
3392 
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 
ahallam 
3410 
d1.addItems(PropertySet(l01,l12,l23,l30)) 
356 
ahallam 
3392 
d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo")) 
357 
ahallam 
3410 
MakeDomain(d1) 
358 
ahallam 
3392 

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 
ahallam 
3410 
MakeDomain(d2) 
364 
ahallam 
3392 
\end{python} 
365 
ahallam 
3410 
Finally, a system call to Gmsh is required to merge and then appropriately 
366 
ahallam 
3392 
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 
ahallam 
3410 
To read in the file the function is called 
385 
ahallam 
3392 
\begin{python} 
386 


domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain 
387 


\end{python} 
388 
ahallam 
3410 
where the integer argument is the number of domain dimensions. 
389 
ahallam 
3392 
% 
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 
ahallam 
3406 
\end{figure} 
404 


\clearpage 
405 



406 


\section{Unbounded problems} 
407 


With a variable elementsize, it is now possible to solve the potential problem 
408 
ahallam 
3410 
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} 