3-D Wave Propagation

In this next example we want to calculate the displacement field $ u\hackscore{i}$ for any time $ t>0$ by solving the wave equation:

$\displaystyle \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0$     (30)

in a three dimensional block of length $ L$ in $ x\hackscore{0}$ and $ x\hackscore{1}$ direction and height $ H$ in $ x\hackscore{2}$ direction. $ \rho$ is the known density which may be a function of its location. $ \sigma\hackscore{ij}$ is the stress field which in case of an isotropic, linear elastic material is given by
$\displaystyle \sigma\hackscore{ij}$ $\displaystyle =$ $\displaystyle \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})$ (31)

where $ \lambda$ and $ \mu$ are the Lame coefficients and $ \delta\hackscore{ij}$ denotes the Kronecker symbol. On the boundary the normal stress is given by
$\displaystyle \sigma\hackscore{ij}n\hackscore{j}=0$     (32)

for all time $ t>0$.

Figure 1.7: Input Displacement at Source Point ( $ \alpha =0.7$, $ t\hackscore {0}=3$, $ U\hackscore {0}=1$).

Figure 1.8: Input Acceleration at Source Point ( $ \alpha =0.7$, $ t\hackscore {0}=3$, $ U\hackscore {0}=1$).

Here we are modelling a point source at the point $ x\hackscore C$ in the $ x\hackscore{0}$-direction which is a negative pulse of amplitude $ U\hackscore{0}$ followed by the same positive pulse. In mathematical terms we use

$\displaystyle u\hackscore{0}(x\hackscore C,t)= U\hackscore{0} \; \sqrt{2} \frac...
...hackscore{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}} \ $     (33)

for all $ t\ge0$ where $ \alpha$ is the width of the pulse and $ t\hackscore{0}$ is the time when the pulse changes from negative to positive. In the simulations we will choose $ \alpha=0.3$ and $ t\hackscore{0}=2$, see Figure 1.7 and will apply the source as a constraint in a in a sphere of small radius around the point $ x\hackscore C$.

We use an explicit time integration scheme to calculate the displacement field $ u$ at 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 with constant time step size $ h$ which is defined by

$\displaystyle u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)}$     (34)

for all $ n=2,3,\ldots$. It is designed to solve a system of equations of the form
$\displaystyle u\hackscore{,tt}=G(u)$     (35)

where one sets $ a^{(n)}=G(u^{(n-1)})$.

In our case $ a^{(n)}$ is given by

$\displaystyle \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}$     (36)

and boundary conditions
$\displaystyle \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0$     (37)

derived from Equation (1.33) where
$\displaystyle \sigma\hackscore{ij}^{(n-1)}$ $\displaystyle =$ $\displaystyle \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}).$ (38)

We also need to apply the constraint
$\displaystyle a^{(n)}\hackscore{0}(x\hackscore C,t)= U\hackscore{0}
\frac{\sqr...
...t-t\hackscore{0}}{\alpha})e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}}$     (39)

which is derived from equation 1.34 by calculating the second order time derivative, see Figure 1.8. Now we have converted our problem for displacement, $ u^{(n)}$, into a problem for acceleration, $ a^{(n)}$, which now depends on the solution at the previous two time steps, $ u^{(n-1)}$ and $ u^{(n-2)}$.

In each time step we have to solve this problem to get the acceleration $ a^{(n)}$, and we will use the LinearPDE class of the esys.escript.linearPDEs to do so. The general form of the PDE defined through the LinearPDE class is discussed in Section 4.1. The form which is relevant here is

$\displaystyle D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .$ (40)

The natural boundary condition

$\displaystyle n\hackscore{j}X\hackscore{ij} =0$ (41)

is used. With $ u=a^{(n)}$ we can identify the values to be assigned to $ D$ and $ X$:

\begin{displaymath}\begin{array}{ll} D\hackscore{ij}=\rho \delta\hackscore{ij}& X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ \end{array}\end{displaymath} (42)

Moreover we need to define the location $ r$ where the constraint 1.41 is applied. We will apply the constraint on a small sphere of radius $ R$ around $ x\hackscore C$ (we will use 3p.c. of the width of the domain):

$\displaystyle q\hackscore{i}(x) = \left\{ \begin{array}{rc} 1 & \Vert x-x_c\Vert\le R \\ 0 & \mbox{otherwise} \end{array} \right.$ (43)

The following script defines a the function wavePropagation which implements the Verlet scheme to solve our wave propagation problem. The argument domain which is a Domain class object defines the domain of the problem. h and tend are the time step size and the end time of the simulation. lam, mu and rho are material properties.
\begin{python}
def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0...
...\
tensor = stress, Ux = u[0] )
\par
return ts, u_pc0, u_pc1, u_pc2
\end{python}
Notice that all coefficients of the PDE which are independent of time $ t$ are set outside the while loop. This is very efficient since it allows the LinearPDE class to reuse information as much as possible when iterating over time.

The statement
\begin{python}
mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
\end{python}
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix formed from the coefficient D. The approximation allows, at the cost of additional error, very fast solution of the PDE. When using lumping the presence of A, B or C will produce wrong results.

There are a few new esys.escript functions in this example: grad(u) returns the gradient $ u\hackscore{i,j}$ of $ u$ (in fact grad(g)[i,j]= $ =u\hackscore{i,j}$). There are restrictions on the argument of the grad function, for instance the statement grad(grad(u)) will raise an exception. trace(g) returns the sum of the main diagonal elements g[k,k] of g and transpose(g) returns the matrix transpose of g (ie. $ \var{transpose(g)[i,j]}=\var{g[j,i]}$).

