In this next example we want to calculate the displacement field
for any time
by solving the wave equation:
Here we are modelling a point source at the point
in the
-direction
which is a negative pulse of amplitude
followed by the same
positive pulse. In mathematical terms we use
We use an explicit time integration scheme to calculate the displacement field at
certain time marks
where
with time step size
. In the following the upper index
refers to values at time
. We use the Verlet scheme with constant time step size
which is defined by
In our case is given by
In each time step we have to solve this problem to get the acceleration , 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
The statement
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
of
(in fact grad(g)[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.
).
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
) from functions such as the displacement vector u. Typically the Locator is used in the following form:
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 .
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
so one should choose
from
![]() ![]() ![]() ![]() ![]() |
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
km (
and
directions ie. l0 and l1).
The ne gives the number of elements in
and
directions.
The depth of the block is aligned with the
-direction.
The depth (l2) of the block in
the
-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
-direction so that the
computation would be faster. Brick(
) is an esys.finley function which creates a rectangular mesh
with
elements over the brick
.
The domain.getSize() return the local element size . 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
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:
The U_pc.out stores 4 columns of data:
respectively, where
are the
components of the displacement vector
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 2will 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:
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