Grid refaktoring

[Kirill Streltsov]

I think that a few bugs have been introduced by the refactoring:

1) Why does Grid need a Simulation AND simWidth, simHeight?
2) SimpleGrid does not really work now, because the interpolation of the current fails. It tries to access the index numCellsX+1, the accessor makes it to numCellsX+2, but the maximum index is numCellsX+1.
3) It would be nice to have an empty grid. If I now use the GUI and enforce the display of fields right from the start, it shows me 3x2 (or something like that...) arrows, but in the simulation constructor a 10x10 grid is introduced:

        grid = GridFactory.createSimpleGrid(this, 10, 10, width, height);

Btw this can not work, because a few lines above width and height are set to 0 therefore cellWidth and cellHeight are set to 0. In the set method we then divide by 0...and maybe in some other places too...

[Andreas Ipp]

It seems you did not take into account the shifted indices in this commit:

https://github.com/kstrel/openpixi/commit/50ee9087d2885c1dae27bc0991be966a4605f38c

The fields that are accessible by accessors always have indices from 0 to numCells - 1, so one should shift all indices by one:

Ex[x][y] -> setEx(x-1, y-1).

Alternatively, one can change the for-loop to run from (0, rows-1) instead of (1, rows).

The new class abstracts the extra row around the grid, so you have to use (-1, -1) to access what previously was (0,0):

https://github.com/openpixi/openpixi/blob/master/pixi/src/main/java/org/openpixi/pixi/physics/grid/Grid.java

[Jan Kis]

1) The future plan of mine was to decouple the simulation from the grid completely, that is why I introduced there already simWidth and simHeight

2) You mean the default grid created in the simulation constructor? Give me please closer details on what should I run to invoke the problem.

3) As discussed yesterday the grid can change to interface and then we will have an empty grid.

I will address all three issues this weekend :) 

[Kirill Streltsov]

@Andreas: I understand that, but I did not want to touch the "extra belt" that Jan has added. As far as I understood it, what has been (0,0) before is (1,1) now. And I think that this is exactly what I use here.

If I would treat (0,0) as a normal grid point we would have always +2 cells now, instead of the amount that actually should be there.

[Kirill Streltsov]

@Jan: go to Particle2DPanel.paintComponent(), set

if(drawFields)
to
if(true)

and start the MainControlApplet. You will gain the same result if you do this with if(drawCurrentGrid).

[Andreas Ipp]

It is difficult to see, because of the many empty lines, but in this commit, Jan changed the loop ranges:

https://github.com/openpixi/openpixi/commit/1ee020a75aa1e02eb90a4f00a7ba84c41803d31e

at the same time when he changed the indices in the accessors. So you should do the same in the loop ranges (and outside the loop you should replace 0 -> -1, etc.)

[Kirill Streltsov]

Jan, can you please clear things up, I think me and Andreas have different interpretations you your changes.

I understand it the following way: I want to have a 10x10 grid, so the values of numCellsX,Y are = 10 but the arrays 12x12. If I want to access the lower left boundary term of the SIMULATION I write getArray(0,0), and the getter shifts my request to the non boundary term of the ARRAY which is (1,1). The ARRAY-Boundaries are not ment to be written to or accessed deliberately since they were introduced with the purpose to enable derivatives at the SIMULATION boundaries, because these need the terms left and right from getArray(0,0).

If I interpret getArray(-1.-1) as the SIMULATION boundary, then I not only get a 12x12 grid even though I wanted to have only 10x10, but also can not go over the boundaries just as it was before the introduction of accessors...

If I want the electric field at the SIMULATION boundary I have to make a derivative of phi, so I write

setE(0,0, getPhi(1,0)-getPhi(-1,0) / (2*cellWidth))

right?

@Andreas: you wrote:
"In the past what was (1,1) is now (0,0). In order to access the old (1,1) you have to write (0,0). In order to access the old (0,0) (which was the extra border around the grid) you now have to access (-1,-1)."

But I do not want to access the old positions....before Jan's changes the grid was two cells smaller!

If I still haven't understood this till tomorrow, I will drop by at your office after 11.

[Andreas Ipp]

Ok I see. I was confused by the many boundary cases.

So where else do you suspect the error to be? If there is an index out of bounds, somewhere the shift of indices was apparently not taken into account properly. Or do you have some other idea where the bug could be? Do you notice anything strange if you step through it with the debugger?

Your unit tests worked before you merged Jan's changes, didn't they? And all unit tests worked in Jan's fork. So the error must have been introduced when merging the two? Or not?

[Jan Kis]

Oki, so here it comes:

When I was merging the code all unit tests passed.

When you create a grid 10 x 10 then numCellsX and numCellsY is 10.

