/[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 4657 - (hide annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years ago) by jfenwick
File MIME type: text/x-python
File size: 4156 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26