/[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 3409 - (show annotations)
Thu Dec 9 07:01:23 2010 UTC (8 years, 2 months ago) by ahallam
File MIME type: text/x-python
File size: 4014 byte(s)
minor changes to examples before holiday.
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)
84 mypde.setValue(q=q,r=0)
85 mypde.setSymmetryOn()
86 sol=mypde.getSolution()
87
88 g_field=grad(sol) #The gravitational acceleration g.
89 g_fieldz=g_field*[0,1] #The vertical component of the g field.
90 gz=length(g_fieldz) #The magnitude of the vertical component.
91 # Save the output to file.
92 saveVTK(os.path.join(save_path,"ex10a.vtu"),\
93 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
94
95 ##################################################REGRIDDING & PLOTTING
96
97
98 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
99 pl.matplotlib.pyplot.autumn()
100 pl.contourf(xi,yi,zi,10)
101 pl.xlabel("Horizontal Displacement (m)")
102 pl.ylabel("Depth (m)")
103 pl.savefig(os.path.join(save_path,"Ucontour.png"))
104 print "Solution has been plotted ..."
105
106 cut=int(len(xi)/2)
107
108 pl.clf()
109
110 r=np.linspace(0,mx/2,100)
111 m=2*pl.pi*10*10*200*-G/(r*r)
112
113 pl.plot(xi,zi[:,cut])
114 #pl.plot(r+2500,m)
115 pl.title("Potential Profile")
116 pl.xlabel("Horizontal Displacement (m)")
117 pl.ylabel("Potential")
118 pl.savefig(os.path.join(save_path,"Upot00.png"))
119
120 out=np.array([xi,zi[:,cut]])
121 pl.savetxt('profile1.asc',out.transpose())
122 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26