DataWindow uses Incorrect Coordinate Ordering

[TedStudley]

DataWindow currently uses matrix-style coordinate ordering of (rows, columns), while the standard finite-difference coordinate ordering is (columns, rows).

Since DataWindow is supposed to be a helper class to make dealing with the underlying array representation easier, it makes very little sense that it use the wrong coordinate ordering. Refactor the class to use the correct coordinate ordering.

[icherkashin]

"DataWindow currently uses matrix-style coordinate ordering of (rows, columns), while the standard finite-difference coordinate ordering is (columns, rows)."

The conventional ordering of indices is (i,j), where i is i-th row, and j is j-th column. The dataWindow class, however, uses (j,i) ordering, which has been a source of great confusion. Considering what you've said, are we at a contradiction here?

[egpuckett]

Ivan is correct. The (i,j) cell has it's lower left hand corner at (i_h, j_h), say. The precise mapping from "index space" (i,j) to 'physical space' (x,y) location is part of the design. However, the first index should be 'i' here and i_h is a location in the x-dimension / on the x-axis, j_h is a location in the y-dimension / on the y-axis. So

uVelocity(i,j)

is the velocity u(x,y) in the x-direction near

(x,y) = (i_h, j_h)

The precise location depends on whether u is edge centered or cell centered. Currently Tau's algorithm for solving Stokes equations works with u centered on left / right edges of the (i,j) cell, so the physical location would be

(x,y) =( i_h, (j+1/2)_h)

ASSUMING the lower left hand corner of the (i,j) cell is at

(x,y) =( i_h, j_h)

which is what I propose.

I believe Harsha has switched the indices in DataWindow and has a running code but it is not tested on the tauBenchmark yet, since we don't have a working convergence test. A working set of convergence tests is very high on my list of priorities, Tau being the first.

One other thing, we've been discussing using boost's multi array, since

1. It is implemented to be fast
2. I believe we can turn bounds checking on / off. (Is this correct.)
3. I believe we can use negative indices so we can have a mapping from index space to X = -h/2 for example

[egpuckett]

We need to include Harsha in this discussion and make no changes to the indexing convention until the four of us have carefully reviewed them. There still seems to be some confusion. I recommend not thinking of the finite difference grid in terms of rows and columns like a matrix, but as a grid that overlies physical space. There may / are issues that must be treated carefully when taking the velocities, etc, and creating a vector U that includes those velocities, pressure, etc. and a matrix A for solving the equation

A U = f

This needs to be handled carefully if we change the index convention.

[TedStudley]

Perfect example of why this is even a problem in the first place. I've always been terrible about remembering the proper ordering for i/j, much less whether (0,0) is in the top/bottom/left/right corner. I think when I originally wrote Data window because Gerry originally told me that the ordering was wrong, I probably 'fixed' the wrong thing.

I think before going any further on this, we should sufficiently document (with diagrams!) what the i/j to x/y relationship is. This can be done with Doxygen after the new build system is merged into master.

[TedStudley]

Moving to using boost::multi_array could be useful, although it won't really provide a significant improvement in the areas you've mentioned.

DataWindow is already likely either as fast or faster than multi_array, being just a simple wrapper on top of the raw memory array. If multi_array does some extra fancy memory tricks to make itself faster (which I doubt it does), it would render itself unusable for our purposes, since we need to keep the memory contiguous in order to be interoperable with Eigen.

We have bounds-checking in DataWindow, and using them in multi_array would be likely to invalidate the third reason for using it.

Being able to use the negative indices to reference the boundary values is doable, but would significantly complicate (and possibly slow down) working with Eigen due to the necessity to use inner and outer strides with Eigen::Map. Working with boundary values, at least, would likely introduce more cache misses since far fewer actual boundary values would fit into cache.

I'll look into a solution for consolidating the interior and boundary values in a single object. This would have the added benefit that we can switch out the implementation of the boundary values far more easily, for example not duplicating the boundary value at every cell on the boundary when using constant single-value boundaries. It's unlikely that there will be a bare-metal efficient way to use negative indices to reference boundary values, since we would either have to use a conditional (adding a few extra operations per access) or use a single contiguous memory region (which would cause cache misses and make calls to Eigen more unweildly). Once again though, efficiency concerns are likely unwarranted here, since the cost of calls to the Eigen solver drastically outweigh pretty much any efficiency gains/losses in DataWindow.

