/[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 5707 - (hide annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 7 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 jfenwick 3981 ##############################################################################
2 jfenwick 2697 #
3 jfenwick 5593 # Copyright (c) 2009-2015 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 jfenwick 2697 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 sshaw 5707 from __future__ import division, print_function
16 jfenwick 2697
17 jfenwick 5593 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 jfenwick 2697 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 caltinay 3346 from esys.weipa import saveVTK
40 jfenwick 2697
41 gross 2904
42 jfenwick 2697 # routine to find consecutive coordinates of a loop in pycad
43     def getLoopCoords(loop):
44 ahallam 2829 # 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 jfenwick 2697
55 jfenwick 3981 ##############################################################################
56 jfenwick 2697 # subroutine: cbphones
57     # Allows us to record the values of a PDE at various
58     # specified locations in the model.
59     # Arguments:
60 ahallam 2829 # 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 jfenwick 3981 ##############################################################################
66 jfenwick 2697 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 jfenwick 3981 ##############################################################################
83 jfenwick 2697 # subroutine: wavesolver2d
84     # Can solve a generic 2D wave propagation problem with a
85     # point source in a homogeneous medium.
86     # Arguments:
87 ahallam 2829 # 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 jfenwick 3981 ##############################################################################
96 ahallam 2829 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath,output="vtk"):
97 jfenwick 2697 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 jfenwick 3892 print("src_radius = ",src_radius)
110 jfenwick 2697
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 ahallam 2829 # while t<tend:
135     while t<1.:
136    
137 jfenwick 2697 # ... get current stress ....
138 ahallam 2829 t=1.
139 jfenwick 2697 ##OLD WAY
140 ahallam 2829 break
141     g=grad(u)
142     stress=lam*trace(g)*kmat+mu*(g+transpose(g))
143 jfenwick 2697 ### ... 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 ahallam 2829 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
150     u_p1 = mypde.getSolution()
151 jfenwick 2697 # ... shift displacements ....
152 ahallam 2829 u_m1=u
153     u=u_p1
154 jfenwick 2697 #stress =
155 ahallam 2829 t+=h
156     n+=1
157 jfenwick 3892 print(n,"-th time step t ",t)
158 ahallam 2829 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
159 jfenwick 2697
160     # print "u at point charge=",u_pc
161 ahallam 2829 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 jfenwick 2697 # save displacements at point source to file for t > 0
169 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))
170 jfenwick 2697
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 ahallam 2829 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 jfenwick 2697
182 ahallam 2829 #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 jfenwick 2697 u_pc_data.close()
201    
202    
203 jfenwick 3981 ##############################################################################
204 jfenwick 2697 # 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 ahallam 2829 # 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 jfenwick 3981 ##############################################################################
217 jfenwick 2697 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 jfenwick 3892 print("src_radius = ",src_radius)
231 jfenwick 2697
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 jfenwick 3892 print(n,"-th time step t ",t)
277 jfenwick 2697 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