1 |
############################################################################## |
2 |
# |
3 |
# Copyright (c) 2003-2015 by The University of Queensland |
4 |
# http://www.uq.edu.au |
5 |
# |
6 |
# Primary Business: Queensland, Australia |
7 |
# Licensed under the Open Software License version 3.0 |
8 |
# http://www.opensource.org/licenses/osl-3.0.php |
9 |
# |
10 |
# Development until 2012 by Earth Systems Science Computational Center (ESSCC) |
11 |
# Development 2012-2013 by School of Earth Sciences |
12 |
# Development from 2014 by Centre for Geoscience Computing (GeoComp) |
13 |
# |
14 |
############################################################################## |
15 |
from __future__ import division, print_function |
16 |
|
17 |
__copyright__="""Copyright (c) 2003-2015 by The University of Queensland |
18 |
http://www.uq.edu.au |
19 |
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 |
import matplotlib |
25 |
matplotlib.use('agg') #For interactive use, you can comment out this line |
26 |
#It's just here to make testing easier |
27 |
import matplotlib.pyplot as plt |
28 |
from numpy import zeros,ones |
29 |
from esys.escript import * |
30 |
from esys.escript.linearPDEs import LinearPDE, SolverOptions |
31 |
from esys.escript.pdetools import Locator |
32 |
try: |
33 |
from esys.dudley import Brick |
34 |
HAVE_DUDLEY = True |
35 |
except ImportError: |
36 |
HAVE_DUDLEY = False |
37 |
from esys.weipa import saveVTK |
38 |
|
39 |
if not HAVE_DUDLEY: |
40 |
print("Dudley module not available") |
41 |
else: |
42 |
ne=32 # number of cells in x_0 and x_1 directions |
43 |
width=10000. # length in x_0 and x_1 directions |
44 |
lam=3.462e9 |
45 |
mu=3.462e9 |
46 |
rho=1154. |
47 |
tend=10. # to ran a full simulation change tend to 60. |
48 |
alpha=0.7 |
49 |
t0=3. |
50 |
|
51 |
|
52 |
U0=1. # maximum displacement |
53 |
mkDir("data") # create directory data if it does not exist already. |
54 |
|
55 |
def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0): |
56 |
# lists to collect displacement at point source |
57 |
ts, u_pc0,u_pc1,u_pc2=[], [], [], [] |
58 |
x=domain.getX() |
59 |
# ... open new PDE ... |
60 |
mypde=LinearPDE(domain) |
61 |
mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING) |
62 |
kron=kronecker(mypde.getDim()) |
63 |
|
64 |
dunit=numpy.array([1.,0.,0.]) # defines direction of point source |
65 |
|
66 |
mypde.setValue(D=kron*rho, q=whereNegative(length(x-xc)-src_radius)*dunit) |
67 |
# ... set initial values .... |
68 |
n=0 |
69 |
# for first two time steps |
70 |
u=Vector(0.,Solution(domain)) |
71 |
u_last=Vector(0.,Solution(domain)) |
72 |
t=0 |
73 |
|
74 |
# define the location of the point source |
75 |
L=Locator(domain,xc) |
76 |
# find potential at point source |
77 |
u_pc=L.getValue(u) |
78 |
print("u at point charge = %s"%u_pc) |
79 |
ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2]) |
80 |
|
81 |
while t<tend: |
82 |
t+=h |
83 |
# ... get current stress .... |
84 |
g=grad(u) |
85 |
stress=lam*trace(g)*kron+mu*(g+transpose(g)) |
86 |
# ... get new acceleration .... |
87 |
amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2) |
88 |
mypde.setValue(X=-stress, r=dunit*amplitude) |
89 |
a=mypde.getSolution() |
90 |
# ... get new displacement ... |
91 |
u_new=2*u-u_last+h**2*a |
92 |
# ... shift displacements .... |
93 |
u_last=u |
94 |
u=u_new |
95 |
n+=1 |
96 |
print("time step %d, t = %s"%(n,t)) |
97 |
u_pc=L.getValue(u) |
98 |
print("u at point charge = %s"%u_pc) |
99 |
ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2]) |
100 |
|
101 |
# ... save current acceleration in units of gravity and displacements |
102 |
if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81, |
103 |
displacement = length(u), tensor = stress, Ux = u[0] ) |
104 |
return ts, u_pc0,u_pc1,u_pc2 |
105 |
|
106 |
# |
107 |
# create domain: |
108 |
# |
109 |
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne) |
110 |
# |
111 |
# sety time step size: |
112 |
# |
113 |
h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize()) |
114 |
print("time step size = %s"%h) |
115 |
# |
116 |
# spherical source at middle of bottom face |
117 |
# |
118 |
xc=[width/2.,width/2.,0.] |
119 |
# define small radius around point xc |
120 |
src_radius = 0.03*width |
121 |
print("src_radius = %s"%src_radius) |
122 |
# |
123 |
# run it |
124 |
# |
125 |
ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0) |
126 |
# |
127 |
# create a plot: |
128 |
# |
129 |
if getMPIRankWorld() == 0: |
130 |
plt.title("Displacement at Point Source") |
131 |
plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1) |
132 |
plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1) |
133 |
plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1) |
134 |
plt.xlabel('time') |
135 |
plt.ylabel('displacement') |
136 |
plt.legend() |
137 |
plt.savefig('u_pc.png', format='png') |
138 |
# or save displacement |
139 |
u_pc_data=FileWriter('./data/U_pc.out') |
140 |
for i in range(len(ts)) : |
141 |
u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i])) |
142 |
u_pc_data.close() |
143 |
|