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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6927 - (hide annotations)
Sun Dec 8 22:26:40 2019 UTC (11 months, 3 weeks ago) by acodd
File MIME type: application/x-tex
File size: 18131 byte(s)
user guide updates

1 ahallam 3153
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland
4 jfenwick 3989 % http://www.uq.edu.au
5 ahallam 3153 %
6     % Primary Business: Queensland, Australia
7 jfenwick 6112 % Licensed under the Apache License, version 2.0
8     % http://www.apache.org/licenses/LICENSE-2.0
9 ahallam 3153 %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 % Development 2012-2013 by School of Earth Sciences
12     % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3989 %
14     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 ahallam 3153
16     \section{Newtonian Potential}
17    
18 ahallam 3232 In this chapter the gravitational potential field is developed for \esc.
19     Gravitational fields are present in many modelling scenarios, including
20 ahallam 3410 geophysical investigations, planetary motion and attraction and micro-particle
21 ahallam 3373 interactions. Gravitational fields also present an opportunity to demonstrate
22 caltinay 5662 the saving and visualisation of vector data for \mayavi, and the construction of
23 ahallam 3410 variable sized meshes.
24 ahallam 3153
25     The gravitational potential $U$ at a point $P$ due to a region with a mass
26 ahallam 3232 distribution of density $\rho(P)$, is given by Poisson's equation
27     \citep{Blakely1995}
28 ahallam 3153 \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 ahallam 3410 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 ahallam 3153 \begin{equation}
36 jfenwick 3308 -\left(A_{jl} u_{,l} \right)_{,j} = Y
37 ahallam 3153 \end{equation}
38 acodd 6927 The LHS side is equivalent to
39 ahallam 3153 \begin{equation} \label{eqn:ex10a}
40 acodd 6927 -\nabla \cdot A \nabla u
41 ahallam 3153 \end{equation}
42 ahallam 3392 and when $A=\delta_{jl}$, \autoref{eqn:ex10a} is equivalent to
43 ahallam 3153 \begin{equation*}
44     -\nabla^2 u
45     \end{equation*}
46 ahallam 3373 Thus Poisson's \autoref{eqn:poisson} satisfies the general form when
47 ahallam 3153 \begin{equation}
48 jfenwick 3308 A=\delta_{jl} \text{ and } Y= 4\pi\gamma\rho
49 ahallam 3153 \end{equation}
50 ahallam 3373 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 ahallam 3410 of the potential is governed by an inverse-square distance law. Thus, in the
57 ahallam 3373 limit as $r$ increases;
58     \begin{equation}
59     \lim_{r\to\infty} U(P) = 0
60     \end{equation}
61 ahallam 3410 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 ahallam 3373
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 ahallam 3410 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 ahallam 3373
75 ahallam 3410 For this example we have set all of the boundaries to zero but only one boundary
76 ahallam 3373 point needs to be set for the problem to be solvable. The normal flux condition
77 ahallam 3410 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 ahallam 3232 which can have an influence upon the solution.
80 ahallam 3153
81     Setting the boundary condition is relatively simple using the \verb!q! and
82 ahallam 3232 \verb!r! variables of the general form. First \verb!q! is defined as a masking
83     function on the boundary using
84 ahallam 3153 \begin{python}
85     q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
86 ahallam 3410 mypde.setValue(q=q,r=0)
87 ahallam 3153 \end{python}
88 ahallam 3232 This identifies the points on the boundary and \verb!r! is simply
89 ahallam 3410 ser to \verb!r=0.0!. This is a Dirichlet boundary condition.
90 ahallam 3153
91 ahallam 3392 \clearpage
92 ahallam 3232 \section{Gravity Pole}
93 ahallam 3153 \sslist{example10a.py}
94 ahallam 3410 A gravity pole is used in this example to demonstrate the vector characteristics
95 ahallam 3392 of gravity, and also to demonstrate how this information can be exported for
96 caltinay 5662 visualisation to \mayavi or an equivalent using the VTK data format.
97 ahallam 3153
98 ahallam 3232 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 ahallam 3373 $\vec{g}$ is a vector and thus, has a a vertical component $g_{z}$ where
108 ahallam 3232 \begin{equation}
109 jfenwick 3308 g_{z}=\vec{g}\cdot\hat{z}
110 ahallam 3232 \end{equation}
111     Finally, there is the magnitude of the vertical component $g$ of
112 jfenwick 3308 $g_{z}$
113 ahallam 3232 \begin{equation}
114 jfenwick 3308 g=|g_{z}|
115 ahallam 3232 \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 ahallam 3410 g_field=grad(sol) #The gravitational acceleration g.
120 ahallam 3232 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 caltinay 5662 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 ahallam 3373 \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 caltinay 5662 object in the \mayavi navigation tree the 4 values saved to the VTK file are
135 ahallam 3373 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 ahallam 3410 object in the explorer tree and selecting the cell to point filter as in
143 ahallam 3373 \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 ahallam 3153 \begin{figure}[ht]
150     \centering
151 ahallam 3373 \includegraphics[width=0.75\textwidth]{figures/mayavi2_openfile.png}
152 caltinay 5662 \caption{Open a file in \mayavi}
153 ahallam 3373 \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 ahallam 3153 \includegraphics[width=0.75\textwidth]{figures/ex10apot.png}
173 ahallam 3392 \caption{Newtonian potential with $\vec{g}$ field directionality. The magnitude
174     of the field is reflected in the size of the vector arrows.}
175 ahallam 3153 \label{fig:ex10pot}
176     \end{figure}
177 ahallam 3373 \clearpage
178 ahallam 3153
179     \section{Gravity Well}
180     \sslist{example10b.py}
181 ahallam 3232 Let us now investigate the effect of gravity in three dimensions. Consider a
182 ahallam 3410 volume which contains a spherical mass anomaly and a gravitational potential
183 ahallam 3373 which decays to zero at the base of the model.
184 ahallam 3153
185 ahallam 3232 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 ahallam 3410 integrated into the \esc solution script.
189 ahallam 3232
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 ahallam 3410 No modifications to the PDE solution section are required. Thus the migration
203     from a 2D to a 3D problem is almost trivial.
204 ahallam 3232
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 ahallam 3410 its strength which is illustrated by the colour. Finally, the streamlines
210     highlight the directional flow of the gravity field in this example.
211 ahallam 3232
212 ahallam 3406 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 ahallam 3410 solution. The profiles in \autoref{fig:ex10p} further demonstrates the affect
219 ahallam 3406 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 ahallam 3410 appropriately dimensioned model for the problem, and and sampling location.
222 ahallam 3392
223 ahallam 3153 \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 ahallam 3392 gravitational potential \textemdash small model.}
228 ahallam 3153 \label{fig:ex10bpot}
229     \end{figure}
230    
231 ahallam 3392 \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 caltinay 4286 \caption{Profile of the gravitational profile along x where $y=0,z=250$ for
243 ahallam 3392 various sized domains.}
244     \label{fig:ex10p}
245     \end{figure}
246     \clearpage
247    
248 ahallam 3448 % \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 ahallam 3232
274 ahallam 3392 \section{Variable mesh-element sizes}
275 ahallam 3406 \sslist{example10m.py}
276     We saw in a previous section that the domain needed to be sufficiently
277 ahallam 3410 large for the boundary conditions to be satisfied. This can be troublesome when
278 ahallam 3406 trying to solve problems that require a dense mesh, either for solution
279 ahallam 3410 resolution of stability reasons. The computational cost of solving a large
280     number of elements can be prohibitive.S
281 ahallam 3406
282 ahallam 3392 With the help of Gmsh, it is possible to create a mesh for \esc, which has a
283 ahallam 3406 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 ahallam 3392
288     To create a variable size mesh, multiple meshing domains are required. The
289 ahallam 3410 domains must share points, boundaries and surfaces so that they are joined; and
290 ahallam 3392 no sub-domains are allowed to overlap. Whilst this initially seems complicated,
291     it is quite simple to implement.
292    
293 ahallam 3410 This example creates a mesh which contains a high resolution sub-domain at its
294 ahallam 3406 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 ahallam 3392 \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 ahallam 3410 The two sub-domains now have a common geometry and no over-laping features as
348 ahallam 3392 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 ahallam 3410 d1.addItems(PropertySet(l01,l12,l23,l30))
358 ahallam 3392 d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo"))
359 ahallam 3410 MakeDomain(d1)
360 ahallam 3392
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 ahallam 3410 MakeDomain(d2)
366 ahallam 3392 \end{python}
367 ahallam 3410 Finally, a system call to Gmsh is required to merge and then appropriately
368 ahallam 3392 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 ahallam 3410 To read in the file the function is called
387 ahallam 3392 \begin{python}
388     domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain
389     \end{python}
390 ahallam 3410 where the integer argument is the number of domain dimensions.
391 ahallam 3392 %
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 ahallam 3406 \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 ahallam 3410 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 ahallam 3448 \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 ahallam 3410
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 ahallam 3448 There is a methodology which can help establish an appropriate zero mass region
448     to a domain.
449 jfenwick 3989 \clearpage

  ViewVC Help
Powered by ViewVC 1.1.26