/[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 3135 - (hide annotations)
Thu Sep 2 03:28:39 2010 UTC (8 years, 5 months ago) by ahallam
File MIME type: text/x-python
File size: 3743 byte(s)
Updates to cookbook scripts. Gravitation Examples.
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     # 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 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     import pylab as pl #Plotting package
42     import numpy as np
43    
44 ahallam 3099 ########################################################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 ahallam 3135 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 ahallam 3099 #PDE related
57 ahallam 3135 rho=200.0
58     rholoc=[2500,-2500]
59 ahallam 3099 G=6.67300*10E-11
60    
61     ################################################ESTABLISHING PARAMETERS
62     #the folder to put our outputs in, leave blank "" for script path
63 ahallam 3135 save_path= os.path.join("data","example10")
64 ahallam 3099 #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 ahallam 3135 mask=wherePositive(10-length(x-rholoc))
71 ahallam 3099 rho=rho*mask
72     kro=kronecker(domain)
73    
74 ahallam 3135 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
75 ahallam 3099 ###############################################ESCRIPT PDE CONSTRUCTION
76    
77     mypde=LinearPDE(domain)
78 ahallam 3135 mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
79     mypde.setSymmetryOn()
80 ahallam 3099 sol=mypde.getSolution()
81 ahallam 3135
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