/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (19 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 7198 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 sshaw 5288 from __future__ import division, print_function
2 jfenwick 3981 ##############################################################################
3 ahallam 3055 #
4 jfenwick 6651 # Copyright (c) 2009-2018 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ahallam 3055 #
7     # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 ahallam 3055 #
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 6651 __copyright__="""Copyright (c) 2009-2018 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3055 Primary Business: Queensland, Australia"""
20 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
21     http://www.apache.org/licenses/LICENSE-2.0"""
22 ahallam 3055 __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     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 sshaw 4821 from esys.escript.linearPDEs import LinearPDE, SolverOptions
43 caltinay 3346 from esys.weipa import saveVTK
44 sshaw 5288 try:
45     # This imports the rectangle domain function
46     from esys.finley import Rectangle, ReadMesh
47     HAVE_FINLEY = True
48     except ImportError:
49     print("Finley module not available")
50     HAVE_FINLEY = False
51 ahallam 3055 ########################################################MPI WORLD CHECK
52     if getMPISizeWorld() > 1:
53 sshaw 4576 import sys
54     print("This example will not run in an MPI world.")
55     sys.exit(0)
56 ahallam 3055
57 sshaw 5288 if HAVE_FINLEY:
58     #################################################ESTABLISHING VARIABLES
59     # where to save output data
60     savepath = "data/example09a"
61     meshpath = "data/example09m"
62     mkDir(savepath)
63     #Geometric and material property related variables.
64     step=10.0 # the element size
65 ahallam 3069
66 sshaw 5288 vel2=1800.; vel1=3000.
67     rho2=2300.; rho1=3100. #density
68     mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus
69     lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant
70 ahallam 3069
71 sshaw 5288 ####################################################TESTING SWITCH
72     testing=True
73     if testing:
74     print('The testing end time is currently selected. This severely limits the number of time iterations.')
75     print("Try changing testing to False for more iterations.")
76     tend=0.001
77     #Model Parameters
78     mx=40.
79     my=40.
80     mz=20.
81     outputs=5
82     else:
83     tend=0.1 # end time
84     #Model Parameters
85     mx=100.0 #x width of model
86     my=100.0 #y width of model
87     mz=50.0 #depth of model
88     outputs=200
89 ahallam 3195
90 sshaw 5288 ####################################################TIME RELATED VARIABLES
91     h=0.00005 # time step
92     # data recording times
93     rtime=0.0 # first time to record
94     rtime_inc=tend/outputs # time increment to record
95     #Check to make sure number of time steps is not too large.
96     print("Time step size= ",h, "Expected number of outputs= ",tend/h)
97 ahallam 3055
98 sshaw 5288 ####################################################CREATING THE SOURCE FUNCTION
99     U0=0.1 # amplitude of point source
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 sshaw 5621
105 sshaw 5288 ls = int(srclength/h)
106     print('source length',ls)
107 sshaw 5621
108 sshaw 5288 source=np.zeros(ls,'float') # source array
109 sshaw 5621 decay1=np.zeros(ls,'float') # decay curve one
110     decay2=np.zeros(ls,'float') # decay curve two
111     time=np.zeros(ls,'float') # time values
112     g=np.log(0.01)/ls
113    
114 sshaw 5288 ampmax=0
115     for it in range(0,ls):
116     t = it*h
117     tt = t-t0
118     dum1 = np.exp(-a * tt * tt)
119     source[it] = -2. * a * tt * dum1
120     if (abs(source[it]) > ampmax):
121     ampmax = abs(source[it])
122 sshaw 5621 time[it]=t*h
123 ahallam 3055
124 sshaw 5288 # will introduce a spherical source at middle left of bottom face
125     xc=[mx/2,my/2,0]
126 ahallam 3055
127 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
128     domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
129     x=domain.getX() # get the locations of the nodes in the domain
130 ahallam 3055
131 sshaw 5288 lam=Scalar(0,Function(domain))
132     lam.setTaggedValue("vintfa",lam1)
133     lam.setTaggedValue("vintfb",lam2)
134     mu=Scalar(0,Function(domain))
135     mu.setTaggedValue("vintfa",mu1)
136     mu.setTaggedValue("vintfb",mu2)
137     rho=Scalar(0,Function(domain))
138     rho.setTaggedValue("vintfa",rho1)
139     rho.setTaggedValue("vintfb",rho2)
140 ahallam 3069
141 sshaw 5288 ##########################################################ESTABLISH PDE
142     mypde=LinearPDE(domain) # create pde
143     mypde.setSymmetryOn() # turn symmetry on
144     # turn lumping on for more efficient solving
145     #mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
146     kmat = kronecker(domain) # create the kronecker delta function of the domain
147     mypde.setValue(D=rho*kmat) #set the general form value D
148 ahallam 3055
149 sshaw 5288 ############################################FIRST TIME STEPS AND SOURCE
150     # define small radius around point xc
151     src_rad = 20; print("sourc radius= ",src_rad)
152     # set initial values for first two time steps with source terms
153     xb=FunctionOnBoundary(domain).getX()
154     yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
155     stop=Scalar(0.0,FunctionOnBoundary(domain))
156     stop.setTaggedValue("stop",1.0)
157 jfenwick 6095 src_dir=np.array([0.,0.,1.0]) # defines direction of point source as down
158 ahallam 3089
159 sshaw 5288 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
160     # initial value of displacement at point source is constant (U0=0.01)
161     # for first two time steps
162     u=[0.0,0.0,0.0]*wherePositive(x)
163     u_m1=u
164 ahallam 3055
165 sshaw 5288 ####################################################ITERATION VARIABLES
166     n=0 # iteration counter
167     t=0 # time counter
168     ##############################################################ITERATION
169     while t<tend:
170     # get current stress
171     g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
172     mypde.setValue(X=-stress) # set PDE values
173     accel = mypde.getSolution() #get PDE solution for accelleration
174     u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
175     u_p1=u_p1#*abc # apply boundary conditions
176     u_m1=u; u=u_p1 # shift values by 1
177     # save current displacement, acceleration and pressure
178     if (t >= rtime):
179     saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
180     acceleration=length(accel),tensor=stress)
181     rtime=rtime+rtime_inc #increment data save time
182     # increment loop values
183     t=t+h; n=n+1
184     if (n < ls):
185     mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
186     print("time step %d, t=%s"%(n,t))

  ViewVC Help
Powered by ViewVC 1.1.26