A big problem with FDS is that many times there is much more calculation than we actually need. It is easer to understand this if we imagine that we are heating a long metal wire in the middle. It will be very hot there and the heat will slowly move to the ends.
The heat distribution near the fire will be changing quickly, while it will be almost constant at the ends. It is obvious that we need more precise calculations near the fire than at the ends.
This is a situation that is quite common in many problems that we simulate with PDEs. Because of this characteristic, the idea of adding an adaptive grid (using different h at different parts of the problem space) to FDS appeared quickly.
This is when it starts to show its real power - due to the large increase in the speed of calculation. There are many articles written about this subject, but most of them are dedicated to mathematical aspects.
The implementation of this kind of grid in a computer program can turnout to be a much bigger problem than expected. When using an adaptive grid there are two important steps that need to be done.
The first is to create a grid that has high resolution at the correct places and that allows easy transition between indexes of high and low resolution areas. The second is to correct the corresponding equations. One method of solving these problems is explained in detail in this article.
With the math out of the way, we can say a couple things about the software implementation. We shall explain it with some pseudo code, and some comments. As expected u will be kept in an array, but we will also need one more to keep a copy of it in the previous step.
double u[MAX], uprev[MAX]
Now we wish to translate equations (7,8) into some code. We will create a function to calculate Diff first, for all the elements in the array.
void Diff(double u,double
diff  = diff [MAX-1]=0;
for(int i =1; i< MAX -1;i++)
diff [i]= K*dt*a*( u[i+1] - 2*u[i] + u[i])/(h*h)
At this point we can say that a is a parameter that holds some physical characteristics and it can be set to 1. Notice that the first and the last elements in the array are set to 0. This is the border and things are a bit different there. Some explanations will follow in a bit. Using Diff we implement the approximation of Runge-Kutta:
void Calculate (double u,double uPrev)
double S1[Max], S2[Max], S3[Max], S4[Max];
uTemp = uPrev;
Diff(uTemp, S1, 1);
uTemp = uPrev + 0.5*S1;
Diff(uTemp, S2, 0.5);
uTemp = uPrev + 0.5*S2;
Diff(uTemp, S3, 0.5);
uTemp = uPrev + 1*S3;
Diff(uTemp, S4, 1);
u = uprev + (S1 + 2*S2 + 2*S3 + S4)/6;
This function looks pretty much as expected, besides FixBorder appearing out of the blue. This is one more problem with FDS - the border. It is obvious that the array has to have some beginning and some end, but at these points we can't calculate Diff. How we solve this problem is very important, because in a number of time steps this will affect all the elements in the array.
Actually what we do at the border can have many different effects on our function which corresponds to different physical models. A lot of research has been done in math about this subject and many different approaches have been used. In our example we shall use a simple one, that isn't very mathematically strict. We shall just keep "going in the same direction".
u = u + (u-u);
u[Max-1] = u[Max-2] + (u[Max-2]-u[Max-3]);
Finally we have the calculation loop.
uprev = u;
The use of the scheme has been put in a class, with which we wish to signify that the use of FDS is practically the same for all possible problems. The only difference would be Diff that is used in the implementation of Runge-Kutha time step.
If an adaptive grid was used there would be some extra lines. First some that would create the grid before the main loop, and some that would recalculate the grid every m step.
There are many processes in nature that have been translated into PDEs, and with this method they could be added to games. This was a presentation of an improved FDS and its implementation, for a very simple problem.
The use of it for more complex problems would in its basics have the same steps. They are finding the correct equations for the problem, creating the scheme, creating the adaptive grid and implementing the functions.
As an example of a problem in 2+1 dimensions you can read this article. This is a powerful method, but solving PDEs is a very hard problem in many cases, and this method might not always give very accurate and stable results.
There are other methods like implicit FDS, finite elements, and the use of Fourier Transform that may be able to solve them, but in some cases it is still best just to use tricks.