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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 18145 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26