/[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 3148 - (hide annotations)
Fri Sep 3 02:09:47 2010 UTC (9 years ago) by jfenwick
File MIME type: text/x-python
File size: 6228 byte(s)
Another attempt to patch the X issue

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 ahallam 3055
44     ########################################################MPI WORLD CHECK
45     if getMPISizeWorld() > 1:
46     import sys
47     print "This example will not run in an MPI world."
48     sys.exit(0)
49    
50     #################################################ESTABLISHING VARIABLES
51     # where to save output data
52 ahallam 3069 savepath = "data/example09"
53 ahallam 3055 mkDir(savepath)
54     #Geometric and material property related variables.
55 ahallam 3089 mx = 200. # model lenght
56     my = 200. # model width
57     mz=100.0
58     step=10.0 # the element size
59 ahallam 3069 ndx = int(mx/step) # steps in x direction
60     ndy = int(my/step) # steps in y direction
61     ndz = int(mz/step)
62    
63     vel2=1800.; vel1=3000.
64     rho2=2300.; rho1=3100. #density
65 ahallam 3089 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus
66     lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant
67 ahallam 3069
68 ahallam 3055 # Time related variables.
69 ahallam 3089 tend=0.1 # end time
70     h=0.0001 # time step
71 ahallam 3055 # data recording times
72     rtime=0.0 # first time to record
73     rtime_inc=tend/50.0 # time increment to record
74     #Check to make sure number of time steps is not too large.
75     print "Time step size= ",h, "Expected number of outputs= ",tend/h
76    
77     U0=0.1 # amplitude of point source
78     ls=500 # length of the source
79     source=np.zeros(ls,'float') # source array
80     decay1=np.zeros(ls,'float') # decay curve one
81     decay2=np.zeros(ls,'float') # decay curve two
82     time=np.zeros(ls,'float') # time values
83     g=np.log(0.01)/ls
84    
85     dfeq=50 #Dominant Frequency
86     a = 2.0 * (np.pi * dfeq)**2.0
87     t0 = 5.0 / (2.0 * np.pi * dfeq)
88     srclength = 5. * t0
89     ls = int(srclength/h)
90     print 'source length',ls
91     source=np.zeros(ls,'float') # source array
92     ampmax=0
93     for it in range(0,ls):
94     t = it*h
95     tt = t-t0
96     dum1 = np.exp(-a * tt * tt)
97     source[it] = -2. * a * tt * dum1
98     if (abs(source[it]) > ampmax):
99     ampmax = abs(source[it])
100     time[t]=t*h
101     #tdecay=decay1*decay2*U0
102     #decay1=decay1*U0; decay2=decay2*U0
103 ahallam 3089 #pl.clf();
104     #pl.plot(source)
105 ahallam 3055 #pl.plot(time,decay1);pl.plot(time,decay2);
106     #pl.plot(time,tdecay)
107 ahallam 3089 #pl.savefig(os.path.join(savepath,'source.png'))
108 ahallam 3055
109     # will introduce a spherical source at middle left of bottom face
110 ahallam 3089 xc=[50,50,0]
111 ahallam 3055
112     ####################################################DOMAIN CONSTRUCTION
113 ahallam 3069 domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain
114     x=domain.getX() # get the locations of the nodes in the domain
115 ahallam 3055
116 ahallam 3069 lam=Scalar(0,Function(domain))
117     lam.setTaggedValue("vintfa",lam1)
118     lam.setTaggedValue("vintfb",lam2)
119     mu=Scalar(0,Function(domain))
120     mu.setTaggedValue("vintfa",mu1)
121     mu.setTaggedValue("vintfb",mu2)
122     rho=Scalar(0,Function(domain))
123     rho.setTaggedValue("vintfa",rho1)
124     rho.setTaggedValue("vintfb",rho2)
125    
126 ahallam 3055 ##########################################################ESTABLISH PDE
127     mypde=LinearPDE(domain) # create pde
128     mypde.setSymmetryOn() # turn symmetry on
129     # turn lumping on for more efficient solving
130 ahallam 3075 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
131 ahallam 3055 kmat = kronecker(domain) # create the kronecker delta function of the domain
132 ahallam 3069 mypde.setValue(D=rho*kmat) #set the general form value D
133 ahallam 3055
134     ############################################FIRST TIME STEPS AND SOURCE
135     # define small radius around point xc
136 ahallam 3089 src_length = 20; print "src_length = ",src_length
137 ahallam 3055 # set initial values for first two time steps with source terms
138 ahallam 3069 xb=FunctionOnBoundary(domain).getX()
139 ahallam 3089 #sy=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
140     y=Vector(0.0,FunctionOnBoundary(domain))
141    
142 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
143 ahallam 3089
144     #sy=sy*src_dir
145     #sy.setTaggedValue("stop")
146     y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1))
147 ahallam 3055 mypde.setValue(y=y) #set the source as a function on the boundary
148     # initial value of displacement at point source is constant (U0=0.01)
149     # for first two time steps
150 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x)
151 ahallam 3055 u_m1=u
152    
153     ####################################################ITERATION VARIABLES
154     n=0 # iteration counter
155     t=0 # time counter
156     ##############################################################ITERATION
157     while t<tend:
158     # get current stress
159 ahallam 3069 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
160 ahallam 3055 mypde.setValue(X=-stress) # set PDE values
161     accel = mypde.getSolution() #get PDE solution for accelleration
162     u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
163 ahallam 3069 u_p1=u_p1#*abc # apply boundary conditions
164 ahallam 3055 u_m1=u; u=u_p1 # shift values by 1
165     # save current displacement, acceleration and pressure
166     if (t >= rtime):
167 ahallam 3069 saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\
168 ahallam 3055 acceleration=length(accel),tensor=stress)
169     rtime=rtime+rtime_inc #increment data save time
170     # increment loop values
171     t=t+h; n=n+1
172     if (n < ls):
173 ahallam 3089 y.setTaggedValue("stop",source[n]*src_dir)
174     mypde.setValue(y=y) #set the source as a function on the boundary
175 ahallam 3055 print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26