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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3232 - (hide annotations)
Fri Oct 1 02:08:38 2010 UTC (8 years, 11 months ago) by ahallam
File MIME type: application/x-tex
File size: 6699 byte(s)
Updates to cookbook, new chapter DC Resis and IP. Using new packages hyperref, natbib
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    

  ViewVC Help
Powered by ViewVC 1.1.26