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

Contents of /trunk/doc/examples/cookbook/example09a.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3387 - (show annotations)
Thu Nov 25 07:09:23 2010 UTC (8 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 6469 byte(s)
Fixed some annoying typos.

1
2 ########################################################
3 #
4 # Copyright (c) 2009-2010 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) 2009-2010 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 ############################################################FILE HEADER
23 # example09.py
24 # Antony Hallam
25 # Seismic Wave Equation Simulation using acceleration solution.
26 # 3D model with multiple layers.
27
28 #######################################################EXTERNAL MODULES
29 from esys.escript import *
30 from esys.finley import Rectangle
31 import os
32 # smoothing operator
33 from esys.escript.pdetools import Projector, Locator
34 from esys.escript.unitsSI import *
35 import numpy as np
36 import matplotlib
37 matplotlib.use('agg') #It's just here for automated testing
38
39 import pylab as pl
40 import matplotlib.cm as cm
41 from esys.escript.linearPDEs import LinearPDE
42 from esys.finley import ReadMesh
43 from esys.weipa import saveVTK
44
45 ########################################################MPI WORLD CHECK
46 if getMPISizeWorld() > 1:
47 import sys
48 print "This example will not run in an MPI world."
49 sys.exit(0)
50
51 #################################################ESTABLISHING VARIABLES
52 # where to save output data
53 savepath = "data/example09"
54 mkDir(savepath)
55 #Geometric and material property related variables.
56 mx = 200. # model lenght
57 my = 200. # model width
58 mz=100.0
59 step=10.0 # the element size
60 ndx = int(mx/step) # steps in x direction
61 ndy = int(my/step) # steps in y direction
62 ndz = int(mz/step)
63
64 vel2=1800.; vel1=3000.
65 rho2=2300.; rho1=3100. #density
66 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus
67 lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant
68
69 # Time related variables.
70 testing=True
71 if testing:
72 print 'The testing end time is currently selected. This severely limits the number of time iterations.'
73 print "Try changing testing to False for more iterations."
74 tend=0.001
75 else:
76 tend=0.1 # end time
77
78 h=0.0001 # time step
79 # data recording times
80 rtime=0.0 # first time to record
81 rtime_inc=tend/50.0 # time increment to record
82 #Check to make sure number of time steps is not too large.
83 print "Time step size= ",h, "Expected number of outputs= ",tend/h
84
85 U0=0.1 # amplitude of point source
86 ls=500 # length of the source
87 source=np.zeros(ls,'float') # source array
88 decay1=np.zeros(ls,'float') # decay curve one
89 decay2=np.zeros(ls,'float') # decay curve two
90 time=np.zeros(ls,'float') # time values
91 g=np.log(0.01)/ls
92
93 dfeq=50 #Dominant Frequency
94 a = 2.0 * (np.pi * dfeq)**2.0
95 t0 = 5.0 / (2.0 * np.pi * dfeq)
96 srclength = 5. * t0
97 ls = int(srclength/h)
98 print 'source length',ls
99 source=np.zeros(ls,'float') # source array
100 ampmax=0
101 for it in range(0,ls):
102 t = it*h
103 tt = t-t0
104 dum1 = np.exp(-a * tt * tt)
105 source[it] = -2. * a * tt * dum1
106 if (abs(source[it]) > ampmax):
107 ampmax = abs(source[it])
108 time[t]=t*h
109 #tdecay=decay1*decay2*U0
110 #decay1=decay1*U0; decay2=decay2*U0
111 #pl.clf();
112 #pl.plot(source)
113 #pl.plot(time,decay1);pl.plot(time,decay2);
114 #pl.plot(time,tdecay)
115 #pl.savefig(os.path.join(savepath,'source.png'))
116
117 # will introduce a spherical source at middle left of bottom face
118 xc=[50,50,0]
119
120 ####################################################DOMAIN CONSTRUCTION
121 domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain
122 x=domain.getX() # get the locations of the nodes in the domain
123
124 lam=Scalar(0,Function(domain))
125 lam.setTaggedValue("vintfa",lam1)
126 lam.setTaggedValue("vintfb",lam2)
127 mu=Scalar(0,Function(domain))
128 mu.setTaggedValue("vintfa",mu1)
129 mu.setTaggedValue("vintfb",mu2)
130 rho=Scalar(0,Function(domain))
131 rho.setTaggedValue("vintfa",rho1)
132 rho.setTaggedValue("vintfb",rho2)
133
134 ##########################################################ESTABLISH PDE
135 mypde=LinearPDE(domain) # create pde
136 mypde.setSymmetryOn() # turn symmetry on
137 # turn lumping on for more efficient solving
138 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
139 kmat = kronecker(domain) # create the kronecker delta function of the domain
140 mypde.setValue(D=rho*kmat) #set the general form value D
141
142 ############################################FIRST TIME STEPS AND SOURCE
143 # define small radius around point xc
144 src_length = 20; print "src_length = ",src_length
145 # set initial values for first two time steps with source terms
146 xb=FunctionOnBoundary(domain).getX()
147 #sy=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
148 y=Vector(0.0,FunctionOnBoundary(domain))
149
150 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
151
152 #sy=sy*src_dir
153 #sy.setTaggedValue("stop")
154 y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1))
155 mypde.setValue(y=y) #set the source as a function on the boundary
156 # initial value of displacement at point source is constant (U0=0.01)
157 # for first two time steps
158 u=[0.0,0.0,0.0]*wherePositive(x)
159 u_m1=u
160
161 ####################################################ITERATION VARIABLES
162 n=0 # iteration counter
163 t=0 # time counter
164 ##############################################################ITERATION
165 while t<tend:
166 # get current stress
167 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
168 mypde.setValue(X=-stress) # set PDE values
169 accel = mypde.getSolution() #get PDE solution for accelleration
170 u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
171 u_p1=u_p1#*abc # apply boundary conditions
172 u_m1=u; u=u_p1 # shift values by 1
173 # save current displacement, acceleration and pressure
174 if (t >= rtime):
175 saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\
176 acceleration=length(accel),tensor=stress)
177 rtime=rtime+rtime_inc #increment data save time
178 # increment loop values
179 t=t+h; n=n+1
180 if (n < ls):
181 y.setTaggedValue("stop",source[n]*src_dir)
182 mypde.setValue(y=y) #set the source as a function on the boundary
183 print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26