/[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 5707 - (hide annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 7 months ago) by sshaw
File MIME type: text/x-python
File size: 4710 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 jfenwick 3981 ##############################################################################
2 ahallam 3099 #
3 jfenwick 5593 # Copyright (c) 2009-2015 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ahallam 3099 #
6     # Primary Business: Queensland, Australia
7     # Licensed under the Open Software License version 3.0
8     # http://www.opensource.org/licenses/osl-3.0.php
9     #
10 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 sshaw 5707 from __future__ import division, print_function
16 ahallam 3099
17 jfenwick 5593 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3099 Primary Business: Queensland, Australia"""
20     __license__="""Licensed under the Open Software License version 3.0
21     http://www.opensource.org/licenses/osl-3.0.php"""
22     __url__="https://launchpad.net/escript-finley"
23    
24     """
25     Author: Antony Hallam antony.hallam@uqconnect.edu.au
26     """
27    
28     ############################################################FILE HEADER
29     # example10a.py
30 ahallam 3195 # Model of gravitational Potential for a gravity POLE.
31 ahallam 3099
32     #######################################################EXTERNAL MODULES
33     # To solve the problem it is necessary to import the modules we require.
34 caltinay 4087 import matplotlib
35     matplotlib.use('agg') #It's just here for automated testing
36    
37 ahallam 3099 from esys.escript import * # This imports everything from the escript library
38     from esys.escript.unitsSI import *
39     from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
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 sshaw 5288 try:
51     from esys.finley import Rectangle
52     HAVE_FINLEY = True
53     except ImportError:
54     print("Finley module not available")
55     HAVE_FINLEY = False
56 ahallam 3099 ########################################################MPI WORLD CHECK
57     if getMPISizeWorld() > 1:
58 jfenwick 3892 print("This example will not run in an MPI world.")
59 ahallam 3099 sys.exit(0)
60    
61 jfenwick 5361 try:
62     from mpl_toolkits.natgrid import _natgrid
63     HAVE_NATGRID=True
64     except ImportError:
65     HAVE_NATGRID=False
66    
67    
68     if HAVE_FINLEY and HAVE_NATGRID:
69 sshaw 5288 #################################################ESTABLISHING VARIABLES
70     #Domain related.
71     mx = 5000*m #meters - model length
72     my = -5000*m #meters - model width
73     ndx = 100 # mesh steps in x direction
74     ndy = 100 # mesh steps in y direction - one dimension means one element
75     #PDE related
76     rho=200.0
77     rholoc=[2500,-2500]
78     G=6.67300*10E-11
79 ahallam 3099
80 sshaw 5288 ################################################ESTABLISHING PARAMETERS
81     #the folder to put our outputs in, leave blank "" for script path
82     save_path= os.path.join("data","example10")
83     #ensure the dir exists
84     mkDir(save_path)
85 ahallam 3099
86 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
87     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
88     x=Solution(domain).getX()
89     mask=wherePositive(10-length(x-rholoc))
90     rho=rho*mask
91     kro=kronecker(domain)
92 ahallam 3099
93 sshaw 5288 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
94     ###############################################ESCRIPT PDE CONSTRUCTION
95 ahallam 3099
96 sshaw 5288 mypde=LinearPDE(domain)
97     mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
98     mypde.setValue(q=q,r=0)
99     mypde.setSymmetryOn()
100     sol=mypde.getSolution()
101 ahallam 3135
102 sshaw 5288 g_field=grad(sol) #The gravitational acceleration g.
103     g_fieldz=g_field*[0,1] #The vertical component of the g field.
104     gz=length(g_fieldz) #The magnitude of the vertical component.
105     # Save the output to file.
106     saveVTK(os.path.join(save_path,"ex10a.vtu"),\
107     grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
108 ahallam 3135
109 sshaw 5288 ##################################################REGRIDDING & PLOTTING
110 ahallam 3135
111    
112 sshaw 5288 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
113     pl.matplotlib.pyplot.autumn()
114     pl.contourf(xi,yi,zi,10)
115     pl.xlabel("Horizontal Displacement (m)")
116     pl.ylabel("Depth (m)")
117     pl.savefig(os.path.join(save_path,"Ucontour.png"))
118     print("Solution has been plotted ...")
119 ahallam 3135
120 sshaw 5288 cut=int(len(xi)//2)
121 ahallam 3135
122 sshaw 5288 pl.clf()
123 ahallam 3135
124 sshaw 5705 r=np.linspace(0.0000001,mx/2,100) # starting point would be 0 but that would cause division by zero later
125 sshaw 5288 m=2*pl.pi*10*10*200*-G/(r*r)
126 ahallam 3135
127 sshaw 5288 pl.plot(xi,zi[:,cut])
128     #pl.plot(r+2500,m)
129     pl.title("Potential Profile")
130     pl.xlabel("Horizontal Displacement (m)")
131     pl.ylabel("Potential")
132     pl.savefig(os.path.join(save_path,"Upot00.png"))
133 ahallam 3135
134 sshaw 5288 out=np.array([xi,zi[:,cut]])
135     pl.savetxt('profile1.asc',out.transpose())
136     pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26