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

  ViewVC Help
Powered by ViewVC 1.1.26