You should work with the accessors as if you had normal grid 10 by 10. Additionaly, to this you can also access positions -1 and 10 but you where should always be 0 (in current scenario). Thus, when making changes to the grid do never call setter on -1 or 10. -1 and 10 are logically not part of the grid they were just introduced so that when you sweep over the grid and calculate a point on the grid from surrounding points you do not have to check for border cases because there is by default 0. 

[Andreas Ipp]

Ok. I also had a close look at the code now and I agree that there are a couple of bugs left from Jan's modification of the code. (or to phrase it differently, the conversion was not fully completed yet).

Unfortunately, we didn't have unit tests to catch those issues.

* class Grid:

 - resetCurrentAndCharge(): The for-loops should loop through all elements, so a "+2" needs to be added here. This fixes the problem that Kirill mentioned that currents at the border grow.

 - same in storeFields().

* class CloudInCell:

 - xCellPosition and yCellPosition are defined from 1 to numCellsX + 1, yet they are used directly in getEx()... probably they should be shifted by -1 before they are used in getEx().

Probably similar issues are present in other classes. It is a good question how to write suitable unit tests to catch these things...

Here is a step-by-step instruction how to reproduce the issue with the 3-2 arrows mentioned by Kirill:

* Launch the interactive app

* Select cell > calculate fields > current  => 10x10 arrows are shown

* Press "reset" -> only 2 arrows are shown. Actually, all 100 arrows are drawn, but their distance is much larger and so they are not shown on the screen, because width and height used in the Grid class don't agree with the width and height used in the Simulation class. I think this is the issue that Kirill wanted to point out that could appear when duplicating values rather than referencing them.

* Press "calculate fields" twice to turn on and off: All 100 arrows are shown.

* Press "reset" again: Only 2 arrows are shown.

So I guess I'll pull Kirill's changes into the master, and then Jan can work on these issues over the weekend.

[Kirill Streltsov]

It was a bit offtopic in the Parallelization thread so I post it here. I have used setRho instead of addRho in the interpolator and therefore the unit test did not pass.

But still some issues remain: the interpolator writes to the -1 and numCells+2 components. The -1 component is only accessed if some random grid size is chosen, like 62x81. The numCells+2 component is (nearly) always accessed. Therefore I have to run the loop from -1 to numCells+2...The problem is here:

            int xCellPosition = (int) (Math.floor(p.getX() / g.getCellWidth()));
            int yCellPosition = (int) (Math.floor(p.getY() / g.getCellHeight()));

My first explanation for the access of numCells+2 was that getX(), getY() are <= width, height respectively. For the case where the = sign applies this would lead to an access of numCells+2. But then I have changed the condition to <= width-0.1 and the componend was still accessed...

An other thing is this: while the test passes with an accuracy of 10^-15 on a "well sized" grid like 100x100 or 10x10 where I do not have to run the loop from -1, I have to reduce the accuracy limit in the other case.

[Jan Kis]

But still some issues remain: the interpolator writes to the -1 and numCells+2 components. The -1 component is only accessed if some random grid size is chosen, like 62x81. The numCells+2 component is (nearly) always accessed. Therefore I have to run the loop from -1 to numCells+2...The problem is here:

            int xCellPosition = (int) (Math.floor(p.getX() / g.getCellWidth()));
            int yCellPosition = (int) (Math.floor(p.getY() / g.getCellHeight()));

My first explanation for the access of numCells+2 was that getX(), getY() are <= width, height respectively. For the case where the = sign applies this would lead to an access of numCells+2. But then I have changed the condition to <= width-0.1 and the componend was still accessed...

I believe we have a serious bug in the grid which is also related to what you are saying.

The logical grid has two properties:

1) it has cells and these cells are separated by 

2) grid points

If we want a grid with 10 by 10 cells we need 11 by 11 grid points.

Currently we only create 10 by 10 grid points (technically it is 10 + 2: 2 for the extra border so that we can read -1 and size() indexes). Nevertheless, we need 11 by 11 grid points. Consequently, we can introduce a new getter in the grid class getNumPointsX() and getNumPointsY(). Furthermore, the entire code which is using getNumCellsX/Y() needs to be revisited and checked whether we really meant number of cells or number of grid points.

I can do the necessary changes in the grid class but I would need Kirill to repair the remaining classes (Interpolators and Field solvers). However, as with the parallelization we would need to change the grid class again it is questionable whether we should do this changes now or later when we come up with the design of grid suitable for parallelization.

[Andreas Ipp]

But still some issues remain: the interpolator writes to the -1 and numCells+2 components. The -1 component is only accessed if some random grid size is chosen, like 62x81. The numCells+2 component is (nearly) always accessed. Therefore I have to run the loop from -1 to numCells+2...The problem is here:

            int xCellPosition = (int) (Math.floor(p.getX() / g.getCellWidth()));
            int yCellPosition = (int) (Math.floor(p.getY() / g.getCellHeight()));

