/[escript]/trunk/doc/examples/cookbook/example10a.py
ViewVC logotype

Annotation of /trunk/doc/examples/cookbook/example10a.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3195 - (hide annotations)
Wed Sep 22 00:28:04 2010 UTC (8 years, 5 months ago) by ahallam
File MIME type: text/x-python
File size: 3922 byte(s)
Shortened runtime of cookbook examples to aid testing.
1 ahallam 3099
2     ########################################################
3     #
4     # Copyright (c) 2009 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
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     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="https://launchpad.net/escript-finley"
21    
22     """
23     Author: Antony Hallam antony.hallam@uqconnect.edu.au
24     """
25    
26     ############################################################FILE HEADER
27     # example10a.py
28 ahallam 3195 # Model of gravitational Potential for a gravity POLE.
29 ahallam 3099
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 ahallam 3135 from math import pi, sqrt, sin, cos
38 ahallam 3099
39 ahallam 3135 from esys.escript.pdetools import Projector
40     from cblib import toRegGrid
41 jfenwick 3148
42     import matplotlib
43     matplotlib.use('agg') #It's just here for automated testing
44    
45 ahallam 3135 import pylab as pl #Plotting package
46     import numpy as np
47    
48 ahallam 3099 ########################################################MPI WORLD CHECK
49     if getMPISizeWorld() > 1:
50     import sys
51     print "This example will not run in an MPI world."
52     sys.exit(0)
53    
54     #################################################ESTABLISHING VARIABLES
55     #Domain related.
56 ahallam 3135 mx = 5000*m #meters - model length
57     my = -5000*m #meters - model width
58     ndx = 100 # mesh steps in x direction
59     ndy = 100 # mesh steps in y direction - one dimension means one element
60 ahallam 3099 #PDE related
61 ahallam 3135 rho=200.0
62     rholoc=[2500,-2500]
63 ahallam 3099 G=6.67300*10E-11
64    
65     ################################################ESTABLISHING PARAMETERS
66     #the folder to put our outputs in, leave blank "" for script path
67 ahallam 3135 save_path= os.path.join("data","example10")
68 ahallam 3099 #ensure the dir exists
69     mkDir(save_path)
70    
71     ####################################################DOMAIN CONSTRUCTION
72     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
73     x=Solution(domain).getX()
74 ahallam 3135 mask=wherePositive(10-length(x-rholoc))
75 ahallam 3099 rho=rho*mask
76     kro=kronecker(domain)
77    
78 ahallam 3135 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
79 ahallam 3099 ###############################################ESCRIPT PDE CONSTRUCTION
80    
81     mypde=LinearPDE(domain)
82 ahallam 3135 mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
83     mypde.setSymmetryOn()
84 ahallam 3099 sol=mypde.getSolution()
85 ahallam 3135
86     g_field=grad(sol) #The graviational accelleration g.
87 ahallam 3195 g_fieldz=g_field*[0,1] #The vertical component of the g field.
88     gz=length(g_fieldz) #The magnitude of the vertical component.
89 ahallam 3135 # Save the output to file.
90     saveVTK(os.path.join(save_path,"ex10a.vtu"),\
91     grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
92    
93     ##################################################REGRIDDING & PLOTTING
94    
95    
96     xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
97     pl.matplotlib.pyplot.autumn()
98     pl.contourf(xi,yi,zi,10)
99     pl.xlabel("Horizontal Displacement (m)")
100     pl.ylabel("Depth (m)")
101     pl.savefig(os.path.join(save_path,"Ucontour.png"))
102     print "Solution has been plotted ..."
103    
104     cut=int(len(xi)/2)
105    
106     pl.clf()
107    
108     r=np.linspace(0,mx/2,100)
109     m=2*pl.pi*10*10*200*-G/(r*r)
110    
111     pl.plot(xi,zi[:,cut])
112     #pl.plot(r+2500,m)
113     pl.title("Potential Profile")
114     pl.xlabel("Horizontal Displacement (m)")
115     pl.ylabel("Potential")
116     pl.savefig(os.path.join(save_path,"Upot00.png"))
117    
118     out=np.array([xi,zi[:,cut]])
119     pl.savetxt('profile1.asc',out.transpose())
120     pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26