/[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 3405 - (hide annotations)
Wed Dec 8 01:09:06 2010 UTC (8 years, 9 months ago) by ahallam
File MIME type: text/x-python
File size: 6444 byte(s)
Changes to example09 for memory footprint. One change too help wording in layer_cake function. example09b removed from testing. Not for impending release.
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 3405 savepath = "data/example09a"
54 ahallam 3389 meshpath = "data/example09m"
55 ahallam 3055 mkDir(savepath)
56     #Geometric and material property related variables.
57 ahallam 3089 step=10.0 # the element size
58 ahallam 3069
59     vel2=1800.; vel1=3000.
60     rho2=2300.; rho1=3100. #density
61 ahallam 3089 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus
62     lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant
63 ahallam 3069
64 ahallam 3405 ####################################################TESTING SWITCH
65 ahallam 3195 testing=True
66     if testing:
67 ahallam 3405 print 'The testing end time is currently selected. This severely limits the number of time iterations.'
68     print "Try changing testing to False for more iterations."
69     tend=0.001
70     #Model Parameters
71     mx=40.
72     my=40.
73     mz=20.
74     outputs=5
75 ahallam 3195 else:
76 ahallam 3405 tend=0.1 # end time
77     #Model Parameters
78     mx=100.0 #x width of model
79     my=100.0 #y width of model
80     mz=50.0 #depth of model
81     outputs=200
82 ahallam 3195
83 ahallam 3405 ####################################################TIME RELATED VARIABLES
84     h=0.00005 # time step
85 ahallam 3055 # data recording times
86     rtime=0.0 # first time to record
87 ahallam 3405 rtime_inc=tend/outputs # time increment to record
88 ahallam 3055 #Check to make sure number of time steps is not too large.
89     print "Time step size= ",h, "Expected number of outputs= ",tend/h
90    
91 ahallam 3405 ####################################################CREATING THE SOURCE FUNCTION
92 ahallam 3055 U0=0.1 # amplitude of point source
93     ls=500 # length of the source
94     source=np.zeros(ls,'float') # source array
95     decay1=np.zeros(ls,'float') # decay curve one
96     decay2=np.zeros(ls,'float') # decay curve two
97     time=np.zeros(ls,'float') # time values
98     g=np.log(0.01)/ls
99    
100     dfeq=50 #Dominant Frequency
101     a = 2.0 * (np.pi * dfeq)**2.0
102     t0 = 5.0 / (2.0 * np.pi * dfeq)
103     srclength = 5. * t0
104     ls = int(srclength/h)
105     print 'source length',ls
106     source=np.zeros(ls,'float') # source array
107     ampmax=0
108     for it in range(0,ls):
109     t = it*h
110     tt = t-t0
111     dum1 = np.exp(-a * tt * tt)
112     source[it] = -2. * a * tt * dum1
113     if (abs(source[it]) > ampmax):
114     ampmax = abs(source[it])
115     time[t]=t*h
116    
117     # will introduce a spherical source at middle left of bottom face
118 ahallam 3405 xc=[mx/2,my/2,0]
119 ahallam 3055
120     ####################################################DOMAIN CONSTRUCTION
121 ahallam 3389 domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
122 ahallam 3069 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 3405 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_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 3405 src_rad = 20; print "sourc radius= ",src_rad
145 ahallam 3055 # set initial values for first two time steps with source terms
146 ahallam 3069 xb=FunctionOnBoundary(domain).getX()
147 ahallam 3405 yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
148     stop=Scalar(0.0,FunctionOnBoundary(domain))
149     stop.setTaggedValue("stop",1.0)
150 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
151 ahallam 3089
152 ahallam 3405 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
153 ahallam 3055 # initial value of displacement at point source is constant (U0=0.01)
154     # for first two time steps
155 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x)
156 ahallam 3055 u_m1=u
157    
158     ####################################################ITERATION VARIABLES
159     n=0 # iteration counter
160     t=0 # time counter
161     ##############################################################ITERATION
162     while t<tend:
163     # get current stress
164 ahallam 3069 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
165 ahallam 3055 mypde.setValue(X=-stress) # set PDE values
166     accel = mypde.getSolution() #get PDE solution for accelleration
167     u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
168 ahallam 3069 u_p1=u_p1#*abc # apply boundary conditions
169 ahallam 3055 u_m1=u; u=u_p1 # shift values by 1
170     # save current displacement, acceleration and pressure
171     if (t >= rtime):
172 ahallam 3405 saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
173 ahallam 3055 acceleration=length(accel),tensor=stress)
174     rtime=rtime+rtime_inc #increment data save time
175     # increment loop values
176     t=t+h; n=n+1
177     if (n < ls):
178 ahallam 3405 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
179 ahallam 3055 print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26