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

  ViewVC Help
Powered by ViewVC 1.1.26