/[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 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (4 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 6547 byte(s)
python3ified things, replaced mixed whitespace and xrange calls
1
2 ##############################################################################
3 #
4 # Copyright (c) 2009-2013 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2009-2013 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 ############################################################FILE HEADER
24 # example09.py
25 # Antony Hallam
26 # Seismic Wave Equation Simulation using acceleration solution.
27 # 3D model with multiple layers.
28
29 #######################################################EXTERNAL MODULES
30 import matplotlib
31 matplotlib.use('agg') #It's just here for automated testing
32 from esys.escript import *
33 from esys.finley import Rectangle
34 import os
35 # smoothing operator
36 from esys.escript.pdetools import Projector, Locator
37 from esys.escript.unitsSI import *
38 import numpy as np
39
40 import pylab as pl
41 import matplotlib.cm as cm
42 from esys.escript.linearPDEs import LinearPDE
43 from esys.finley import ReadMesh
44 from esys.weipa import saveVTK
45
46 ########################################################MPI WORLD CHECK
47 if getMPISizeWorld() > 1:
48 import sys
49 print("This example will not run in an MPI world.")
50 sys.exit(0)
51
52 #################################################ESTABLISHING VARIABLES
53 # where to save output data
54 savepath = "data/example09a"
55 meshpath = "data/example09m"
56 mkDir(savepath)
57 #Geometric and material property related variables.
58 step=10.0 # the element size
59
60 vel2=1800.; vel1=3000.
61 rho2=2300.; rho1=3100. #density
62 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus
63 lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant
64
65 ####################################################TESTING SWITCH
66 testing=True
67 if testing:
68 print('The testing end time is currently selected. This severely limits the number of time iterations.')
69 print("Try changing testing to False for more iterations.")
70 tend=0.001
71 #Model Parameters
72 mx=40.
73 my=40.
74 mz=20.
75 outputs=5
76 else:
77 tend=0.1 # end time
78 #Model Parameters
79 mx=100.0 #x width of model
80 my=100.0 #y width of model
81 mz=50.0 #depth of model
82 outputs=200
83
84 ####################################################TIME RELATED VARIABLES
85 h=0.00005 # time step
86 # data recording times
87 rtime=0.0 # first time to record
88 rtime_inc=tend/outputs # time increment to record
89 #Check to make sure number of time steps is not too large.
90 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
91
92 ####################################################CREATING THE SOURCE FUNCTION
93 U0=0.1 # amplitude of point source
94 ls=500 # length of the source
95 source=np.zeros(ls,'float') # source array
96 decay1=np.zeros(ls,'float') # decay curve one
97 decay2=np.zeros(ls,'float') # decay curve two
98 time=np.zeros(ls,'float') # time values
99 g=np.log(0.01)/ls
100
101 dfeq=50 #Dominant Frequency
102 a = 2.0 * (np.pi * dfeq)**2.0
103 t0 = 5.0 / (2.0 * np.pi * dfeq)
104 srclength = 5. * t0
105 ls = int(srclength/h)
106 print('source length',ls)
107 source=np.zeros(ls,'float') # source array
108 ampmax=0
109 for it in range(0,ls):
110 t = it*h
111 tt = t-t0
112 dum1 = np.exp(-a * tt * tt)
113 source[it] = -2. * a * tt * dum1
114 if (abs(source[it]) > ampmax):
115 ampmax = abs(source[it])
116 time[t]=t*h
117
118 # will introduce a spherical source at middle left of bottom face
119 xc=[mx/2,my/2,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().HRZ_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_rad = 20; print("sourc radius= ",src_rad)
146 # set initial values for first two time steps with source terms
147 xb=FunctionOnBoundary(domain).getX()
148 yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
149 stop=Scalar(0.0,FunctionOnBoundary(domain))
150 stop.setTaggedValue("stop",1.0)
151 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
152
153 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
154 # initial value of displacement at point source is constant (U0=0.01)
155 # for first two time steps
156 u=[0.0,0.0,0.0]*wherePositive(x)
157 u_m1=u
158
159 ####################################################ITERATION VARIABLES
160 n=0 # iteration counter
161 t=0 # time counter
162 ##############################################################ITERATION
163 while t<tend:
164 # get current stress
165 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
166 mypde.setValue(X=-stress) # set PDE values
167 accel = mypde.getSolution() #get PDE solution for accelleration
168 u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
169 u_p1=u_p1#*abc # apply boundary conditions
170 u_m1=u; u=u_p1 # shift values by 1
171 # save current displacement, acceleration and pressure
172 if (t >= rtime):
173 saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
174 acceleration=length(accel),tensor=stress)
175 rtime=rtime+rtime_inc #increment data save time
176 # increment loop values
177 t=t+h; n=n+1
178 if (n < ls):
179 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
180 print("time step %d, t=%s"%(n,t))

  ViewVC Help
Powered by ViewVC 1.1.26