The next example is more interesting than parallelPi. Our objective is to solve a two-dimensional diffusion PDE:

where u(x,y,t) is the temperature of the surface being modeled, and D is the so-called diffusion coefficient.

We will choose as the surface a rectangular plate:

It will be comforting as we test our code to know an exact solution of the diffusion equation.

By examination of the diffusion equation, one may readily write down an exact solution:

which is consistent with Dirichlet boundary conditions of u = 0 along the edges of the surface, and with the initial condition:

However, it is easier to interpret numerical results for a uniform initial condition. The entire plate is at the same initial temperature, and at time t = 0 is immersed into a zero-temperature bath. To obtain an exact solution consistent with a uniform initial condition we use the method of separation of variables. Again, the Dirichlet boundary conditions we choose is temperature u equals 0 along the edges of the surface.

We write the solution as a product of three single-variable functions: u(x,y,t) = X(x) Y(y) T(t). The diffusion equation then looks like this:

If this relationship is to hold for arbitrary (x,y,z), then each term must be constant. This leads to three single-variable ordinary differential equations plus a relationship for the constants:

We chose the sign and form of these constants to make it easy to satisfy the thermodynamic law of entropy coupled with our Dirichlet boundary conditions. The temperature of the surface must vanish at the edges and decrease with time.

The solution of the ODE for T is a decreasing exponential function. To satisfy the boundary conditions the solutions of the X and Y ODEs are sine functions:

Therefore, the general solution of the diffusion equation for the rectangular plate given the zero-temperature boundary conditions is:

where the expansion coefficients A are to be determined by the initial conditions, the temperature of the rectangular plate at t = 0, u(x,y,0) = f(x,y).

Using the orthogonality of the sine function:

We obtain an expression for the expansion coefficients:

Now let's assume the simplist initial condition, a uniform temperature on the plate (except for the boundary of course), f(x,y) = Uo. Then it is easy to show that:

And at last the exact solution for u = 0 boundary conditions and given a uniform initial temperature Uo is:

This is of course an infinite sum, but with an exponential term that decreases with increasing indices m,n. In choosing a cutoff for m,n be aware that for small time t many terms may be needed to achieve satisfactory accuracy. In the section below, we will look at a finite difference based solution to our diffusion equation.

We superimpose a computational grid upon the plate we are modeling:

which shows the numbering of the x and y grid indices i and j, respectively. The number of interior grid points in the x-direction is nni, in the y-direction nnj. The boundaries are i = 0, i = nni + 1, j = 0, and j = nnj + 1.

Clearly the x-direction and y-direction step sizes are:

We'll use a friendly notation for the temperature:

where we've introduced the time index n and the time step size dt.

We'll use the Forward Time Centered Space (FTCS) finite-difference approximation to the diffusion equation:

where i = 1, 2, ... , nni and j = 1, 2, ... , nnj. Also n = 0, 1, 2, ... .

In practice you will not want to use the unsophisticated FTCS representation because of its instability. The instability is avoided, however, if the time step size dt obeys the so-called Courant condition:

The disadvantage of this condition is that for complicated geometries for which dx and dy must be small, dt becomes very small, and many iterations are required to evolve the system though time.

The finite-difference equation may be rearranged:

Using the temperature at time step n, the temperature at time step n + 1 is computed.

Later when considering parallelization we will be especially interested in the stencil. The stencil is the set of grid points needed to compute the single ij grid point at the next time step. For the above finite-difference representation of the two-dimensional diffusion equation, we have a 5-point stencil that must be computed to obtain u_ij at the next time step:

For computing u_ij adjacent to a boundary, one of these stencil points will be a boundary point.