/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years ago) by jfenwick
File MIME type: text/x-python
File size: 4156 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



1
2 ##############################################################################
3 #
4 # Copyright (c) 2009-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2009-2014 by 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.finley import Rectangle # This imports the rectangle domain function from finley
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 ########################################################MPI WORLD CHECK
52 if getMPISizeWorld() > 1:
53 import sys
54 print("This example will not run in an MPI world.")
55 sys.exit(0)
56
57 #################################################ESTABLISHING VARIABLES
58 #Domain related.
59 mx = 5000*m #meters - model length
60 my = -5000*m #meters - model width
61 ndx = 100 # mesh steps in x direction
62 ndy = 100 # mesh steps in y direction - one dimension means one element
63 #PDE related
64 rho=200.0
65 rholoc=[2500,-2500]
66 G=6.67300*10E-11
67
68 ################################################ESTABLISHING PARAMETERS
69 #the folder to put our outputs in, leave blank "" for script path
70 save_path= os.path.join("data","example10")
71 #ensure the dir exists
72 mkDir(save_path)
73
74 ####################################################DOMAIN CONSTRUCTION
75 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
76 x=Solution(domain).getX()
77 mask=wherePositive(10-length(x-rholoc))
78 rho=rho*mask
79 kro=kronecker(domain)
80
81 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
82 ###############################################ESCRIPT PDE CONSTRUCTION
83
84 mypde=LinearPDE(domain)
85 mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
86 mypde.setValue(q=q,r=0)
87 mypde.setSymmetryOn()
88 sol=mypde.getSolution()
89
90 g_field=grad(sol) #The gravitational acceleration g.
91 g_fieldz=g_field*[0,1] #The vertical component of the g field.
92 gz=length(g_fieldz) #The magnitude of the vertical component.
93 # Save the output to file.
94 saveVTK(os.path.join(save_path,"ex10a.vtu"),\
95 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
96
97 ##################################################REGRIDDING & PLOTTING
98
99
100 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
101 pl.matplotlib.pyplot.autumn()
102 pl.contourf(xi,yi,zi,10)
103 pl.xlabel("Horizontal Displacement (m)")
104 pl.ylabel("Depth (m)")
105 pl.savefig(os.path.join(save_path,"Ucontour.png"))
106 print("Solution has been plotted ...")
107
108 cut=int(len(xi)/2)
109
110 pl.clf()
111
112 r=np.linspace(0,mx/2,100)
113 m=2*pl.pi*10*10*200*-G/(r*r)
114
115 pl.plot(xi,zi[:,cut])
116 #pl.plot(r+2500,m)
117 pl.title("Potential Profile")
118 pl.xlabel("Horizontal Displacement (m)")
119 pl.ylabel("Potential")
120 pl.savefig(os.path.join(save_path,"Upot00.png"))
121
122 out=np.array([xi,zi[:,cut]])
123 pl.savetxt('profile1.asc',out.transpose())
124 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26