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

  ViewVC Help
Powered by ViewVC 1.1.26