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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26