/[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 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 9620 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



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

  ViewVC Help
Powered by ViewVC 1.1.26