Computational Physics Lab 5

Boundary Value Problems

This lab has two goals:

Learn to overcome python execution speed bottlenecks by writing python modules in C++.
Explore basic PDE solving techniques for boundary value problems.

Part 1. Combining Python and C++

In the practice of solving computational problems, you will often come to a point where the algorithm execution speed becomes a bottleneck of your method. In fact, reaching this point in Python is not that difficult: iterating over large arrays and performing arithmetic operations on array members could be couple orders of magnitude slower in pure Python than in compiled languages such as C/C++ or Fortran.

So, why not just program everything in C++ or Fotran? There are many reasons. Perhaps, the most compelling one is that time-critical parts of your code usually comprise only a small fraction of the complete program. Virtually all program developers find that they are significantly more productive when they use high-level languages such as Python or Ruby. Many features of well-designed interpreted languages contribute to this phenomenon: clean syntax, consistent exception handling policies, automatic garbage collection, built-in heterogeneous containers, named function arguments with reasonable defaults, rich standard library, simple interfaces to graphics and GUI tools, high-level networking and database access, etc. Absence of the compilation step and interactive debugging capabilities are also helpful. Python simply makes a more efficient development platform than, let say, C++, and it usually does much more per line of code.

At times, however, the overall payoff from fast code execution exceeds the payoff from efficient program development. Solving differential equations on large multidimensional grids is one such example. Mixed language programming can be a very effective approach for this type of problem -- code the time-critical parts of the numerical algorithm in C++ or some other compiled language, but use Python to steer the computations and to build the complete solution of the problem. Besides the number crunching, construction of the complete solution usually involves setting up the computational framework, inputting the data, and presenting the results in a graphical form. Python is perfectly suitable for programming these tasks. Commercial numerical modeling and data analysis environments (such as Matlab and SAS) follow the same approach by allowing their users to write C/C++ extensions.

I am going to assume that you have basic understanding of C++. If this is not true please notify me, and we will figure something out (a few tutorial sessions can be arranged).

In this part of the lab, you will interface two simple C++ classes to Python using the SWIG tool (SWIG stands for "Simplified Wrapper and Interface Generator"). This interfacing method is quite well documented in chapters 5 and 10 of the H.P. Langtangen's book Python Scripting for Computational Science (one of the books recommended in the course syllabus). The book is available here in the electronic form. You can figure out the download/pdf password if you purchase or borrow the printed version of the book. Other notable tools for Python/C++ interfacing are Boost.Python, SIP, and weave. A somewhat different approach is taken by PyCXX which provides a C++ wrapper to Python C API and aids in writing Python modules in C++ from scratch.

