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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3387 - (hide 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 ahallam 3055
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 jfenwick 3148 import matplotlib
37     matplotlib.use('agg') #It's just here for automated testing
38    
39 ahallam 3055 import pylab as pl
40     import matplotlib.cm as cm
41     from esys.escript.linearPDEs import LinearPDE
42 ahallam 3069 from esys.finley import ReadMesh
43 caltinay 3346 from esys.weipa import saveVTK
44 ahallam 3055
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 ahallam 3069 savepath = "data/example09"
54 ahallam 3055 mkDir(savepath)
55     #Geometric and material property related variables.
56 ahallam 3089 mx = 200. # model lenght
57     my = 200. # model width
58     mz=100.0
59     step=10.0 # the element size
60 ahallam 3069 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 ahallam 3089 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 ahallam 3069
69 ahallam 3055 # Time related variables.
70 ahallam 3195 testing=True
71     if testing:
72 caltinay 3387 print 'The testing end time is currently selected. This severely limits the number of time iterations.'
73 ahallam 3195 print "Try changing testing to False for more iterations."
74     tend=0.001
75     else:
76     tend=0.1 # end time
77    
78 ahallam 3089 h=0.0001 # time step
79 ahallam 3055 # 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 ahallam 3089 #pl.clf();
112     #pl.plot(source)
113 ahallam 3055 #pl.plot(time,decay1);pl.plot(time,decay2);
114     #pl.plot(time,tdecay)
115 ahallam 3089 #pl.savefig(os.path.join(savepath,'source.png'))
116 ahallam 3055
117     # will introduce a spherical source at middle left of bottom face
118 ahallam 3089 xc=[50,50,0]
119 ahallam 3055
120     ####################################################DOMAIN CONSTRUCTION
121 ahallam 3069 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 ahallam 3055
124 ahallam 3069 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 ahallam 3055 ##########################################################ESTABLISH PDE
135     mypde=LinearPDE(domain) # create pde
136     mypde.setSymmetryOn() # turn symmetry on
137     # turn lumping on for more efficient solving
138 ahallam 3075 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
139 ahallam 3055 kmat = kronecker(domain) # create the kronecker delta function of the domain
140 ahallam 3069 mypde.setValue(D=rho*kmat) #set the general form value D
141 ahallam 3055
142     ############################################FIRST TIME STEPS AND SOURCE
143     # define small radius around point xc
144 ahallam 3089 src_length = 20; print "src_length = ",src_length
145 ahallam 3055 # set initial values for first two time steps with source terms
146 ahallam 3069 xb=FunctionOnBoundary(domain).getX()
147 ahallam 3089 #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 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
151 ahallam 3089
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 ahallam 3055 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 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x)
159 ahallam 3055 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 ahallam 3069 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
168 ahallam 3055 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 ahallam 3069 u_p1=u_p1#*abc # apply boundary conditions
172 ahallam 3055 u_m1=u; u=u_p1 # shift values by 1
173     # save current displacement, acceleration and pressure
174     if (t >= rtime):
175 ahallam 3069 saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\
176 ahallam 3055 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 ahallam 3089 y.setTaggedValue("stop",source[n]*src_dir)
182     mypde.setValue(y=y) #set the source as a function on the boundary
183 ahallam 3055 print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26