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

Revision 3135 - (show annotations)
Thu Sep 2 03:28:39 2010 UTC (11 years ago) by ahallam
File MIME type: text/x-python
File size: 3743 byte(s)
```Updates to cookbook scripts. Gravitation Examples.
```
 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. 29 30 #######################################################EXTERNAL MODULES 31 # To solve the problem it is necessary to import the modules we require. 32 from esys.escript import * # This imports everything from the escript library 33 from esys.escript.unitsSI import * 34 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE 35 from esys.finley import Rectangle # This imports the rectangle domain function from finley 36 import os, sys #This package is necessary to handle saving our data. 37 from math import pi, sqrt, sin, cos 38 39 from esys.escript.pdetools import Projector 40 from cblib import toRegGrid 41 import pylab as pl #Plotting package 42 import numpy as np 43 44 ########################################################MPI WORLD CHECK 45 if getMPISizeWorld() > 1: 46 import sys 47 print "This example will not run in an MPI world." 48 sys.exit(0) 49 50 #################################################ESTABLISHING VARIABLES 51 #Domain related. 52 mx = 5000*m #meters - model length 53 my = -5000*m #meters - model width 54 ndx = 100 # mesh steps in x direction 55 ndy = 100 # mesh steps in y direction - one dimension means one element 56 #PDE related 57 rho=200.0 58 rholoc=[2500,-2500] 59 G=6.67300*10E-11 60 61 ################################################ESTABLISHING PARAMETERS 62 #the folder to put our outputs in, leave blank "" for script path 63 save_path= os.path.join("data","example10") 64 #ensure the dir exists 65 mkDir(save_path) 66 67 ####################################################DOMAIN CONSTRUCTION 68 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 69 x=Solution(domain).getX() 70 mask=wherePositive(10-length(x-rholoc)) 71 rho=rho*mask 72 kro=kronecker(domain) 73 74 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx) 75 ###############################################ESCRIPT PDE CONSTRUCTION 76 77 mypde=LinearPDE(domain) 78 mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0) 79 mypde.setSymmetryOn() 80 sol=mypde.getSolution() 81 82 g_field=grad(sol) #The graviational accelleration g. 83 g_fieldz=g_field*[0,1] # 84 gz=length(g_fieldz) 85 # Save the output to file. 86 saveVTK(os.path.join(save_path,"ex10a.vtu"),\ 87 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz) 88 89 ##################################################REGRIDDING & PLOTTING 90 91 92 xi, yi, zi = toRegGrid(sol, nx=50, ny=50) 93 pl.matplotlib.pyplot.autumn() 94 pl.contourf(xi,yi,zi,10) 95 pl.xlabel("Horizontal Displacement (m)") 96 pl.ylabel("Depth (m)") 97 pl.savefig(os.path.join(save_path,"Ucontour.png")) 98 print "Solution has been plotted ..." 99 100 cut=int(len(xi)/2) 101 102 pl.clf() 103 104 r=np.linspace(0,mx/2,100) 105 m=2*pl.pi*10*10*200*-G/(r*r) 106 107 pl.plot(xi,zi[:,cut]) 108 #pl.plot(r+2500,m) 109 pl.title("Potential Profile") 110 pl.xlabel("Horizontal Displacement (m)") 111 pl.ylabel("Potential") 112 pl.savefig(os.path.join(save_path,"Upot00.png")) 113 114 out=np.array([xi,zi[:,cut]]) 115 pl.savetxt('profile1.asc',out.transpose()) 116 pl.clf()

 ViewVC Help Powered by ViewVC 1.1.26