/[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 4087 - (hide 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 ahallam 3099
2 jfenwick 3981 ##############################################################################
3 ahallam 3099 #
4 jfenwick 3911 # Copyright (c) 2009-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ahallam 3099 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 ahallam 3099
16 jfenwick 3911 __copyright__="""Copyright (c) 2009-2012 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 ahallam 3099 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 ahallam 3195 # Model of gravitational Potential for a gravity POLE.
30 ahallam 3099
31     #######################################################EXTERNAL MODULES
32     # To solve the problem it is necessary to import the modules we require.
33 caltinay 4087 import matplotlib
34     matplotlib.use('agg') #It's just here for automated testing
35    
36 ahallam 3099 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 caltinay 3346 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
41 ahallam 3099 import os, sys #This package is necessary to handle saving our data.
42 ahallam 3135 from math import pi, sqrt, sin, cos
43 ahallam 3099
44 ahallam 3135 from esys.escript.pdetools import Projector
45 jfenwick 3148
46 gross 3430 from cblib import toRegGrid
47 ahallam 3135 import pylab as pl #Plotting package
48     import numpy as np
49    
50 ahallam 3099 ########################################################MPI WORLD CHECK
51     if getMPISizeWorld() > 1:
52     import sys
53 jfenwick 3892 print("This example will not run in an MPI world.")
54 ahallam 3099 sys.exit(0)
55    
56     #################################################ESTABLISHING VARIABLES
57     #Domain related.
58 ahallam 3135 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 ahallam 3099 #PDE related
63 ahallam 3135 rho=200.0
64     rholoc=[2500,-2500]
65 ahallam 3099 G=6.67300*10E-11
66    
67     ################################################ESTABLISHING PARAMETERS
68     #the folder to put our outputs in, leave blank "" for script path
69 ahallam 3135 save_path= os.path.join("data","example10")
70 ahallam 3099 #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 ahallam 3135 mask=wherePositive(10-length(x-rholoc))
77 ahallam 3099 rho=rho*mask
78     kro=kronecker(domain)
79    
80 ahallam 3135 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
81 ahallam 3099 ###############################################ESCRIPT PDE CONSTRUCTION
82    
83     mypde=LinearPDE(domain)
84 ahallam 3409 mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
85     mypde.setValue(q=q,r=0)
86 ahallam 3135 mypde.setSymmetryOn()
87 ahallam 3099 sol=mypde.getSolution()
88 ahallam 3135
89 ahallam 3409 g_field=grad(sol) #The gravitational acceleration g.
90 ahallam 3195 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 ahallam 3135 # 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 jfenwick 3892 print("Solution has been plotted ...")
106 ahallam 3135
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