[icherkashin]

DataWindow is already likely either as fast or faster than multi_array, being just a simple wrapper on top of the raw memory array. If multi_array does some extra fancy memory tricks to make itself faster (which I doubt it does), it would render itself unusable for our purposes,

I don't know the source that claimed it would be faster: it is Dr. Puckett's interpretation of my words, perhaps. In any case, its performance will suffice for our purposes.

since we need to keep the memory contiguous in order to be interoperable with Eigen.

_Boost::multiarray is simply, in Stroustrup's terminology, a resource handle - it provides for effective access and management of data without modifying it. The contiguity and storage order of the data is the property of the data itself, not of the handle thereto: _multiarray will preserve the column-major ordering of the array, as well as its contiguity.

We have bounds-checking in DataWindow, and using them in multi_array would be likely to invalidate the third reason for using it.

_Multiarray has indices range checking by default and they can be turned off. Why would it invalidate the third reason - it has no relation to using negative indices?

Being able to use the negative indices to reference the boundary values is doable, but would significantly complicate (and possibly slow down) working with Eigen due to the necessity to use inner and outer strides with Eigen::Map. Working with boundary values, at least, would likely introduce more cache misses since far fewer actual boundary values would fit into cache.

This is completely untrue: Eigen is given only the data that it is needs to be supplied with. For example, we will simply supply an array of velocities, with or without boundary values, depending on whether they are needed in the matrix - this is what is done currently in the process of solving Stokes equations. In other words, which resource handle to the data is used is irrelevant to Eigen - what matters is the actual data that it is given. We shouldn't be solving a problem that does not actually exist. Furthermore, negative indices are convenient in the advection algorithms, but this part of the code is not even working properly at the moment: the Stokes part will essentially remained untouched by the transition to _multiarray.

I'll look into a solution for consolidating the interior and boundary values in a single object. This would have the added benefit that we can switch out the implementation of the boundary values far more easily, for example not duplicating the boundary value at every cell on the boundary when using constant single-value boundaries. It's unlikely that there will be a bare-metal efficient way to use negative indices to reference boundary values, since we would either have to use a conditional (adding a few extra operations per access) or use a single contiguous memory region (which would cause cache misses and make calls to Eigen more unweildly). Once again though, efficiency concerns are likely unwarranted here, since the cost of calls to the Eigen solver drastically outweigh pretty much any efficiency gains/losses in DataWindow.

I think you are trying to address some abstract problem of implementing negative indices, as if it is a fundamental capability that is needed everywhere in the code. In fact, negative indices are only necessary (or convenient, to say the least) in the thermal advection part of the code. There are no implications of using negative indices that go beyond the dataWindow class: using _multiarray will eliminate the necessity to even modify this class since _multiarray is already capable of using arbitrary index bases in the arrays. Perhaps, if you have an example of using negative indices where the concerns you have listed are relevant, you should present it.

I don't wish to sound adverse, but I am convinced that the issues I discussed above would be false targets for this project that would simply waste our time and effort. Moreover, based on my analysis of the code and its requirements, it wasn't even necessary to create the dataWindow - it simply duplicates existing capabilities of libraries like _multiarray. Unfortunately, it took a long time to come to this conclusion. We certainly do not want to repeat the same mistakes by addressing problems that do not really exist.

I'm looking forward to hearing your response to my claims, especially from Professor Puckett who de facto establishes the requirements for this project. What is your third-party perspective on these issues?

[icherkashin]

I think before going any further on this, we should sufficiently document (with diagrams!) what the i/j to x/y relationship is. This can be done with Doxygen after the new build system is merged into master.

Please see my email with the subject Example grid illustrating indexing convention. Unfortunately, I could not attach a Postscript file here, and I could not convert it successfully to formats accepted here.

[TedStudley]

