# Diff of /trunk/doc/examples/cookbook/cblib.py

revision 2666 by ahallam, Thu Sep 10 02:58:44 2009 UTC revision 2667 by jfenwick, Thu Sep 17 01:49:11 2009 UTC
# Line 1  Line 1
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2009 by University of Queensland
4  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
5  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
6  #  #
# Line 10  Line 10
10  #  #
11  ########################################################  ########################################################
12
13  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2009 by University of Queensland
14  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
15  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
16  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 23  Author: Antony Hallam antony.hallam@uqco Line 23  Author: Antony Hallam antony.hallam@uqco
23  """  """
24
25  from esys.pycad import CurveLoop  from esys.pycad import CurveLoop
26  import numpy as np
# To run the function it is necessary to import the modules we
# require.
# This imports everything from the escript library
#from esys.escript import Locator
27  # numpy for array handling  # numpy for array handling
28  import numpy as np  import numpy as np
# pylab for matplotlib and plotting
import pylab as pl
29  # tools for dealing with PDEs - contains locator  # tools for dealing with PDEs - contains locator
30  from esys.escript.pdetools import Locator  from esys.escript.pdetools import Locator
31
32    from esys.escript import *
33    from esys.escript.linearPDEs import LinearPDE
34
35  # routine to find consecutive coordinates of a loop in pycad  # routine to find consecutive coordinates of a loop in pycad
36  def getLoopCoords(loop):  def getLoopCoords(loop):
37      # return all construction points of input      # return all construction points of input
# Line 76  def toXYTuple(coords): Line 73  def toXYTuple(coords):
73      coordX = coords[:,0]; coordY = coords[:,1] #X and Y components.      coordX = coords[:,0]; coordY = coords[:,1] #X and Y components.
74      return coordX,coordY      return coordX,coordY
75
76
77    ########################################################
78    # subroutine: cbphones
79    # Allows us to record the values of a PDE at various
80    # specified locations in the model.
81    # Arguments:
82    #   domain  : domain of model
83    #   U       : Current time state displacement solution.
84    #   phones  : Geophone Locations
85    #   dim     : model dimesions
86    #   savepath: where to output the data files local is default
87    ########################################################
88    def cbphones(domain,U,phones,dim,savepath=""):
89       #find the number of geophones
90       nphones = len(phones)
91       u_pot = np.zeros([nphones,dim],float)
92
93       for i in range(0,nphones):
94         # define the location of the phone source
95         L=Locator(domain,numpy.array(phones[i]))
96         # find potential at point source.
97         temp = L.getValue(U)
98         for j in range(0,dim):
99           u_pot[i,j]=temp[j]
100
101       # open file to save displacement at point source
102       return u_pot
103
104    ########################################################
105    # subroutine: wavesolver2d
106    # Can solve a generic 2D wave propagation problem with a
107    # point source in a homogeneous medium.
108    # Arguments:
109    #   domain  : domain to solve over
110    #   h       : time step
111    #   tend    : end time
112    #   lam, mu : lames constants for domain
113    #   rho : density of domain
114    #   U0  : magnitude of source
115    #   xc  : source location in domain (Vector)
116    #   savepath: where to output the data files
117    ########################################################
118    def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
119       from esys.escript.linearPDEs import LinearPDE
120       x=domain.getX()
121       # ... open new PDE ...
122       mypde=LinearPDE(domain)
123       #mypde.setSolverMethod(LinearPDE.LUMPING)
124       mypde.setSymmetryOn()
125       kmat = kronecker(domain)
126       mypde.setValue(D=kmat*rho)
127
128       # define small radius around point xc
129       # Lsup(x) returns the maximum value of the argument x
130       src_radius = 50#2*Lsup(domain.getSize())
131       print "src_radius = ",src_radius
132
133       dunit=numpy.array([0.,1.]) # defines direction of point source
134
135
136       # ... set initial values ....
137       n=0
138       # initial value of displacement at point source is constant (U0=0.01)
139       # for first two time steps
140       u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
141       u_m1=u
142       t=0
143
144       u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
145       u_pc_x1 = u_pot[0,0]
146       u_pc_y1 = u_pot[0,1]
147       u_pc_x2 = u_pot[1,0]
148       u_pc_y2 = u_pot[1,1]
149       u_pc_x3 = u_pot[2,0]
150       u_pc_y3 = u_pot[2,1]
151
152       # open file to save displacement at point source
153       u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
154       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))
155
156       while t<tend:
157         # ... get current stress ....
158
159    ##OLD WAY
160         g=grad(u)
161         stress=lam*trace(g)*kmat+mu*(g+transpose(g))
162         ### ... get new acceleration ....
163         #mypde.setValue(X=-stress)
164         #a=mypde.getSolution()
165         ### ... get new displacement ...
166         #u_p1=2*u-u_m1+h*h*a
167    ###NEW WAY
168         mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
169         u_p1 = mypde.getSolution()
170         # ... shift displacements ....
171         u_m1=u
172         u=u_p1
173         #stress =
174         t+=h
175         n+=1
176         print n,"-th time step t ",t
177         u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
178
179    #     print "u at point charge=",u_pc
180         u_pc_x1 = u_pot[0,0]
181         u_pc_y1 = u_pot[0,1]
182         u_pc_x2 = u_pot[1,0]
183         u_pc_y2 = u_pot[1,1]
184         u_pc_x3 = u_pot[2,0]
185         u_pc_y3 = u_pot[2,1]
186
187         # save displacements at point source to file for t > 0
188         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))
189
190         # ... save current acceleration in units of gravity and displacements
191         #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
192         #displacement = length(u), tensor = stress, Ux = u[0] )
193         saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
194
195       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    #   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    ########################################################
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  # joel wrote this to create directories for you
293  import os  import os
294  def needdirs(dirlist):  def needdirs(dirlist):

Legend:
 Removed from v.2666 changed lines Added in v.2667

 ViewVC Help Powered by ViewVC 1.1.26