My first explanation for the access of numCells+2 was that getX(), getY() are <= width, height respectively. For the case where the = sign applies this would lead to an access of numCells+2. But then I have changed the condition to <= width-0.1 and the componend was still accessed...

Do you check (maybe with an assert) that the values of p.getX() and p.getY() are where you expect them to be? (like 0 <= p.getX() < width, etc). It may be that the particle pusher pushes the particles outside of the allowed region and then you would see what you see, but this is not a problem of the interpolator, but of the particle pusher.

 

I believe we have a serious bug in the grid which is also related to what you are saying.

The logical grid has two properties:

1) it has cells and these cells are separated by 

2) grid points

If we want a grid with 10 by 10 cells we need 11 by 11 grid points.

Depends on what those grid points mean. If the meaning of the field is, e.g. the charge in the center of each cell then you need as many grid points as cells. I think charges were defined in the center of a cell, and electric fields at the border between two cells (but they are also labelled by simple integers). So the index of a field may mean different things, and this can also depend on the interpolator used.

I can do the necessary changes in the grid class but I would need Kirill to repair the remaining classes (Interpolators and Field solvers). However, as with the parallelization we would need to change the grid class again it is questionable whether we should do this changes now or later when we come up with the design of grid suitable for parallelization.

If there is a problem that we don't fully understand, we should first fix the problem before we move on.

[Jan Kis]

Depends on what those grid points mean. If the meaning of the field is, e.g. the charge in the center of each cell then you need as many grid points as cells. I think charges were defined in the center of a cell, and electric fields at the border between two cells (but they are also labelled by simple integers). So the index of a field may mean different things, and this can also depend on the interpolator used.

Thank you very much Andreas, this is very valuable explanation. If the index of the field can really mean with one interpolator center of the cell and in another interpolator border of the cell then the bug remains because we only allocate space for the cell center points. 

I suggest to accommodate the code according to your explanation so that it is clear from the code that we can use the fields in two different ways. In order to do so I suggest the following changes:

1) create a field of size number of border points

2) have two types of getters 

  - getCenterPointsX() == getNumCellsX()

  - getBorderPointsX() == getNumCellsX() + 1

3) redo the existing code to use these so that it is explicitly visible which interpolator uses the center of the cell and which uses the borders of the cell

@(Kirill and Andreas) Does this sound resonable? If so I can do the changes 1 and 2 and help Kirill to do the 3. 

[Kirill Streltsov]

No I haven't checked it with asserts, one should definitely do that. Although this should not be a problem since the particle pusher is not accesed before the interpolation is performed. If the particle pusher was to blame, then we would have the same problem on the left boundary, but we dont. Not even a prepare or complete is called before the interpolation, but that would only change the velocity anyway...

Yes, the indices may mean different positions. Look at the FDTD wikipedia article. Currently I have assumed that the chargedensity is defined at the four surrounding grid points of a cell. The E fields, in case of the YeeGrid, are defined in the middle of each cell edge and are perpendicular to them.

I thought that in the algorithms we always go from 0 to <numCells and if these algorithms try to access -1 or whatever, the accessors take care of it. I can not think of a clever way to implement this other than doing the boundary checking in the accessors...

Regarding your suggestion Jan, I am not sure if thats a good idea. We would introduce even more accessors that do exactly the same thing...No matter where the value is defined I still need a "belt" around the simulation. The only difference is that those quantities that are defined on the grid points need one grid cell more, but one could take care of this transparently...I would prefer an explanation via comment.

I think that the boundary problem as well as the parallelization are interconnected and should be dealt with simultaniously. I will try to help out, but I am currently very busy with sutying for an exam that I have next friday...

[Andreas Ipp]

I also don't think there is a necessity to introduce additional getters. Someone needs to carefully go through the existing code and fix the bugs. Maybe not much change is necessary to fix the problems. The only problem is to find the time to do it.

Since it seems this is more of a physicist's task, it shouldn't keep Jan from starting to plan the parallelization. Maybe this could happen in a separate branch, so that there is no time pressure to finalize version 0.4, but work on 0.5 can already start.

The only thing I could imagine could be generalized in the current approach is the size of the "safety belt". Currently it is hard-coded to 1 cell ("magic number" ;-) ). If this could be set to an arbitrary number (0, 1, 2, ..) depending on the algorithm, the Grid classes wouldn't / shouldn't need to know what the values that they store mean (whether it is cell center, border center, or grid points). The caller is reposible for reserving the right amount of cells and "safety" cells around the grid, and what they store in those fields is completely up to the calling routines.