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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3392 - (show annotations)
Thu Dec 2 03:26:29 2010 UTC (8 years, 9 months ago) by ahallam
File MIME type: application/x-tex
File size: 15063 byte(s)
Updates to cookbook. Includes new section on variable meshing. To be reviewed.
1
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 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 present an 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 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 \begin{equation}
33 -\left(A_{jl} u_{,l} \right)_{,j} = Y
34 \end{equation}
35 One recognises that the LHS side is equivalent to
36 \begin{equation} \label{eqn:ex10a}
37 -\nabla A \nabla u
38 \end{equation}
39 and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to
40 \begin{equation*}
41 -\nabla^2 u
42 \end{equation*}
43 Thus Poisson's \autoref{eqn:poisson} satisfies the general form when
44 \begin{equation}
45 A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho
46 \end{equation}
47 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 necessary to give careful consideration to the boundary conditions of the model,
75 which can have an influence upon the solution.
76
77 Setting the boundary condition is relatively simple using the \verb!q! and
78 \verb!r! variables of the general form. First \verb!q! is defined as a masking
79 function on the boundary using
80 \begin{python}
81 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
82 \end{python}
83 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
86 \clearpage
87 \section{Gravity Pole}
88 \sslist{example10a.py}
89 A gravity pole is used in this example to demonstrate the radial directionality
90 of gravity, and also to demonstrate how this information can be exported for
91 visualisation to Mayavi or an equivalent using the VTK data format.
92
93 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 $\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where
103 \begin{equation}
104 g_{z}=\vec{g}\cdot\hat{z}
105 \end{equation}
106 Finally, there is the magnitude of the vertical component $g$ of
107 $g_{z}$
108 \begin{equation}
109 g=|g_{z}|
110 \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 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 \begin{figure}[ht]
145 \centering
146 \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 \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
168 \caption{Newtonian potential with $\vec{g}$ field directionality. The magnitude
169 of the field is reflected in the size of the vector arrows.}
170 \label{fig:ex10pot}
171 \end{figure}
172 \clearpage
173
174 \section{Gravity Well}
175 \sslist{example10b.py}
176 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 which decays to zero at the base of the model.
179
180 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 migration of 2D to a 3D problem almost trivial.
199
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 \autoref{fig:ex10bpot2} and \autoref{fig:ex10p} highlight the importance of the
208 boundary conditions and the size of the model. The results are clearly very
209 different between the two model sizes and the profile shows the effect of model
210 size for various models.
211
212 \begin{figure}[htp]
213 \centering
214 \includegraphics[width=0.75\textwidth]{figures/ex10bpot.png}
215 \caption{Gravity well with iso surfaces and streamlines of the vector
216 gravitational potential \textemdash small model.}
217 \label{fig:ex10bpot}
218 \end{figure}
219
220 \begin{figure}[htp]
221 \centering
222 \includegraphics[width=0.75\textwidth]{figures/ex10bpot2.png}
223 \caption{Gravity well with iso surfaces and streamlines of the vector
224 gravitational potential \textemdash large model.}
225 \label{fig:ex10bpot2}
226 \end{figure}
227
228 \begin{figure}[htp]
229 \centering
230 \includegraphics[width=0.85\textwidth]{figures/ex10p_boundeff.pdf}
231 \caption{Profile of the graviational provile along x where $y=0,z=250$ for
232 various sized domains.}
233 \label{fig:ex10p}
234 \end{figure}
235 \clearpage
236
237 \section{Gravity Surface over a fault model.}
238 \sslist{example10c.py,example10m.py}
239 This model demonstrates the gravity result for a more complicated domain which
240 contains a fault. Additional information will be added when geophysical boundary
241 conditions for a gravity scenario have been established.
242
243 \begin{figure}[htp]
244 \centering
245 \subfigure[The geometry of the fault model in example10c.py.]
246 {\label{fig:ex10cgeo}
247 \includegraphics[width=0.8\textwidth]{figures/ex10potfaultgeo.png}} \\
248 \subfigure[The fault of interest from the fault model in
249 example10c.py.]
250 {\label{fig:ex10cmsh}
251 \includegraphics[width=0.8\textwidth]{figures/ex10potfaultmsh.png}}
252 \end{figure}
253
254 \begin{figure}[htp]
255 \centering
256 \includegraphics[width=0.8\textwidth]{figures/ex10cpot.png}
257 \caption{Gravitational potential of the fault model with primary layers and
258 faults identified as isosurfaces.}
259 \label{fig:ex10cpot}
260 \end{figure}
261 \clearpage
262
263 \section{Variable mesh-element sizes}
264 \sslist{example10m.py and example10d.py}
265 With the help of Gmsh, it is possible to create a mesh for \esc, which has a
266 variable element size. This will be extremely useful for the gravitational
267 potential, and trying to impose appropriate boundary conditions.
268
269 A large domain can be created that has sufficient resultion close to the source
270 and extends to distances large enough that the infinity clause is satisfied,
271 without the requirement for large number of cells.
272
273 To create a variable size mesh, multiple meshing domains are required. The
274 domains mush share points, boundaries and surfaces so that they are joined; and
275 no sub-domains are allowed to overlap. Whilst this initially seems complicated,
276 it is quite simple to implement.
277
278 This example creates a mesh which contains a high resolution sub-domain in its
279 center. We begin by defining two curve loops which describe the large of
280 big-subdomain and the smaller sub-domain which is to contain the high resolution
281 portion of the mesh.
282 \begin{python}
283 ################################################BIG DOMAIN
284 #ESTABLISHING PARAMETERS
285 width=10000. #width of model
286 depth=10000. #depth of model
287 bele_size=500. #big element size
288 #DOMAIN CONSTRUCTION
289 p0=Point(0.0, 0.0)
290 p1=Point(width, 0.0)
291 p2=Point(width, depth)
292 p3=Point(0.0, depth)
293 # Join corners in anti-clockwise manner.
294 l01=Line(p0, p1)
295 l12=Line(p1, p2)
296 l23=Line(p2, p3)
297 l30=Line(p3, p0)
298
299 cbig=CurveLoop(l01,l12,l23,l30)
300
301 ################################################SMALL DOMAIN
302 #ESTABLISHING PARAMETERS
303 xwidth=2000.0 #x width of model
304 zdepth=2000.0 #y width of model
305 sele_size=10. #small element size
306 #TRANSFORM
307 xshift=width/2-xwidth/2
308 zshift=depth/2-zdepth/2
309 #DOMAIN CONSTRUCTION
310 p4=Point(xshift, zshift)
311 p5=Point(xwidth+xshift, zshift)
312 p6=Point(xwidth+xshift, zdepth+zshift)
313 p7=Point(xshift, zdepth+zshift)
314 # Join corners in anti-clockwise manner.
315 l45=Line(p4, p5)
316 l56=Line(p5, p6)
317 l67=Line(p6, p7)
318 l74=Line(p7, p4)
319
320 csmall=CurveLoop(l45,l56,l67,l74)
321 \end{python}
322 The small sub-domain curve can then be used to create a surface.
323 \begin{python}
324 ssmall=PlaneSurface(csmall)
325 \end{python}
326 However, so that the two domains do not overlap, when the big sub-domain
327 curveloop is used to create a surface it must contain a hole. The hole is
328 defined by the small sub-domain curveloop.
329 \begin{python}
330 sbig=PlaneSurface(cbig,holes=[csmall])
331 \end{python}
332 The two sub-domains now have a common geometry and no overlaping features as
333 per \autoref{fig:ex10mgeo}. Notice, that both domains have a normal in the
334 same direction.
335
336 The next step, is exporting each sub-domain individually, with an appropriate
337 element size. This is carried out using the \pycad Design command.
338 \begin{python}
339 # Design the geometry for the big mesh.
340 d1=Design(dim=2, element_size=bele_size, order=1)
341 d1.addItems(sbig)
342 d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo"))
343
344 # Design the geometry for the small mesh.
345 d2=Design(dim=2, element_size=sele_size, order=1)
346 d2.addItems(ssmall)
347 d2.setScriptFileName(os.path.join(save_path,"example10m_small.geo"))
348 \end{python}
349 Finally, a system call to Gmsh is required to merge and then appropriate
350 mesh the two sub-domains together.
351 \begin{python}
352 # Join the two meshes using Gmsh and then apply a 2D meshing algorithm.
353 # The small mesh must come before the big mesh in the merging call!!@!!@!
354 sp.call("gmsh -2 "+
355 os.path.join(save_path,"example10m_small.geo")+" "+
356 os.path.join(save_path,"example10m_big.geo")+" -o "+
357 os.path.join(save_path,"example10m.msh"),shell=True)
358 \end{python}
359 The ``-2'' option is responsible for the 2D meshing, and the ``-o'' option
360 provides the output path. The resulting mesh is depicted in
361 \autoref{fig:ex10mmsh}
362
363 To use the Gmsh ``*.msh'' file in the solution script, the mesh reading function
364 ``ReadGmsh'' is required. It can be imported via;
365 \begin{python}
366 from esys.finley import ReadGmsh
367 \end{python}
368 To read in the file the function is called like
369 \begin{python}
370 domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain
371 \end{python}
372 %
373 \begin{figure}[ht]
374 \centering
375 \includegraphics[width=0.8\textwidth]{figures/ex10m_geo.png}
376 \caption{Geometry of two surfaces for a single domain.}
377 \label{fig:ex10mgeo}
378 \end{figure}
379
380 \begin{figure}[ht]
381 \centering
382 \includegraphics[width=0.8\textwidth]{figures/ex10m_msh.png}
383 \caption{Mesh of merged surfaces, showing variable element size. Elements
384 range from 10m in the centroid to 500m at the boundary.}
385 \label{fig:ex10mmsh}
386 \end{figure}

  ViewVC Help
Powered by ViewVC 1.1.26