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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4821 - (hide annotations)
Tue Apr 1 04:58:33 2014 UTC (5 years, 5 months ago) by sshaw
File MIME type: text/x-python
File size: 6573 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 3389
2 jfenwick 3981 ##############################################################################
3 ahallam 3389 #
4 jfenwick 4657 # Copyright (c) 2009-2014 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ahallam 3389 #
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 3389
17 jfenwick 4657 __copyright__="""Copyright (c) 2009-2014 by University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3389 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 3389 from esys.escript import *
34     from esys.finley import Rectangle
35     from esys.weipa import saveVTK
36     import os
37     # smoothing operator
38     from esys.escript.pdetools import Projector, Locator
39     from esys.escript.unitsSI import *
40     import numpy as np
41    
42     import pylab as pl
43     import matplotlib.cm as cm
44 sshaw 4821 from esys.escript.linearPDEs import LinearPDE, SolverOptions
45 ahallam 3389 from esys.finley import ReadMesh
46    
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 3389
53     #################################################ESTABLISHING VARIABLES
54     # where to save output data
55     savepath = "data/example09c"
56     meshpath = "data/example09n"
57     mkDir(savepath)
58     #Geometric and material property related variables.
59     domain=ReadMesh(os.path.join(savepath,'example09n.fly')) # create the domain
60     x=Solution(domain).getX()
61     #parameters layers 1,2,3,4 and fault
62     prho=np.array([2200.,2500.,3200.,4500.,5500.]) #density
63     pvel=np.array([1500.,2200.,3000.,3200.,5000.]) #velocity
64     pmu=pvel**2.*prho/4. #bulk modulus
65     plam=pvel**2.*prho/2. #lames constant
66     nlayers=4
67     width=300.0
68     rho=Scalar(0,Function(domain))
69     vel=Scalar(0,Function(domain))
70     mu=Scalar(0,Function(domain))
71     lam=Scalar(0,Function(domain))
72    
73 jfenwick 3892 print(0.5*np.sqrt(prho/(plam+2*pmu))*0.5)
74 ahallam 3389
75     for i in range(0,nlayers):
76     rho.setTaggedValue('lblock%d'%i,prho[i])
77     rho.setTaggedValue('rblock%d'%i,prho[i])
78     vel.setTaggedValue('lblock%d'%i,pvel[i])
79     vel.setTaggedValue('rblock%d'%i,pvel[i])
80     mu.setTaggedValue('lblock%d'%i,pmu[i])
81     mu.setTaggedValue('rblock%d'%i,pmu[i])
82     lam.setTaggedValue('lblock%d'%i,plam[i])
83     lam.setTaggedValue('rblock%d'%i,plam[i])
84     i=nlayers
85     rho.setTaggedValue('fault',prho[i])
86     vel.setTaggedValue('fault',pvel[i])
87     mu.setTaggedValue('fault',pmu[i])
88     lam.setTaggedValue('fault',plam[i])
89    
90    
91     # Time related variables.
92     testing=False
93     if testing:
94 sshaw 4576 print('The testing end time is currently selected. This severely limits the number of time iterations.')
95     print("Try changing testing to False for more iterations.")
96     tend=0.1
97 ahallam 3389 else:
98 sshaw 4576 tend=0.1 # end time
99 ahallam 3389
100     h=0.00001 # time step
101     # data recording times
102     rtime=0.0 # first time to record
103     rtime_inc=tend/750.0 # time increment to record
104     #Check to make sure number of time steps is not too large.
105 jfenwick 3892 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
106 ahallam 3389
107     U0=0.1 # amplitude of point source
108     ls=500 # length of the source
109     source=np.zeros(ls,'float') # source array
110     decay1=np.zeros(ls,'float') # decay curve one
111     decay2=np.zeros(ls,'float') # decay curve two
112     time=np.zeros(ls,'float') # time values
113     g=np.log(0.01)/ls
114    
115     dfeq=50 #Dominant Frequency
116     a = 2.0 * (np.pi * dfeq)**2.0
117     t0 = 5.0 / (2.0 * np.pi * dfeq)
118     srclength = 5. * t0
119     ls = int(srclength/h)
120 jfenwick 3892 print('source length',ls)
121 ahallam 3389 source=np.zeros(ls,'float') # source array
122     ampmax=0
123     for it in range(0,ls):
124     t = it*h
125     tt = t-t0
126     dum1 = np.exp(-a * tt * tt)
127     source[it] = -2. * a * tt * dum1
128     if (abs(source[it]) > ampmax):
129     ampmax = abs(source[it])
130     time[t]=t*h
131    
132     # will introduce a spherical source at middle left of bottom face
133     xc=[150,0]
134    
135     ##########################################################ESTABLISH PDE
136     mypde=LinearPDE(domain) # create pde
137     mypde.setSymmetryOn() # turn symmetry on
138     # turn lumping on for more efficient solving
139 sshaw 4821 mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
140 ahallam 3389 kmat = kronecker(domain) # create the kronecker delta function of the domain
141     mypde.setValue(D=rho*kmat) #set the general form value D
142    
143     ############################################FIRST TIME STEPS AND SOURCE
144     # define small radius around point xc
145 jfenwick 3892 src_length = 10; print("src_length = ",src_length)
146 ahallam 3389 # set initial values for first two time steps with source terms
147     xb=FunctionOnBoundary(domain).getX()
148     yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)
149     stop=Scalar(0.0,FunctionOnBoundary(domain))
150     stop.setTaggedValue("top",1.0)
151     src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
152    
153     mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
154    
155     # initial value of displacement at point source is constant (U0=0.01)
156     # for first two time steps
157     u=[0.0,0.0]*x
158     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 3389 g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
167     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 3389 u_m1=u; u=u_p1 # shift values by 1
172     # save current displacement, acceleration and pressure
173     if (t >= rtime):
174     saveVTK(os.path.join(savepath,"ex09c.%05d.vtu"%n),displacement=length(u),\
175     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     mypde.setValue(y=source[n]*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