/[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 3346 - (hide annotations)
Fri Nov 12 01:19:02 2010 UTC (8 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 9562 byte(s)
Replaced usage of esys.escript.util.saveVTK by weipa.saveVTK in all python
scripts.

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

  ViewVC Help
Powered by ViewVC 1.1.26