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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3410 - (hide annotations)
Thu Dec 9 07:02:52 2010 UTC (8 years, 9 months ago) by ahallam
File MIME type: application/x-tex
File size: 17831 byte(s)
Minor changes to cookbook, gravity example.
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 ahallam 3410 geophysical investigations, planetary motion and attraction and micro-particle
19 ahallam 3373 interactions. Gravitational fields also present an opportunity to demonstrate
20 ahallam 3410 the saving and visualisation of vector data for Mayavi, and the construction of
21     variable sized meshes.
22 ahallam 3153
23     The gravitational potential $U$ at a point $P$ due to a region with a mass
24 ahallam 3232 distribution of density $\rho(P)$, is given by Poisson's equation
25     \citep{Blakely1995}
26 ahallam 3153 \begin{equation} \label{eqn:poisson}
27     \nabla^2 U(P) = -4\pi\gamma\rho(P)
28     \end{equation}
29     where $\gamma$ is the gravitational constant.
30 ahallam 3410 Consider now the \esc general form,
31     \autoref{eqn:poisson} requires only two coefficients,
32     $A$ and $Y$, thus the relevant terms of the general form are reduced to;
33 ahallam 3153 \begin{equation}
34 jfenwick 3308 -\left(A_{jl} u_{,l} \right)_{,j} = Y
35 ahallam 3153 \end{equation}
36 ahallam 3373 One recognises that the LHS side is equivalent to
37 ahallam 3153 \begin{equation} \label{eqn:ex10a}
38     -\nabla A \nabla u
39     \end{equation}
40 ahallam 3392 and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to
41 ahallam 3153 \begin{equation*}
42     -\nabla^2 u
43     \end{equation*}
44 ahallam 3373 Thus Poisson's \autoref{eqn:poisson} satisfies the general form when
45 ahallam 3153 \begin{equation}
46 jfenwick 3308 A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho
47 ahallam 3153 \end{equation}
48 ahallam 3373 The boundary condition is the last parameter that requires consideration. The
49     potential $U$ is related to the mass of a sphere by
50     \begin{equation}
51     U(P)=-\gamma \frac{m}{r^2}
52     \end{equation} where $m$ is the mass of the sphere and $r$ is the distance from
53     the center of the mass to the observation point $P$. Plainly, the magnitude
54 ahallam 3410 of the potential is governed by an inverse-square distance law. Thus, in the
55 ahallam 3373 limit as $r$ increases;
56     \begin{equation}
57     \lim_{r\to\infty} U(P) = 0
58     \end{equation}
59 ahallam 3410 Provided that the domain being solved is large enough, and the source mass is
60     contained within a central region of the domain, the potential will decay to
61     zero. This is a dirichlet boundary condition where $U=0$.
62 ahallam 3373
63     This boundary condition can be satisfied when there is some mass suspended in a
64     free-space. For geophysical models where the mass of interest may be an anomaly
65     inside a host rock, the anomaly can be isolated by subtracting the density of the
66 ahallam 3410 host rock from the model. This creates a fictitious free space model that will
67     satisfy the analytic boundary conditions. The result is that
68     $Y=4\pi\gamma\Delta\rho$, where $\Delta\rho=\rho-\rho_0$ and $\rho_0$ is the
69     baseline or host density. This of course means that the true gravity response is
70     not being modelled, but the response due to an anomalous mass with a
71     density contrast $\Delta\rho$.
72 ahallam 3373
73 ahallam 3410 For this example we have set all of the boundaries to zero but only one boundary
74 ahallam 3373 point needs to be set for the problem to be solvable. The normal flux condition
75 ahallam 3410 is also zero by default. Note, that for a more realistic and complicated models
76     it may be necessary to give careful consideration to the boundary conditions of the model,
77 ahallam 3232 which can have an influence upon the solution.
78 ahallam 3153
79     Setting the boundary condition is relatively simple using the \verb!q! and
80 ahallam 3232 \verb!r! variables of the general form. First \verb!q! is defined as a masking
81     function on the boundary using
82 ahallam 3153 \begin{python}
83     q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
84 ahallam 3410 mypde.setValue(q=q,r=0)
85 ahallam 3153 \end{python}
86 ahallam 3232 This identifies the points on the boundary and \verb!r! is simply
87 ahallam 3410 ser to \verb!r=0.0!. This is a Dirichlet boundary condition.
88 ahallam 3153
89 ahallam 3392 \clearpage
90 ahallam 3232 \section{Gravity Pole}
91 ahallam 3153 \sslist{example10a.py}
92 ahallam 3410 A gravity pole is used in this example to demonstrate the vector characteristics
93 ahallam 3392 of gravity, and also to demonstrate how this information can be exported for
94 ahallam 3410 visualisation to Mayavi or an equivalent using the VTK data format.
95 ahallam 3153
96 ahallam 3232 The solution script for this section is very simple. First the domain is
97     constructed, then the parameters of the model are set, and finally the steady
98     state solution is found. There are quite a few values that can now be derived
99     from the solution and saved to file for visualisation.
100    
101     The potential $U$ is related to the gravitational response $\vec{g}$ via
102     \begin{equation}
103     \vec{g} = \nabla U
104     \end{equation}
105 ahallam 3373 $\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where
106 ahallam 3232 \begin{equation}
107 jfenwick 3308 g_{z}=\vec{g}\cdot\hat{z}
108 ahallam 3232 \end{equation}
109     Finally, there is the magnitude of the vertical component $g$ of
110 jfenwick 3308 $g_{z}$
111 ahallam 3232 \begin{equation}
112 jfenwick 3308 g=|g_{z}|
113 ahallam 3232 \end{equation}
114     These values are derived from the \esc solution \verb!sol! to the potential $U$
115     using the following commands
116     \begin{python}
117 ahallam 3410 g_field=grad(sol) #The gravitational acceleration g.
118 ahallam 3232 g_fieldz=g_field*[0,1] #The vertical component of the g field.
119     gz=length(g_fieldz) #The magnitude of the vertical component.
120     \end{python}
121     This data can now be simply exported to a VTK file via
122     \begin{python}
123     # Save the output to file.
124     saveVTK(os.path.join(save_path,"ex10a.vtu"),\
125     grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
126     \end{python}
127    
128 ahallam 3373 It is quite simple to visualise the data from the gravity solution in Mayavi2.
129     With Mayavi2 open go to File, Load data, Open file \ldots as in
130     \autoref{fig:mayavi2openfile} and select the saved data file. The data will
131     have then been loaded and is ready for visualisation. Notice that under the data
132     object in the Mayavi2 navigation tree the 4 values saved to the VTK file are
133     available (\autoref{fig:mayavi2data}). There are two vector values,
134     \verb|gfield| and \verb|gfieldz|. Note that to plot both of these on the same
135     chart requires that the data object be imported twice.
136    
137     The point scalar data \verb|grav_pot| is the gravitational potential and it is
138     easily plotted using a surface module. To visualise the cell data a filter is
139     required that converts to point data. This is done by right clicking on the data
140 ahallam 3410 object in the explorer tree and selecting the cell to point filter as in
141 ahallam 3373 \autoref{fig:mayavi2cell2point}.
142    
143     The settings can then be modified to suit personal taste. An example of the
144     potential and gravitational field vectors is illustrated in
145     \autoref{fig:ex10pot}.
146    
147 ahallam 3153 \begin{figure}[ht]
148     \centering
149 ahallam 3373 \includegraphics[width=0.75\textwidth]{figures/mayavi2_openfile.png}
150     \caption{Open a file in Mayavi2}
151     \label{fig:mayavi2openfile}
152     \end{figure}
153    
154     \begin{figure}[ht]
155     \centering
156     \includegraphics[width=0.75\textwidth]{figures/mayavi2_data.png}
157     \caption{The 4 types of data in the imported VTK file.}
158     \label{fig:mayavi2data}
159     \end{figure}
160    
161     \begin{figure}[ht]
162     \centering
163     \includegraphics[width=0.75\textwidth]{figures/mayavi2_cell2point.png}
164     \caption{Converting cell data to point data.}
165     \label{fig:mayavi2cell2point}
166     \end{figure}
167    
168     \begin{figure}[ht]
169     \centering
170 ahallam 3153 \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
171 ahallam 3392 \caption{Newtonian potential with $\vec{g}$ field directionality. The magnitude
172     of the field is reflected in the size of the vector arrows.}
173 ahallam 3153 \label{fig:ex10pot}
174     \end{figure}
175 ahallam 3373 \clearpage
176 ahallam 3153
177     \section{Gravity Well}
178     \sslist{example10b.py}
179 ahallam 3232 Let us now investigate the effect of gravity in three dimensions. Consider a
180 ahallam 3410 volume which contains a spherical mass anomaly and a gravitational potential
181 ahallam 3373 which decays to zero at the base of the model.
182 ahallam 3153
183 ahallam 3232 The script used to solve this model is very similar to that used for the gravity
184     pole in the previous section, but with an extra spatial dimension. As for all
185     the 3D problems examined in this cookbook, the extra dimension is easily
186 ahallam 3410 integrated into the \esc solution script.
187 ahallam 3232
188     The domain is redefined from a rectangle to a Brick;
189     \begin{python}
190     domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz)
191     \end{python}
192     the source is relocated along $z$;
193     \begin{python}
194     x=x-[mx/2,my/2,mz/2]
195     \end{python}
196     and, the boundary conditions are updated.
197     \begin{python}
198     q=whereZero(x[2]-inf(x[2]))
199     \end{python}
200 ahallam 3410 No modifications to the PDE solution section are required. Thus the migration
201     from a 2D to a 3D problem is almost trivial.
202 ahallam 3232
203     \autoref{fig:ex10bpot} illustrates the strength of a PDE solution. Three
204     different visualisation types help define and illuminate properties of the data.
205     The cut surfaces of the potential are similar to a 2D section for a given x or y
206     and z. The iso-surfaces illuminate the 3D shape of the gravity field, as well as
207 ahallam 3410 its strength which is illustrated by the colour. Finally, the streamlines
208     highlight the directional flow of the gravity field in this example.
209 ahallam 3232
210 ahallam 3406 The boundary conditions were discussed previously, but not thoroughly
211     investigated. It was stated, that in the limit as the boundary becomes more
212     remote from the source, the potential will reduce to zero.
213     \autoref{fig:ex10bpot2} is the solution to the same gravity problem
214     but with a slightly larger domain. It is obvious in this case that
215     the previous domain size was too small to accurately represent the
216 ahallam 3410 solution. The profiles in \autoref{fig:ex10p} further demonstrates the affect
217 ahallam 3406 the domain size has upon the solution. As the domain size increases, the
218     profiles begin to converge. This convergence is a good indicator of an
219 ahallam 3410 appropriately dimensioned model for the problem, and and sampling location.
220 ahallam 3392
221 ahallam 3153 \begin{figure}[htp]
222     \centering
223     \includegraphics[width=0.75\textwidth]{figures/ex10bpot.png}
224     \caption{Gravity well with iso surfaces and streamlines of the vector
225 ahallam 3392 gravitational potential \textemdash small model.}
226 ahallam 3153 \label{fig:ex10bpot}
227     \end{figure}
228    
229 ahallam 3392 \begin{figure}[htp]
230     \centering
231     \includegraphics[width=0.75\textwidth]{figures/ex10bpot2.png}
232     \caption{Gravity well with iso surfaces and streamlines of the vector
233     gravitational potential \textemdash large model.}
234     \label{fig:ex10bpot2}
235     \end{figure}
236    
237     \begin{figure}[htp]
238     \centering
239     \includegraphics[width=0.85\textwidth]{figures/ex10p_boundeff.pdf}
240     \caption{Profile of the graviational provile along x where $y=0,z=250$ for
241     various sized domains.}
242     \label{fig:ex10p}
243     \end{figure}
244     \clearpage
245    
246 ahallam 3153 \section{Gravity Surface over a fault model.}
247     \sslist{example10c.py,example10m.py}
248 ahallam 3232 This model demonstrates the gravity result for a more complicated domain which
249     contains a fault. Additional information will be added when geophysical boundary
250     conditions for a gravity scenario have been established.
251    
252 ahallam 3153 \begin{figure}[htp]
253     \centering
254     \subfigure[The geometry of the fault model in example10c.py.]
255     {\label{fig:ex10cgeo}
256     \includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\
257     \subfigure[The fault of interest from the fault model in
258     example10c.py.]
259     {\label{fig:ex10cmsh}
260     \includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}}
261     \end{figure}
262    
263     \begin{figure}[htp]
264     \centering
265     \includegraphics[width=0.8\textwidth]{figures/ex10cpot.png}
266     \caption{Gravitational potential of the fault model with primary layers and
267     faults identified as isosurfaces.}
268     \label{fig:ex10cpot}
269     \end{figure}
270 ahallam 3392 \clearpage
271 ahallam 3153
272 ahallam 3392 \section{Variable mesh-element sizes}
273 ahallam 3406 \sslist{example10m.py}
274     We saw in a previous section that the domain needed to be sufficiently
275 ahallam 3410 large for the boundary conditions to be satisfied. This can be troublesome when
276 ahallam 3406 trying to solve problems that require a dense mesh, either for solution
277 ahallam 3410 resolution of stability reasons. The computational cost of solving a large
278     number of elements can be prohibitive.S
279 ahallam 3406
280 ahallam 3392 With the help of Gmsh, it is possible to create a mesh for \esc, which has a
281 ahallam 3406 variable element size. Such an approach can significantly reduce the number of
282     elements that need to be solved, and a large domain can be created that has
283     sufficient resolution close to the source and extends to distances large enough
284     that the infinity clause is satisfied.
285 ahallam 3392
286     To create a variable size mesh, multiple meshing domains are required. The
287 ahallam 3410 domains must share points, boundaries and surfaces so that they are joined; and
288 ahallam 3392 no sub-domains are allowed to overlap. Whilst this initially seems complicated,
289     it is quite simple to implement.
290    
291 ahallam 3410 This example creates a mesh which contains a high resolution sub-domain at its
292 ahallam 3406 center. We begin by defining two curve loops which describe the large or
293     big sub-domain and the smaller sub-domain which is to contain the high
294     resolution portion of the mesh.
295 ahallam 3392 \begin{python}
296     ################################################BIG DOMAIN
297     #ESTABLISHING PARAMETERS
298     width=10000. #width of model
299     depth=10000. #depth of model
300     bele_size=500. #big element size
301     #DOMAIN CONSTRUCTION
302     p0=Point(0.0, 0.0)
303     p1=Point(width, 0.0)
304     p2=Point(width, depth)
305     p3=Point(0.0, depth)
306     # Join corners in anti-clockwise manner.
307     l01=Line(p0, p1)
308     l12=Line(p1, p2)
309     l23=Line(p2, p3)
310     l30=Line(p3, p0)
311    
312     cbig=CurveLoop(l01,l12,l23,l30)
313    
314     ################################################SMALL DOMAIN
315     #ESTABLISHING PARAMETERS
316     xwidth=2000.0 #x width of model
317     zdepth=2000.0 #y width of model
318     sele_size=10. #small element size
319     #TRANSFORM
320     xshift=width/2-xwidth/2
321     zshift=depth/2-zdepth/2
322     #DOMAIN CONSTRUCTION
323     p4=Point(xshift, zshift)
324     p5=Point(xwidth+xshift, zshift)
325     p6=Point(xwidth+xshift, zdepth+zshift)
326     p7=Point(xshift, zdepth+zshift)
327     # Join corners in anti-clockwise manner.
328     l45=Line(p4, p5)
329     l56=Line(p5, p6)
330     l67=Line(p6, p7)
331     l74=Line(p7, p4)
332    
333     csmall=CurveLoop(l45,l56,l67,l74)
334     \end{python}
335     The small sub-domain curve can then be used to create a surface.
336     \begin{python}
337     ssmall=PlaneSurface(csmall)
338     \end{python}
339     However, so that the two domains do not overlap, when the big sub-domain
340     curveloop is used to create a surface it must contain a hole. The hole is
341     defined by the small sub-domain curveloop.
342     \begin{python}
343     sbig=PlaneSurface(cbig,holes=[csmall])
344     \end{python}
345 ahallam 3410 The two sub-domains now have a common geometry and no over-laping features as
346 ahallam 3392 per \autoref{fig:ex10mgeo}. Notice, that both domains have a normal in the
347     same direction.
348    
349     The next step, is exporting each sub-domain individually, with an appropriate
350     element size. This is carried out using the \pycad Design command.
351     \begin{python}
352     # Design the geometry for the big mesh.
353     d1=Design(dim=2, element_size=bele_size, order=1)
354     d1.addItems(sbig)
355 ahallam 3410 d1.addItems(PropertySet(l01,l12,l23,l30))
356 ahallam 3392 d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo"))
357 ahallam 3410 MakeDomain(d1)
358 ahallam 3392
359     # Design the geometry for the small mesh.
360     d2=Design(dim=2, element_size=sele_size, order=1)
361     d2.addItems(ssmall)
362     d2.setScriptFileName(os.path.join(save_path,"example10m_small.geo"))
363 ahallam 3410 MakeDomain(d2)
364 ahallam 3392 \end{python}
365 ahallam 3410 Finally, a system call to Gmsh is required to merge and then appropriately
366 ahallam 3392 mesh the two sub-domains together.
367     \begin{python}
368     # Join the two meshes using Gmsh and then apply a 2D meshing algorithm.
369     # The small mesh must come before the big mesh in the merging call!!@!!@!
370     sp.call("gmsh -2 "+
371     os.path.join(save_path,"example10m_small.geo")+" "+
372     os.path.join(save_path,"example10m_big.geo")+" -o "+
373     os.path.join(save_path,"example10m.msh"),shell=True)
374     \end{python}
375     The ``-2'' option is responsible for the 2D meshing, and the ``-o'' option
376     provides the output path. The resulting mesh is depicted in
377     \autoref{fig:ex10mmsh}
378    
379     To use the Gmsh ``*.msh'' file in the solution script, the mesh reading function
380     ``ReadGmsh'' is required. It can be imported via;
381     \begin{python}
382     from esys.finley import ReadGmsh
383     \end{python}
384 ahallam 3410 To read in the file the function is called
385 ahallam 3392 \begin{python}
386     domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain
387     \end{python}
388 ahallam 3410 where the integer argument is the number of domain dimensions.
389 ahallam 3392 %
390     \begin{figure}[ht]
391     \centering
392     \includegraphics[width=0.8\textwidth]{figures/ex10m_geo.png}
393     \caption{Geometry of two surfaces for a single domain.}
394     \label{fig:ex10mgeo}
395     \end{figure}
396    
397     \begin{figure}[ht]
398     \centering
399     \includegraphics[width=0.8\textwidth]{figures/ex10m_msh.png}
400     \caption{Mesh of merged surfaces, showing variable element size. Elements
401     range from 10m in the centroid to 500m at the boundary.}
402     \label{fig:ex10mmsh}
403 ahallam 3406 \end{figure}
404     \clearpage
405    
406     \section{Unbounded problems}
407     With a variable element-size, it is now possible to solve the potential problem
408 ahallam 3410 over a very large mesh. To test the accuracy of the solution, we will compare
409     the \esc result with the analytic solution for the vertical gravitational
410     acceleration $g_z$ of an infinite horizontal cylinder.
411    
412     For a horizontal cylinder with a circular cross-section with infinite strike,
413     the analytic solution is give by
414     \begin{equation}
415     g_z = 2\gamma\pi R^2 \Delta\rho \frac{z}{(x^2+z^2)}
416     \end{equation}
417     where $\gamma$ is the gravitational constant (as defined previously), $R$ is the
418     radius of the cylinder, $\Delta\rho$ is the density contrast and $x$ and $z$ are
419     geometric factors, relating the observation point to the center of the source
420     via the horizontal and vertical displacements respectively.
421    
422     The accuracy of the solution was tested using a square domain. For each test the
423     dimensions of the domain were modified, being set to 5, 10, 20 and 40 Km. The
424     results are compared with the analytic solution and are depicted in
425     \autoref{fig:ex10q boundeff} and \autoref{fig:ex10q boundeff zoom}. Clearly,
426     as the domain size increases, the results are valid at greater
427     distances from the source. The same is true at the anomaly peak, where the
428     variation around the source diminishes with an increasing domain size.
429    
430     \begin{figure}[ht]
431     \centering
432     \includegraphics[width=0.8\textwidth]{figures/ex10q_boundeff.pdf}
433     \caption{Solution profile 1000.0 meters from the source as the domain size
434     increases.}
435     \label{fig:ex10q boundeff}
436     \end{figure}
437    
438     \begin{figure}[ht]
439     \centering
440     \includegraphics[width=0.8\textwidth]{figures/ex10q_boundeff_zoom.pdf}
441     \caption{Magnification of \autoref{fig:ex10q boundeff}.}
442     \label{fig:ex10q boundeff zoom}
443     \end{figure}
444     \clearpage
445    
446     \subsection{Inversion using scipy}

  ViewVC Help
Powered by ViewVC 1.1.26