doc/cookbook/example10.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2003-2010 by University of Queensland
% Earth Systems Science Computational Center (ESSCC)
% http://www.uq.edu.au/esscc
%
% Primary Business: Queensland, Australia
% Licensed under the Open Software License version 3.0
% http://www.opensource.org/licenses/osl-3.0.php
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Newtonian Potential}

In this chapter the gravitational potential field is developed for \esc.
Gravitational fields are present in many modelling scenarios, including
geophysical investigations, planetary motion and attraction and microparticle
interactions. Gravitational fields also presents the opportunity to demonstrate
the saving and visualisation of vector data for Mayavi.

The gravitational potential $U$ at a point $P$ due to a region with a mass
distribution of density $\rho(P)$, is given by Poisson's equation
\citep{Blakely1995}
\begin{equation} \label{eqn:poisson}
\nabla^2 U(P) = -4\pi\gamma\rho(P)
\end{equation}
where $\gamma$ is the gravitational constant.
Consider now the \esc general form, which can simply be related to
\auoref{eqn:poisson} using two coefficients. The result is
\begin{equation}
-\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y
\end{equation}
one recognises that the LHS side is equivalent to 
\begin{equation} \label{eqn:ex10a}
-\nabla A \nabla u
\end{equation}
If $A=\delta\hackscore{jl}$ then \autoref{eqn:ex10a} is equivalent to
\begin{equation*}
-\nabla^2 u
\end{equation*}
and thus Poisson's \autoref{eqn:poisson} satisfies the general form when
\begin{equation}
A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho
\end{equation}
At least one boundary point must be set for the problem to be solvable. For this
example we have set all of the boundaries to zero. The normal flux condition is
also zero by default. For a more realistic and complicated models it may be
necessary to give careful consideration to the boundary conditions of the model,
which can have an influence upon the solution.

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

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

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