Please download the code for the "HelloWorld" and "HelloWorld2" classes (this example is discussed in detail in section 5.2.4 of the Langtangen's book): HelloWorld.h, HelloWorld.cpp, HelloWorld2.h, and HelloWorld2.cpp. The class "HelloWorld2" inherits from "HelloWorld" and extends it in a rather trivial way by adding the "gets" function. In order to provide an access to these classes from Python, SWIG needs an interface definition file: hw.i. Here are the file contents:

/* file: hw.i */
%module hw
%{
/* include C++ header files necessary to compile the interface */
#include "HelloWorld.h"
#include "HelloWorld2.h"
%}

%include "HelloWorld.h"

%include "typemaps.i"
%apply double* OUTPUT { double& s_ }
%include "HelloWorld2.h"

%extend HelloWorld {
void print_() { self->message(std::cout); }
}

Note the following features of this file:

Everything inside /* */ is a comment, just like in the C language. // at the beginning of a line can be used as well, as in C++.
% indicates that the following statement is a SWIG directive. There are about two dozen different directives that can be used to give SWIG hints, provide annotation, and change SWIG's behavior in some way or another.
The file starts with the %module directive which declares the module name. This directive must precede all other statements in the file.
The purpose of SWIG is to generate a "wrap" file which contains everything that is needed to build a working module for the target scripting language. Everything inside the %{ %} block is copied verbatim into this file. In particular, it is usually necessary to include the header files of the wrapped classes into the wrap file.
The SWIG %include "HelloWorld.h" directive (do not confuse it with the CPP #include) tells SWIG to extract declarations of the functions to be wrapped from the corresponding header file. Note that, in general, this file does not have to be the same as the header file included in the %{ %} block. In particular, one might want to restrict a set of functions exposed in Python if one does not need an access to the full class functionality -- this results in a slimmer interface library.
The %include directive can be followed by any file name. In particular, one can use the name of another interface file.
SWIG offers so-called "typemaps" for dealing with pointers and references that represent output arguments from a function. The file "typemaps.i" which comes with the SWIG distribution contains some ready-made declarations for specifying pointers and references as input, output, or input/output arguments to functions. This is necessary because most scripting languages, including Python, do not support pointer programming at all -- so that C/C++ pointer operations must be mapped into the supported features of the interpreted language.
The statement %apply double* OUTPUT { double& s_ } marks all subsequent appearances of "double& s" in function argument lists as output (this argument appears later in the "gets" function of the "HelloWorld2" class).
SWIG does not know how to wrap C++ streams correctly, so operator << declared in the "HelloWorld.h" header will not be mapped into Python. However, one can define additional class functions as a part of the interface file using the %extend directive. This directive is used here to invoke the "message" method on the standard output. Note the use of self-> to reference the method inside the %extend directive and the use of print_ (with underscore) to name the function -- use of print (without underscore) would be in conflict with the reserved Python keyword "print".

Generate the "hw" module wrapper files as follows:

swig -python -c++ hw.i

This will create a file named "hw_wrap.cxx" and a file named "hw.py". These files contain the wrapper code for the new module. Compile the C++ wrapper together with the original C++ code and make a shared library for the new module like this:

set python_inc = `python -c "import sys; print sys.prefix"`
set python_ver = `python -c "import sys; print sys.version[:3]"`
g++ -fPIC -O -c -I{$python_inc}/include/python{$python_ver} \
HelloWorld.cpp HelloWorld2.cpp hw_wrap.cxx
g++ -shared -o _hw.so HelloWorld.o HelloWorld2.o hw_wrap.o

The "_hw.so" (with leading underscore) should be the name of the shared library file when "hw" is the name of the module. You can test your new module by copying/pasting the following into the Python prompt and observing the output:

from hw import HelloWorld, HelloWorld2

hw = HelloWorld()
r1 = 2
r2 = 3
hw.set(r1, r2)
s = hw.get()
print "Hello, World! sin(%g + %g)=%g" % (r1, r2, s)
hw.print_()

hw2 = HelloWorld2()
hw2.set(r1, r2)
s = hw2.gets()
print "Hello, World2! sin(%g + %g)=%g" % (r1, r2, s)

SWIG supports most C++ features but not all of them. If you intend to mix C++ with python regularly, I suggest writing a SWIG interface file and wrapping your code at an early stage, before writing the complete C++ implementation of your algorithm (you may have to write stubs for the functions you want to expose in Python). In this way you can make sure that you will not encounter difficulties in constructing Python API for your code.

Part 2. Numerical Analysis of a Boundary Value Problem

Using the Python/C++ interfacing technique from Part 1, complete Exercise 5.7 in Giordano & Nakanishi (page 143). Write your own C++ class for solving this problem. Of course, you don't really want to write two different programs as G&N suggest -- just use slightly different schemes for updating grid values from one iteration to the next depending on the relaxation method used (Jacobi or SOR).

Your class should include a method for converting the final result into a numpy array so that Python tools can be used for subsequent visualization. The following simple class implements this conversion (this class can be used as a basic skeleton for your 2d Laplace equation solver): Arr2d.hh, Arr2d.cc. The method of generating the "Arr2d.py" file and the shared library from the SWIG interface file Arr2d.i remains unchanged from Part 1. To check your wrapped module, adapt and use the Arr2d_example.py script.

I suggest using the matshow and/or the contour matplotlib functions for visualizing the results. Here is an example code snippet which illustrates simultaneous usage of these two functions:

from numpy import *
from pylab import *

def fcn(x, y):
    dx = (x-100)/100.0
    dy = (y-100)/100.0
    return exp(-(dx*dx + 2*dy*dy))

a = fromfunction(fcn, (200,200))

matshow(a)
contour(a, 20, colors = 'k')

show()

Submit your lab report by email before 2 pm 02/26/2008. Include the complete code necessary to solve the problem posed in Part 2 of the lab, the plots of the results, and a brief description of the steps needed to compile and run your code.