[siamak]

Hi There!
I tried solving pressure correction equations using python and c++, I implemented the same equation but I got totally different answer after some iterations. I have checked for several inputs.
Python:

for k in range(1, press_iter+1):
    for i in range(1, GRID+1):
        for j in range(1, GRID+1):
            pp[i, j] = pp[i, j] + (1.7/aPp[i, j])*(aE[i, j]*pp[i+1, j] + aW[i, j]*pp[i-1, j] + aN[i, j]*pp[i, j+1] + aS[i, j]*pp[i, j-1] - source[i, j] - pp[i, j]*aPp[i, j])

C++:

for(k=1;k<=press_iter;k++){
            for(j=1;j<=grid;j++){
                for(i=1;i<=grid;i++){
                    pp[i][j] = pp[i][j] + (1.7/app[i][j])*(ae[i][j]*pp[i+1][j] + aw[i][j]*pp[i-1][j] + an[i][j]*pp[i][j+1] + as[i][j]*pp[i][j-1] - source[i][j] - pp[i][j]*app[i][j]);
                }
            }
        }

Can anybody help why it happens?

[FMDenaro]

But the difference is in terms of number of iterations ? I see the "for" on i,j differently coded, that can affect the performance (not the converged solution).

[sbaffini]

While things between the two codes could differ for reasons well beyond the coding itself (like, for example, which precision you used in C++ for all the variables, vs the one used in python, just to mention the most obvious possibility), let us assume, for the sake of it, that it is something in the code (but we don't have them all to make a call).

We can then probably approach this like it is usually done in refactoring. You have two different pieces of code that, however, after each step are meant to produce the same exact output (note, this is at the bit level in actual code refactoring) if provided with the exact same input. Then the way to proceed is to check in the codes, step by step, the input and the output (e.g., print them in a file and make a diff/meld) until you find the place where the two differ.

Let's take a step further and imagine that, just before the loops you posted, the two codes are still the same in the value of each variable. Now, let me first say that your ranges for the integers in the two codes are different (EDIT: I actually later realized this is not the case two posts below), so I don't even proceed in reading the supposedly identical line, but I'll tell you what you should have done: just like before, printing to file for every loop until you (actually meld or diff) spot the difference.

[sbaffini]

Originally Posted by FMDenaro

But the difference is in terms of number of iterations ? I see the "for" on i,j differently nested, that can affect the performance (not the converged solution).

Well, at the bit level they produce different solutions, something that an unstable or overrelaxed algorithm would be susceptible to. But my call is on the wrong C++ index ranges, which do not seem to include the grid+1 and press_iter+1 values (EDIT: And it is a wrong call because, as I realize below, python range() doesn't include the last element).

[sbaffini]

My bad, my bad, apparently range() in python doesn't include the last element, which I didn't know. So, the Filippo answer is the most relevant here.

[FMDenaro]

Originally Posted by sbaffini

My bad, my bad, apparently range() in python doesn't include the last element, which I didn't know. So, the Filippo answer is the most relevant here.

I had the same doubt, let see if the range is the same …

[sbaffini]

Originally Posted by FMDenaro

I had the same doubt, let see if the range is the same …

I'm bad at both C++ and python, but actually tested the C++ loop to be sure before answering. Yet, I kind of assumed that a range() function would work more reasonably.

for i in [1, 2, 3, 4]

just works as expected, but:

for i in range(1,4)

stops at 3, which I didn't know.

[FMDenaro]

What is more, I am not sure how the array is stored. In case of Fortran the difference in the for cycles would affect the cache access, what happens in this case ?

[sbaffini]

I'm clearly not the one that can answer this question anymore but, the C++ code appears to loop in a column major fashion (in contrast to the C++ row major memory layout), while the python one does the opposite (but I have no idea which memory layout is used in python).

If the two codes have identical values just before the loop then it might be a matter of stability related to the order of the loops.

[siamak]

Originally Posted by FMDenaro

But the difference is in terms of number of iterations ? I see the "for" on i,j differently coded, that can affect the performance (not the converged solution).

thanks in advance
During testing the implemented code in python I switched the order of i and j and I got a different answer. As Other mentioned I even checked the type of the variable but the answer is still different. I have used the same different inputs and got different answers respectively for all inputs. It seems the code snippets produce a different result.

[siamak]

Originally Posted by sbaffini

I'm bad at both C++ and python, but actually tested the C++ loop to be sure before answering. Yet, I kind of assumed that a range() function would work more reasonably.

for i in [1, 2, 3, 4]

just works as expected, but:

for i in range(1,4)

stops at 3, which I didn't know.

I am sure about the ranges and the number of the iteration and the inputs. I really wonder why I am getting a different answer.

[FMDenaro]

Originally Posted by siamak

I am sure about the ranges and the number of the iteration and the inputs. I really wonder why I am getting a different answer.

Again, what is this difference?

[sbaffini]

Originally Posted by siamak

thanks in advance
During testing the implemented code in python I switched the order of i and j and I got a different answer. As Other mentioned I even checked the type of the variable but the answer is still different. I have used the same different inputs and got different answers respectively for all inputs. It seems the code snippets produce a different result.

As my mistake might have clarified, the importan thing here is to not assume anything about anything so, have you checked that, before that loop, all the variables are indeed the same? We left unspecified what we mean by "the same" but, as writing them in the same binary form might be complicated, I guess, when working from different codes and languages, we can maybe assume that you compare them by writing them in ASCII using the underlying precision of the variables. This is the first mandatory check before actually going into that piece of code you posted. Note that "all the variables" here really means all of them, including "grid" and "press_iter", the initial values in pp and, obviously, the matrix coefficients (apparently "1.7" defaults to double precision in both python and C++, but maybe check that as well).

If the previous check passes, you have indeed verified that using a different i-j order changes the result (I guess, within the same criteria you have used to say that the two previous pieces of code produce different results). Is this the case both in python and C++? If it is for both, then you first need to put both codes in a form that does the same loop. That is, given the same input, the two codes will process it in the same order. At this point, I'm not going to tell if this is so or not in your codes, but you can check it.

If they have the same input and process it in the same order but still give you different outputs, then I suggest you first try different initial inputs, like all 10 for pressure, and see what happens.

In any case, this matter is deterministic and the code is fully procedural. There is a specific point where the two codes necessarily start to diverge but were identical before. Simply find that point by checking step after step, that's how it is done.

Once you have find that point, the exact line that produces different results, you can then try to find out why and, finally, judge about if it is a reasonable cause for the different results or not.

To be more explicit, you literally have to move a bunch of prints from top to bottom of both codes until you catch the difference (yet, hopefully, you are familiar with binary search and know how to avoid doing this literally line by line).

I'm writing this because you don't usually spot things in this way, by naked eye, so none of us can actually do this for you. Also, as mentioned by Filippo, we don't even know what the differences are in the first place.

[calim_cfd]

Originally Posted by siamak

I am sure about the ranges and the number of the iteration and the inputs. I really wonder why I am getting a different answer.

yeah.. python coding is a pain here.. it begins the couting at 0 .... so.if you dont say anything... a loop with a size N will stop at N-1...which actually is more accurate since the 10 iteration is the beginning of the next ten

[sbaffini]

Originally Posted by calim_cfd

yeah.. python coding is a pain here.. it begins the couting at 0 .... so.if you dont say anything... a loop with a size N will stop at N-1...which actually is more accurate since the 10 iteration is the beginning of the next ten

So, I see, the thing is mostly connected to C indexing, lists, splitting lists and slicing, all of which remain consistent if range is defined in this way. Probably makes sense, yet painful af to digest

[siamak]

I really appreciate all your answers
It was better to post the code from the first point.
I have changed all the iteration variables and noticed that at the first iteration, the results are close but the differences in the result are increasing.

[sbaffini]

Originally Posted by siamak

I really appreciate all your answers
It was better to post the code from the first point.
I have changed all the iteration variables and noticed that at the first iteration, the results are close but the differences in the result are increasing.

Just to clarify some points:

- with reference to the python code, you have checked that the two codes produce the same results up to the last *CHECKED before the pressure but differ just after it, where you write #Differences in pressure solver?

- such differences are already present in the first iteration?

- couldn't avoid noting that you are using 5 points only for each direction, but you otherwise have issues? Sounds suspicious to me, yet not directly related to the differences

- out of curiosity, what's the hardcoded 1.7?

Now we can't really look and fix any of the two codes, honestly (not me at least). But you could help us help you by providing additional details. For example, what's the amount of the difference?

Now it comes a problem. If the difference is beyond machine precision, it can mean one of two things to me (considering that the pressure iteration line appears identical in the two codes):

1) There is an offending operation in that line that would put the code to work outside of the standard, at least for one of the codes. This might either be due already to the initial values of any of the variables involved or be a result of the operations in that line. You should activate all the controls and flags in C++ to check this (I don't know exactly how this works in python)

2) You didn't really check everything

There is really no escape from this. There must be a point in the code where it is not really doing what you think it is doing and you have to find it.

If doing this produces inconsistent results for the actual code (i.e., that in the end there is no error and the two just produce different results for no apparent reason), you might want help yourself by reducing the complexity of the two codes until you spot the modification that restores equality.

However, all this time might be invested also in a different, more proficuous manner, by actually debugging each code separately, because having issues beyond a 5x5 grid is not ok. There is certainly some problem that a small grid helps to hide. When you have both codes working and verified for accuracy with grid refinement then you can more effectively compare them.

[siamak]

I commented out all other parts of the code and used lots of the test cases for comparing the output of the codes. I have used 5 grids just for speeding the test cases (using pytest and doctest). I did not want to clutter the question with all the test cases I tried.

[sbaffini]

Originally Posted by siamak

I commented out all other parts of the code and used lots of the test cases for comparing the output of the codes. I have used 5 grids just for speeding the test cases (using pytest and doctest). I did not want to clutter the question with all the test cases I tried.

Ok, but if you only answer randomly at the questions, we're not going anywhere here.

Do you have an L1, L2 or Linf norm of the difference when it first appears?

Can you actually plot the difference on a grid? That usually helps identifying issues.

What's that 1.7 in the code?

Have you actually verified both the codes? Do you have convergence plots for the errors?

[siamak]

@sbaffini
Thanks for your answer
I think You got the point, I tried to check the L1 and L2, but it seems after some iterations it goes behind my machine's precision. I have to check it more carefully to be sure about it.