/[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 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (6 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 4026 byte(s)
Copyright changes
1 ahallam 3099
2     ########################################################
3     #
4 jfenwick 3911 # Copyright (c) 2009-2012 by University of Queensland
5 ahallam 3099 # 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 jfenwick 3911 __copyright__="""Copyright (c) 2009-2012 by University of Queensland
15 ahallam 3099 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 ahallam 3195 # Model of gravitational Potential for a gravity POLE.
29 ahallam 3099
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 caltinay 3346 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
37 ahallam 3099 import os, sys #This package is necessary to handle saving our data.
38 ahallam 3135 from math import pi, sqrt, sin, cos
39 ahallam 3099
40 ahallam 3135 from esys.escript.pdetools import Projector
41 jfenwick 3148
42     import matplotlib
43     matplotlib.use('agg') #It's just here for automated testing
44    
45 gross 3430 from cblib import toRegGrid
46 ahallam 3135 import pylab as pl #Plotting package
47     import numpy as np
48    
49 ahallam 3099 ########################################################MPI WORLD CHECK
50     if getMPISizeWorld() > 1:
51     import sys
52 jfenwick 3892 print("This example will not run in an MPI world.")
53 ahallam 3099 sys.exit(0)
54    
55     #################################################ESTABLISHING VARIABLES
56     #Domain related.
57 ahallam 3135 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 ahallam 3099 #PDE related
62 ahallam 3135 rho=200.0
63     rholoc=[2500,-2500]
64 ahallam 3099 G=6.67300*10E-11
65    
66     ################################################ESTABLISHING PARAMETERS
67     #the folder to put our outputs in, leave blank "" for script path
68 ahallam 3135 save_path= os.path.join("data","example10")
69 ahallam 3099 #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 ahallam 3135 mask=wherePositive(10-length(x-rholoc))
76 ahallam 3099 rho=rho*mask
77     kro=kronecker(domain)
78    
79 ahallam 3135 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
80 ahallam 3099 ###############################################ESCRIPT PDE CONSTRUCTION
81    
82     mypde=LinearPDE(domain)
83 ahallam 3409 mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
84     mypde.setValue(q=q,r=0)
85 ahallam 3135 mypde.setSymmetryOn()
86 ahallam 3099 sol=mypde.getSolution()
87 ahallam 3135
88 ahallam 3409 g_field=grad(sol) #The gravitational acceleration g.
89 ahallam 3195 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 ahallam 3135 # 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 jfenwick 3892 print("Solution has been plotted ...")
105 ahallam 3135
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