We initialize the values of u and u_last to be zero. It is important to initialize both with the solution FunctionSpace FunctionSpace as they have to be seen as solutions of PDEs from previous time steps. In fact, the grad does not accept arguments with a certain FunctionSpace, for more details see Section 3.1.3.

The Locator is designed to extract values at a given location (in this case $ x\hackscore C$) from functions such as the displacement vector u. Typically the Locator is used in the following form:
\begin{python}
L=Locator(domain,xc)
u=...
u_pc=L.getValue(u)
\end{python}
The return value u_pc is the value of u at the location xc[*]. The values are collected in the lists u_pc0, u_pc1 and u_pc2 together with the corresponding time marker in ts. The values are handed back to the caller. Later we will show to ways to access these data.

One of the big advantages of the Verlet scheme is the fact that the problem to be solved 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). The problem becomes so simple because we use the stress from the last time step rather then the stress which is actually present at the current time step. Schemes using this approach are called an explicit time integration schemes . The backward Euler scheme we have used in Chapter 1.3 is an example of an implicit scheme . In this case one uses the actual status of each variable at a particular time rather then values from previous time steps. This will lead to a problem which is more expensive to solve, in particular for non-linear problems. Although explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time steps then implicit schemes and can suffer from stability problems. Therefore they need a very careful selection of the time step size $ h$.

An easy, heuristic way of choosing an appropriate time step size is the Courant condition which says that within a time step a information should not travel further than a cell used in the discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as $ \sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $ h$ from

$\displaystyle h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x$     (44)

where $ \Delta x$ is the cell diameter. The factor $ \frac{1}{5}$ is a safety factor considering the heuristics of the formula.

Figure: Selected time steps ( $ n = 11, 22, 32, 36$) of a wave propagation over a $ 10$km$ \times 10$km$ \times 3.125$km block from a point source initially at $ (5$km$ ,5$km$ ,0)$ with time step size $ h=0.02083$. Color represents the displacement. Here the view is oriented onto the bottom face.
\includegraphics[width=2in]{figures/Wave11} \includegraphics[width=2in]{figures/Wave22} \includegraphics[width=2in]{figures/Wave28} \includegraphics[width=2in]{figures/Wave32} \includegraphics[width=2in]{figures/Wave36}

The following script uses the wavePropagation function to solve the wave equation for a point source located at the bottom face of a block. The width of the block in each direction on the bottom face is $ 10$km ( $ x\hackscore 0$ and $ x\hackscore 1$ directions ie. l0 and l1). The ne gives the number of elements in $ x\hackscore{0}$ and $ x\hackscore 1$ directions. The depth of the block is aligned with the $ x\hackscore{2}$-direction. The depth (l2) of the block in the $ x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the of the depth is chosen such that the elements become cubic. We chose 10 for the number of elements in $ x\hackscore{2}$-direction so that the computation would be faster. Brick( $ n\hackscore 0, n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$) is an esys.finley function which creates a rectangular mesh 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]$.
\begin{python}
from esys.finley import Brick
ne=32  ...
The domain.getSize() return the local element size $ \Delta x$. The inf makes sure that the Courant condition 1.46 olds everywhere in the domain.

The script is available as wave.py in the example directory . To visualize the results from the data directory:

 
mayavi2 -d usoln.1.vtu -m SurfaceMap &
You can rotate this figure by clicking on it with the mouse and moving it around. Again use Configure Data to move backwards and forward in time, and also to choose the results (acceleration, displacement or $ u\hackscore x$) by using Select Scalar. Figure 1.9 shows the results for the displacement at various time steps.

Figure 1.10: Amplitude at Point source from the Simulation.

It remains to show some possibilities to inspect the collected data u_pc0, u_pc1 and u_pc2. One way is to write the data to a file and then use an external package such as gnuplot [26], excel or OpenOffice to read the data for further analysis. The following code shows one possible form to write the data to the file ./data/U_pc.out:
\begin{python}
u_pc_data=FileWriter('./data/U_pc.out')
for i in xrange(len(ts)) ...
... %f %f %f\n''%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
u_pc_data.close()
\end{python}
The U_pc.out stores 4 columns of data: $ t,u\hackscore x,u\hackscore y,u\hackscore z$ respectively, where $ u\hackscore x,u\hackscore y,u\hackscore z$ are the $ x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of the displacement vector $ u$ at the point source. These can be plotted easily using any plotting package. In gnuplot [26]the command:

 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, 
'U_pc.out' u 1:4 title 'U_z' w l lw 2
will reproduce Figure 1.10 (As expected this is identical to the input signal shown in Figure 1.7)2. It is pointed out that we are not using the standard python open to write to the file U_pc.out as it is not safe when running esys.escript under MPI, see chapter 2 for more details.

Alternatively, one can implement plotting the results at run time rather than in a post-processing step. This avoids the generation of an intermediate data file. In escript the preferred way of creating 2D plots of time dependent data is matplotlib. The following script creates the plot and writes it into the file u_pc.png in the PNG image format:
\begin{python}
import matplotlib.pyplot as plt
if getMPIRankWorld() == 0:
plt.t...
...displacement')
plt.legend()
plt.savefig('u_pc.png', format='png')
\end{python}
You can add the plt.show() to create a interactive browser window. Please not that through the getMPIRankWorld() == 0 statement the plot is generated on one processor only (in this case the rank 0 processor) when run under MPI.

Both options for processing the point source data are include in the example file wave.py. There other options available to process these data in particular through the SciPy[7] package , eg Fourier transformations. It is beyond the scope of this users guide to document the usage of SciPy[7] for time series analysis but is highly recommended that users use SciPy[7] to any further data processing.

esys@esscc.uq.edu.au