/[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 4576 - (hide annotations)
Mon Dec 9 23:35:30 2013 UTC (5 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 6547 byte(s)
python3ified things, replaced mixed whitespace and xrange calls
1 ahallam 3055
2 jfenwick 3981 ##############################################################################
3 ahallam 3055 #
4 jfenwick 4154 # Copyright (c) 2009-2013 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ahallam 3055 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 ahallam 3055
16 jfenwick 4154 __copyright__="""Copyright (c) 2009-2013 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 ahallam 3055 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 caltinay 4087 import matplotlib
31     matplotlib.use('agg') #It's just here for automated testing
32 ahallam 3055 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 jfenwick 3148
40 ahallam 3055 import pylab as pl
41     import matplotlib.cm as cm
42     from esys.escript.linearPDEs import LinearPDE
43 ahallam 3069 from esys.finley import ReadMesh
44 caltinay 3346 from esys.weipa import saveVTK
45 ahallam 3055
46     ########################################################MPI WORLD CHECK
47     if getMPISizeWorld() > 1:
48 sshaw 4576 import sys
49     print("This example will not run in an MPI world.")
50     sys.exit(0)
51 ahallam 3055
52     #################################################ESTABLISHING VARIABLES
53     # where to save output data
54 ahallam 3405 savepath = "data/example09a"
55 ahallam 3389 meshpath = "data/example09m"
56 ahallam 3055 mkDir(savepath)
57     #Geometric and material property related variables.
58 ahallam 3089 step=10.0 # the element size
59 ahallam 3069
60     vel2=1800.; vel1=3000.
61     rho2=2300.; rho1=3100. #density
62 ahallam 3089 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 ahallam 3069
65 ahallam 3405 ####################################################TESTING SWITCH
66 ahallam 3195 testing=True
67     if testing:
68 jfenwick 3892 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 ahallam 3405 tend=0.001
71     #Model Parameters
72     mx=40.
73     my=40.
74     mz=20.
75     outputs=5
76 ahallam 3195 else:
77 ahallam 3405 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 ahallam 3195
84 ahallam 3405 ####################################################TIME RELATED VARIABLES
85     h=0.00005 # time step
86 ahallam 3055 # data recording times
87     rtime=0.0 # first time to record
88 ahallam 3405 rtime_inc=tend/outputs # time increment to record
89 ahallam 3055 #Check to make sure number of time steps is not too large.
90 jfenwick 3892 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
91 ahallam 3055
92 ahallam 3405 ####################################################CREATING THE SOURCE FUNCTION
93 ahallam 3055 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 jfenwick 3892 print('source length',ls)
107 ahallam 3055 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 ahallam 3405 xc=[mx/2,my/2,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 3405 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_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 jfenwick 3892 src_rad = 20; print("sourc radius= ",src_rad)
146 ahallam 3055 # set initial values for first two time steps with source terms
147 ahallam 3069 xb=FunctionOnBoundary(domain).getX()
148 ahallam 3405 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 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
152 ahallam 3089
153 ahallam 3405 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
154 ahallam 3055 # initial value of displacement at point source is constant (U0=0.01)
155     # for first two time steps
156 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x)
157 ahallam 3055 u_m1=u
158    
159     ####################################################ITERATION VARIABLES
160     n=0 # iteration counter
161     t=0 # time counter
162     ##############################################################ITERATION
163     while t<tend:
164 sshaw 4576 # get current stress
165 ahallam 3069 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
166 ahallam 3055 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 sshaw 4576 u_p1=u_p1#*abc # apply boundary conditions
170 ahallam 3055 u_m1=u; u=u_p1 # shift values by 1
171     # save current displacement, acceleration and pressure
172     if (t >= rtime):
173 ahallam 3405 saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
174 ahallam 3055 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 ahallam 3405 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
180 caltinay 4005 print("time step %d, t=%s"%(n,t))

  ViewVC Help
Powered by ViewVC 1.1.26