Key Exercise 1 (Deadline: Mar 10)

Implement a general Python program for solving the 2D wave equation with damping, variable wave velocity, and du/dn=0 boundary conditions, on a rectangle-shaped domain with uniform mesh partition. Various extensions by different groups can explore how to migrate CPU-time intensive loops to compiled code (Fortran/C).

Apply the program to solve a wave equation problem where the wave velocity is very small in a square-shaped subdomain in the center of the domain, while the wave velocity is of order unity outside this subdomain. As initial condition, choose a u(x,y,0)=cos(pi*x) for x<1/2 and u=0 otherwise, and make sure the length of the domain in x direction much larger than the extent of this initial disturbance (say 0 <= x <=L, with L=20). This initial condition leads to an incoming plane left from the left (x=0), which is disturbed by the subdomain with low wave velocity in the center of the domain. Stop the simulation before the wave hits the right boundary (x=L) and gets reflected there by the boundary condition du/dn=0. Present results with and without damping.

How can you verify that the program works? Or equivalently, if the program doesn't work, how can you find errors? Suggestions: 1) create a test problem with u=const as solution; 2) choose a constant wave velocity, Courant number 1, no damping, wave propagating in the x direction, and use the property that 1D waves are exactly reproduced by the scheme, regardless of the spatial and temporal mesh partition; 3) argue that the solution observed in the real application is physically reasonable (and point out numerical artifacts if you can see such).

Extensions of Key Exercise 1

Goal: learn how to speed up slow loops in Python by migrating code to compiled languages (Fortran, C, C++).

Information on combining Python with Fortran or C++ can be found in Chapters 9 and 10 of a book that can be downloaded by

wget http://simula.no/~hpl/scripting/Langtangen_TCSE3_3rd.pdf

There is also a paper on the effect of using Cython to speed up a wave equation solver. The new tool for treating Fortran code, fwrap, builds on Cython.

Tentative suggestions:

Group 1 and 3: Implement an open boundary condition at x=L, which can let the wave propagate out of the domain (without being reflected). This condition reads u_t + c*u_x = 0. For constant c in the whole domain and a 1D wave, the condition is exact (demonstrate that!), but for 2D waves it will give some reflections.
Group 2: Implement the core loop over spatial points using vectorized expressions in Python.
Group 4: Visualize the waves using MayaVi2 or some other fancy visualization software. Try to visualize both the waves (u) and the variable wave velocity (c) in the same scene. Nice extension: 3D wave equation!
Group 5: Implement the core loop over spatial points in Fortran 77 using fwrap.
Group 6: Implement the core loop over spatial points in C++ using Instant.
Group 7: Implement the core loop over spatial points in Fortran 77 using f2py.

Troubleshooting

Strange solutions

Make sure you test the program in a stepwise fashion as outlined in the lecture on March 10 (start with a constant solution, use a 1D program in parallel with the 2D program, introduce a piecewise constant wave velocity, use a smooth Gaussian function rather than the suggested cosine-"hat" as suggested, etc.).

Make sure you use a correct stability criterion - the max timestep in 2D is a factor 1/sqrt(2) smaller than in 1D.

Problems with plotting

If you use the Easyviz interface to Matplotlib, Matlab, Gnuplot, VTK, etc., it can be wise to work with the latest version. Download and install it by the following commands:

hg clone https://scitools.googlecode.com/hg/ scitools
cd scitools
python setup.py install

When scitools is updated, just go to the scitools directory and run the python setup.py install command.

By default, scitools applies matplotlib for plotting, but matplotlib is not very suited for animating elevated surfaces u(x,y,t) (it runs very slowly). Gnuplot is a better choice. Go into the file

lib/scitools/scitools.cfg

and change the line

backend = matplotlib

backend = gnuplot

Reinstall. Now, gnuplot is the default plotting engine, and it is well suited for animating 2D scalar fields.