/[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 4859 - (hide 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 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     from esys.finley import Rectangle # This imports the rectangle domain function from finley
42 caltinay 3346 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
43 ahallam 3099 import os, sys #This package is necessary to handle saving our data.
44 ahallam 3135 from math import pi, sqrt, sin, cos
45 ahallam 3099
46 ahallam 3135 from esys.escript.pdetools import Projector
47 jfenwick 3148
48 gross 3430 from cblib import toRegGrid
49 ahallam 3135 import pylab as pl #Plotting package
50     import numpy as np
51    
52 ahallam 3099 ########################################################MPI WORLD CHECK
53     if getMPISizeWorld() > 1:
54     import sys
55 jfenwick 3892 print("This example will not run in an MPI world.")
56 ahallam 3099 sys.exit(0)
57    
58     #################################################ESTABLISHING VARIABLES
59     #Domain related.
60 ahallam 3135 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 ahallam 3099 #PDE related
65 ahallam 3135 rho=200.0
66     rholoc=[2500,-2500]
67 ahallam 3099 G=6.67300*10E-11
68    
69     ################################################ESTABLISHING PARAMETERS
70     #the folder to put our outputs in, leave blank "" for script path
71 ahallam 3135 save_path= os.path.join("data","example10")
72 ahallam 3099 #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 ahallam 3135 mask=wherePositive(10-length(x-rholoc))
79 ahallam 3099 rho=rho*mask
80     kro=kronecker(domain)
81    
82 ahallam 3135 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
83 ahallam 3099 ###############################################ESCRIPT PDE CONSTRUCTION
84    
85     mypde=LinearPDE(domain)
86 ahallam 3409 mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
87     mypde.setValue(q=q,r=0)
88 ahallam 3135 mypde.setSymmetryOn()
89 ahallam 3099 sol=mypde.getSolution()
90 ahallam 3135
91 ahallam 3409 g_field=grad(sol) #The gravitational acceleration g.
92 ahallam 3195 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 ahallam 3135 # 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 jfenwick 3892 print("Solution has been plotted ...")
108 ahallam 3135
109 jfenwick 4853 cut=int(len(xi)//2)
110 ahallam 3135
111     pl.clf()
112    
113 jfenwick 4859 r=np.linspace(0.0000001,mx/2,100) # starting point would be 0 but that would cause division by zero later
114 ahallam 3135 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