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

  ViewVC Help
Powered by ViewVC 1.1.26