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

Revision 3406 - (hide annotations)
Wed Dec 8 01:09:37 2010 UTC (8 years, 9 months ago) by ahallam
File MIME type: application/x-tex
File size: 15863 byte(s)
Updates to cookbook documentation for release.

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