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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 10001 byte(s)
Don't panic.
Updating copyright stamps

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

  ViewVC Help
Powered by ViewVC 1.1.26