# Contents of /trunk/doc/examples/cookbook/example10d.py

Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4861 byte(s)
```Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2009 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 __copyright__="""Copyright (c) 2009 by University of Queensland 15 Earth Systems Science Computational Center (ESSCC) 16 http://www.uq.edu.au/esscc 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 __url__= 21 22 """ 23 Author: Antony Hallam antony.hallam@uqconnect.edu.au 24 """ 25 26 ############################################################FILE HEADER 27 # example10a.py 28 # Model of gravitational Potential for an infinite strike horizontal 29 # cylinder. This example tests the boundary condition criterion. 30 31 #######################################################EXTERNAL MODULES 32 # To solve the problem it is necessary to import the modules we require. 33 from esys.escript import * # This imports everything from the escript library 34 from esys.escript.unitsSI import * 35 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE 36 from esys.finley import Rectangle # This imports the rectangle domain function from finley 37 from esys.weipa import saveVTK # This imports the VTK file saver from weipa 38 import os, sys #This package is necessary to handle saving our data. 39 from math import pi, sqrt, sin, cos 40 41 from esys.escript.pdetools import Projector, Locator 42 from cblib import toRegGrid 43 from esys.finley import ReadGmsh 44 45 import matplotlib 46 matplotlib.use('agg') #It's just here for automated testing 47 48 import pylab as pl #Plotting package 49 import numpy as np 50 51 ########################################################MPI WORLD CHECK 52 if getMPISizeWorld() > 1: 53 import sys 54 print("This example will not run in an MPI world.") 55 sys.exit(0) 56 57 #################################################ESTABLISHING VARIABLES 58 #PDE related 59 rho=200.0 #the density contrast of the source 60 G=6.67300*10E-11 #gravitational constant 61 R=100. #radius of the source body 62 ################################################ESTABLISHING PARAMETERS 63 #the folder to put our outputs in, leave blank "" for script path 64 save_path= os.path.join("data","example10") 65 mesh_path = "data/example10m" 66 #ensure the dir exists 67 mkDir(save_path) 68 69 ####################################################DOMAIN CONSTRUCTION 70 #Geometric and material property related variables. 71 domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain 72 x=domain.getX() 73 #Domain related dimensions from the imported file. 74 mx = Lsup(x)*m #meters - model length 75 my = Lsup(x)*m #meters - model width 76 rholoc=[mx/2,my/2] 77 78 mask=wherePositive(R-length(x-rholoc)) 79 rhoe=rho*mask 80 kro=kronecker(domain) 81 82 #Define the analytic solution. 83 def analytic_gz(x,z,R,drho): 84 G=6.67300*10E-11 85 return G*2*np.pi*R*R*drho*(z/(x*x+z*z)) 86 87 #Define the boundary conditions. 88 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx) 89 ###############################################ESCRIPT PDE CONSTRUCTION 90 mypde=LinearPDE(domain) 91 mypde.setValue(A=kro,Y=4.*3.1415*G*rhoe,q=q,r=0.) 92 mypde.setSymmetryOn() 93 sol=mypde.getSolution() 94 95 g_field=grad(sol) #The gravitational acceleration g. 96 g_fieldz=g_field*[0,1] #The vertical component of the g field. 97 gz=length(g_fieldz) #The magnitude of the vertical component. 98 # Save the output to file. 99 saveVTK(os.path.join(save_path,"ex10d.vtu"),\ 100 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz,rho=rhoe) 101 102 ################################################MODEL SIZE SAMPLING 103 smoother = Projector(domain) #Function smoother. 104 z=1000. #Distance of profile from source. 105 sol_angz=[] #Array for analytic gz 106 sol_anx=[] #Array for x 107 #calculate analytic gz and x location. 108 for x in range(int(-mx/2.),int(mx/2.),10): 109 sol_angz.append(analytic_gz(x,z,R,rho)) 110 sol_anx.append(x+mx/2) 111 #save analytic solution 112 pl.savetxt(os.path.join(save_path,"ex10d_as.asc"),sol_angz) 113 114 sol_escgz=[] #Array for escript solution for gz 115 #calculate the location of the profile in the domain 116 for i in range(0,len(sol_anx)): 117 sol_escgz.append([sol_anx[i],my/2.-z]) 118 119 rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record 120 psol=rec.getValue(smoother(gz)) #extract the solution 121 #save X and escript solution profile to file 122 pl.savetxt(os.path.join(save_path,"ex10d_%05d.asc"%mx),psol) 123 pl.savetxt(os.path.join(save_path,"ex10d_x%05d.asc"%mx),sol_anx) 124 125 #plot the pofile and analytic solution using example10q.py 126

 ViewVC Help Powered by ViewVC 1.1.26