/[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 3389 - (hide annotations)
Wed Dec 1 23:03:35 2010 UTC (8 years, 9 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 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 3389 meshpath = "data/example09m"
55 ahallam 3055 mkDir(savepath)
56     #Geometric and material property related variables.
57 ahallam 3089 mx = 200. # model lenght
58     my = 200. # model width
59     mz=100.0
60     step=10.0 # the element size
61 ahallam 3069 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 ahallam 3089 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 ahallam 3069
70 ahallam 3055 # Time related variables.
71 ahallam 3195 testing=True
72     if testing:
73 caltinay 3387 print 'The testing end time is currently selected. This severely limits the number of time iterations.'
74 ahallam 3195 print "Try changing testing to False for more iterations."
75     tend=0.001
76     else:
77     tend=0.1 # end time
78    
79 ahallam 3089 h=0.0001 # time step
80 ahallam 3055 # 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 ahallam 3089 #pl.clf();
113     #pl.plot(source)
114 ahallam 3055 #pl.plot(time,decay1);pl.plot(time,decay2);
115     #pl.plot(time,tdecay)
116 ahallam 3089 #pl.savefig(os.path.join(savepath,'source.png'))
117 ahallam 3055
118     # will introduce a spherical source at middle left of bottom face
119 ahallam 3089 xc=[50,50,0]
120 ahallam 3055
121     ####################################################DOMAIN CONSTRUCTION
122 ahallam 3389 domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
123 ahallam 3069 x=domain.getX() # get the locations of the nodes in the domain
124 ahallam 3055
125 ahallam 3069 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 ahallam 3055 ##########################################################ESTABLISH PDE
136     mypde=LinearPDE(domain) # create pde
137     mypde.setSymmetryOn() # turn symmetry on
138     # turn lumping on for more efficient solving
139 ahallam 3075 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
140 ahallam 3055 kmat = kronecker(domain) # create the kronecker delta function of the domain
141 ahallam 3069 mypde.setValue(D=rho*kmat) #set the general form value D
142 ahallam 3055
143     ############################################FIRST TIME STEPS AND SOURCE
144     # define small radius around point xc
145 ahallam 3089 src_length = 20; print "src_length = ",src_length
146 ahallam 3055 # set initial values for first two time steps with source terms
147 ahallam 3069 xb=FunctionOnBoundary(domain).getX()
148 ahallam 3089 #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 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
152 ahallam 3089
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 ahallam 3055 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 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x)
160 ahallam 3055 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 ahallam 3069 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
169 ahallam 3055 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 ahallam 3069 u_p1=u_p1#*abc # apply boundary conditions
173 ahallam 3055 u_m1=u; u=u_p1 # shift values by 1
174     # save current displacement, acceleration and pressure
175     if (t >= rtime):
176 ahallam 3069 saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\
177 ahallam 3055 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 ahallam 3089 y.setTaggedValue("stop",source[n]*src_dir)
183     mypde.setValue(y=y) #set the source as a function on the boundary
184 ahallam 3055 print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26