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. 
21 

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

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

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

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

86 
\section{Gravity Pole} 
87 
\sslist{example10a.py} 
88 
A gravity pole is used in this example to demonstrate the radial directionality 
89 
of gravity, and also to show how this information can be exported for 
90 
visualisation to Mayavi or an equivalent using the VTK data format. 
91 

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

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

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

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

139 
The settings can then be modified to suit personal taste. An example of the 
140 
potential and gravitational field vectors is illustrated in 
141 
\autoref{fig:ex10pot}. 
142 

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

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

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

164 
\begin{figure}[ht] 
165 
\centering 
166 
\includegraphics[width=0.75\textwidth]{figures/ex10apot.png} 
167 
\caption{Newtonian potential with g field directionality.} 
168 
\label{fig:ex10pot} 
169 
\end{figure} 
170 
\clearpage 
171 

172 
\section{Gravity Well} 
173 
\sslist{example10b.py} 
174 
Let us now investigate the effect of gravity in three dimensions. Consider a 
175 
volume which contains a sphericle mass anomaly and a gravitational potential 
176 
which decays to zero at the base of the model. 
177 

178 
The script used to solve this model is very similar to that used for the gravity 
179 
pole in the previous section, but with an extra spatial dimension. As for all 
180 
the 3D problems examined in this cookbook, the extra dimension is easily 
181 
integrated into the \esc solultion script. 
182 

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

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

205 
\begin{figure}[htp] 
206 
\centering 
207 
\includegraphics[width=0.75\textwidth]{figures/ex10bpot.png} 
208 
\caption{Gravity well with iso surfaces and streamlines of the vector 
209 
gravitational potential.} 
210 
\label{fig:ex10bpot} 
211 
\end{figure} 
212 

213 
\section{Gravity Surface over a fault model.} 
214 
\sslist{example10c.py,example10m.py} 
215 
This model demonstrates the gravity result for a more complicated domain which 
216 
contains a fault. Additional information will be added when geophysical boundary 
217 
conditions for a gravity scenario have been established. 
218 

219 
\begin{figure}[htp] 
220 
\centering 
221 
\subfigure[The geometry of the fault model in example10c.py.] 
222 
{\label{fig:ex10cgeo} 
223 
\includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\ 
224 
\subfigure[The fault of interest from the fault model in 
225 
example10c.py.] 
226 
{\label{fig:ex10cmsh} 
227 
\includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}} 
228 
\end{figure} 
229 

230 
\begin{figure}[htp] 
231 
\centering 
232 
\includegraphics[width=0.8\textwidth]{figures/ex10cpot.png} 
233 
\caption{Gravitational potential of the fault model with primary layers and 
234 
faults identified as isosurfaces.} 
235 
\label{fig:ex10cpot} 
236 
\end{figure} 
237 
