# Diff of /trunk/doc/cookbook/example10.tex

revision 3233 by jfenwick, Mon Oct 4 00:54:26 2010 UTC revision 3373 by ahallam, Tue Nov 23 00:29:07 2010 UTC
# Line 16  Line 16
16  In this chapter the gravitational potential field is developed for \esc.  In this chapter the gravitational potential field is developed for \esc.
17  Gravitational fields are present in many modelling scenarios, including  Gravitational fields are present in many modelling scenarios, including
18  geophysical investigations, planetary motion and attraction and microparticle  geophysical investigations, planetary motion and attraction and microparticle
19  interactions. Gravitational fields also presents the opportunity to demonstrate  interactions. Gravitational fields also present an opportunity to demonstrate
20  the saving and visualisation of vector data for Mayavi.  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  The gravitational potential $U$ at a point $P$ due to a region with a mass
# Line 26  distribution of density $\rho(P)$, is gi Line 26  distribution of density $\rho(P)$, is gi
26  \nabla^2 U(P) = -4\pi\gamma\rho(P)  \nabla^2 U(P) = -4\pi\gamma\rho(P)
27
28  where $\gamma$ is the gravitational constant.  where $\gamma$ is the gravitational constant.
29  Consider now the \esc general form, which can simply be related to  Consider now the \esc general form, which has a simple relationship to
30  \autoref{eqn:poisson} using two coefficients. The result is  \autoref{eqn:poisson}. There are only two non-zero coefficients, $A$ and $Y$,
31    thus the relevant terms of the general form are reduced to;
32
33  -\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y  -\left(A_{jl} u_{,l} \right)_{,j} = Y
34
35  one recognises that the LHS side is equivalent to  One recognises that the LHS side is equivalent to
36   \label{eqn:ex10a}   \label{eqn:ex10a}
37  -\nabla A \nabla u  -\nabla A \nabla u
38
39  If $A=\delta\hackscore{jl}$ then \autoref{eqn:ex10a} is equivalent to  and when a $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to
40  \begin{equation*}  \begin{equation*}
41  -\nabla^2 u  -\nabla^2 u
42  \end{equation*}  \end{equation*}
43  and thus Poisson's \autoref{eqn:poisson} satisfies the general form when  Thus Poisson's \autoref{eqn:poisson} satisfies the general form when
44
45  A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho  A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho
46
47  At least one boundary point must be set for the problem to be solvable. For this  The boundary condition is the last parameter that requires consideration. The
48  example we have set all of the boundaries to zero. The normal flux condition is  potential $U$ is related to the mass of a sphere by
49  also zero by default. For a more realistic and complicated models it may be
50    U(P)=-\gamma \frac{m}{r^2}
51     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 inverse-square distance law. Thus, in the
54    limit as $r$ increases;
55
56    \lim_{r\to\infty} U(P) = 0
57
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    free-space. 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,  necessary to give careful consideration to the boundary conditions of the model,
75  which can have an influence upon the solution.  which can have an influence upon the solution.
76
# Line 73  The potential $U$ is related to the grav Line 98  The potential $U$ is related to the grav
98
99  \vec{g} = \nabla U  \vec{g} = \nabla U
100
101  This for example as a vertical component $g\hackscore{z}$ where  $\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where
102
103  g\hackscore{z}=\vec{g}\cdot\hat{z}  g_{z}=\vec{g}\cdot\hat{z}
104
105  Finally, there is the magnitude of the vertical component $g$ of  Finally, there is the magnitude of the vertical component $g$ of
106  $g\hackscore{z}$  $g_{z}$
107
108  g=|g\hackscore{z}|  g=|g_{z}|
109
110  These values are derived from the \esc solution \verb!sol! to the potential $U$  These values are derived from the \esc solution \verb!sol! to the potential $U$
111  using the following commands  using the following commands
# Line 96  saveVTK(os.path.join(save_path,"ex10a.vt Line 121  saveVTK(os.path.join(save_path,"ex10a.vt
121          grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)          grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
122  \end{python}  \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    \verb|gfield| and \verb|gfieldz|. 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 \verb|grav_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]  \begin{figure}[ht]
165  \centering  \centering
166  \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}  \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
167  \caption{Newtonian potential with g field directionality.}  \caption{Newtonian potential with g field directionality.}
168  \label{fig:ex10pot}  \label{fig:ex10pot}
169  \end{figure}  \end{figure}
170    \clearpage
171
172  \section{Gravity Well}  \section{Gravity Well}
173  \sslist{example10b.py}  \sslist{example10b.py}
174  Let us now investigate the effect of gravity in three dimensions. Consider a  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  volume which contains a sphericle mass anomaly and a gravitational potential
176  which decays to zero at the base of the anomaly.  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  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  pole in the previous section, but with an extra spatial dimension. As for all

Legend:
 Removed from v.3233 changed lines Added in v.3373