/[escript]/trunk/doc/examples/cookbook/example08c.py
ViewVC logotype

Annotation of /trunk/doc/examples/cookbook/example08c.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4576 - (hide annotations)
Mon Dec 9 23:35:30 2013 UTC (5 years, 3 months ago) by sshaw
File MIME type: text/x-python
File size: 9690 byte(s)
python3ified things, replaced mixed whitespace and xrange calls
1 ahallam 3075
2 jfenwick 3981 ##############################################################################
3 ahallam 3075 #
4 jfenwick 4154 # Copyright (c) 2009-2013 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ahallam 3075 #
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     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 ahallam 3075
16 jfenwick 4154 __copyright__="""Copyright (c) 2009-2013 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 ahallam 3075 Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21     __url__="https://launchpad.net/escript-finley"
22    
23     """
24     Author: Antony Hallam antony.hallam@uqconnect.edu.au
25     """
26    
27     ############################################################FILE HEADER
28     # example08c.py
29     # Create either a 2D syncline or anticline model using pycad meshing
30     # tools. Wave equation solution.
31    
32     #######################################################EXTERNAL MODULES
33     import matplotlib
34     matplotlib.use('agg') #It's just here for automated testing
35     from esys.pycad import * #domain constructor
36     from esys.pycad.gmsh import Design #Finite Element meshing package
37     from esys.finley import MakeDomain #Converter for escript
38     import os #file path tool
39     from math import * # math package
40     from esys.escript import *
41     from esys.escript.unitsSI import *
42     from esys.escript.linearPDEs import LinearPDE
43     from esys.escript.pdetools import Projector
44 caltinay 3346 from esys.weipa import saveVTK
45 ahallam 3075 from cblib import toRegGrid, subsample
46 jfenwick 3148 import matplotlib
47     matplotlib.use('agg') #It's just here for automated testing
48    
49 ahallam 3075 import pylab as pl #Plotting package
50     import numpy as np
51    
52     ########################################################MPI WORLD CHECK
53     if getMPISizeWorld() > 1:
54 sshaw 4576 import sys
55     print("This example will not run in an MPI world.")
56     sys.exit(0)
57 ahallam 3075
58     #################################################ESTABLISHING VARIABLES
59     #set modal to 1 for a syncline or -1 for an anticline structural
60 sshaw 4576 #configuration
61 ahallam 3075 modal=-1
62    
63     # the folder to put our outputs in, leave blank "" for script path -
64     # note this folder path must exist to work
65     save_path= os.path.join("data","example08c")
66     mkDir(save_path)
67    
68     ################################################ESTABLISHING PARAMETERS
69     #Model Parameters
70     width=1000.0*m #width of model
71     depth=1000.0*m #depth of model
72     dx=5
73     xstep=dx # calculate the size of delta x
74     ystep=dx # calculate the size of delta y
75    
76     sspl=51 #number of discrete points in spline
77     dsp=width/(sspl-1) #dx of spline steps for width
78     dep_sp=500.0*m #avg depth of spline
79     h_sp=250.*m #heigh of spline
80     orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
81    
82     vel2=1800.; vel1=3000.
83     rho2=2300.; rho1=3100. #density
84     mu2=(vel2**2.)*rho2/8.; mu1=(vel1**2.)*rho1/8. #bulk modulus
85     lam2=mu2*6.; lam1=mu1*6. #lames constant
86    
87    
88     # Time related variables.
89 ahallam 3195 testing=True
90     if testing:
91 sshaw 4576 print('The testing end time is currently selected. This severely limits the number of time iterations.')
92     print("Try changing testing to False for more iterations.")
93     tend=0.001
94 ahallam 3195 else:
95 sshaw 4576 tend=0.5 # end time
96 ahallam 3195
97 ahallam 3075 h=0.0001 # time step
98     # data recording times
99     rtime=0.0 # first time to record
100     rtime_inc=tend/50.0 # time increment to record
101     # will introduce a spherical source at middle left of bottom face
102     xc=[width/2,0]
103     #Check to make sure number of time steps is not too large.
104 jfenwick 3892 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
105 ahallam 3075
106     U0=0.1 # amplitude of point source
107     ls=500 # length of the source
108     source=np.zeros(ls,'float') # source array
109     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     dfeq=50 #Dominant Frequency
115     a = 2.0 * (np.pi * dfeq)**2.0
116     t0 = 5.0 / (2.0 * np.pi * dfeq)
117     srclength = 5. * t0
118     ls = int(srclength/h)
119 jfenwick 3892 print('source length',ls)
120 ahallam 3075 source=np.zeros(ls,'float') # source array
121     ampmax=0
122     for it in range(0,ls):
123     t = it*h
124     tt = t-t0
125     dum1 = np.exp(-a * tt * tt)
126     source[it] = -2. * a * tt * dum1
127     if (abs(source[it]) > ampmax):
128     ampmax = abs(source[it])
129     time[t]=t*h
130    
131     ####################################################DOMAIN CONSTRUCTION
132     # Domain Corners
133     p0=Point(0.0, 0.0, 0.0)
134     p1=Point(0.0, depth, 0.0)
135     p2=Point(width, depth, 0.0)
136     p3=Point(width, 0.0, 0.0)
137    
138     # Generate Material Boundary
139     x=[ Point(i*dsp\
140     ,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
141     for i in range(0,sspl)\
142     ]
143     mysp = Spline(*tuple(x))
144     # Start and end of material boundary.
145     x1=mysp.getStartPoint()
146     x2=mysp.getEndPoint()
147 sshaw 4576
148 ahallam 3075 # Create TOP BLOCK
149     # lines
150     tbl1=Line(p0,x1)
151     tbl2=mysp
152     tbl3=Line(x2,p3)
153     l30=Line(p3, p0)
154     # curve
155     tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
156     # surface
157     tblock = PlaneSurface(tblockloop)
158     # Create BOTTOM BLOCK
159     # lines
160     bbl1=Line(x1,p1)
161     bbl3=Line(p2,x2)
162     bbl4=-mysp
163     l12=Line(p1, p2)
164     # curve
165     bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
166    
167     # surface
168     bblock = PlaneSurface(bblockloop)
169    
170     #clockwise check as splines must be set as polygons in the point order
171     #they were created. Otherwise get a line across plot.
172     bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
173    
174     ################################################CREATE MESH FOR ESCRIPT
175     # Create a Design which can make the mesh
176     d=Design(dim=2, element_size=dx, order=2)
177     # Add the subdomains and flux boundaries.
178     d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
179     PropertySet("linetop",l30))
180     # Create the geometry, mesh and Escript domain
181     d.setScriptFileName(os.path.join(save_path,"example08c.geo"))
182     d.setMeshFileName(os.path.join(save_path,"example08c.msh"))
183     domain=MakeDomain(d, optimizeLabeling=True)
184     x=domain.getX()
185 jfenwick 3892 print("Domain has been generated ...")
186 ahallam 3075
187     lam=Scalar(0,Function(domain))
188     lam.setTaggedValue("top",lam1)
189     lam.setTaggedValue("bottom",lam2)
190     mu=Scalar(0,Function(domain))
191     mu.setTaggedValue("top",mu1)
192     mu.setTaggedValue("bottom",mu2)
193     rho=Scalar(0,Function(domain))
194     rho.setTaggedValue("top",rho1)
195     rho.setTaggedValue("bottom",rho2)
196    
197     ##########################################################ESTABLISH PDE
198     mypde=LinearPDE(domain) # create pde
199     mypde.setSymmetryOn() # turn symmetry on
200     # turn lumping on for more efficient solving
201     #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
202     kmat = kronecker(domain) # create the kronecker delta function of the domain
203     mypde.setValue(D=rho*kmat) #set the general form value D
204    
205     ##########################################################ESTABLISH ABC
206     # Define where the boundary decay will be applied.
207 ahallam 3089 bn=20.
208 ahallam 3075 bleft=xstep*bn; bright=width-(xstep*bn); bbot=depth-(ystep*bn)
209     # btop=ystep*bn # don't apply to force boundary!!!
210    
211     # locate these points in the domain
212     left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot
213    
214     tgamma=0.85 # decay value for exponential function
215     def calc_gamma(G,npts):
216     func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))
217     return func
218    
219     gleft = calc_gamma(tgamma,bleft)
220     gright = calc_gamma(tgamma,bleft)
221     gbottom= calc_gamma(tgamma,ystep*bn)
222    
223 jfenwick 3892 print('gamma', gleft,gright,gbottom)
224 ahallam 3075
225     # calculate decay functions
226     def abc_bfunc(gamma,loc,x,G):
227     func=exp(-1.*(gamma*abs(loc-x))**2.)
228     return func
229    
230     fleft=abc_bfunc(gleft,bleft,x[0],tgamma)
231     fright=abc_bfunc(gright,bright,x[0],tgamma)
232     fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)
233     # apply these functions only where relevant
234     abcleft=fleft*whereNegative(left)
235     abcright=fright*wherePositive(right)
236     abcbottom=fbottom*wherePositive(bottom)
237     # make sure the inside of the abc is value 1
238     abcleft=abcleft+whereZero(abcleft)
239     abcright=abcright+whereZero(abcright)
240     abcbottom=abcbottom+whereZero(abcbottom)
241     # multiply the conditions together to get a smooth result
242     abc=abcleft*abcright*abcbottom
243    
244     ############################################FIRST TIME STEPS AND SOURCE
245     # define small radius around point xc
246 jfenwick 3892 src_length = 40; print("src_length = ",src_length)
247 ahallam 3075 # set initial values for first two time steps with source terms
248     xb=FunctionOnBoundary(domain).getX()
249     y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
250     src_dir=numpy.array([0.,1.]) # defines direction of point source as down
251     y=y*src_dir
252     mypde.setValue(y=y) #set the source as a function on the boundary
253     # initial value of displacement at point source is constant (U0=0.01)
254     # for first two time steps
255     u=[0.0,0.0]*wherePositive(x)
256     u_m1=u
257    
258     ####################################################ITERATION VARIABLES
259     n=0 # iteration counter
260     t=0 # time counter
261     ##############################################################ITERATION
262     while t<tend:
263 sshaw 4576 # get current stress
264 ahallam 3075 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))
265     mypde.setValue(X=-stress*abc) # set PDE values
266     accel = mypde.getSolution() #get PDE solution for accelleration
267     u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
268 sshaw 4576 u_p1=u_p1*abc # apply boundary conditions
269 ahallam 3075 u_m1=u; u=u_p1 # shift values by 1
270     # save current displacement, acceleration and pressure
271     if (t >= rtime):
272     saveVTK(os.path.join(save_path,"ex08c.%05d.vtu"%n),\
273     vector_displacement=u,displacement=length(u),\
274     vector_acceleration=accel,acceleration=length(accel),\
275     tensor=stress)
276     rtime=rtime+rtime_inc #increment data save time
277     # increment loop values
278     t=t+h; n=n+1
279     if (n < ls):
280     y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
281     y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the boundary
282 caltinay 4005 print("time step %d, t=%s"%(n,t))

  ViewVC Help
Powered by ViewVC 1.1.26