/[escript]/release/4.1/doc/examples/cookbook/example10a.py
ViewVC logotype

Annotation of /release/4.1/doc/examples/cookbook/example10a.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5767 - (hide annotations)
Wed Jul 29 01:33:52 2015 UTC (3 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 4703 byte(s)
more examples using lib check for natgrid
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 sshaw 5767 from cblib import toRegGrid, HAVE_NATGRID
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 sshaw 5767 if not HAVE_NATGRID:
62     print("This example requires that natgrid is available to matplotlib")
63 jfenwick 5361
64     if HAVE_FINLEY and HAVE_NATGRID:
65 sshaw 5288 #################################################ESTABLISHING VARIABLES
66     #Domain related.
67     mx = 5000*m #meters - model length
68     my = -5000*m #meters - model width
69     ndx = 100 # mesh steps in x direction
70     ndy = 100 # mesh steps in y direction - one dimension means one element
71     #PDE related
72     rho=200.0
73     rholoc=[2500,-2500]
74     G=6.67300*10E-11
75 ahallam 3099
76 sshaw 5288 ################################################ESTABLISHING PARAMETERS
77     #the folder to put our outputs in, leave blank "" for script path
78     save_path= os.path.join("data","example10")
79     #ensure the dir exists
80     mkDir(save_path)
81 ahallam 3099
82 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
83     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
84     x=Solution(domain).getX()
85     mask=wherePositive(10-length(x-rholoc))
86     rho=rho*mask
87     kro=kronecker(domain)
88 ahallam 3099
89 sshaw 5288 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
90     ###############################################ESCRIPT PDE CONSTRUCTION
91 ahallam 3099
92 sshaw 5288 mypde=LinearPDE(domain)
93     mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
94     mypde.setValue(q=q,r=0)
95     mypde.setSymmetryOn()
96     sol=mypde.getSolution()
97 ahallam 3135
98 sshaw 5288 g_field=grad(sol) #The gravitational acceleration g.
99     g_fieldz=g_field*[0,1] #The vertical component of the g field.
100     gz=length(g_fieldz) #The magnitude of the vertical component.
101     # Save the output to file.
102     saveVTK(os.path.join(save_path,"ex10a.vtu"),\
103     grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
104 ahallam 3135
105 sshaw 5288 ##################################################REGRIDDING & PLOTTING
106 ahallam 3135
107    
108 sshaw 5288 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
109     pl.matplotlib.pyplot.autumn()
110     pl.contourf(xi,yi,zi,10)
111     pl.xlabel("Horizontal Displacement (m)")
112     pl.ylabel("Depth (m)")
113     pl.savefig(os.path.join(save_path,"Ucontour.png"))
114     print("Solution has been plotted ...")
115 ahallam 3135
116 sshaw 5288 cut=int(len(xi)//2)
117 ahallam 3135
118 sshaw 5288 pl.clf()
119 ahallam 3135
120 sshaw 5705 r=np.linspace(0.0000001,mx/2,100) # starting point would be 0 but that would cause division by zero later
121 sshaw 5288 m=2*pl.pi*10*10*200*-G/(r*r)
122 ahallam 3135
123 sshaw 5288 pl.plot(xi,zi[:,cut])
124     #pl.plot(r+2500,m)
125     pl.title("Potential Profile")
126     pl.xlabel("Horizontal Displacement (m)")
127     pl.ylabel("Potential")
128     pl.savefig(os.path.join(save_path,"Upot00.png"))
129 ahallam 3135
130 sshaw 5288 out=np.array([xi,zi[:,cut]])
131     pl.savetxt('profile1.asc',out.transpose())
132     pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26