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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3195 - (show 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
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 # Model of gravitational Potential for a gravity POLE.
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
42 import matplotlib
43 matplotlib.use('agg') #It's just here for automated testing
44
45 import pylab as pl #Plotting package
46 import numpy as np
47
48 ########################################################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 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 #PDE related
61 rho=200.0
62 rholoc=[2500,-2500]
63 G=6.67300*10E-11
64
65 ################################################ESTABLISHING PARAMETERS
66 #the folder to put our outputs in, leave blank "" for script path
67 save_path= os.path.join("data","example10")
68 #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 mask=wherePositive(10-length(x-rholoc))
75 rho=rho*mask
76 kro=kronecker(domain)
77
78 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
79 ###############################################ESCRIPT PDE CONSTRUCTION
80
81 mypde=LinearPDE(domain)
82 mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
83 mypde.setSymmetryOn()
84 sol=mypde.getSolution()
85
86 g_field=grad(sol) #The graviational accelleration g.
87 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 # 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