This was a 'challenge' assignment (at least that's how I perceived it) for us by our Fortran Simulation Lab professor. We were given a dam impounding 100m deep water. We had to solve the Laplace's equation in 2D using the 'implicit scheme' - basically solving a system of linear equations consisting of Laplace's finite difference equation for each unknown node. For the finite differencing domain, 4 grids along the horizontal and 3 grids along the vertical (each grid being a square of unspecified dimensions) were provided. We were also given Dirichlet boundary conditions on the nodes at the downstream end. When I first set out on manually solving the problem as a control for when actually developing code for it, I thought that the heel of the dam should have a head equal to the impounded water's depth, as all points upstream on the ground surface should - and still do think so. But that is how the problem was posed and so I ran with it.
The problem
I wanted to make the code 'general' in that it could work for other scenarios as well - not for just this specific combination of parameters. So, I decided to allow for 6 inputs: contourInterval, numDys, numDxs, DxsTillDam, reservoirHead, downstreamBoundaryHead. reservoirHead and contourInterval should be pretty obvious. numDxs and numDys are just that - total grids/boxes to work with along horizontal and vertical axes. DxsTillDam is the number of grid boxes to the left of dam's heel. downstreamBoundaryHead is the head at the downstream boundary. Based on these parameters, I had to think hard on what the coefficient matrix for the system of linear equations would look like. Since the matrix has a lot of structure in it, it was not that big of a problem. I just looked at the matrix I had prepared while working out the problem by pen and paper and coded in the logic to populate a 2D array with the exact same rule that the elements of the original problem's matrix followed. Another 2D array that functioned as a single column matrix i.e. vector representing the Dirichlet boundary conditions was created side by side while populating the coefficient matrix.
Matrix for reference and Manually solved grid for reference
Mazumder and Wen Shen were indispensable for understanding the matrix formulation.
The next step was to find the points on the grid lines that had the same head for each head. After some thinking, I settled on the idea of sweeping from the left boundary towards the right upto the right boundary(-1 to be precise), one column at a time, looking for a specific head to the immediate right and to the immediate above standing at each node moving up the column to top(-1 to be precise). Since the head being searched for rarely exist on the grid nodes themselves, I interpolated between the current node and the east node, and the current node and the north node to get the exact coordinates of the head - provided it existed between the nodes. I stored the interpolated coordinates off to two matrices. One matrix for the X coordinate, another for Y. In each matrix, each column corresponded to a head, viz. a column for 10m, another for 30m, etc.
Now that I had the equipotential coordinates/points for each head in question, all that was left was to join them. I first tried joining them in the order of their discovery i.e. according to their positions in the coordinate matrices. That didn't work very well as anyone might've guessed. Order matters when joining a bunch of equipotential points. You don't want jigjagging lines going back and forth for what should be a smooth isobaric curve. I then tried joining the points from the top down - basically sorting them in descending order of their Y coordinate. But that didn't work very well either. Since the isobars tend to be U-shaped, a lot of times, this method led to left-right-left jigjags. I then thought of joining the points from left to right. This worked for majority of the isobars. But it almost always messed up one isobar - right in the middle, when the curves shift from left facing to right facing. The isobars to the left and to the right are 'proper' functions i.e. unique Y for any given X. However, at the transition zone, between the two kinds - which happens to be right at the center of the base of dam(not here though, for reason described towards the top of this article) or more precisely at the ground level right between nodes of reservoirHead and 0 - the isobaric curve is somewhat vertical and hence exhibits 'ill' behavior as a function i.e. two or more Y's at some X. It is precisely due to the existence of such isobars that the left to right joining method fails towards the center of the grid. I finally ended up adopting a hybrid technique: joining from top to bottom but with shortest distance to the immediately next bottom point consideration. This involved sorting the equipotential points for a head in descending order of Y coordinate and then sorting by shortest distance between any two consecutive points. This gave better results but for some unknown reason, which may or may not be related to this technique, bottom most isobars starting making jump between left and right. I have not researched too much into this issue but it may just be how flownets are: loop-like and when a loop gets clipped due to our domain boundary, there's bound to be equipotential points to the right and to the left for a head with the bottom points clipped. In such a case, this 'jumpy' behavior is exactly what would be expected. Anyway, I just wrote a couple of anti-jump lines to prevent joining points farther apart than a specific X and Y distance.
My scribbles: 1 2 3
The viewport width and starting position(top left coordinates) are by default set to 600, 200, 8. Sometimes a larger viewport width may be desirable. The user may freely change these without fear of messing up the grid, irrespective of problem-centric parameters listed in the second paragraph.
Code
The problem
I wanted to make the code 'general' in that it could work for other scenarios as well - not for just this specific combination of parameters. So, I decided to allow for 6 inputs: contourInterval, numDys, numDxs, DxsTillDam, reservoirHead, downstreamBoundaryHead. reservoirHead and contourInterval should be pretty obvious. numDxs and numDys are just that - total grids/boxes to work with along horizontal and vertical axes. DxsTillDam is the number of grid boxes to the left of dam's heel. downstreamBoundaryHead is the head at the downstream boundary. Based on these parameters, I had to think hard on what the coefficient matrix for the system of linear equations would look like. Since the matrix has a lot of structure in it, it was not that big of a problem. I just looked at the matrix I had prepared while working out the problem by pen and paper and coded in the logic to populate a 2D array with the exact same rule that the elements of the original problem's matrix followed. Another 2D array that functioned as a single column matrix i.e. vector representing the Dirichlet boundary conditions was created side by side while populating the coefficient matrix.
Matrix for reference and Manually solved grid for reference
Mazumder and Wen Shen were indispensable for understanding the matrix formulation.
The next step was to find the points on the grid lines that had the same head for each head. After some thinking, I settled on the idea of sweeping from the left boundary towards the right upto the right boundary(-1 to be precise), one column at a time, looking for a specific head to the immediate right and to the immediate above standing at each node moving up the column to top(-1 to be precise). Since the head being searched for rarely exist on the grid nodes themselves, I interpolated between the current node and the east node, and the current node and the north node to get the exact coordinates of the head - provided it existed between the nodes. I stored the interpolated coordinates off to two matrices. One matrix for the X coordinate, another for Y. In each matrix, each column corresponded to a head, viz. a column for 10m, another for 30m, etc.
Now that I had the equipotential coordinates/points for each head in question, all that was left was to join them. I first tried joining them in the order of their discovery i.e. according to their positions in the coordinate matrices. That didn't work very well as anyone might've guessed. Order matters when joining a bunch of equipotential points. You don't want jigjagging lines going back and forth for what should be a smooth isobaric curve. I then tried joining the points from the top down - basically sorting them in descending order of their Y coordinate. But that didn't work very well either. Since the isobars tend to be U-shaped, a lot of times, this method led to left-right-left jigjags. I then thought of joining the points from left to right. This worked for majority of the isobars. But it almost always messed up one isobar - right in the middle, when the curves shift from left facing to right facing. The isobars to the left and to the right are 'proper' functions i.e. unique Y for any given X. However, at the transition zone, between the two kinds - which happens to be right at the center of the base of dam(not here though, for reason described towards the top of this article) or more precisely at the ground level right between nodes of reservoirHead and 0 - the isobaric curve is somewhat vertical and hence exhibits 'ill' behavior as a function i.e. two or more Y's at some X. It is precisely due to the existence of such isobars that the left to right joining method fails towards the center of the grid. I finally ended up adopting a hybrid technique: joining from top to bottom but with shortest distance to the immediately next bottom point consideration. This involved sorting the equipotential points for a head in descending order of Y coordinate and then sorting by shortest distance between any two consecutive points. This gave better results but for some unknown reason, which may or may not be related to this technique, bottom most isobars starting making jump between left and right. I have not researched too much into this issue but it may just be how flownets are: loop-like and when a loop gets clipped due to our domain boundary, there's bound to be equipotential points to the right and to the left for a head with the bottom points clipped. In such a case, this 'jumpy' behavior is exactly what would be expected. Anyway, I just wrote a couple of anti-jump lines to prevent joining points farther apart than a specific X and Y distance.
My scribbles: 1 2 3
The viewport width and starting position(top left coordinates) are by default set to 600, 200, 8. Sometimes a larger viewport width may be desirable. The user may freely change these without fear of messing up the grid, irrespective of problem-centric parameters listed in the second paragraph.
Code
 |
| Flownet with reservoirHead=100, numDxs=4, numDys=3, contourInterval=10, downstreamBoundaryHead=50 and DxsTillDam=2 |
 |
| Flownet with reservoirHead=100, numDxs=20, numDys=20, contourInterval=10, downstreamBoundaryHead=50 and DxsTillDam=4 |
No comments:
Post a Comment