/[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 5288 - (show 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 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.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 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
74 ################################################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
80 ####################################################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
87 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
88 ###############################################ESCRIPT PDE CONSTRUCTION
89
90 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
96 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
103 ##################################################REGRIDDING & PLOTTING
104
105
106 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
114 cut=int(len(xi)//2)
115
116 pl.clf()
117
118 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
121 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
128 out=np.array([xi,zi[:,cut]])
129 pl.savetxt('profile1.asc',out.transpose())
130 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26