It's completely possible that I'm missing something about the implementation of multi_array that invalidates my objections, but so far you haven't really told me anything that I don't already know. It seems as though you may have misinterpreted what I'm actually objecting to in several cases, notably the necessity to keep the interior and boundary values each separated and individually contiguous for performance reasons. It's possible this can be done while still allowing multi_array to efficiently and concisely represent them as though they were contiguous, although I haven't seen this yet.

From my 'third-party' perspective, when I originally authored the code, DataWindow was necessary to address a concern that Gerry had about proper index/coordinate ordering. It was improperly implemented due to an unfortunate miscommunication about exactly what ordering was incorrect at the time. At no point in the process of implementing the wrapper class did the requirement of consolidating the interior and boundary values in a single consistently-addressable object ever arise. Using multi_array for this purpose was wholly unnecessary, since all that was needed to fulfill the requirements was a simple wrapper class around a pointer with a custom 'indexing' operator.

Why don't you implement a proof-of-concept using multi_array. If my concerns are misfounded, then It'll be a lot easier to see. If they aren't, then it'll be a lot easier to communicate them with a concrete example.

[icherkashin]

From my 'third-party' perspective, when I originally authored the code, DataWindow was necessary to address a concern that Gerry had about proper index/coordinate ordering.

By "third-party" I meant Professor Puckett, in the sense that he could consider both of our opinions from an independent perspective. I hope you took no offence at this - perhaps it wasn't clear whom I was addressing.

At no point in the process of implementing the wrapper class did the requirement of consolidating the interior and boundary values in a single consistently-addressable object ever arise.

It is interesting, I didn't realize this was your point. Neither was this what I was talking about or arguing for.

Considering this,

Why don't you implement a proof-of-concept using multi_array. If my concerns are misfounded, then It'll be a lot easier to see. If they aren't, then it'll be a lot easier to communicate them with a concrete example.

I agree that this would be the best approach. We've simply lost track of our discussion due to the lack of concrete examples.

Still, to summarize my points:

• _multiarray would provide capabilities equivalent to dataWindow's. In other words, it would simply replace dataWindow in the existing working parts of the code without any functional changes to other parts of the code. For example,

uVelocityWindow(j,i)

would simply become

uVelocity_multi_array(i,j)

with everything else remaining the same.

• _multiarray would provide extra capabilities in those parts of the code that are still in progress. For example, in the advection part, it would allow us to use negative indices to implement ghost cells in case they are necessary. We would still be able to choose which array, with or without the boundary data, to pass to Eigen if it becomes necessary.

I hope this summary above communicates my understanding of the problem in general. Of course, other issues may arise which I may not see at this point but you may. Once I create an example of replacing dataWindow with _multiarray in some part of the code, let's take a look at it together!

[icherkashin]

What do you think of changing the interface to dataWindow

operator(unsigned int col, unsigned int row)

to

operator(unsigned int row, unsigned int col)

i.e. from (j,i) to (i,j)?

I've done it before but we've decided not to implement it at that moment, but now Harsha, I've heard, has also done it. This will not change the underlying column-major storage ordering yet will allow us to use conventional indexing. The only issue I see is to carefully reflect this change in the interface everywhere where it is used in the code.

Using _multiarray is something I consider for near future but, certainly, not for now. I still think it is useful to discuss it.

[TedStudley]

I'm sorry; I incorrectly parsed the comment as though it were directed towards myself, and it would be difficult not to take offense, given the situation. I apologise for the misunderstanding.

