/[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 3135 - (show 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
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 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