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

Contents of /trunk/doc/examples/cookbook/example10a.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 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 ##############################################################################
2 #
3 # Copyright (c) 2009-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 http://www.uq.edu.au
19 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 # Model of gravitational Potential for a gravity POLE.
31
32 #######################################################EXTERNAL MODULES
33 # To solve the problem it is necessary to import the modules we require.
34 import matplotlib
35 matplotlib.use('agg') #It's just here for automated testing
36
37 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.weipa import saveVTK # This imports the VTK file saver from weipa
41 import os, sys #This package is necessary to handle saving our data.
42 from math import pi, sqrt, sin, cos
43
44 from esys.escript.pdetools import Projector
45
46 from cblib import toRegGrid
47 import pylab as pl #Plotting package
48 import numpy as np
49
50 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 ########################################################MPI WORLD CHECK
57 if getMPISizeWorld() > 1:
58 print("This example will not run in an MPI world.")
59 sys.exit(0)
60
61 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 #################################################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
80 ################################################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
86 ####################################################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
93 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
94 ###############################################ESCRIPT PDE CONSTRUCTION
95
96 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
102 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
109 ##################################################REGRIDDING & PLOTTING
110
111
112 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
120 cut=int(len(xi)//2)
121
122 pl.clf()
123
124 r=np.linspace(0.0000001,mx/2,100) # starting point would be 0 but that would cause division by zero later
125 m=2*pl.pi*10*10*200*-G/(r*r)
126
127 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
134 out=np.array([xi,zi[:,cut]])
135 pl.savetxt('profile1.asc',out.transpose())
136 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26