test/python/seismic_wave.py

"""
seismic wave propagation
                                                                                                                                                                                                     
@var __author__: name of author
@var __licence__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""
                                                                                                                                                                                                     
__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__author__="Lutz Gross, l.gross@uq.edu.au"
__url__="http://www.iservo.edu.au/esys/escript"
__version__="$Revision$"
__date__="$Date$"

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
import time

output=True
n_end=10000

resolution=1000.  # number of elements per m in the finest region
o=1              # element order

l=100000.           # width and length m (without obsorber)
h=30000.            # height in m        (without obsorber)
d_absorber=l*0.10   # thickness of absorbing layer

l_sand=20000.          # thickness of sand region on surface
h_sand=5000.           # thickness of sand layer under the water

l_x_water=10000.       # length of water in x
l_y_water=10000.       # length of water in y
h_water=2000.          # depth of water region

x_sand=l/2-l_x_water/2-l_sand # x coordinate of location of sand region (without obsorber)
y_sand=l/2-l_y_water/2-l_sand # y coordinate of location of sand region (without obsorber)


# origin
origin={"x": -d_absorber, "y" : -d_absorber , "z" : -h-d_absorber }
# location and geometrical size of event reltive to origin:
xc=[l*0.2,l*0.3,-h*0.7]
src_radius  = 2*resolution
# direction of event:
event=numarray.array([0.,0.,1.])*1.e6
# time and length of the event
tc=2.
tc_length=0.5

# material properties:
bedrock=0
absorber=1
water=2
sand=3

rho_tab={}
rho_tab[bedrock]=8e3
rho_tab[absorber]=rho_tab[bedrock]
rho_tab[water]=1e3
rho_tab[sand]=5e3

mu_tab={}
mu_tab[bedrock]=1.7e11
mu_tab[absorber]=mu_tab[bedrock]
mu_tab[water]=0.
mu_tab[sand]=1.5e10

lmbd_tab={}
lmbd_tab[bedrock]=1.7e11
lmbd_tab[absorber]=lmbd_tab[bedrock]
lmbd_tab[water]=1.e9
lmbd_tab[sand]=1.5e10

eta_tab={}
eta_tab[absorber]=-log(0.05)*sqrt(rho_tab[absorber]*(lmbd_tab[absorber]+2*mu_tab[absorber]))/d_absorber
eta_tab[sand]=eta_tab[absorber]/40.
eta_tab[water]=eta_tab[absorber]/40.
eta_tab[bedrock]=eta_tab[absorber]/40.


# material properties:
bedrock=0
absorber=1
water=2
sand=3

rho={}
rho[bedrock]=8e3
rho[absorber]=rho[bedrock]
rho[water]=1e3
rho[sand]=5e3

mu={}
mu[bedrock]=1.7e11
mu[absorber]=mu[bedrock]
mu[water]=0.
mu[sand]=1.5e10

lmbd={}
lmbd[bedrock]=1.7e11
lmbd_absorber=lmbd[bedrock]
lmbd[water]=1.e9
lmbd[sand]=1.5e10

eta={}
eta[absorber]=-log(0.05)*sqrt(rho[absorber]*(lmbd_absorber+2*mu[absorber]))/d_absorber
eta[sand]=eta[absorber]/40.
eta[water]=eta[absorber]/40.
eta[bedrock]=eta[absorber]/40.

if output:
   print "event location = ",xc
   print "radius of event = ",src_radius
   print "time of event = ",tc
   print "length of event = ",tc_length
   print "direction = ",event

t_end=30.
dt_write=0.1


def getDomain():
    """
    this defines a dom as a brick of length and width l and hight h

      
    """
    global netotal
    
    v_p={}
    for tag in rho_tab.keys():
       v_p[tag]=sqrt((2*mu_tab[tag]+lmbd_tab[tag])/rho_tab[tag])
    v_p_ref=min(v_p.values())
    print "velocities: bedrock = %s, sand = %s, water =%s, absorber =%s, reference =%s"%(v_p[bedrock],v_p[sand],v_p[water],v_p[absorber],v_p_ref)

    sections={}
    sections["x"]=[d_absorber, x_sand, l_sand, l_x_water, l_sand, l-x_sand-2*l_sand-l_x_water, d_absorber]
    sections["y"]=[d_absorber, y_sand, l_sand, l_y_water, l_sand, l-y_sand-2*l_sand-l_y_water, d_absorber]
    sections["z"]=[d_absorber,h-h_water-h_sand,h_sand,h_water]
    if output:
      print "sections x = ",sections["x"]
      print "sections y = ",sections["y"]
      print "sections z = ",sections["z"]

    mats= [ 
            [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],

            [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],

            [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
              [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
              [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],

            [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
              [absorber, bedrock , sand    , water   , sand    , bedrock , absorber],
              [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
              [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
              [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ] ]
    
    num_elem={}
    for d in sections:
       num_elem[d]=[]
       for i in range(len(sections[d])):
           if d=="x":
              v_p_min=v_p[mats[0][0][i]]
              for q in range(len(sections["y"])):
                 for r in range(len(sections["z"])):
                    v_p_min=min(v_p[mats[r][q][i]],v_p_min)
           elif d=="y":
              v_p_min=v_p[mats[0][i][0]]
              for q in range(len(sections["x"])):
                 for r in range(len(sections["z"])):
                    v_p_min=min(v_p[mats[r][i][q]],v_p_min)
           elif d=="z":
              v_p_min=v_p[mats[i][0][0]]
              for q in range(len(sections["x"])):
                 for r in range(len(sections["y"])):
                    v_p_min=min(v_p[mats[i][r][q]],v_p_min)
           num_elem[d].append(max(1,int(sections[d][i] * v_p_ref/v_p_min /resolution+0.5)))
       
    ne_x=sum(num_elem["x"])
    ne_y=sum(num_elem["y"])
    ne_z=sum(num_elem["z"])
    netotal=ne_x*ne_y*ne_z
    if output: print "grid : %s x %s x %s (%s elements)"%(ne_x,ne_y,ne_z,netotal)
    dom=Brick(ne_x,ne_y,ne_z,l0=o*ne_x,l1=o*ne_y,l2=o*ne_z,order=o)
    x_old=dom.getX()
    x_new=0

    for d in sections:
        if d=="x": 
            i=0
            f=[1,0,0]
        if d=="y": 
            i=1
            f=[0,1,0]
        if d=="z": 
            i=2
            f=[0,0,1]
        x=x_old[i]

        p=origin[d]
        ne=0
        s=0.
 
        for i in range(len(sections[d])-1):
            msk=whereNonPositive(x-o*ne+0.5)
            s=s*msk + (sections[d][i]/(o*num_elem[d][i])*(x-o*ne)+p)*(1.-msk)
            ne+=num_elem[d][i]
            p+=sections[d][i]
        x_new=x_new + s * f
    dom.setX(x_new)

    fs=Function(dom)
    x=Function(dom).getX()
    x0=x[0]
    x1=x[1]
    x2=x[2]
    p_z=origin["z"]
    for i in range(len(mats)):
       f_z=wherePositive(x2-p_z)*wherePositive(x2-p_z+sections["z"][i])
       p_y=origin["y"]
       for j in range(len(mats[i])):
         f_y=wherePositive(x1-p_y)*wherePositive(x1-p_z+sections["y"][j])
         p_x=origin["x"]
         for k in range(len(mats[i][j])):
             f_x=wherePositive(x0-p_x)*wherePositive(x0-p_x+sections["x"][k]) 
             fs.setTags(mats[i][j][k],f_x*f_y*f_z)
             p_x+=sections["x"][k]
         p_y+=sections["y"][j]
       p_z+=sections["z"][i]
    return dom

def getMaterialProperties(dom):
   rho =Scalar(rho_tab[bedrock],Function(dom))
   eta =Scalar(eta_tab[bedrock],Function(dom))
   mu  =Scalar(mu_tab[bedrock],Function(dom))
   lmbd=Scalar(lmbd_tab[bedrock],Function(dom))
   tags=Scalar(bedrock,Function(dom))
   
   for tag in rho_tab.keys():
      rho.setTaggedValue(tag,rho_tab[tag])
      eta.setTaggedValue(tag,eta_tab[tag])
      mu.setTaggedValue(tag,mu_tab[tag])
      lmbd.setTaggedValue(tag,lmbd_tab[tag])
      tags.setTaggedValue(tag,tag)
   return rho,mu,lmbd,eta

def wavePropagation(dom,rho,mu,lmbd,eta):
   x=Function(dom).getX()
   # ... open new PDE ...
   mypde=LinearPDE(dom)
   mypde.setSolverMethod(LinearPDE.LUMPING)
   k=kronecker(Function(dom))
   mypde.setValue(D=k*rho)

   dt=(1./5.)*inf(dom.getSize()/sqrt((2*mu+lmbd)/rho))
   if output: print "time step size = ",dt
   # ... set initial values ....
   n=0
   t=0
   t_write=0.
   n_write=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=Vector(0.,Solution(dom))
   v=Vector(0.,Solution(dom))
   a=Vector(0.,Solution(dom))
   a2=Vector(0.,Solution(dom))
   v=Vector(0.,Solution(dom))

   starttime = time.clock()
   while t<t_end and n<n_end:
     if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),
     # prediction:
     u_pr=u+dt*v+(dt**2/2)*a+(dt**3/6)*a2
     v_pr=v+dt*a+(dt**2/2)*a2
     a_pr=a+dt*a2
     # ... get current stress ....
     eps=symmetric(grad(u_pr))
     stress=lmbd*trace(eps)*k+2*mu*eps
     # ... force due to event:
     if abs(t-tc)<5*tc_length:
        F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event
        if output: print Lsup(F)
     else:
        if output: print 0.
     # ... get new acceleration ....
     mypde.setValue(X=-stress,Y=F-eta*v_pr)
     a=mypde.getSolution()
     # ... get new displacement ...
     da=a-a_pr
     u=u_pr+(dt**2/12.)*da
     v=v_pr+(5*dt/12.)*da
     a2+=da/dt
     # ... save current acceleration in units of gravity and displacements 
     if output:
          if t>=t_write: 
             saveVTK("disp.%i.vtu"%n,displacement=u, amplitude=length(u))
             t_write+=dt_write
             n_write+=1
     t+=dt
     n+=1

   endtime = time.clock()
   totaltime = endtime-starttime
   global netotal
   print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)
if __name__ =="__main__":
   dom=getDomain()
   rho,mu,lmbd,eta=getMaterialProperties(dom)
   wavePropagation(dom,rho,mu,lmbd,eta)
1	gross	842	"""
2			seismic wave propagation
3
4			@var __author__: name of author
5			@var __licence__: licence agreement
6			@var __url__: url entry point on documentation
7			@var __version__: version
8			@var __date__: date of the version
9			"""
10
11			__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
12			http://www.access.edu.au
13			Primary Business: Queensland, Australia"""
14			__license__="""Licensed under the Open Software License version 3.0
15			http://www.opensource.org/licenses/osl-3.0.php"""
16			__author__="Lutz Gross, l.gross@uq.edu.au"
17			__url__="http://www.iservo.edu.au/esys/escript"
18			__version__="$Revision$"
19			__date__="$Date$"
20
21			from esys.escript import *
22			from esys.escript.linearPDEs import LinearPDE
23			from esys.finley import Brick
24	gross	843	import time
25	gross	842
26	gross	844	output=True
27			n_end=10000
28	gross	843
29	gross	865	resolution=1000. # number of elements per m in the finest region
30			o=1 # element order
31	gross	842
32	gross	865	l=100000. # width and length m (without obsorber)
33			h=30000. # height in m (without obsorber)
34			d_absorber=l*0.10 # thickness of absorbing layer
35	gross	842
36	gross	866	l_sand=20000. # thickness of sand region on surface
37			h_sand=5000. # thickness of sand layer under the water
38	gross	842
39	gross	865	l_x_water=10000. # length of water in x
40			l_y_water=10000. # length of water in y
41			h_water=2000. # depth of water region
42	gross	842
43	gross	865	x_sand=l/2-l_x_water/2-l_sand # x coordinate of location of sand region (without obsorber)
44			y_sand=l/2-l_y_water/2-l_sand # y coordinate of location of sand region (without obsorber)
45	gross	842
46
47	gross	865	# origin
48			origin={"x": -d_absorber, "y" : -d_absorber , "z" : -h-d_absorber }
49			# location and geometrical size of event reltive to origin:
50			xc=[l0.2,l0.3,-h*0.7]
51	gross	862	src_radius = 2*resolution
52	gross	842	# direction of event:
53	gross	862	event=numarray.array([0.,0.,1.])*1.e6
54	gross	842	# time and length of the event
55	gross	865	tc=2.
56			tc_length=0.5
57
58			# material properties:
59			bedrock=0
60			absorber=1
61			water=2
62			sand=3
63
64			rho_tab={}
65			rho_tab[bedrock]=8e3
66			rho_tab[absorber]=rho_tab[bedrock]
67			rho_tab[water]=1e3
68			rho_tab[sand]=5e3
69
70			mu_tab={}
71			mu_tab[bedrock]=1.7e11
72			mu_tab[absorber]=mu_tab[bedrock]
73			mu_tab[water]=0.
74			mu_tab[sand]=1.5e10
75
76			lmbd_tab={}
77			lmbd_tab[bedrock]=1.7e11
78			lmbd_tab[absorber]=lmbd_tab[bedrock]
79			lmbd_tab[water]=1.e9
80			lmbd_tab[sand]=1.5e10
81
82			eta_tab={}
83			eta_tab[absorber]=-log(0.05)sqrt(rho_tab[absorber](lmbd_tab[absorber]+2*mu_tab[absorber]))/d_absorber
84			eta_tab[sand]=eta_tab[absorber]/40.
85			eta_tab[water]=eta_tab[absorber]/40.
86			eta_tab[bedrock]=eta_tab[absorber]/40.
87
88	gross	866
89			# material properties:
90			bedrock=0
91			absorber=1
92			water=2
93			sand=3
94
95			rho={}
96			rho[bedrock]=8e3
97			rho[absorber]=rho[bedrock]
98			rho[water]=1e3
99			rho[sand]=5e3
100
101			mu={}
102			mu[bedrock]=1.7e11
103			mu[absorber]=mu[bedrock]
104			mu[water]=0.
105			mu[sand]=1.5e10
106
107			lmbd={}
108			lmbd[bedrock]=1.7e11
109			lmbd_absorber=lmbd[bedrock]
110			lmbd[water]=1.e9
111			lmbd[sand]=1.5e10
112
113			eta={}
114			eta[absorber]=-log(0.05)sqrt(rho[absorber](lmbd_absorber+2*mu[absorber]))/d_absorber
115			eta[sand]=eta[absorber]/40.
116			eta[water]=eta[absorber]/40.
117			eta[bedrock]=eta[absorber]/40.
118
119	gross	843	if output:
120			print "event location = ",xc
121			print "radius of event = ",src_radius
122			print "time of event = ",tc
123			print "length of event = ",tc_length
124			print "direction = ",event
125	gross	842
126	gross	862	t_end=30.
127	gross	865	dt_write=0.1
128	gross	842
129	gross	865
130	gross	842	def getDomain():
131			"""
132			this defines a dom as a brick of length and width l and hight h
133
134
135			"""
136	gross	843	global netotal
137	gross	844
138	gross	865	v_p={}
139			for tag in rho_tab.keys():
140			v_p[tag]=sqrt((2*mu_tab[tag]+lmbd_tab[tag])/rho_tab[tag])
141			v_p_ref=min(v_p.values())
142			print "velocities: bedrock = %s, sand = %s, water =%s, absorber =%s, reference =%s"%(v_p[bedrock],v_p[sand],v_p[water],v_p[absorber],v_p_ref)
143	gross	844
144	gross	865	sections={}
145			sections["x"]=[d_absorber, x_sand, l_sand, l_x_water, l_sand, l-x_sand-2*l_sand-l_x_water, d_absorber]
146			sections["y"]=[d_absorber, y_sand, l_sand, l_y_water, l_sand, l-y_sand-2*l_sand-l_y_water, d_absorber]
147			sections["z"]=[d_absorber,h-h_water-h_sand,h_sand,h_water]
148			if output:
149			print "sections x = ",sections["x"]
150			print "sections y = ",sections["y"]
151			print "sections z = ",sections["z"]
152	gross	844
153	gross	865	mats= [
154			[ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
155			[absorber, absorber, absorber, absorber, absorber, absorber, absorber],
156			[absorber, absorber, absorber, absorber, absorber, absorber, absorber],
157			[absorber, absorber, absorber, absorber, absorber, absorber, absorber],
158			[absorber, absorber, absorber, absorber, absorber, absorber, absorber],
159			[absorber, absorber, absorber, absorber, absorber, absorber, absorber],
160			[absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
161	gross	844
162	gross	865	[ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
163			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
164			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
165			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
166			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
167			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
168			[absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
169	gross	844
170	gross	865	[ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
171			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
172			[absorber, bedrock , sand , sand , sand , bedrock , absorber],
173			[absorber, bedrock , sand , sand , sand , bedrock , absorber],
174			[absorber, bedrock , sand , sand , sand , bedrock , absorber],
175			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
176			[absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
177
178			[ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
179			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
180			[absorber, bedrock , sand , sand , sand , bedrock , absorber],
181			[absorber, bedrock , sand , water , sand , bedrock , absorber],
182			[absorber, bedrock , sand , sand , sand , bedrock , absorber],
183			[absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
184			[absorber, absorber, absorber, absorber, absorber, absorber, absorber] ] ]
185	gross	844
186	gross	865	num_elem={}
187			for d in sections:
188			num_elem[d]=[]
189			for i in range(len(sections[d])):
190			if d=="x":
191			v_p_min=v_p[mats[0][0][i]]
192			for q in range(len(sections["y"])):
193			for r in range(len(sections["z"])):
194			v_p_min=min(v_p[mats[r][q][i]],v_p_min)
195			elif d=="y":
196			v_p_min=v_p[mats[0][i][0]]
197			for q in range(len(sections["x"])):
198			for r in range(len(sections["z"])):
199			v_p_min=min(v_p[mats[r][i][q]],v_p_min)
200			elif d=="z":
201			v_p_min=v_p[mats[i][0][0]]
202			for q in range(len(sections["x"])):
203			for r in range(len(sections["y"])):
204			v_p_min=min(v_p[mats[i][r][q]],v_p_min)
205			num_elem[d].append(max(1,int(sections[d][i] * v_p_ref/v_p_min /resolution+0.5)))
206
207			ne_x=sum(num_elem["x"])
208			ne_y=sum(num_elem["y"])
209			ne_z=sum(num_elem["z"])
210			netotal=ne_xne_yne_z
211			if output: print "grid : %s x %s x %s (%s elements)"%(ne_x,ne_y,ne_z,netotal)
212			dom=Brick(ne_x,ne_y,ne_z,l0=one_x,l1=one_y,l2=o*ne_z,order=o)
213			x_old=dom.getX()
214			x_new=0
215	gross	844
216	gross	865	for d in sections:
217			if d=="x":
218			i=0
219			f=[1,0,0]
220			if d=="y":
221			i=1
222			f=[0,1,0]
223			if d=="z":
224			i=2
225			f=[0,0,1]
226			x=x_old[i]
227	gross	844
228	gross	865	p=origin[d]
229			ne=0
230			s=0.
231
232			for i in range(len(sections[d])-1):
233			msk=whereNonPositive(x-o*ne+0.5)
234			s=smsk + (sections[d][i]/(onum_elem[d][i])(x-one)+p)*(1.-msk)
235			ne+=num_elem[d][i]
236			p+=sections[d][i]
237			x_new=x_new + s * f
238			dom.setX(x_new)
239	gross	844
240	gross	842	fs=Function(dom)
241			x=Function(dom).getX()
242	gross	865	x0=x[0]
243			x1=x[1]
244			x2=x[2]
245			p_z=origin["z"]
246			for i in range(len(mats)):
247			f_z=wherePositive(x2-p_z)*wherePositive(x2-p_z+sections["z"][i])
248			p_y=origin["y"]
249			for j in range(len(mats[i])):
250			f_y=wherePositive(x1-p_y)*wherePositive(x1-p_z+sections["y"][j])
251			p_x=origin["x"]
252			for k in range(len(mats[i][j])):
253			f_x=wherePositive(x0-p_x)*wherePositive(x0-p_x+sections["x"][k])
254			fs.setTags(mats[i][j][k],f_xf_yf_z)
255			p_x+=sections["x"][k]
256			p_y+=sections["y"][j]
257			p_z+=sections["z"][i]
258	gross	842	return dom
259
260			def getMaterialProperties(dom):
261	gross	865	rho =Scalar(rho_tab[bedrock],Function(dom))
262			eta =Scalar(eta_tab[bedrock],Function(dom))
263			mu =Scalar(mu_tab[bedrock],Function(dom))
264			lmbd=Scalar(lmbd_tab[bedrock],Function(dom))
265			tags=Scalar(bedrock,Function(dom))
266
267			for tag in rho_tab.keys():
268			rho.setTaggedValue(tag,rho_tab[tag])
269			eta.setTaggedValue(tag,eta_tab[tag])
270			mu.setTaggedValue(tag,mu_tab[tag])
271			lmbd.setTaggedValue(tag,lmbd_tab[tag])
272			tags.setTaggedValue(tag,tag)
273			return rho,mu,lmbd,eta
274	gross	842
275	gross	865	def wavePropagation(dom,rho,mu,lmbd,eta):
276	gross	842	x=Function(dom).getX()
277			# ... open new PDE ...
278			mypde=LinearPDE(dom)
279			mypde.setSolverMethod(LinearPDE.LUMPING)
280			k=kronecker(Function(dom))
281			mypde.setValue(D=k*rho)
282
283	gross	865	dt=(1./5.)inf(dom.getSize()/sqrt((2mu+lmbd)/rho))
284	gross	843	if output: print "time step size = ",dt
285	gross	842	# ... set initial values ....
286			n=0
287			t=0
288	gross	865	t_write=0.
289			n_write=0
290	gross	842	# initial value of displacement at point source is constant (U0=0.01)
291			# for first two time steps
292	gross	866	u=Vector(0.,Solution(dom))
293	gross	865	v=Vector(0.,Solution(dom))
294	gross	866	a=Vector(0.,Solution(dom))
295			a2=Vector(0.,Solution(dom))
296			v=Vector(0.,Solution(dom))
297	gross	842
298	gross	843	starttime = time.clock()
299			while t<t_end and n<n_end:
300			if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),
301	gross	866	# prediction:
302			u_pr=u+dtv+(dt2/2)a+(dt*3/6)a2
303			v_pr=v+dta+(dt2/2)a2
304			a_pr=a+dt*a2
305	gross	842	# ... get current stress ....
306	gross	866	eps=symmetric(grad(u_pr))
307	gross	865	stress=lmbdtrace(eps)k+2mueps
308	gross	842	# ... force due to event:
309	gross	866	if abs(t-tc)<5*tc_length:
310			F=exp(-((t-tc)/tc_length)*2)exp(-(length(x-xc)/src_radius)*2)event
311			if output: print Lsup(F)
312			else:
313			if output: print 0.
314	gross	842	# ... get new acceleration ....
315	gross	866	mypde.setValue(X=-stress,Y=F-eta*v_pr)
316	gross	842	a=mypde.getSolution()
317			# ... get new displacement ...
318	gross	866	da=a-a_pr
319			u=u_pr+(dt*2/12.)da
320			v=v_pr+(5dt/12.)da
321			a2+=da/dt
322	gross	842	# ... save current acceleration in units of gravity and displacements
323	gross	843	if output:
324	gross	862	if t>=t_write:
325	gross	865	saveVTK("disp.%i.vtu"%n,displacement=u, amplitude=length(u))
326			t_write+=dt_write
327			n_write+=1
328	gross	842	t+=dt
329			n+=1
330
331	gross	843	endtime = time.clock()
332			totaltime = endtime-starttime
333			global netotal
334			print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)
335	gross	842	if __name__ =="__main__":
336			dom=getDomain()
337	gross	865	rho,mu,lmbd,eta=getMaterialProperties(dom)
338			wavePropagation(dom,rho,mu,lmbd,eta)