Implementation of the Reaction Diffusion Simulation

Update – This post was written before the release of Silverlight 3.0b, which provides a number of enhancements relevant to this implementation, such as a WriteableBitmap and Pixel Shaders

Rendering

The first obstacle to implementing the RD simulation is that Silverlight 2.0 does not by default provide a means of generating dynamic images. WPF has a WriteableBitmap, but no equivalent exists in Silverlight. However, it does support PNG streams so we can dynamically update a bitmap by encoding it to PNG on the fly. For this I have used Joe Stegman’s PNGEncoder class, which I have modified slightly to deal with RGB data and to reduce memory usage.

Performance notes

Performance of the simulator is dominated by the time it takes to update the arrays containing the concentration values of the chemicals. Actual rendering time is negligible, even though the image has to be dynamically encoded to a PNG file in order for Silverlight to display it. There are two ways of implementing 2D arrays in Silverlight with C#, rectangular arrays or C style jagged arrays. If this were a desktop application, the most performant method would be to use rectangular arrays and unsafe access, however that option is not available to us in Silverlight so it turns out that, for safe access, jagged arrays are faster.

For example, originally the simulation step looked like this

for (int j = 0; j < H; j++)
{
 int jm = j != 0 ? j - 1 : H - 1;
 int jp = j != H - 1 ? j + 1 : 0;
 for (int i = 0; i < W; i++)
 {
  int im = i != 0 ? i - 1 : W - 1;
  int ip = i != W - 1 ? i + 1 : 0;
  double d2u = U[j,im]+U[j,ip]+U[jm,i]+U[jp,i]-4*U[j,i];
  double d2v = V[j,im]+V[j,ip]+V[jm,i]+V[jp,i]-4*V[j,i];
  double uv2 = U[j,i]*V[j,i]*V[j,i];
  Uu[j,i] = du*d2u-uv2+F+U[j,i]*one_minus_F;
  Vv[j,i] = dv*d2v+uv2+V[j,i]*one_minus_F_minus_k;
 }
}

On my machine this took on the order of 600ms to perform 100 iterations. Using jagged arrays instead, we have the following

for (int j = 0; j < H; j++)
{
 int jm = j != 0 ? j - 1 : H - 1;
 int jp = j != H - 1 ? j + 1 : 0;
 for (int i = 0; i < W; i++)
 {
  int im = i != 0 ? i - 1 : W - 1;
  int ip = i != W - 1 ? i + 1 : 0;
  double d2u = U[j][im]+U[j][ip]+U[jm][i]+U[jp][i]-4*U[j][i];
  double d2v = V[j][im]+V[j][ip]+V[jm][i]+V[jp][i]-4*V[j][i];
  double uv2 = U[j][i]*V[j][i]*V[j][i];
  Uu[j][i] = du*d2u-uv2+ F+U[j][i]*one_minus_F;
  Vv[j][i] = dv*d2v+uv2+V[j][i]*one_minus_F_minus_k;
 }
}

Just replacing the array type and making no other changes this now takes about 500ms, an improvement of about 13%. However, now that we are using jagged arrays we can do better and move the row accesses out of the inner loop.

for (int j = 0; j < H; j++)
{
 int jm = j != 0 ? j - 1 : H - 1;
 int jp = j != H - 1 ? j + 1 : 0;
 var u_m = U[jm];
 var u = U[j];
 var u_p = U[jp];
 var v_m = V[jm];
 var v = V[j];
 var v_p = V[jp];
 var uu = Uu[j];
 var vv = Vv[j];
 for (int i = 0; i < W; i++)
 {
  int im = i != 0 ? i - 1 : W - 1;
  int ip = i != W - 1 ? i + 1 : 0;
  double d2u = u[im]+u[ip]+u_m[i]+u_p[i]-4*u[i];
  double d2v = v[im]+v[ip]+v_m[i]+v_p[i]-4*v[i];
  double uv2 = u[i]*v[i] *v[i];
  uu[i] = du*d2u-uv2+F+u[i]*one_minus_F;
  vv[i] = dv*d2v+uv2+v[i]*one_minus_F_minus_k;
 }
}

This brings the count down to about 330ms, an improvement of about 40% over the original implementation.

We’ve taken the sequential code about as far as we can go. The next step is to try to take advantage of muliple cores that may be available on the user’s machine.