To explain a bit of this: this program calculates fluid flow around a cylinder. In the original version, it was doing only symmetric-flow cases, so it only solved the equations for half of the flow, using a mathematical "mirror" through the cylinder center to do the other half. However, I'm going to want to use it for some non-symmetric cases, and so it needed to be altered to calculate the full flow rather than half of it.
Now, the bug was that things were going unstable (pressures trying to get to negative infinity very quickly) at one point that happened to be next to the mirror. So, after poking at it and not having any luck, I decided to see if it was something with the mirror, and put in the code to make it go full-circle without a mirror. And the problem disappeared.... Then I went home.
Today started with a bit of a digression. Another thing I've been needing to do is make the registration brochure for a conference that my advisor is hosting -- I've got earlier versions of that, so it's just an update with new material, but it still ended up taking all afternoon because it needed to be completely rearranged to fit on a 4-fold rather than a 3-fold brochure layout to include all the info.
So, after supper, there was more debugging. I had left the program in an almost-running state, except that it was producing asymmetric results for symmetric input. This isn't one of the situations where that's supposed to happen, so clearly there was a bug somewhere.
Now, this program has in it a bit of cleverness called "multigrid", which is probably worth describing to make things easier to explain....
First, the reason we need cleverness: We want to solve this problem on a fairly fine grid of datapoints, so that the solution we get will accurately approximate the true solution (which is effectively on a continuum, rather than discrete points). However, the amount of work that it takes to get a solution goes as N squared, where N is the number of points, so for "sufficiently fine" grids it amounts to several hours for a fairly simple problem.
The cleverness, then, is as follows: we take the problem as posed on a fine grid, and make a coarse-grid approximation of it. Then, we move towards the solution on the coarse grid, and interpolate those results back to the fine grid. This goes much faster than just moving on the fine grid, but because we keep going back to the fine grid (and doing a step towards the solution there each time we return to it), we get the fine-grid solution. But we can do even better -- we approximate the coarse grid with an even coarser grid, and so forth up to something that's four points wide. The end result of all this is that it takes much, much less time to get a solution -- a minute, instead of several hours, at least in theory.
So, anyhow, the asymmetry was only showing up when I went up a level to a coarser grid. Turns out that it was basically an off-by-one error on moving back to the fine grid, so that it effectively twisted the solution slightly before putting it back down.
I fixed that, and one level worked. So, I tried two levels of coarsening, and the asymmetry was back with a vengance.
After much poking, I discovered that one of the arrays was somehow changing values between when I set it in a subroutine, and when that subroutine returned. To make a long story short -- what the subroutine called "vrfr" was not what the main program called "vrfr" -- I'd somehow mixed up some variable-passing statements, and so this particular subroutine was getting its own vrfr to play with that no other part of the program ever saw. Argh.
So it seems to be working. And still printing out large piles of debug information, because I haven't stripped that out yet....
Next up: While it may be working, it's not really solving the right problem -- it's solving for flow around a cylinder with extra viscosity added to make things behave better when it's not close to the final solution. I need to make it not need that; probably what will happen is having that in at the start of each run, and then slowly taking it out once things settle down to something close to a solution.
And then I need to try some more interesting sorts of flow coming in and see if I can get data, and then I need to generate some piles of data.
Then, phase two of this project....