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

  ViewVC Help
Powered by ViewVC 1.1.26