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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26