1 |
|
2 |
######################################################## |
3 |
# |
4 |
# Copyright (c) 2009 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 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 |
# You can shorten the execution time by reducing variable tend from 60 to 0.5 |
23 |
|
24 |
# Importing all the necessary modules required. |
25 |
from esys.escript import * |
26 |
from esys.finley import Rectangle |
27 |
import sys |
28 |
import os |
29 |
from cblib1 import needdirs, wavesolver2d |
30 |
|
31 |
# Establish a save path. |
32 |
savepath = "data/wavesolver2d001nwtest" |
33 |
needdirs([savepath]) |
34 |
|
35 |
|
36 |
#Geometric and material property related variables. |
37 |
mx = 1000 # model lenght |
38 |
my = 200 # model width |
39 |
ndx = 200 # steps in x direction |
40 |
ndy = 40 # steps in y direction |
41 |
lam=3.462e9 #lames constant |
42 |
mu=3.462e9 #bulk modulus |
43 |
rho=1154. #density |
44 |
# Time related variables. |
45 |
tend=0.5 #end time |
46 |
#calculating )the timestep |
47 |
h=(1./5.)*sqrt(rho/(lam+2*mu))*(mx/ndx) |
48 |
#Check to make sure number of time steps is not too large. |
49 |
print "Time step size= ",h, "Expected number of outputs= ",tend/h |
50 |
|
51 |
#uncomment the following lines to give the user a chance to stop |
52 |
#proceeder = raw_input("Is this ok?(y/n)") |
53 |
#Exit if user thinks too many outputs. |
54 |
#if proceeder == "n": |
55 |
# sys.exit() |
56 |
|
57 |
U0=0.01 # amplitude of point source |
58 |
# spherical source at middle of bottom face |
59 |
|
60 |
xc=[300,200] |
61 |
|
62 |
mydomain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) |
63 |
wavesolver2d(mydomain,h,tend,lam,mu,rho,U0,xc,savepath) |
64 |
|