| 1 |
% |
| 2 |
% $Id$ |
| 3 |
% |
| 4 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 5 |
% |
| 6 |
% Copyright 2003-2007 by ACceSS MNRF |
| 7 |
% Copyright 2007 by University of Queensland |
| 8 |
% |
| 9 |
% http://esscc.uq.edu.au |
| 10 |
% Primary Business: Queensland, Australia |
| 11 |
% Licensed under the Open Software License version 3.0 |
| 12 |
% http://www.opensource.org/licenses/osl-3.0.php |
| 13 |
% |
| 14 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 15 |
% |
| 16 |
|
| 17 |
\section{3-D Wave Propagation} |
| 18 |
\label{WAVE CHAP} |
| 19 |
|
| 20 |
In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation |
| 21 |
\index{wave equation}: |
| 22 |
\begin{eqnarray}\label{WAVE general problem} |
| 23 |
\rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0 |
| 24 |
\end{eqnarray} |
| 25 |
in a three dimensional block of length $L$ in $x\hackscore{0}$ |
| 26 |
and $x\hackscore{1}$ direction and height $H$ |
| 27 |
in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location. |
| 28 |
$\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by |
| 29 |
\begin{eqnarray} \label{WAVE stress} |
| 30 |
\sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) |
| 31 |
\end{eqnarray} |
| 32 |
where $\lambda$ and $\mu$ are the Lame coefficients |
| 33 |
\index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. |
| 34 |
On the boundary the normal stress is given by |
| 35 |
\begin{eqnarray} \label{WAVE natural} |
| 36 |
\sigma\hackscore{ij}n\hackscore{j}=0 |
| 37 |
\end{eqnarray} |
| 38 |
for all time $t>0$. |
| 39 |
|
| 40 |
At initial time $t=0$ the displacement |
| 41 |
$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: |
| 42 |
\begin{eqnarray} \label{WAVE initial} |
| 43 |
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 \; |
| 44 |
\end{eqnarray} |
| 45 |
for all $x$ in the domain. |
| 46 |
|
| 47 |
Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we |
| 48 |
set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point |
| 49 |
$x\hackscore C$. |
| 50 |
|
| 51 |
We use an explicit time integration scheme to calculate the displacement field $u$ at |
| 52 |
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$ |
| 53 |
which is defined by |
| 54 |
\begin{eqnarray} \label{WAVE dyn 2} |
| 55 |
u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ |
| 56 |
\end{eqnarray} |
| 57 |
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form |
| 58 |
\begin{eqnarray} \label{WAVE dyn 2b} |
| 59 |
u\hackscore{,tt}=G(u) |
| 60 |
\end{eqnarray} |
| 61 |
where one sets $a^{(n)}=G(u^{(n-1)})$. |
| 62 |
|
| 63 |
In our case $a^{(n)}$ is given by |
| 64 |
\begin{eqnarray}\label{WAVE dyn 3} |
| 65 |
\rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j} |
| 66 |
\end{eqnarray} |
| 67 |
and boundary conditions |
| 68 |
\begin{eqnarray} \label{WAVE natural at n} |
| 69 |
\sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0 |
| 70 |
\end{eqnarray} |
| 71 |
derived from \eqn{WAVE natural} where |
| 72 |
\begin{eqnarray} \label{WAVE dyn 3a} |
| 73 |
\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}). |
| 74 |
\end{eqnarray} |
| 75 |
Now we have converted our problem for displacement, $u^{(n)}$, into a problem for |
| 76 |
acceleration, $a^(n)$, which now depends |
| 77 |
on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. |
| 78 |
|
| 79 |
In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will |
| 80 |
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through |
| 81 |
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is |
| 82 |
\begin{equation}\label{WAVE dyn 100} |
| 83 |
D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . |
| 84 |
\end{equation} |
| 85 |
The natural boundary condition |
| 86 |
\begin{equation}\label{WAVE dyn 101} |
| 87 |
n\hackscore{j}X\hackscore{ij} =0 |
| 88 |
\end{equation} |
| 89 |
is used. |
| 90 |
|
| 91 |
With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: |
| 92 |
\begin{equation}\label{WAVE dyn 6} |
| 93 |
\begin{array}{ll} |
| 94 |
D\hackscore{ij}=\rho \delta\hackscore{ij}& |
| 95 |
X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ |
| 96 |
\end{array} |
| 97 |
\end{equation} |
| 98 |
|
| 99 |
|
| 100 |
%Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero |
| 101 |
%(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. |
| 102 |
|
| 103 |
The following script defines a the function \function{wavePropagation} which |
| 104 |
implements the Verlet scheme to solve our wave propagation problem. |
| 105 |
The argument \var{domain} which is a \Domain class object |
| 106 |
defines the domain of the problem. \var{h} and \var{tend} are the time step size |
| 107 |
and the end time of the simulation. \var{lam}, \var{mu} and |
| 108 |
\var{rho} are material properties. |
| 109 |
\begin{python} |
| 110 |
from esys.linearPDEs import LinearPDE |
| 111 |
from numarray import identity,zeros,ones |
| 112 |
from esys.escript import * |
| 113 |
from esys.escript.pdetools import Locator |
| 114 |
def wavePropagation(domain,h,tend,lam,mu,rho,U0): |
| 115 |
x=domain.getX() |
| 116 |
# ... open new PDE ... |
| 117 |
mypde=LinearPDE(domain) |
| 118 |
mypde.setSolverMethod(LinearPDE.LUMPING) |
| 119 |
kronecker=identity(mypde.getDim()) |
| 120 |
# spherical source at middle of bottom face |
| 121 |
xc=[width/2.,width/2.,0.] |
| 122 |
# define small radius around point xc |
| 123 |
# Lsup(x) returns the maximum value of the argument x |
| 124 |
src_radius = 0.1*Lsup(domain.getSize()) |
| 125 |
print "src_radius = ",src_radius |
| 126 |
dunit=numarray.array([1.,0.,0.]) # defines direction of point source |
| 127 |
mypde.setValue(D=kronecker*rho) |
| 128 |
# ... set initial values .... |
| 129 |
n=0 |
| 130 |
# initial value of displacement at point source is constant (U0=0.01) |
| 131 |
# for first two time steps |
| 132 |
u=U0*whereNegative(length(x-xc)-src_radius)*dunit |
| 133 |
u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit |
| 134 |
t=0 |
| 135 |
# define the location of the point source |
| 136 |
L=Locator(domain,xc) |
| 137 |
# find potential at point source |
| 138 |
u_pc=L.getValue(u) |
| 139 |
print "u at point charge=",u_pc |
| 140 |
u_pc_x = u_pc[0] |
| 141 |
u_pc_y = u_pc[1] |
| 142 |
u_pc_z = u_pc[2] |
| 143 |
# open file to save displacement at point source |
| 144 |
u_pc_data=open('./data/U_pc.out','w') |
| 145 |
u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) |
| 146 |
while t<tend: |
| 147 |
# ... get current stress .... |
| 148 |
g=grad(u) |
| 149 |
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) |
| 150 |
# ... get new acceleration .... |
| 151 |
mypde.setValue(X=-stress) |
| 152 |
a=mypde.getSolution() |
| 153 |
# ... get new displacement ... |
| 154 |
u_new=2*u-u_last+h**2*a |
| 155 |
# ... shift displacements .... |
| 156 |
u_last=u |
| 157 |
u=u_new |
| 158 |
t+=h |
| 159 |
n+=1 |
| 160 |
print n,"-th time step t ",t |
| 161 |
L=Locator(domain,xc) |
| 162 |
u_pc=L.getValue(u) |
| 163 |
print "u at point charge=",u_pc |
| 164 |
u_pc_x=u_pc[0] |
| 165 |
u_pc_y=u_pc[1] |
| 166 |
u_pc_z=u_pc[2] |
| 167 |
# save displacements at point source to file for t > 0 |
| 168 |
u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) |
| 169 |
# ... save current acceleration in units of gravity and displacements |
| 170 |
if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81, |
| 171 |
displacement = length(u), Ux = u[0] ) |
| 172 |
u_pc_data.close() |
| 173 |
\end{python} |
| 174 |
Notice that |
| 175 |
all coefficients of the PDE which are independent of time $t$ are set outside the \code{while} |
| 176 |
loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible |
| 177 |
when iterating over time. |
| 178 |
|
| 179 |
there are a few new \escript functions in this example: |
| 180 |
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). |
| 181 |
There are restrictions on the argument of the \function{grad} function, for instance |
| 182 |
the statement \code{grad(grad(u))} will raise an exception. |
| 183 |
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} |
| 184 |
and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$). |
| 185 |
|
| 186 |
We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small |
| 187 |
sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$. |
| 188 |
These satisfy the initial conditions defined in \eqn{WAVE initial}. |
| 189 |
|
| 190 |
The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}. |
| 191 |
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ |
| 192 |
respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of |
| 193 |
the displacement vector $u$ at the point source. These can be |
| 194 |
plotted easily using any plotting package. In gnuplot the command: |
| 195 |
\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, |
| 196 |
'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}. |
| 197 |
%\begin{python} |
| 198 |
%u=Vector(0.,ContinuousFunction(domain)) |
| 199 |
%\end{python} |
| 200 |
%declares \var{u} as a vector-valued function of its coordinates |
| 201 |
%in the domain (i.e. with two or three components in the case of |
| 202 |
%a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal |
| 203 |
%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} |
| 204 |
%specifies the `smoothness` of the function. In our case we want the displacement field to be |
| 205 |
%a continuous function over the \Domain \var{domain}. On other option would be |
| 206 |
%\var{Function(domain)} in which case continuity is not guaranteed. In fact the |
| 207 |
%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while |
| 208 |
%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. |
| 209 |
|
| 210 |
The statement \code{myPDE.setLumpingOn()} |
| 211 |
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix |
| 212 |
formed from the coefficient \var{D}. |
| 213 |
The approximation allows, at the cost of |
| 214 |
additional error, very fast |
| 215 |
solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var |
| 216 |
{C} will produce wrong results. |
| 217 |
|
| 218 |
|
| 219 |
One of the big advantages of the Verlet scheme is the fact that the problem to be solved |
| 220 |
in each time step is very simple and does not involve any spatial derivatives (which is what allows us to use lumping in this simulation). |
| 221 |
The problem becomes so simple because we use the stress from the last time step rather then the stress which is |
| 222 |
actually present at the current time step. Schemes using this approach are called an explicit time integration |
| 223 |
schemes \index{explicit scheme} \index{time integration!explicit}. The |
| 224 |
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is |
| 225 |
an example of an implicit scheme |
| 226 |
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of |
| 227 |
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is |
| 228 |
more expensive to solve, in particular for non-linear problems. |
| 229 |
Although |
| 230 |
explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time |
| 231 |
steps then implicit schemes and can suffer from stability problems. Therefore they need a |
| 232 |
very careful selection of the time step size $h$. |
| 233 |
|
| 234 |
An easy, heuristic way of choosing an appropriate |
| 235 |
time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition} |
| 236 |
which says that within a time step a information should not travel further than a cell used in the |
| 237 |
discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as |
| 238 |
$\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from |
| 239 |
\begin{eqnarray}\label{WAVE dyn 66} |
| 240 |
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x |
| 241 |
\end{eqnarray} |
| 242 |
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of |
| 243 |
the formula. |
| 244 |
|
| 245 |
The following script uses the \code{wavePropagation} function to solve the |
| 246 |
wave equation for a point source located at the bottom face of a block. The width of the block in |
| 247 |
each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}). |
| 248 |
\code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions. |
| 249 |
The depth of the block is aligned with the $x\hackscore{2}$-direction. |
| 250 |
The depth (\code{l2}) of the block in |
| 251 |
the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the |
| 252 |
of the depth is chosen such that the elements |
| 253 |
become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the |
| 254 |
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 |
| 255 |
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]$. |
| 256 |
\begin{python} |
| 257 |
from esys.finley import Brick |
| 258 |
ne=32 # number of cells in x_0 and x_1 directions |
| 259 |
width=10000. # length in x_0 and x_1 directions |
| 260 |
lam=3.462e9 |
| 261 |
mu=3.462e9 |
| 262 |
rho=1154. |
| 263 |
tend=60 |
| 264 |
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne) |
| 265 |
print "time step size = ",h |
| 266 |
U0=0.01 # amplitude of point source |
| 267 |
|
| 268 |
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.) |
| 269 |
wavePropagation(mydomain,h,tend,lam,mu,rho,U0) |
| 270 |
\end{python} |
| 271 |
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 |
| 272 |
working directory. |
| 273 |
To visualize the results from the data directory: |
| 274 |
\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. |
| 275 |
Again use \code{Configure Data} to move backwards and forward in time, and |
| 276 |
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. |
| 277 |
|
| 278 |
\begin{figure}{t!} |
| 279 |
\centerline{\includegraphics[width=3.in,angle=270]{figures/WavePC.eps}} |
| 280 |
\caption{Amplitude at Point source} |
| 281 |
\label{WAVE FIG 1} |
| 282 |
\end{figure} |
| 283 |
|
| 284 |
\begin{figure}{t} |
| 285 |
\begin{center} |
| 286 |
\includegraphics[width=2in,angle=270]{figures/Wavet0.eps} |
| 287 |
\includegraphics[width=2in,angle=270]{figures/Wavet1.eps} |
| 288 |
\includegraphics[width=2in,angle=270]{figures/Wavet3.eps} |
| 289 |
\includegraphics[width=2in,angle=270]{figures/Wavet10.eps} |
| 290 |
\includegraphics[width=2in,angle=270]{figures/Wavet30.eps} |
| 291 |
\includegraphics[width=2in,angle=270]{figures/Wavet288.eps} |
| 292 |
\end{center} |
| 293 |
\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 |
| 294 |
from a point source initially at $(5\mbox{km},5\mbox{km},0)$ |
| 295 |
with time step size $h=0.02083$. Color represents the displacement. |
| 296 |
Here the view is oriented onto the bottom face. |
| 297 |
\label{WAVE FIG 2}} |
| 298 |
\end{figure} |
Parent Directory
| 
Revision Log