/[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 5288 - (hide annotations)
Tue Dec 2 23:18:40 2014 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 4587 byte(s)
fixing tests for cases where required domains not built
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 sshaw 5288 if HAVE_FINLEY:
63     #################################################ESTABLISHING VARIABLES
64     #Domain related.
65     mx = 5000*m #meters - model length
66     my = -5000*m #meters - model width
67     ndx = 100 # mesh steps in x direction
68     ndy = 100 # mesh steps in y direction - one dimension means one element
69     #PDE related
70     rho=200.0
71     rholoc=[2500,-2500]
72     G=6.67300*10E-11
73 ahallam 3099
74 sshaw 5288 ################################################ESTABLISHING PARAMETERS
75     #the folder to put our outputs in, leave blank "" for script path
76     save_path= os.path.join("data","example10")
77     #ensure the dir exists
78     mkDir(save_path)
79 ahallam 3099
80 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
81     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
82     x=Solution(domain).getX()
83     mask=wherePositive(10-length(x-rholoc))
84     rho=rho*mask
85     kro=kronecker(domain)
86 ahallam 3099
87 sshaw 5288 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
88     ###############################################ESCRIPT PDE CONSTRUCTION
89 ahallam 3099
90 sshaw 5288 mypde=LinearPDE(domain)
91     mypde.setValue(A=kro,Y=4.*3.1415*G*rho)
92     mypde.setValue(q=q,r=0)
93     mypde.setSymmetryOn()
94     sol=mypde.getSolution()
95 ahallam 3135
96 sshaw 5288 g_field=grad(sol) #The gravitational acceleration g.
97     g_fieldz=g_field*[0,1] #The vertical component of the g field.
98     gz=length(g_fieldz) #The magnitude of the vertical component.
99     # Save the output to file.
100     saveVTK(os.path.join(save_path,"ex10a.vtu"),\
101     grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz)
102 ahallam 3135
103 sshaw 5288 ##################################################REGRIDDING & PLOTTING
104 ahallam 3135
105    
106 sshaw 5288 xi, yi, zi = toRegGrid(sol, nx=50, ny=50)
107     pl.matplotlib.pyplot.autumn()
108     pl.contourf(xi,yi,zi,10)
109     pl.xlabel("Horizontal Displacement (m)")
110     pl.ylabel("Depth (m)")
111     pl.savefig(os.path.join(save_path,"Ucontour.png"))
112     print("Solution has been plotted ...")
113 ahallam 3135
114 sshaw 5288 cut=int(len(xi)//2)
115 ahallam 3135
116 sshaw 5288 pl.clf()
117 ahallam 3135
118 sshaw 5288 r=np.linspace(0.0000001,mx/2,100) # starting point would be 0 but that would cause division by zero later
119     m=2*pl.pi*10*10*200*-G/(r*r)
120 ahallam 3135
121 sshaw 5288 pl.plot(xi,zi[:,cut])
122     #pl.plot(r+2500,m)
123     pl.title("Potential Profile")
124     pl.xlabel("Horizontal Displacement (m)")
125     pl.ylabel("Potential")
126     pl.savefig(os.path.join(save_path,"Upot00.png"))
127 ahallam 3135
128 sshaw 5288 out=np.array([xi,zi[:,cut]])
129     pl.savetxt('profile1.asc',out.transpose())
130     pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26