/[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 3346 - (show annotations)
Fri Nov 12 01:19:02 2010 UTC (8 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 3998 byte(s)
Replaced usage of esys.escript.util.saveVTK by weipa.saveVTK in all python
scripts.

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 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
37 import os, sys #This package is necessary to handle saving our data.
38 from math import pi, sqrt, sin, cos
39
40 from esys.escript.pdetools import Projector
41 from cblib import toRegGrid
42
43 import matplotlib
44 matplotlib.use('agg') #It's just here for automated testing
45
46 import pylab as pl #Plotting package
47 import numpy as np
48
49 ########################################################MPI WORLD CHECK
50 if getMPISizeWorld() > 1:
51 import sys
52 print "This example will not run in an MPI world."
53 sys.exit(0)
54
55 #################################################ESTABLISHING VARIABLES
56 #Domain related.
57 mx = 5000*m #meters - model length
58 my = -5000*m #meters - model width
59 ndx = 100 # mesh steps in x direction
60 ndy = 100 # mesh steps in y direction - one dimension means one element
61 #PDE related
62 rho=200.0
63 rholoc=[2500,-2500]
64 G=6.67300*10E-11
65
66 ################################################ESTABLISHING PARAMETERS
67 #the folder to put our outputs in, leave blank "" for script path
68 save_path= os.path.join("data","example10")
69 #ensure the dir exists
70 mkDir(save_path)
71
72 ####################################################DOMAIN CONSTRUCTION
73 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
74 x=Solution(domain).getX()
75 mask=wherePositive(10-length(x-rholoc))
76 rho=rho*mask
77 kro=kronecker(domain)
78
79 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
80 ###############################################ESCRIPT PDE CONSTRUCTION
81
82 mypde=LinearPDE(domain)
83 mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
84 mypde.setSymmetryOn()
85 sol=mypde.getSolution()
86
87 g_field=grad(sol) #The graviational accelleration g.
88 g_fieldz=g_field*[0,1] #The vertical component of the g field.
89 gz=length(g_fieldz) #The magnitude of the vertical component.
90 # Save the output to file.
91 saveVTK(os.path.join(save_path,"ex10a.vtu"),\
92 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
93
94 ##################################################REGRIDDING & PLOTTING
95
96
97 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
98 pl.matplotlib.pyplot.autumn()
99 pl.contourf(xi,yi,zi,10)
100 pl.xlabel("Horizontal Displacement (m)")
101 pl.ylabel("Depth (m)")
102 pl.savefig(os.path.join(save_path,"Ucontour.png"))
103 print "Solution has been plotted ..."
104
105 cut=int(len(xi)/2)
106
107 pl.clf()
108
109 r=np.linspace(0,mx/2,100)
110 m=2*pl.pi*10*10*200*-G/(r*r)
111
112 pl.plot(xi,zi[:,cut])
113 #pl.plot(r+2500,m)
114 pl.title("Potential Profile")
115 pl.xlabel("Horizontal Displacement (m)")
116 pl.ylabel("Potential")
117 pl.savefig(os.path.join(save_path,"Upot00.png"))
118
119 out=np.array([xi,zi[:,cut]])
120 pl.savetxt('profile1.asc',out.transpose())
121 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26