 # Contents of /branches/doubleplusgood/doc/cookbook/example10.tex

Revision 4345 - (show annotations)
Fri Mar 29 07:09:41 2013 UTC (5 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 18081 byte(s)
Spelling fixes

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