/[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 4859 - (show annotations)
Thu Apr 10 00:23:29 2014 UTC (4 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4306 byte(s)
Avoid division by zero

1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2009-2014 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-2014 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.finley import Rectangle # This imports the rectangle domain function from finley
42 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
43 import os, sys #This package is necessary to handle saving our data.
44 from math import pi, sqrt, sin, cos
45
46 from esys.escript.pdetools import Projector
47
48 from cblib import toRegGrid
49 import pylab as pl #Plotting package
50 import numpy as np
51
52 ########################################################MPI WORLD CHECK
53 if getMPISizeWorld() > 1:
54 import sys
55 print("This example will not run in an MPI world.")
56 sys.exit(0)
57
58 #################################################ESTABLISHING VARIABLES
59 #Domain related.
60 mx = 5000*m #meters - model length
61 my = -5000*m #meters - model width
62 ndx = 100 # mesh steps in x direction
63 ndy = 100 # mesh steps in y direction - one dimension means one element
64 #PDE related
65 rho=200.0
66 rholoc=[2500,-2500]
67 G=6.67300*10E-11
68
69 ################################################ESTABLISHING PARAMETERS
70 #the folder to put our outputs in, leave blank "" for script path
71 save_path= os.path.join("data","example10")
72 #ensure the dir exists
73 mkDir(save_path)
74
75 ####################################################DOMAIN CONSTRUCTION
76 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
77 x=Solution(domain).getX()
78 mask=wherePositive(10-length(x-rholoc))
79 rho=rho*mask
80 kro=kronecker(domain)
81
82 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
83 ###############################################ESCRIPT PDE CONSTRUCTION
84
85 mypde=LinearPDE(domain)
86 mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
87 mypde.setValue(q=q,r=0)
88 mypde.setSymmetryOn()
89 sol=mypde.getSolution()
90
91 g_field=grad(sol) #The gravitational acceleration g.
92 g_fieldz=g_field*[0,1] #The vertical component of the g field.
93 gz=length(g_fieldz) #The magnitude of the vertical component.
94 # Save the output to file.
95 saveVTK(os.path.join(save_path,"ex10a.vtu"),\
96 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
97
98 ##################################################REGRIDDING & PLOTTING
99
100
101 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
102 pl.matplotlib.pyplot.autumn()
103 pl.contourf(xi,yi,zi,10)
104 pl.xlabel("Horizontal Displacement (m)")
105 pl.ylabel("Depth (m)")
106 pl.savefig(os.path.join(save_path,"Ucontour.png"))
107 print("Solution has been plotted ...")
108
109 cut=int(len(xi)//2)
110
111 pl.clf()
112
113 r=np.linspace(0.0000001,mx/2,100) # starting point would be 0 but that would cause division by zero later
114 m=2*pl.pi*10*10*200*-G/(r*r)
115
116 pl.plot(xi,zi[:,cut])
117 #pl.plot(r+2500,m)
118 pl.title("Potential Profile")
119 pl.xlabel("Horizontal Displacement (m)")
120 pl.ylabel("Potential")
121 pl.savefig(os.path.join(save_path,"Upot00.png"))
122
123 out=np.array([xi,zi[:,cut]])
124 pl.savetxt('profile1.asc',out.transpose())
125 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26