/[escript]/trunk/doc/cookbook/example10.tex
ViewVC logotype

Annotation of /trunk/doc/cookbook/example10.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3373 - (hide annotations)
Tue Nov 23 00:29:07 2010 UTC (10 years ago) by ahallam
File MIME type: application/x-tex
File size: 9629 byte(s)
New figures.
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 ahallam 3373 interactions. Gravitational fields also present an opportunity to demonstrate
20 ahallam 3232 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 3373 Consider now the \esc general form, which has a simple relationship to
30     \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 ahallam 3153 \begin{equation}
33 jfenwick 3308 -\left(A_{jl} u_{,l} \right)_{,j} = Y
34 ahallam 3153 \end{equation}
35 ahallam 3373 One recognises that the LHS side is equivalent to
36 ahallam 3153 \begin{equation} \label{eqn:ex10a}
37     -\nabla A \nabla u
38     \end{equation}
39 ahallam 3373 and when a $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to
40 ahallam 3153 \begin{equation*}
41     -\nabla^2 u
42     \end{equation*}
43 ahallam 3373 Thus Poisson's \autoref{eqn:poisson} satisfies the general form when
44 ahallam 3153 \begin{equation}
45 jfenwick 3308 A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho
46 ahallam 3153 \end{equation}
47 ahallam 3373 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 inverse-square 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     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 ahallam 3232 necessary to give careful consideration to the boundary conditions of the model,
75     which can have an influence upon the solution.
76 ahallam 3153
77     Setting the boundary condition is relatively simple using the \verb!q! and
78 ahallam 3232 \verb!r! variables of the general form. First \verb!q! is defined as a masking
79     function on the boundary using
80 ahallam 3153 \begin{python}
81     q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
82     \end{python}
83 ahallam 3232 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 ahallam 3153
86 ahallam 3232 \section{Gravity Pole}
87 ahallam 3153 \sslist{example10a.py}
88 ahallam 3232 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 ahallam 3153
92 ahallam 3232 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 ahallam 3373 $\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where
102 ahallam 3232 \begin{equation}
103 jfenwick 3308 g_{z}=\vec{g}\cdot\hat{z}
104 ahallam 3232 \end{equation}
105     Finally, there is the magnitude of the vertical component $g$ of
106 jfenwick 3308 $g_{z}$
107 ahallam 3232 \begin{equation}
108 jfenwick 3308 g=|g_{z}|
109 ahallam 3232 \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 ahallam 3373 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 ahallam 3153 \begin{figure}[ht]
144     \centering
145 ahallam 3373 \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 ahallam 3153 \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
167     \caption{Newtonian potential with g field directionality.}
168     \label{fig:ex10pot}
169     \end{figure}
170 ahallam 3373 \clearpage
171 ahallam 3153
172     \section{Gravity Well}
173     \sslist{example10b.py}
174 ahallam 3232 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 ahallam 3373 which decays to zero at the base of the model.
177 ahallam 3153
178 ahallam 3232 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 iso-surfaces 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 ahallam 3153 \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 ahallam 3232 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 ahallam 3153 \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    

  ViewVC Help
Powered by ViewVC 1.1.26