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

Revision 3410 - (show annotations)
Thu Dec 9 07:02:52 2010 UTC (8 years, 8 months ago) by ahallam
File MIME type: application/x-tex
File size: 17831 byte(s)
Minor changes to cookbook, gravity example.

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