/[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 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 6451 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1
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 import matplotlib
37 matplotlib.use('agg') #It's just here for automated testing
38
39 import pylab as pl
40 import matplotlib.cm as cm
41 from esys.escript.linearPDEs import LinearPDE
42 from esys.finley import ReadMesh
43 from esys.weipa import saveVTK
44
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 savepath = "data/example09a"
54 meshpath = "data/example09m"
55 mkDir(savepath)
56 #Geometric and material property related variables.
57 step=10.0 # the element size
58
59 vel2=1800.; vel1=3000.
60 rho2=2300.; rho1=3100. #density
61 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
64 ####################################################TESTING SWITCH
65 testing=True
66 if testing:
67 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 else:
76 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
83 ####################################################TIME RELATED VARIABLES
84 h=0.00005 # time step
85 # data recording times
86 rtime=0.0 # first time to record
87 rtime_inc=tend/outputs # time increment to record
88 #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 ####################################################CREATING THE SOURCE FUNCTION
92 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 xc=[mx/2,my/2,0]
119
120 ####################################################DOMAIN CONSTRUCTION
121 domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
122 x=domain.getX() # get the locations of the nodes in the domain
123
124 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 ##########################################################ESTABLISH PDE
135 mypde=LinearPDE(domain) # create pde
136 mypde.setSymmetryOn() # turn symmetry on
137 # turn lumping on for more efficient solving
138 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
139 kmat = kronecker(domain) # create the kronecker delta function of the domain
140 mypde.setValue(D=rho*kmat) #set the general form value D
141
142 ############################################FIRST TIME STEPS AND SOURCE
143 # define small radius around point xc
144 src_rad = 20; print("sourc radius= ",src_rad)
145 # set initial values for first two time steps with source terms
146 xb=FunctionOnBoundary(domain).getX()
147 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 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
151
152 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
153 # initial value of displacement at point source is constant (U0=0.01)
154 # for first two time steps
155 u=[0.0,0.0,0.0]*wherePositive(x)
156 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 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
165 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 u_p1=u_p1#*abc # apply boundary conditions
169 u_m1=u; u=u_p1 # shift values by 1
170 # save current displacement, acceleration and pressure
171 if (t >= rtime):
172 saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
173 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 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
179 print(n,"-th time step t ",t)

  ViewVC Help
Powered by ViewVC 1.1.26