/[escript]/trunk/doc/examples/cookbook/cblib1.py
ViewVC logotype

Contents of /trunk/doc/examples/cookbook/cblib1.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 9814 byte(s)
Round 1 of copyright fixes
1 ##############################################################################
2 #
3 # Copyright (c) 2009-2013 by University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development since 2012 by School of Earth Sciences
12 #
13 ##############################################################################
14
15 __copyright__="""Copyright (c) 2009-2013 by University of Queensland
16 http://www.uq.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 A collection of routines to use in cookbook examples.
24 Author: Antony Hallam antony.hallam@uqconnect.edu.au
25 """
26
27 from esys.pycad import CurveLoop
28
29 # numpy for array handling
30 import numpy as np
31
32 # tools for dealing with PDEs - contains locator
33 from esys.escript.pdetools import Locator, Projector
34
35 from esys.escript import *
36 from esys.escript.linearPDEs import LinearPDE
37 from esys.weipa import saveVTK
38
39
40 # routine to find consecutive coordinates of a loop in pycad
41 def getLoopCoords(loop):
42 # return all construction points of input
43 temp = loop.getConstructionPoints()
44 #create a numpy array for xyz components or construction points
45 coords = np.zeros([len(temp),3],float)
46 #place construction points in array
47 for i in range(0,len(temp)):
48 coords[i,:]=temp[i].getCoordinates()
49 #return a numpy array
50 return coords
51
52
53 ##############################################################################
54 # subroutine: cbphones
55 # Allows us to record the values of a PDE at various
56 # specified locations in the model.
57 # Arguments:
58 # domain : domain of model
59 # U : Current time state displacement solution.
60 # phones : Geophone Locations
61 # dim : model dimesions
62 # savepath: where to output the data files local is default
63 ##############################################################################
64 def cbphones(domain,U,phones,dim,savepath=""):
65 #find the number of geophones
66 nphones = len(phones)
67 u_pot = np.zeros([nphones,dim],float)
68
69 for i in range(0,nphones):
70 # define the location of the phone source
71 L=Locator(domain,numpy.array(phones[i]))
72 # find potential at point source.
73 temp = L.getValue(U)
74 for j in range(0,dim):
75 u_pot[i,j]=temp[j]
76
77 # open file to save displacement at point source
78 return u_pot
79
80 ##############################################################################
81 # subroutine: wavesolver2d
82 # Can solve a generic 2D wave propagation problem with a
83 # point source in a homogeneous medium.
84 # Arguments:
85 # domain : domain to solve over
86 # h : time step
87 # tend : end time
88 # lam, mu : lames constants for domain
89 # rho : density of domain
90 # U0 : magnitude of source
91 # xc : source location in domain (Vector)
92 # savepath: where to output the data files
93 ##############################################################################
94 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath,output="vtk"):
95 from esys.escript.linearPDEs import LinearPDE
96 x=domain.getX()
97 # ... open new PDE ...
98 mypde=LinearPDE(domain)
99 #mypde.setSolverMethod(LinearPDE.LUMPING)
100 mypde.setSymmetryOn()
101 kmat = kronecker(domain)
102 mypde.setValue(D=kmat*rho)
103
104 # define small radius around point xc
105 # Lsup(x) returns the maximum value of the argument x
106 src_radius = 50#2*Lsup(domain.getSize())
107 print("src_radius = ",src_radius)
108
109 dunit=numpy.array([0.,1.]) # defines direction of point source
110
111
112 # ... set initial values ....
113 n=0
114 # initial value of displacement at point source is constant (U0=0.01)
115 # for first two time steps
116 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
117 u_m1=u
118 t=0
119
120 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
121 u_pc_x1 = u_pot[0,0]
122 u_pc_y1 = u_pot[0,1]
123 u_pc_x2 = u_pot[1,0]
124 u_pc_y2 = u_pot[1,1]
125 u_pc_x3 = u_pot[2,0]
126 u_pc_y3 = u_pot[2,1]
127
128 # open file to save displacement at point source
129 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
130 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
131
132 # while t<tend:
133 while t<1.:
134
135 # ... get current stress ....
136 t=1.
137 ##OLD WAY
138 break
139 g=grad(u)
140 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
141 ### ... get new acceleration ....
142 #mypde.setValue(X=-stress)
143 #a=mypde.getSolution()
144 ### ... get new displacement ...
145 #u_p1=2*u-u_m1+h*h*a
146 ###NEW WAY
147 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
148 u_p1 = mypde.getSolution()
149 # ... shift displacements ....
150 u_m1=u
151 u=u_p1
152 #stress =
153 t+=h
154 n+=1
155 print(n,"-th time step t ",t)
156 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
157
158 # print "u at point charge=",u_pc
159 u_pc_x1 = u_pot[0,0]
160 u_pc_y1 = u_pot[0,1]
161 u_pc_x2 = u_pot[1,0]
162 u_pc_y2 = u_pot[1,1]
163 u_pc_x3 = u_pot[2,0]
164 u_pc_y3 = u_pot[2,1]
165
166 # save displacements at point source to file for t > 0
167 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
168
169 # ... save current acceleration in units of gravity and displacements
170 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
171 #displacement = length(u), tensor = stress, Ux = u[0] )
172 if output == "vtk":
173 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
174 else:
175 quT=qu.toListOfTuples()
176 #Projector is used to smooth the data.
177 proj=Projector(mymesh)
178 smthT=proj(T)
179
180 #move data to a regular grid for plotting
181 xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
182
183 # contour the gridded data,
184 # select colour
185 pl.matplotlib.pyplot.autumn()
186 pl.clf()
187 # contour temperature
188 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
189 # labels and formatting
190 pl.clabel(CS, inline=1, fontsize=8)
191 pl.title("Heat Refraction across a clinal structure.")
192 pl.xlabel("Horizontal Displacement (m)")
193 pl.ylabel("Depth (m)")
194 pl.legend()
195 if getMPIRankWorld() == 0: #check for MPI processing
196 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
197
198 u_pc_data.close()
199
200
201 ##############################################################################
202 # subroutine: wavesolver2d
203 # Can solve a generic 2D wave propagation problem with a
204 # point source in a homogeneous medium with friction.
205 # Arguments:
206 # domain : domain to solve over
207 # h : time step
208 # tend : end time
209 # lam, mu : lames constants for domain
210 # rho : density of domain
211 # U0 : magnitude of source
212 # xc : source location in domain (Vector)
213 # savepath: where to output the data files
214 ##############################################################################
215 def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):
216 x=domain.getX()
217 # ... open new PDE ...
218 mypde=LinearPDE(domain)
219 #mypde.setSolverMethod(LinearPDE.LUMPING)
220 mypde.setSymmetryOn()
221 kmat = kronecker(domain)
222 mypde.setValue(D=kmat)
223 b=0.9
224
225 # define small radius around point xc
226 # Lsup(x) returns the maximum value of the argument x
227 src_radius = 50#2*Lsup(domain.getSize())
228 print("src_radius = ",src_radius)
229
230 dunit=numpy.array([0.,1.]) # defines direction of point source
231
232
233 # ... set initial values ....
234 n=0
235 # initial value of displacement at point source is constant (U0=0.01)
236 # for first two time steps
237 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
238 u_m1=u
239 t=0
240
241 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
242 u_pc_x1 = u_pot[0,0]
243 u_pc_y1 = u_pot[0,1]
244 u_pc_x2 = u_pot[1,0]
245 u_pc_y2 = u_pot[1,1]
246 u_pc_x3 = u_pot[2,0]
247 u_pc_y3 = u_pot[2,1]
248
249 # open file to save displacement at point source
250 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
251 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
252
253 while t<tend:
254 # ... get current stress ....
255
256 ##OLD WAY
257 g=grad(u)
258 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
259 ### ... get new acceleration ....
260 #mypde.setValue(X=-stress)
261 #a=mypde.getSolution()
262 ### ... get new displacement ...
263 #u_p1=2*u-u_m1+h*h*a
264 ###NEW WAY
265 y = ((rho/(-rho-b*h))*(u_m1-2*u))+(((b*h)/(-rho-(b*h)))*-u)
266 mypde.setValue(X=-stress*((h*h)/(-rho-h*b)),Y=y)
267 u_p1 = mypde.getSolution()
268 # ... shift displacements ....
269 u_m1=u
270 u=u_p1
271 #stress =
272 t+=h
273 n+=1
274 print(n,"-th time step t ",t)
275 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
276
277 # print "u at point charge=",u_pc
278 u_pc_x1 = u_pot[0,0]
279 u_pc_y1 = u_pot[0,1]
280 u_pc_x2 = u_pot[1,0]
281 u_pc_y2 = u_pot[1,1]
282 u_pc_x3 = u_pot[2,0]
283 u_pc_y3 = u_pot[2,1]
284
285 # save displacements at point source to file for t > 0
286 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
287
288 # ... save current acceleration in units of gravity and displacements
289 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
290 #displacement = length(u), tensor = stress, Ux = u[0] )
291 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
292
293 u_pc_data.close()
294

  ViewVC Help
Powered by ViewVC 1.1.26