/[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 5448 - (show annotations)
Fri Feb 6 05:31:37 2015 UTC (4 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 4722 byte(s)
Updating all the dates
1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2009-2015 by University of Queensland
6 # http://www.uq.edu.au
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 __copyright__="""Copyright (c) 2009-2015 by University of Queensland
19 http://www.uq.edu.au
20 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 # Model of gravitational Potential for a gravity POLE.
32
33 #######################################################EXTERNAL MODULES
34 # To solve the problem it is necessary to import the modules we require.
35 import matplotlib
36 matplotlib.use('agg') #It's just here for automated testing
37
38 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 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
42 import os, sys #This package is necessary to handle saving our data.
43 from math import pi, sqrt, sin, cos
44
45 from esys.escript.pdetools import Projector
46
47 from cblib import toRegGrid
48 import pylab as pl #Plotting package
49 import numpy as np
50
51 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 ########################################################MPI WORLD CHECK
58 if getMPISizeWorld() > 1:
59 print("This example will not run in an MPI world.")
60 sys.exit(0)
61
62 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 #################################################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
81 ################################################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
87 ####################################################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
94 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
95 ###############################################ESCRIPT PDE CONSTRUCTION
96
97 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
103 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
110 ##################################################REGRIDDING & PLOTTING
111
112
113 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
121 cut=int(len(xi)//2)
122
123 pl.clf()
124
125 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
128 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
135 out=np.array([xi,zi[:,cut]])
136 pl.savetxt('profile1.asc',out.transpose())
137 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26