As such, I decided I would try to model one of the simplest thermodynamic systems: a hot plate in a room. The problem is incredibly simple: There is a room of finite size and low temperature. At one end of the room is a hot plate that has constant high temperature. The idea is to model the flow of heat throughout the system. For this example, we can assume the room is 2-dimensional.
The most straightforward thing to do is obviously to divide up the room into a grid of small squares. If these squares are small, it can be assumed that the temperature across the area of the square is constant, so each square can have a single temperature associated with it. The problem is then reduced to trying to figure out how the temperatures of these little squares changes over time.
The broad thermodynamic solution is fairly simple. The idea is that the rate of change of a particular point is proportional to the temperature difference between that particular point and the neighboring points. This makes intuitive sense – something really really cold heats up faster than something room temperature. As we are dealing with 2-dimensional squares, the neighboring squares can be the ones above, below, and on either side of a given square. Finding a single “surrounding” temperature can be done by simply taking the average of these four adjacent squares.
The formula can then be written dT/dt = alpha*deltaT. deltaT can therefore be written as (Tup + Tdown + Tright + Tleft)/4 – Tlocal. (little t is time, big T is temperature). This is a differential equation, and may or may not be able to be solved with analytical methods. I have in front of me, however, a computer, so perhaps I could use it to try to solve this problem. In a computer simulation, time can be simulated by computation time, so each dt can be approximated by each iteration of an outer loop. Performing the calculation on the right side creates the rate of change of the temperature. This must be multiplied by that dt to create the temperature delta that needs to be updated. If we assume that dt is a constant for each loop iteration, we can assume that it has the value of 1, because using any other value of dt will simply scale the time dimension of the final result. Therefore, dividing by one can be neglected, and the equation can simply be reduced to temperature += alpha*deltaT. This calculation must be performed for each point in the space, or for each little square.
Because this calculation must be performed for each of the little squares, it can be done in parallel (as there is no dependency on the squares). There is a trick, however. The temperature of neighboring squares is used in the calculation of the temperature of the current square. This means that the output of the calculation is the same type and semantics as the input. It would be easy to straightforwardly update the current temperature inside a thread. However, the neighboring squares may or may not have been calculated. If they have been calculated already, they would have used the old value for the current temperature. If they have yet to be calculated, they will use the new value for the current temperature, because you just overwrote the old one. The way to get around this is to use two buffers, an input buffer and an output buffer (or a read buffer and a write buffer; call them what you like). Reading in the values from the read buffer and writing to another buffer means that there won’t be any inconsistency in the system. After doing the computation for each of the little squares, simply swap the buffers for the next iteration.
I decided to write the program in Haskell, because it is my favorite language. The program requires an iteration over all the little squares in the requisite area, as discussed before. For each one of those little squares, the neighbors need to be accessed to find the surrounding temperature of the current square. Haskell’s lists are linked lists, however, so lookup is not O(1). The iteration could be done by using 5 iterators: one for the current square and four for the surrounding squares. This would be rather complicated to do, so I elected instead to use Haskell’s Vector library. Vectors have O(1) lookup time, so I could simply just query for the neighboring pixels.
Coding up the iteration was fairly simple. More difficult, though, was displaying the result. Viewing the temperatures as text is pretty visually unappealing, so I decided I would try again to use libGD. This turned out to be another headache (see this for previous headaches). Since that previous post, a Snow Leopard release of GHC has been released in the Ports library, so I have migrated to that build. I used Cabal to install Graphics.GD without any hitch. However, the library didn’t appear to work very well. I tested out the library in GHCi by creating an image, writing to it, and saving it to disk. However, GHCi segfaulted whenever an image fell out of scope. I’m assuming that this has to do with a bug in the image closing code that the garbage collector would invoke. I found that the code would work just fine if I used a text file and invoked runghc on it, however, so I decided to just do that instead of use the interactive mode. However, only about half of the drawing commands appeared to be showing up with this solution. Out of desperation, I tried to compile the program into an executable with ghc –make. This appeared to solve the problem completely, and all the drawing commands appeared in the resulting image. This solution works, but it is a bit of a hack. Oh well – it appears the GHC maintainers aren’t so concerned with certain specific libraries.
After running the program, I observed something quite interesting. It has to do with the walls of the area I am interested in. Finding the neighboring squares is pretty straightforward when you’re on the border of the arena: simply disregard the pixels that are beyond the boundary. This means that the denominator in the average calculation found above isn’t always a “4.” It is the number of neighboring squares that are within the bounds of the area in question. This led to a very straightforward looking heat transfer from one side of the area to the other. More interestingly, I tried hardcoding this denominator to 4 to see what would happen. Functionally, this means that all four surrounding squares are inspected, but the ones outside the area in question have a temperature of 0. This actually led to a very interesting looking heat transfer pattern – a sort of rounding of the border away from the edges of the area.
I also generated a couple thousand frames and found that after a while, the temperature doesn’t appear to be changing much if the area in question is large enough. This makes me wonder about convergence of the equation. It is true that the smaller the temperature difference, the slower the temperature change. This means that the temperature changes less and less as time goes on. However, after generating a few thousand frames, the temperature doesn’t seem to be changing at all. Intuitively, I know that the stable point of the system is when the entire area is the same temperature as the hot plate. I wonder if this is the kind of system where the analog, mathematical formula does not converge (outside of the stable point) but the digital, binary approximation does converge to something that is not the stable point. Food for thought.
I’m hoping to code this same program again in C and in CUDA to compare runtimes. CUDA should be an interesting build, as most of this problem is not computation time but communication time. Communicating between blocks in CUDA is an interesting problem, and I’m curious as to how to do it in the most optimal way. More on that later!
My Haskell code can be viewed here.
No comments:
Post a Comment