The potential $U$ is related to the gravitational response $\vec{g}$ via
\begin{equation}
\vec{g} = \nabla U
\end{equation}
This for example as a vertical component $g\hackscore{z}$ where
\begin{equation}
g\hackscore{z}=\vec{g}\cdot\hat{z}
\end{equation}
Finally, there is the magnitude of the vertical component $g$ of
$g\hackscore{z}$
\begin{equation}
g=|g\hackscore{z}|
\end{equation}
These values are derived from the \esc solution \verb!sol! to the potential $U$
using the following commands
\begin{python}
g_field=grad(sol) #The graviational accelleration g.
g_fieldz=g_field*[0,1] #The vertical component of the g field.
gz=length(g_fieldz) #The magnitude of the vertical component.
\end{python}
This data can now be simply exported to a VTK file via
\begin{python}
# Save the output to file.
saveVTK(os.path.join(save_path,"ex10a.vtu"),\
        grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
\end{python}

\begin{figure}[ht]
\centering
\includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
\caption{Newtonian potential with g field directionality.}
\label{fig:ex10pot}
\end{figure}

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

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

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

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

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

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

\begin{figure}[htp]
\centering
\subfigure[The geometry of the fault model in example10c.py.]
    {\label{fig:ex10cgeo}
     \includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\
\subfigure[The fault of interest from the fault model in
    example10c.py.] 
    {\label{fig:ex10cmsh}
     \includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}}
\end{figure}

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

1	ahallam	3153
2			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3			%
4			% Copyright (c) 2003-2010 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/osl-3.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			geophysical investigations, planetary motion and attraction and microparticle
19			interactions. Gravitational fields also presents the opportunity to demonstrate
20			the saving and visualisation of vector data for Mayavi.
21	ahallam	3153
22			The gravitational potential $U$ at a point $P$ due to a region with a mass
23	ahallam	3232	distribution of density $\rho(P)$, is given by Poisson's equation
24			\citep{Blakely1995}
25	ahallam	3153	\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	ahallam	3232	Consider now the \esc general form, which can simply be related to
30			\auoref{eqn:poisson} using two coefficients. The result is
31	ahallam	3153	\begin{equation}
32			-\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y
33			\end{equation}
34	ahallam	3232	one recognises that the LHS side is equivalent to
35	ahallam	3153	\begin{equation} \label{eqn:ex10a}
36			-\nabla A \nabla u
37			\end{equation}
38	ahallam	3232	If $A=\delta\hackscore{jl}$ then \autoref{eqn:ex10a} is equivalent to
39	ahallam	3153	\begin{equation*}
40			-\nabla^2 u
41			\end{equation*}
42	ahallam	3232	and thus Poisson's \autoref{eqn:poisson} satisfies the general form when
43	ahallam	3153	\begin{equation}
44			A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho
45			\end{equation}
46			At least one boundary point must be set for the problem to be solvable. For this
47			example we have set all of the boundaries to zero. The normal flux condition is
48	ahallam	3232	also zero by default. For a more realistic and complicated models it may be
49			necessary to give careful consideration to the boundary conditions of the model,
50			which can have an influence upon the solution.
51	ahallam	3153
52			Setting the boundary condition is relatively simple using the \verb!q! and
53	ahallam	3232	\verb!r! variables of the general form. First \verb!q! is defined as a masking
54			function on the boundary using
55	ahallam	3153	\begin{python}
56			q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
57			\end{python}
58	ahallam	3232	This identifies the points on the boundary and \verb!r! is simply
59			ser to \verb!r=0.0!. This is a dirichlet boundary condition.
60	ahallam	3153
61	ahallam	3232	\section{Gravity Pole}
62	ahallam	3153	\sslist{example10a.py}
63	ahallam	3232	A gravity pole is used in this example to demonstrate the radial directionality
64			of gravity, and also to show how this information can be exported for
65			visualisation to Mayavi or an equivalent using the VTK data format.
66	ahallam	3153
67	ahallam	3232	The solution script for this section is very simple. First the domain is
68			constructed, then the parameters of the model are set, and finally the steady
69			state solution is found. There are quite a few values that can now be derived
70			from the solution and saved to file for visualisation.
71
72			The potential $U$ is related to the gravitational response $\vec{g}$ via
73			\begin{equation}
74			\vec{g} = \nabla U
75			\end{equation}
76			This for example as a vertical component $g\hackscore{z}$ where
77			\begin{equation}
78			g\hackscore{z}=\vec{g}\cdot\hat{z}
79			\end{equation}
80			Finally, there is the magnitude of the vertical component $g$ of
81			$g\hackscore{z}$
82			\begin{equation}
83			g=\|g\hackscore{z}\|
84			\end{equation}
85			These values are derived from the \esc solution \verb!sol! to the potential $U$
86			using the following commands
87			\begin{python}
88			g_field=grad(sol) #The graviational accelleration g.
89			g_fieldz=g_field*[0,1] #The vertical component of the g field.
90			gz=length(g_fieldz) #The magnitude of the vertical component.
91			\end{python}
92			This data can now be simply exported to a VTK file via
93			\begin{python}
94			# Save the output to file.
95			saveVTK(os.path.join(save_path,"ex10a.vtu"),\
96			grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
97			\end{python}
98
99	ahallam	3153	\begin{figure}[ht]
100			\centering
101			\includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
102			\caption{Newtonian potential with g field directionality.}
103			\label{fig:ex10pot}
104			\end{figure}
105
106			\section{Gravity Well}
107			\sslist{example10b.py}
108	ahallam	3232	Let us now investigate the effect of gravity in three dimensions. Consider a
109			volume which contains a sphericle mass anomaly and a gravitational potential
110			which decays to zero at the base of the anomaly.
111	ahallam	3153
112	ahallam	3232	The script used to solve this model is very similar to that used for the gravity
113			pole in the previous section, but with an extra spatial dimension. As for all
114			the 3D problems examined in this cookbook, the extra dimension is easily
115			integrated into the \esc solultion script.
116
117			The domain is redefined from a rectangle to a Brick;
118			\begin{python}
119			domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz)
120			\end{python}
121			the source is relocated along $z$;
122			\begin{python}
123			x=x-[mx/2,my/2,mz/2]
124			\end{python}
125			and, the boundary conditions are updated.
126			\begin{python}
127			q=whereZero(x[2]-inf(x[2]))
128			\end{python}
129			No modifications to the PDE solution section are required. This make the
130			migration of 2D to a 3D problem almost trivial.
131
132			\autoref{fig:ex10bpot} illustrates the strength of a PDE solution. Three
133			different visualisation types help define and illuminate properties of the data.
134			The cut surfaces of the potential are similar to a 2D section for a given x or y
135			and z. The iso-surfaces illuminate the 3D shape of the gravity field, as well as
136			its strength which is give by the colour. Finally, the streamlines highlight the
137			directional flow of the gravity field in this example.
138
139	ahallam	3153	\begin{figure}[htp]
140			\centering
141			\includegraphics[width=0.75\textwidth]{figures/ex10bpot.png}
142			\caption{Gravity well with iso surfaces and streamlines of the vector
143			gravitational potential.}
144			\label{fig:ex10bpot}
145			\end{figure}
146
147			\section{Gravity Surface over a fault model.}
148			\sslist{example10c.py,example10m.py}
149	ahallam	3232	This model demonstrates the gravity result for a more complicated domain which
150			contains a fault. Additional information will be added when geophysical boundary
151			conditions for a gravity scenario have been established.
152
153	ahallam	3153	\begin{figure}[htp]
154			\centering
155			\subfigure[The geometry of the fault model in example10c.py.]
156			{\label{fig:ex10cgeo}
157			\includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\
158			\subfigure[The fault of interest from the fault model in
159			example10c.py.]
160			{\label{fig:ex10cmsh}
161			\includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}}
162			\end{figure}
163
164			\begin{figure}[htp]
165			\centering
166			\includegraphics[width=0.8\textwidth]{figures/ex10cpot.png}
167			\caption{Gravitational potential of the fault model with primary layers and
168			faults identified as isosurfaces.}
169			\label{fig:ex10cpot}
170			\end{figure}
171