20 |
\begin{eqnarray} \label{WAVE natural} |
\begin{eqnarray} \label{WAVE natural} |
21 |
\sigma\hackscore{ij}n\hackscore{j}=0 |
\sigma\hackscore{ij}n\hackscore{j}=0 |
22 |
\end{eqnarray} |
\end{eqnarray} |
23 |
for all time $t>0$. At initial time $t=0$ the displacement |
for all time $t>0$. |
24 |
$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are give: |
|
25 |
|
At initial time $t=0$ the displacement |
26 |
|
$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: |
27 |
\begin{eqnarray} \label{WAVE initial} |
\begin{eqnarray} \label{WAVE initial} |
28 |
u\hackscore{i}(0,x)=0 & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; |
u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array} \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; |
29 |
\end{eqnarray} |
\end{eqnarray} |
30 |
for all $x$ in the domain. |
for all $x$ in the domain. |
31 |
|
|
32 |
The points on the left face of the domain, ie. points $x=(x\hackscore{i})$ with $x\hackscore{0}=0$, |
Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we |
33 |
are moved into the $x\hackscore{2}$ direction by $u^{max}$. This is expressed by the constraint |
set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point |
34 |
\begin{eqnarray} \label{WAVE impact} |
$x\hackscore C$. |
|
u\hackscore{2}(x,t)=s(t) \mbox{ for } x=(x\hackscore{i}) with x\hackscore{0}=0 |
|
|
\end{eqnarray} |
|
|
where |
|
|
\begin{eqnarray} \label{WAVE impact s} |
|
|
s(t)=u^{max} (1-e^{-(\tau^{-1}t)^3}) |
|
|
\end{eqnarray} |
|
|
$\tau$ measures the time needed to reach the maximum displacement $u^{max}$ ( |
|
|
actually $\tau (ln 2)^{\frac{1}{3}}) \approx 0.9 \tau$ is the time until $\frac{u^{max}}{2}$ is reached). |
|
|
The smaller $\tau$ the faster $u^{max}$ is reached. |
|
35 |
|
|
36 |
We use an explicit time-integration scheme to calculate the displacement field $u$ at |
We use an explicit time-integration scheme to calculate the displacement field $u$ at |
37 |
certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$ |
certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$ |
39 |
\begin{eqnarray} \label{WAVE dyn 2} |
\begin{eqnarray} \label{WAVE dyn 2} |
40 |
u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ |
u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ |
41 |
\end{eqnarray} |
\end{eqnarray} |
42 |
for all $n=1,2,\ldots$. It is designed to solve a system of equations of the form |
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form |
43 |
\begin{eqnarray} \label{WAVE dyn 2b} |
\begin{eqnarray} \label{WAVE dyn 2b} |
44 |
u\hackscore{,tt}=G(u) |
u\hackscore{,tt}=G(u) |
45 |
\end{eqnarray} |
\end{eqnarray} |
55 |
\end{eqnarray} |
\end{eqnarray} |
56 |
derived from \eqn{WAVE natural} where |
derived from \eqn{WAVE natural} where |
57 |
\begin{eqnarray} \label{WAVE dyn 3a} |
\begin{eqnarray} \label{WAVE dyn 3a} |
58 |
\sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}) |
\sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}). |
|
\end{eqnarray}. |
|
|
Additionally $a^{(n)}$ has to fullfill the constraint |
|
|
\begin{eqnarray}\label{WAVE dyn 4} |
|
|
a^{(n)}\hackscore{2}(x)= s\hackscore{,tt}(t^{(n)})\mbox{ where } s\hackscore{,tt}(t)=\frac{u^{max}}{\tau^2} |
|
|
(6\tau^{-1}t- 9(\tau^{-1}t)^4) e^{-(\tau^{-1}t)^3} |
|
59 |
\end{eqnarray} |
\end{eqnarray} |
60 |
derived from \eqn{WAVE impact}. |
Now we have converted our problem for displacement, $u^{(n)}$, into a problem for |
61 |
|
acceleration, $a^(n)$, which now depends |
62 |
|
on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. |
63 |
|
|
64 |
In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will |
In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will |
65 |
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through |
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through |
66 |
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are |
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The forms which are relevant here are |
67 |
\begin{equation}\label{WAVE dyn 100} |
\begin{equation}\label{WAVE dyn 100} |
68 |
D\hackscore{ij} u\hackscore{j} = - X\hackscore{ij,j}\; . |
D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . |
69 |
\end{equation} |
\end{equation} |
70 |
Implicitly, the natural boundary condition |
Implicitly, the natural boundary condition |
71 |
\begin{equation}\label{WAVE dyn 101} |
\begin{equation}\label{WAVE dyn 101} |
75 |
\begin{equation}\label{WAVE dyn 101} |
\begin{equation}\label{WAVE dyn 101} |
76 |
u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0 |
u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0 |
77 |
\end{equation} |
\end{equation} |
78 |
where $r$ and $q$ are given functions. |
where $r$ and $q$ are given functions. However in this problem we don't need to define |
79 |
|
any constraints for our problem. |
80 |
|
|
81 |
With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$, $r$ and $q$: |
With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: |
82 |
\begin{equation}\label{WAVE dyn 6} |
\begin{equation}\label{WAVE dyn 6} |
83 |
\begin{array}{ll} |
\begin{array}{ll} |
84 |
D\hackscore{ij}=\rho \delta\hackscore{ij}& |
D\hackscore{ij}=\rho \delta\hackscore{ij}& |
85 |
X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ |
X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ |
|
r\hackscore{i}=\delta\hackscore{i2} \cdot s\hackscore{,tt}(t^{(n)}) & q\hackscore{i}=\delta\hackscore{i2} \cdot x\hackscore{0}.\texttt{whereZero()} |
|
86 |
\end{array} |
\end{array} |
87 |
\end{equation} |
\end{equation} |
88 |
Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero |
|
89 |
(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. |
|
90 |
|
%Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero |
91 |
|
%(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. |
92 |
|
|
93 |
In the following script defines a the function \function{wavePropagation} which |
The following script defines a the function \function{wavePropagation} which |
94 |
implements the Verlet scheme to solve our wave propagation problem. |
implements the Verlet scheme to solve our wave propagation problem. |
95 |
The argument \var{domain} which is a \Domain class object |
The argument \var{domain} which is a \Domain class object |
96 |
defines the domain of the problem. \var{h} and \var{tend} is the step size |
defines the domain of the problem. \var{h} and \var{tend} are the step size |
97 |
and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and |
and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and |
98 |
\var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time: |
\var{rho} are material properties. |
99 |
\begin{python} |
\begin{python} |
100 |
from esys.linearPDEs import LinearPDE |
from esys.linearPDEs import LinearPDE |
101 |
from numarray import identity |
from numarray import identity,zeros,ones |
102 |
from esys.escript import * |
from esys.escript import * |
103 |
def wavePropagation(domain,h,tend,lam,mu,rho,s_tt): |
from esys.escript.pdetools import Locator |
104 |
|
def wavePropagation(domain,h,tend,lam,mu,rho,U0): |
105 |
x=domain.getX() |
x=domain.getX() |
106 |
# ... open new PDE ... |
# ... open new PDE ... |
107 |
mypde=LinearPDE(domain) |
mypde=LinearPDE(domain) |
108 |
mypde.setLumpingOn() |
mypde.setSolverMethod(LinearPDE.LUMPING) |
109 |
kronecker=identity(mypde.getDim()) |
kronecker=identity(mypde.getDim()) |
110 |
mypde.setValue(D=kronecker*rho, \ |
# spherical source at middle of bottom face |
111 |
q=x[0].whereZero()*kronecker[1,:]) |
xc=[width/2.,width/2.,0.] |
112 |
|
# define small radius around point xc |
113 |
|
# Lsup(x) returns the maximum value of the argument x |
114 |
|
src_radius = 0.1*Lsup(domain.getSize()) |
115 |
|
print "src_radius = ",src_radius |
116 |
|
dunit=numarray.array([1.,0.,0.]) # defines direction of point source |
117 |
|
mypde.setValue(D=kronecker*rho) |
118 |
# ... set initial values .... |
# ... set initial values .... |
119 |
n=0 |
n=0 |
120 |
u=Vector(0,ContinuousFunction(domain)) |
# initial value of displacement at point source is constant (U0=0.01) |
121 |
u_last=Vector(0,ContinuousFunction(domain)) |
# for first two time steps |
122 |
|
u=U0*whereNegative(length(x-xc)-src_radius)*dunit |
123 |
|
u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit |
124 |
t=0 |
t=0 |
125 |
|
# define the location of the point source |
126 |
|
L=Locator(domain,xc) |
127 |
|
# find potential at point source |
128 |
|
u_pc=L.getValue(u) |
129 |
|
print "u at point charge=",u_pc |
130 |
|
u_pc_x = u_pc[0] |
131 |
|
u_pc_y = u_pc[1] |
132 |
|
u_pc_z = u_pc[2] |
133 |
|
# open file to save displacement at point source |
134 |
|
u_pc_data=open('./data/U_pc.out','w') |
135 |
|
u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) |
136 |
while t<tend: |
while t<tend: |
137 |
# ... get current stress .... |
# ... get current stress .... |
138 |
g=grad(u) |
g=grad(u) |
139 |
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) |
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) |
140 |
# ... get new acceleration .... |
# ... get new acceleration .... |
141 |
mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:]) |
mypde.setValue(X=-stress) |
142 |
a=mypde.getSolution() |
a=mypde.getSolution() |
143 |
# ... get new displacement ... |
# ... get new displacement ... |
144 |
u_new=2*u-u_last+h**2*a |
u_new=2*u-u_last+h**2*a |
148 |
t+=h |
t+=h |
149 |
n+=1 |
n+=1 |
150 |
print n,"-th time step t ",t |
print n,"-th time step t ",t |
151 |
print "a=",inf(a),sup(a) |
L=Locator(domain,xc) |
152 |
print "u=",inf(u),sup(u) |
u_pc=L.getValue(u) |
153 |
# ... save current acceleration in units of gravity |
print "u at point charge=",u_pc |
154 |
if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10)) |
u_pc_x=u_pc[0] |
155 |
|
u_pc_y=u_pc[1] |
156 |
|
u_pc_z=u_pc[2] |
157 |
|
# save displacements at point source to file for t > 0 |
158 |
|
u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) |
159 |
|
# ... save current acceleration in units of gravity and displacements |
160 |
|
if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81, |
161 |
|
displacement = length(u), Ux = u[0] ) |
162 |
|
u_pc_data.close() |
163 |
\end{python} |
\end{python} |
164 |
Notice that |
Notice that |
165 |
all coefficients of the PDE which are independent from time $t$ are set outside the \code{while} |
all coefficients of the PDE which are independent from time $t$ are set outside the \code{while} |
166 |
loop. This allows the \LinearPDE class to ruse information as much as possible |
loop. This allows the \LinearPDE class to reuse information as much as possible |
167 |
when iterating over time. |
when iterating over time. |
168 |
|
|
169 |
We have seen most of the functions in previous examples but there are some new functions here: |
We have seen most of the functions in previous examples but there are some new functions here: |
174 |
and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie. |
and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie. |
175 |
$\var{transpose(g)[i,j]}=\var{g[j,i]}$). |
$\var{transpose(g)[i,j]}=\var{g[j,i]}$). |
176 |
|
|
177 |
The initialization of \var{u} and \var{u_last} needs some explanation. The statement |
We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small |
178 |
\begin{python} |
sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$. |
179 |
u=Vector(0.,ContinuousFunction(domain)) |
These satisfy the initial conditions defined in \eqn{WAVE initial}. |
180 |
\end{python} |
|
181 |
declares \var{u} as a vector-valued function of its coordinates |
The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}. |
182 |
in the domain (i.e. with two or three components in the case of |
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ |
183 |
a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal |
respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of |
184 |
to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} |
the displacement vector $u$ at the point source. These can be |
185 |
specifies the `smoothness` of the function. In our case we want the displacement field to be |
plotted easily using any plotting package. In gnuplot the command: |
186 |
a continuous function over the \Domain \var{domain}. On other option would be |
\code{plot 'U_pc.out' u 1:2 title 'U_x' w l lw 2, 'U_pc.out' u 1:3 title 'U_y' w l lw 2, |
187 |
\var{Function(domain)} in which case continuity is not garanteed. In fact the |
'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}. |
188 |
argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while |
%\begin{python} |
189 |
returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. |
%u=Vector(0.,ContinuousFunction(domain)) |
190 |
|
%\end{python} |
191 |
|
%declares \var{u} as a vector-valued function of its coordinates |
192 |
|
%in the domain (i.e. with two or three components in the case of |
193 |
|
%a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal |
194 |
|
%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} |
195 |
|
%specifies the `smoothness` of the function. In our case we want the displacement field to be |
196 |
|
%a continuous function over the \Domain \var{domain}. On other option would be |
197 |
|
%\var{Function(domain)} in which case continuity is not garanteed. In fact the |
198 |
|
%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while |
199 |
|
%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. |
200 |
|
|
201 |
The statement \code{myPDE.setLumpingOn()} |
The statement \code{myPDE.setLumpingOn()} |
202 |
switches on the usage of an aggressive approximation of the PDE operator, in this case |
switches on the usage of an aggressive approximation of the PDE operator, in this case |
203 |
formed by the coefficient \var{D} (actually the discrete |
formed by the coefficient \var{D} (actually the discrete |
204 |
version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of |
version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of |
205 |
an additional error, very fast |
an additional error, very fast |
206 |
solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var |
solving of the PDE when the solution is requested. In the general case, the presence of \var{A}, \var{B} or \var |
207 |
{C} \code{setLumpingOn} will produce wrong results. |
{C} \code{setLumpingOn} will produce wrong results. |
208 |
|
|
|
It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent |
|
|
with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time |
|
|
direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the |
|
|
constraint defined by \eqn{WAVE impact} fulfills |
|
|
the initial conditions in \eqn{WAVE initial}. |
|
|
|
|
209 |
|
|
210 |
One of the big advantages of the Verlet scheme is the fact that the problem to be solved |
One of the big advantages of the Verlet scheme is the fact that the problem to be solved |
211 |
in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping). |
in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping). |
212 |
The problem becomes so simple because we use the stress from the last time step rather then the stress which |
The problem becomes so simple because we use the stress from the last time step rather then the stress which is |
213 |
actually is present at the current time step. Schemes using this approach are called an explicit time integration |
actually present at the current time step. Schemes using this approach are called an explicit time integration |
214 |
scheme \index{explicit scheme} \index{time integration!explicit}. The |
scheme \index{explicit scheme} \index{time integration!explicit}. The |
215 |
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is |
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is |
216 |
an example of an implicit schemes |
an example of an implicit schemes |
217 |
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of |
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of |
218 |
each variable at a particular time rather then values from previous time steps. This will lead to problem which is |
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is |
219 |
more expensive to solve, in particular for non-linear problems. |
more expensive to solve, in particular for non-linear problems. |
220 |
Although |
Although |
221 |
explicit time integration schemes are cheep to finalize a single time step, they need a significant smaller time |
explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time |
222 |
step then implicit schemes and can suffer from stability problems. Therefor they need a |
steps then implicit schemes and can suffer from stability problems. Therefore they need a |
223 |
very careful selection of the time step size $h$. |
very careful selection of the time step size $h$. |
224 |
|
|
225 |
An easy, heuristic way of choosing an appropriate |
An easy, heuristic way of choosing an appropriate |
230 |
\begin{eqnarray}\label{WAVE dyn 66} |
\begin{eqnarray}\label{WAVE dyn 66} |
231 |
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x |
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x |
232 |
\end{eqnarray} |
\end{eqnarray} |
233 |
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristic of |
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of |
234 |
the formula. |
the formula. |
235 |
|
|
236 |
The following script uses the \code{wavePropagtion} function to solve the |
The following script uses the \code{wavePropagation} function to solve the |
237 |
wave equation over a block in crust of the earth. The depth of the block is |
wave equation for a point source located at the bottom face of a block. The width of the block in |
238 |
$10km$ and the width in each direction is $100km$. The depth of the block is alligned |
each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}). |
239 |
with the $x\hackscore{0}$-direction. \var{ne} gives the number of elements in $x\hackscore{0}$-direction. The number |
\code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions. |
240 |
of elements in $x\hackscore{1}$- and $x\hackscore{2}$-direction is chosen such that the elements |
The depth of the block is aligned with the $x\hackscore{2}$-direction. |
241 |
become cubic. The impact is $2m$ acting within a time of about $10sec$: |
The depth (\code{l2}) of the block in |
242 |
|
the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the |
243 |
|
of the depth is chosen such that the elements |
244 |
|
become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the |
245 |
|
computation would be faster. \code{Brick($n\hackscore 0, n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$)} is an \finley function which creates a rectangular mesh |
246 |
|
with $n\hackscore 0 \times n\hackscore 1 \times n\hackscore2$ elements over the brick $[0,l\hackscore 0] \times [0,l\hackscore 1] \times [0,l\hackscore 2]$. |
247 |
\begin{python} |
\begin{python} |
248 |
ne=10 # number of cells in x_0-direction |
from esys.finley import Brick |
249 |
depth=10000. # length in x_0-direction |
ne=32 # number of cells in x_0 and x_1 directions |
250 |
width=100000. # length in x_1 and x_2 direction |
width=10000. # length in x_0 and x_1 directions |
251 |
lam=3.462e9 |
lam=3.462e9 |
252 |
mu=3.462e9 |
mu=3.462e9 |
253 |
rho=1154. |
rho=1154. |
|
tau=10. |
|
|
umax=2. |
|
254 |
tend=60 |
tend=60 |
255 |
h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne) |
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne) |
|
def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3) |
|
|
|
|
256 |
print "time step size = ",h |
print "time step size = ",h |
257 |
mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width) |
U0=0.01 # amplitude of point source |
|
wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt) |
|
|
\end{python} |
|
|
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. |
|
258 |
|
|
259 |
|
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.) |
260 |
% \begin{figure} |
wavePropagation(mydomain,h,tend,lam,mu,rho,U0) |
261 |
% \centerline{\includegraphics[width=\figwidth]{DiffusionDomain}} |
\end{python} |
262 |
% \caption{Temperature Diffusion Problem with Circular Heat Source} |
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. Before running the code make sure you have created a directory called \code{data} in the current |
263 |
% \label{WAVE FIG 1} |
working directory. |
264 |
% \end{figure} |
To visualize the results from the data directory: |
265 |
|
\begin{verbatim} mayavi -d usoln.1.vtu -m SurfaceMap &\end{verbatim} You can rotate this figure by clicking on it with the mouse and moving it around. |
266 |
\begin{figure} |
Again use \code{Configure Data} to move backwards and forward in time, and |
267 |
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} |
also to choose the results (acceleration, displacement or $u\hackscore x$) by using \code{Select Scalar}. Figure \ref{WAVE FIG 2} shows the results for the displacement at various time steps. |
268 |
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} |
|
269 |
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} |
\begin{figure}{t!} |
270 |
%\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} |
\centerline{\includegraphics[width=3.in,angle=270]{WavePC.ps}} |
271 |
\caption{Selected time steps of a wave propagation over a $10km \times 100km \times 100km$ |
\caption{Amplitude at Point source} |
272 |
with time step size $h=0.0666$. Color represents the acceleration.} |
\label{WAVE FIG 1} |
273 |
|
\end{figure} |
274 |
|
|
275 |
|
\begin{figure}{t} |
276 |
|
\begin{center} |
277 |
|
\includegraphics[width=2in,angle=270]{Wavet0.ps} |
278 |
|
\includegraphics[width=2in,angle=270]{Wavet1.ps} |
279 |
|
\includegraphics[width=2in,angle=270]{Wavet3.ps} |
280 |
|
\includegraphics[width=2in,angle=270]{Wavet10.ps} |
281 |
|
\includegraphics[width=2in,angle=270]{Wavet30.ps} |
282 |
|
\includegraphics[width=2in,angle=270]{Wavet288.ps} |
283 |
|
\end{center} |
284 |
|
\caption{Selected time steps ($n = 0,1,30,100,300$ and 2880) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block |
285 |
|
from a point source initially at $(5\mbox{km},5\mbox{km},0)$ |
286 |
|
with time step size $h=0.02083$. Color represents the displacement. |
287 |
|
Here the view is oriented onto the bottom face.} |
288 |
\label{WAVE FIG 2} |
\label{WAVE FIG 2} |
289 |
\end{figure} |
\end{figure} |