/[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 4821 - (hide annotations)
Tue Apr 1 04:58:33 2014 UTC (5 years, 6 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 ahallam 3055
2 jfenwick 3981 ##############################################################################
3 ahallam 3055 #
4 jfenwick 4657 # Copyright (c) 2009-2014 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 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 #
15     ##############################################################################
16 ahallam 3055
17 jfenwick 4657 __copyright__="""Copyright (c) 2009-2014 by University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3055 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 caltinay 4087 import matplotlib
32     matplotlib.use('agg') #It's just here for automated testing
33 ahallam 3055 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 jfenwick 3148
41 ahallam 3055 import pylab as pl
42     import matplotlib.cm as cm
43 sshaw 4821 from esys.escript.linearPDEs import LinearPDE, SolverOptions
44 ahallam 3069 from esys.finley import ReadMesh
45 caltinay 3346 from esys.weipa import saveVTK
46 ahallam 3055
47     ########################################################MPI WORLD CHECK
48     if getMPISizeWorld() > 1:
49 sshaw 4576 import sys
50     print("This example will not run in an MPI world.")
51     sys.exit(0)
52 ahallam 3055
53     #################################################ESTABLISHING VARIABLES
54     # where to save output data
55 ahallam 3405 savepath = "data/example09a"
56 ahallam 3389 meshpath = "data/example09m"
57 ahallam 3055 mkDir(savepath)
58     #Geometric and material property related variables.
59 ahallam 3089 step=10.0 # the element size
60 ahallam 3069
61     vel2=1800.; vel1=3000.
62     rho2=2300.; rho1=3100. #density
63 ahallam 3089 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 ahallam 3069
66 ahallam 3405 ####################################################TESTING SWITCH
67 ahallam 3195 testing=True
68     if testing:
69 jfenwick 3892 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 ahallam 3405 tend=0.001
72     #Model Parameters
73     mx=40.
74     my=40.
75     mz=20.
76     outputs=5
77 ahallam 3195 else:
78 ahallam 3405 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 ahallam 3195
85 ahallam 3405 ####################################################TIME RELATED VARIABLES
86     h=0.00005 # time step
87 ahallam 3055 # data recording times
88     rtime=0.0 # first time to record
89 ahallam 3405 rtime_inc=tend/outputs # time increment to record
90 ahallam 3055 #Check to make sure number of time steps is not too large.
91 jfenwick 3892 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
92 ahallam 3055
93 ahallam 3405 ####################################################CREATING THE SOURCE FUNCTION
94 ahallam 3055 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 jfenwick 3892 print('source length',ls)
108 ahallam 3055 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 ahallam 3405 xc=[mx/2,my/2,0]
121 ahallam 3055
122     ####################################################DOMAIN CONSTRUCTION
123 ahallam 3389 domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
124 ahallam 3069 x=domain.getX() # get the locations of the nodes in the domain
125 ahallam 3055
126 ahallam 3069 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 ahallam 3055 ##########################################################ESTABLISH PDE
137     mypde=LinearPDE(domain) # create pde
138     mypde.setSymmetryOn() # turn symmetry on
139     # turn lumping on for more efficient solving
140 sshaw 4821 #mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
141 ahallam 3055 kmat = kronecker(domain) # create the kronecker delta function of the domain
142 ahallam 3069 mypde.setValue(D=rho*kmat) #set the general form value D
143 ahallam 3055
144     ############################################FIRST TIME STEPS AND SOURCE
145     # define small radius around point xc
146 jfenwick 3892 src_rad = 20; print("sourc radius= ",src_rad)
147 ahallam 3055 # set initial values for first two time steps with source terms
148 ahallam 3069 xb=FunctionOnBoundary(domain).getX()
149 ahallam 3405 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 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
153 ahallam 3089
154 ahallam 3405 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
155 ahallam 3055 # initial value of displacement at point source is constant (U0=0.01)
156     # for first two time steps
157 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x)
158 ahallam 3055 u_m1=u
159    
160     ####################################################ITERATION VARIABLES
161     n=0 # iteration counter
162     t=0 # time counter
163     ##############################################################ITERATION
164     while t<tend:
165 sshaw 4576 # get current stress
166 ahallam 3069 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
167 ahallam 3055 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 sshaw 4576 u_p1=u_p1#*abc # apply boundary conditions
171 ahallam 3055 u_m1=u; u=u_p1 # shift values by 1
172     # save current displacement, acceleration and pressure
173     if (t >= rtime):
174 ahallam 3405 saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
175 ahallam 3055 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 ahallam 3405 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
181 caltinay 4005 print("time step %d, t=%s"%(n,t))

  ViewVC Help
Powered by ViewVC 1.1.26