/[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 3389 - (show annotations)
Wed Dec 1 23:03:35 2010 UTC (8 years, 8 months ago) by ahallam
File MIME type: text/x-python
File size: 6498 byte(s)
Some changes to the cookbook examples. Starting Inversion experimentation.
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 meshpath = "data/example09m"
55 mkDir(savepath)
56 #Geometric and material property related variables.
57 mx = 200. # model lenght
58 my = 200. # model width
59 mz=100.0
60 step=10.0 # the element size
61 ndx = int(mx/step) # steps in x direction
62 ndy = int(my/step) # steps in y direction
63 ndz = int(mz/step)
64
65 vel2=1800.; vel1=3000.
66 rho2=2300.; rho1=3100. #density
67 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus
68 lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant
69
70 # Time related variables.
71 testing=True
72 if testing:
73 print 'The testing end time is currently selected. This severely limits the number of time iterations.'
74 print "Try changing testing to False for more iterations."
75 tend=0.001
76 else:
77 tend=0.1 # end time
78
79 h=0.0001 # time step
80 # data recording times
81 rtime=0.0 # first time to record
82 rtime_inc=tend/50.0 # time increment to record
83 #Check to make sure number of time steps is not too large.
84 print "Time step size= ",h, "Expected number of outputs= ",tend/h
85
86 U0=0.1 # amplitude of point source
87 ls=500 # length of the source
88 source=np.zeros(ls,'float') # source array
89 decay1=np.zeros(ls,'float') # decay curve one
90 decay2=np.zeros(ls,'float') # decay curve two
91 time=np.zeros(ls,'float') # time values
92 g=np.log(0.01)/ls
93
94 dfeq=50 #Dominant Frequency
95 a = 2.0 * (np.pi * dfeq)**2.0
96 t0 = 5.0 / (2.0 * np.pi * dfeq)
97 srclength = 5. * t0
98 ls = int(srclength/h)
99 print 'source length',ls
100 source=np.zeros(ls,'float') # source array
101 ampmax=0
102 for it in range(0,ls):
103 t = it*h
104 tt = t-t0
105 dum1 = np.exp(-a * tt * tt)
106 source[it] = -2. * a * tt * dum1
107 if (abs(source[it]) > ampmax):
108 ampmax = abs(source[it])
109 time[t]=t*h
110 #tdecay=decay1*decay2*U0
111 #decay1=decay1*U0; decay2=decay2*U0
112 #pl.clf();
113 #pl.plot(source)
114 #pl.plot(time,decay1);pl.plot(time,decay2);
115 #pl.plot(time,tdecay)
116 #pl.savefig(os.path.join(savepath,'source.png'))
117
118 # will introduce a spherical source at middle left of bottom face
119 xc=[50,50,0]
120
121 ####################################################DOMAIN CONSTRUCTION
122 domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
123 x=domain.getX() # get the locations of the nodes in the domain
124
125 lam=Scalar(0,Function(domain))
126 lam.setTaggedValue("vintfa",lam1)
127 lam.setTaggedValue("vintfb",lam2)
128 mu=Scalar(0,Function(domain))
129 mu.setTaggedValue("vintfa",mu1)
130 mu.setTaggedValue("vintfb",mu2)
131 rho=Scalar(0,Function(domain))
132 rho.setTaggedValue("vintfa",rho1)
133 rho.setTaggedValue("vintfb",rho2)
134
135 ##########################################################ESTABLISH PDE
136 mypde=LinearPDE(domain) # create pde
137 mypde.setSymmetryOn() # turn symmetry on
138 # turn lumping on for more efficient solving
139 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
140 kmat = kronecker(domain) # create the kronecker delta function of the domain
141 mypde.setValue(D=rho*kmat) #set the general form value D
142
143 ############################################FIRST TIME STEPS AND SOURCE
144 # define small radius around point xc
145 src_length = 20; print "src_length = ",src_length
146 # set initial values for first two time steps with source terms
147 xb=FunctionOnBoundary(domain).getX()
148 #sy=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
149 y=Vector(0.0,FunctionOnBoundary(domain))
150
151 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
152
153 #sy=sy*src_dir
154 #sy.setTaggedValue("stop")
155 y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1))
156 mypde.setValue(y=y) #set the source as a function on the boundary
157 # initial value of displacement at point source is constant (U0=0.01)
158 # for first two time steps
159 u=[0.0,0.0,0.0]*wherePositive(x)
160 u_m1=u
161
162 ####################################################ITERATION VARIABLES
163 n=0 # iteration counter
164 t=0 # time counter
165 ##############################################################ITERATION
166 while t<tend:
167 # get current stress
168 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
169 mypde.setValue(X=-stress) # set PDE values
170 accel = mypde.getSolution() #get PDE solution for accelleration
171 u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
172 u_p1=u_p1#*abc # apply boundary conditions
173 u_m1=u; u=u_p1 # shift values by 1
174 # save current displacement, acceleration and pressure
175 if (t >= rtime):
176 saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\
177 acceleration=length(accel),tensor=stress)
178 rtime=rtime+rtime_inc #increment data save time
179 # increment loop values
180 t=t+h; n=n+1
181 if (n < ls):
182 y.setTaggedValue("stop",source[n]*src_dir)
183 mypde.setValue(y=y) #set the source as a function on the boundary
184 print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26