/[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 4087 - (show annotations)
Thu Nov 22 22:28:01 2012 UTC (6 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 4088 byte(s)
Moved matplotlib imports in test scripts before escript since there is an
import chain which pulls it so the use() function has no effect.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2009-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2009-2012 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 Author: Antony Hallam antony.hallam@uqconnect.edu.au
25 """
26
27 ############################################################FILE HEADER
28 # example10a.py
29 # Model of gravitational Potential for a gravity POLE.
30
31 #######################################################EXTERNAL MODULES
32 # To solve the problem it is necessary to import the modules we require.
33 import matplotlib
34 matplotlib.use('agg') #It's just here for automated testing
35
36 from esys.escript import * # This imports everything from the escript library
37 from esys.escript.unitsSI import *
38 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
39 from esys.finley import Rectangle # This imports the rectangle domain function from finley
40 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
41 import os, sys #This package is necessary to handle saving our data.
42 from math import pi, sqrt, sin, cos
43
44 from esys.escript.pdetools import Projector
45
46 from cblib import toRegGrid
47 import pylab as pl #Plotting package
48 import numpy as np
49
50 ########################################################MPI WORLD CHECK
51 if getMPISizeWorld() > 1:
52 import sys
53 print("This example will not run in an MPI world.")
54 sys.exit(0)
55
56 #################################################ESTABLISHING VARIABLES
57 #Domain related.
58 mx = 5000*m #meters - model length
59 my = -5000*m #meters - model width
60 ndx = 100 # mesh steps in x direction
61 ndy = 100 # mesh steps in y direction - one dimension means one element
62 #PDE related
63 rho=200.0
64 rholoc=[2500,-2500]
65 G=6.67300*10E-11
66
67 ################################################ESTABLISHING PARAMETERS
68 #the folder to put our outputs in, leave blank "" for script path
69 save_path= os.path.join("data","example10")
70 #ensure the dir exists
71 mkDir(save_path)
72
73 ####################################################DOMAIN CONSTRUCTION
74 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
75 x=Solution(domain).getX()
76 mask=wherePositive(10-length(x-rholoc))
77 rho=rho*mask
78 kro=kronecker(domain)
79
80 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
81 ###############################################ESCRIPT PDE CONSTRUCTION
82
83 mypde=LinearPDE(domain)
84 mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
85 mypde.setValue(q=q,r=0)
86 mypde.setSymmetryOn()
87 sol=mypde.getSolution()
88
89 g_field=grad(sol) #The gravitational acceleration g.
90 g_fieldz=g_field*[0,1] #The vertical component of the g field.
91 gz=length(g_fieldz) #The magnitude of the vertical component.
92 # Save the output to file.
93 saveVTK(os.path.join(save_path,"ex10a.vtu"),\
94 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
95
96 ##################################################REGRIDDING & PLOTTING
97
98
99 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
100 pl.matplotlib.pyplot.autumn()
101 pl.contourf(xi,yi,zi,10)
102 pl.xlabel("Horizontal Displacement (m)")
103 pl.ylabel("Depth (m)")
104 pl.savefig(os.path.join(save_path,"Ucontour.png"))
105 print("Solution has been plotted ...")
106
107 cut=int(len(xi)/2)
108
109 pl.clf()
110
111 r=np.linspace(0,mx/2,100)
112 m=2*pl.pi*10*10*200*-G/(r*r)
113
114 pl.plot(xi,zi[:,cut])
115 #pl.plot(r+2500,m)
116 pl.title("Potential Profile")
117 pl.xlabel("Horizontal Displacement (m)")
118 pl.ylabel("Potential")
119 pl.savefig(os.path.join(save_path,"Upot00.png"))
120
121 out=np.array([xi,zi[:,cut]])
122 pl.savetxt('profile1.asc',out.transpose())
123 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26