I believe that the plan for fixing the argument order of DataWindow::operator() is to merge in @Vandemar's changes, which should have fixed the interface as well as all usages of the DataWindow class in the code. Whenever they're ready, just open a pull request and somebody can review them (probably not me, seeing as my tendency to get the ordering mixed up is the reason we're making this pull request to begin with...) I'll assign Harsha to the issue and he can make a pull request to close it whenever the changes are ready.

[icherkashin]

I'm sorry; I incorrectly parsed the comment as though it were directed towards myself, and it would be difficult not to take offense, given the situation. I apologise for the misunderstanding.

Of course, Ted, everything is fine, you shouldn't apologize! Such situations are inevitable, and I'm glad we are able to resolve problems like this!

I believe that the plan for fixing the argument order of DataWindow::operator() is to merge in @Vandemar's changes, which should have fixed the interface as well as all usages of the DataWindow class in the code. Whenever they're ready, just open a pull request and somebody can review them (probably not me, seeing as my tendency to get the ordering mixed up is the reason we're making this pull request to begin with...) I'll assign Harsha to the issue and he can make a pull request to close it whenever the changes are ready.

I don't see the interface changed there. However, there is some other class that Harsha has added https://github.com/Vandemar/mc-mini/commit/8ea051b9e2553ade757c98a842396639f520078c. I don't understand Harsha's changes - he should explain what he did.

[TedStudley]

I think that the swap is in this branch, although I haven't had a chance to look it over yet.

[hlokavarapu]

The code modifications related to the domain wrapper class was meant as a means of temporarily swapping the indicies and to eventually port over to the boost multi_array library. For now though, I am going to remove that particular branch in the near future.

[icherkashin]

An indexing issue arises when estimating convergence for the Tau problem with the known exact solution. When the conventional (i.e. ij) indexing is used in the exact solution functions, there is no convergence. However, when the indices (or, to be even more precise here, the coordinates) are swapped, everything converges with the second order accuracy.

# Exact solution for uVelocity restricted to the cell-centers of the grid by interpolation (averaging).
 def v_x_exact_interpolated(i,j):
   # NOTE: This is the correct index ordering, i.e. "ij" <=> "xy"
  #return 0.5 * ( cos(i*h_x)*sin( (j+0.5)*h_y) + cos( (i+1)*h_x )*sin( (j+0.5)*h_y))
  # NOTE: This is the incorrect index ordering, i.e. "ji" <=> "xy"
  # Convergence takes place for this (incorrect) index ordering.
  return 0.5 * ( cos((j+0.5)*h_y)*sin( i*h_x) + cos( (j+0.5)*h_y )*sin((i+1)*h_x ))

# Exact solution for vVelocity restricted to the cell-centers of the grid by interpolation (averaging).
 def v_y_exact_interpolated(i,j):
  # NOTE: This is the correct index ordering, i.e. "ij" <=> "xy"
  #return 0.5 * (-sin( (i+0.5)*h_x )*cos(j*h_y) - sin( (i+0.5)*h_x )*cos( (j+1)*h_y))
  # NOTE: This is the incorrect index ordering, i.e. "ji" <=> "xy"
  # Convergence takes place for this (incorrect) index ordering.
  return 0.5 * (-sin( j*h_y )*cos((i+0.5)*h_x) - sin( (j+1)*h_y )*cos( (i+0.5)*h_x)) 

Possible source of this problem: incorrect indexing in the initialization of the velocity boundary values - https://github.com/TedStudley/mc-mini/blob/master/source/problem/initialization.cpp#L151-L157:

  if (boundaryModel == "tauBenchmark") {
    for (int i = 0; i < M; ++i)
      for (int j = 0; j < 2; ++j)
        uVelocityBoundaryWindow (j, i) = cos (j * N * h) * sin ((i + 0.5) * h);
    for (int i = 0; i < 2; ++i)
      for (int j = 0; j < N; ++j)
        vVelocityBoundaryWindow (j, i) = -sin ((j + 0.5) * h) * cos (i * M * h);

[egpuckett]

One suggestion: Please double check to be certain that the ordering is correct in the exact solution. I know that it is very likely you have the exact solution correct, but a second pair of eyes, say Harsha's, may spot something.

[egpuckett]

Also, would you please put your branch of mc-mini in which you obtain convergence with the indices switched on your GitHub account. Please call it something obvious like 'tauBenchmark-devel'.

[hlokavarapu]

if (boundaryModel == "tauBenchmark") { for (int i = 0; i < M; ++i) for (int j = 0; j < 2; ++j) uVelocityBoundaryWindow (j, i) = cos (j * N * h) * sin ((i + 0.5) * h); for (int i = 0; i < 2; ++i) for (int j = 0; j < N; ++j) vVelocityBoundaryWindow (j, i) = -sin ((j + 0.5) * h) * cos (i * M * h);

Comparing u and v velocity equations with the Tau paper in the domain and as well as the boundaries, I find no discrepancies.