/[escript]/trunk/doc/cookbook/cblib/wavesolver2d.py
ViewVC logotype

Contents of /trunk/doc/cookbook/cblib/wavesolver2d.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2434 - (show annotations)
Thu May 21 07:59:53 2009 UTC (11 years, 4 months ago) by ahallam
File MIME type: text/x-python
File size: 3967 byte(s)
Week 8: Still working on seismic wave propagation. Need to come up with a smooth source function next week. All TEXT files up to date.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 # You can shorten the execution time by reducing variable tend from 60 to 0.5
23
24 # Importing all the necessary modules required.
25 from esys.escript import *
26 from esys.escript.pdetools import Locator
27 from esys.escript.linearPDEs import LinearPDE
28 from esys.finley import Rectangle
29 from numarray import identity,zeros,ones
30 import os
31
32 ########################################################
33 # subroutine: wavesolver2d
34 # Can solve a generic 2D wave propagation problem with a
35 # point source in a homogeneous medium.
36 # Arguments:
37 # domain : domain to solve over
38 # h : time step
39 # tend : end time
40 # lam, mu : lames constants for domain
41 # rho : density of domain
42 # U0 : magnitude of source
43 # xc : source location in domain (Vector)
44 # savepath: where to output the data files
45 ########################################################
46 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
47 x=domain.getX()
48 # ... open new PDE ...
49 mypde=LinearPDE(domain)
50 # mypde.setSolverMethod(LinearPDE.LUMPING)
51 mypde.setSymmetryOn()
52 kmat = kronecker(domain)
53 mypde.setValue(D=kmat*rho)
54
55 # define small radius around point xc
56 # Lsup(x) returns the maximum value of the argument x
57 src_radius = 50#2*Lsup(domain.getSize())
58 print "src_radius = ",src_radius
59
60 dunit=numarray.array([1.,0.]) # defines direction of point source
61
62
63 # ... set initial values ....
64 n=0
65 # initial value of displacement at point source is constant (U0=0.01)
66 # for first two time steps
67 u=U0*whereNegative(length(x-xc)-src_radius)*dunit
68 print u
69 u_m1=U0*whereNegative(length(x-xc)-src_radius)*dunit
70 t=0
71
72 # define the location of the point source
73 L=Locator(domain,numarray.array(xc))
74 # find potential at point source
75 u_pc=L.getValue(u)
76 print "u at point charge=",u_pc
77
78 u_pc_x = u_pc[0]
79 u_pc_y = u_pc[1]
80
81 # open file to save displacement at point source
82 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
83 u_pc_data.write("%f %f %f\n"%(t,u_pc_x,u_pc_y))
84
85 #NEW WAY
86 #g=grad(u)
87 #stress=lam*trace(g)*kmat+mu*(g+transpose(g))
88
89 while t<tend:
90 # ... get current stress ....
91
92 ##OLD WAY
93 g=grad(u)
94 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
95 ### ... get new acceleration ....
96 mypde.setValue(X=-stress)
97 a=mypde.getSolution()
98 ### ... get new displacement ...
99 u_p1=2*u-u_m1+h*h*a
100 ###NEW WAY
101 #mypde.setValue(X=-stress*h**2)
102 #mypde.setValue(Y=rho*(2*u-u_m1))
103 #u_p1 = mypde.getSolution()
104
105 # ... shift displacements ....
106 u_m1=u
107 u=u_p1
108 #stress =
109 t+=h
110 n+=1
111 print n,"-th time step t ",t
112 u_pc=L.getValue(u)
113 # print "u at point charge=",u_pc
114
115 u_pc_x=u_pc[0]
116 u_pc_y=u_pc[1]
117
118 # save displacements at point source to file for t > 0
119 u_pc_data.write("%f %f %f\n"%(t,u_pc_x,u_pc_y))
120
121 # ... save current acceleration in units of gravity and displacements
122 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
123 #displacement = length(u), tensor = stress, Ux = u[0] )
124 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
125 u_pc_data.close()

  ViewVC Help
Powered by ViewVC 1.1.26