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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3029 - (hide annotations)
Fri May 21 02:01:37 2010 UTC (10 years, 2 months ago) by ahallam
File MIME type: text/x-python
File size: 4282 byte(s)
small updates and work on seismic code
1 ahallam 3025
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     ############################################################FILE HEADER
23 ahallam 3029 # example08a.py
24 ahallam 3025 # Antony Hallam
25     # Seismic Wave Equation Simulation using acceleration solution.
26    
27     #######################################################EXTERNAL MODULES
28     from esys.escript import *
29     from esys.finley import Rectangle
30     import sys
31     import os
32     # smoothing operator
33     from esys.escript.pdetools import Projector, Locator
34     from esys.escript.unitsSI import *
35     import numpy as np
36     import pylab as pl
37     import matplotlib.cm as cm
38     from esys.escript.linearPDEs import LinearPDE
39    
40     #################################################ESTABLISHING VARIABLES
41     # where to save output data
42 ahallam 3029 savepath = "data/example08a"
43 ahallam 3025 mkDir(savepath)
44     #Geometric and material property related variables.
45     mx = 1000. # model lenght
46     my = -1000. # model width
47     ndx = 200 # steps in x direction
48     ndy = 200 # steps in y direction
49     xstep=mx/ndx # calculate the size of delta x
50     ystep=abs(my/ndy) # calculate the size of delta y
51     lam=3.462e9 #lames constant
52     mu=3.462e9 #bulk modulus
53     rho=1154. #density
54     # Time related variables.
55     tend=1.5 # end time
56     h=0.001 # time step
57 ahallam 3029 # data recording times
58     rtime=0.0 # first time to record
59     rtime_inc=tend/20.0 # time increment to record
60 ahallam 3025 #Check to make sure number of time steps is not too large.
61     print "Time step size= ",h, "Expected number of outputs= ",tend/h
62    
63     U0=0.01 # amplitude of point source
64     # will introduce a spherical source at middle left of bottom face
65     xc=[ndx/2-ndx/4,0]
66    
67     ####################################################DOMAIN CONSTRUCTION
68     domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) # create the domain
69     x=domain.getX() # get the locations of the nodes in the domani
70    
71     ##########################################################ESTABLISH PDE
72     mypde=LinearPDE(domain) # create pde
73     mypde.setSymmetryOn() # turn symmetry on
74     # turn lumping on for more efficient solving
75     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
76     kmat = kronecker(domain) # create the kronecker delta function of the domain
77     mypde.setValue(D=kmat*rho) #set the general form value D
78    
79     ############################################FIRST TIME STEPS AND SOURCE
80     # define small radius around point xc
81     src_length = 20; print "src_length = ",src_length
82     # set initial values for first two time steps with source terms
83 ahallam 3029 y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
84 ahallam 3025 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
85 ahallam 3029 y=y*src_dir
86 ahallam 3025 mypde.setValue(y=y) #set the source as a function on the boundary
87     # initial value of displacement at point source is constant (U0=0.01)
88     # for first two time steps
89     u=[0.0,0.0]*whereNegative(x)
90     u_m1=u
91    
92     ####################################################ITERATION VARIABLES
93     n=0 # iteration counter
94     t=0 # time counter
95     ##############################################################ITERATION
96     while t<tend:
97     # get current stress
98     g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))
99     mypde.setValue(X=-stress) # set PDE values
100     accel = mypde.getSolution() #get PDE solution for accelleration
101     u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
102     u_m1=u; u=u_p1 # shift values by 1
103     # save current displacement, acceleration and pressure
104     if (t >= rtime):
105 ahallam 3029 saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\
106     acceleration=length(accel),tensor=stress)
107 ahallam 3025 rtime=rtime+rtime_inc #increment data save time
108     # increment loop values
109     t=t+h; n=n+1
110     print n,"-th time step t ",t

  ViewVC Help
Powered by ViewVC 1.1.26