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

Revision 3308 - (show annotations)
Tue Oct 26 03:24:54 2010 UTC (10 years, 1 month ago) by jfenwick
File MIME type: application/x-tex
File size: 6619 byte(s)
Moved listing and class up a level.
Removed hackscore from cookbook.
Made cookbook use the new stuff - removed old stuff

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 % 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 14 \section{Newtonian Potential} 15 16 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 22 The gravitational potential $U$ at a point $P$ due to a region with a mass 23 distribution of density $\rho(P)$, is given by Poisson's equation 24 \citep{Blakely1995} 25 \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 Consider now the \esc general form, which can simply be related to 30 \autoref{eqn:poisson} using two coefficients. The result is 31 \begin{equation} 32 -\left(A_{jl} u_{,l} \right)_{,j} = Y 33 \end{equation} 34 one recognises that the LHS side is equivalent to 35 \begin{equation} \label{eqn:ex10a} 36 -\nabla A \nabla u 37 \end{equation} 38 If $A=\delta_{jl}$ then \autoref{eqn:ex10a} is equivalent to 39 \begin{equation*} 40 -\nabla^2 u 41 \end{equation*} 42 and thus Poisson's \autoref{eqn:poisson} satisfies the general form when 43 \begin{equation} 44 A=\delta_{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 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 53 \verb!r! variables of the general form. First \verb!q! is defined as a masking 54 function on the boundary using 55 \begin{python} 56 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx) 57 \end{python} 58 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 61 \section{Gravity Pole} 62 \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_{z}$ where 77 \begin{equation} 78 g_{z}=\vec{g}\cdot\hat{z} 79 \end{equation} 80 Finally, there is the magnitude of the vertical component $g$ of 81 $g_{z}$ 82 \begin{equation} 83 g=|g_{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] 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 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[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 \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 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] 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