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

revision 3153 by ahallam, Sun Sep 5 01:31:49 2010 UTC revision 3232 by ahallam, Fri Oct 1 02:08:38 2010 UTC
# Line 13  Line 13
13
14  \section{Newtonian Potential}  \section{Newtonian Potential}
15
16  In this chapter we examine the gravitational potential field due to source  In this chapter the gravitational potential field is developed for \esc.
17  bodies and geological structure. The use of vector data and visualisation in  Gravitational fields are present in many modelling scenarios, including
18  Mayavi is investigated.  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
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
23  distribution $\rho(P)$ is given by Poisson's equation  distribution of density $\rho(P)$, is given by Poisson's equation
24    \citep{Blakely1995}
25  \begin{equation} \label{eqn:poisson}  \begin{equation} \label{eqn:poisson}
26  \nabla^2 U(P) = -4\pi\gamma\rho(P)  \nabla^2 U(P) = -4\pi\gamma\rho(P)
27  \end{equation}  \end{equation}
28  where $\gamma$ is the gravitational constant.  where $\gamma$ is the gravitational constant.
29  If this is now related to the simplified general form  Consider now the \esc general form, which can simply be related to
30    \auoref{eqn:poisson} using two coefficients. The result is
31  \begin{equation}  \begin{equation}
32  -\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y  -\left(A\hackscore{jl} u\hackscore{,l} \right)\hackscore{,j} = Y
33  \end{equation}  \end{equation}
34  where the LHS side is equivalent to  one recognises that the LHS side is equivalent to
35  \begin{equation} \label{eqn:ex10a}  \begin{equation} \label{eqn:ex10a}
36  -\nabla A \nabla u  -\nabla A \nabla u
37  \end{equation}  \end{equation}
38  If $A=\delta\hackscore{jl}$ then equation \ref{eqn:ex10a} is equivalent to  If $A=\delta\hackscore{jl}$ then \autoref{eqn:ex10a} is equivalent to
39  \begin{equation*}  \begin{equation*}
40  -\nabla^2 u  -\nabla^2 u
41  \end{equation*}  \end{equation*}
42  and thus Poisson's Equation \ref{eqn:poisson} satisfies the general form when  and thus Poisson's \autoref{eqn:poisson} satisfies the general form when
43  \begin{equation}  \begin{equation}
44  A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho  A=\delta\hackscore{jl} \text{ and } Y= 4\pi\gamma\rho
45  \end{equation}  \end{equation}
46  At least one boundary point must be set for the problem to be solvable. For this  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  example we have set all of the boundaries to zero. The normal flux condition is
48  also zero by default.  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
52  Setting the boundary condition is relatively simple using the \verb!q! and  Setting the boundary condition is relatively simple using the \verb!q! and
53  \verb!r! variables of the general form, where  \verb!r! variables of the general form. First \verb!q! is defined as a masking
54    function on the boundary using
55  \begin{python}  \begin{python}
56  q=whereZero(x-my)+whereZero(x)+whereZero(x)+whereZero(x-mx)  q=whereZero(x-my)+whereZero(x)+whereZero(x)+whereZero(x-mx)
57  \end{python}  \end{python}
58  identifies the points on the boundary and \verb!r=0.0! sets the dirichlet  This identifies the points on the boundary and \verb!r! is simply
59  condition.  ser to \verb!r=0.0!. This is a dirichlet boundary condition.
60
61  \section{Gravity Profile}  \section{Gravity Pole}
62  \sslist{example10a.py}  \sslist{example10a.py}
63    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
67    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  \begin{figure}[ht]  \begin{figure}[ht]
100  \centering  \centering
# Line 63  condition. Line 105  condition.
105
106  \section{Gravity Well}  \section{Gravity Well}
107  \sslist{example10b.py}  \sslist{example10b.py}
108    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
112    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-inf(x))
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  \begin{figure}[htp]  \begin{figure}[htp]
140  \centering  \centering
# Line 74  gravitational potential.} Line 146  gravitational potential.}
146
147  \section{Gravity Surface over a fault model.}  \section{Gravity Surface over a fault model.}
148  \sslist{example10c.py,example10m.py}  \sslist{example10c.py,example10m.py}
149    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  \begin{figure}[htp]  \begin{figure}[htp]
154  \centering  \centering
155  \subfigure[The geometry of the fault model in example10c.py.]  \subfigure[The geometry of the fault model in example10c.py.]

Legend:
 Removed from v.3153 changed lines Added in v.3232