Poisson Solver

[Kirill Streltsov]

I have been working on the PoissonSolver. You can check my progress in my PoissonSolver branch.

It now calculates the electric field and saves it in the Ex and Ey arrays of the grid. I am currently testing it...Here you can see the potential and electric field of a point charge:
http://htu.at/~kirill/bak/potential.jpg
http://htu.at/~kirill/bak/efield.png

The transformation assumes that periodic boundaries are present
http://htu.at/~kirill/bak/potential-linecharge.png
therefore I have made the class non static and renamed it to PoissonSolverPediodic.

As we discussed a while ago there is still a problem with the poisson equation in fourier space, which corresponds to these lines

        //solve Poisson equation in Fourier space
        for(int i = 0; i < columns; i++) {
            for(int j = 0; j < rows; j++) {
                double d = (4 - 2 * Math.cos((2 * Math.PI * i) / g.numCellsX) - 2 * Math.cos((2 * Math.PI * j) / g.numCellsY));
                if (d != 0) {
                phi[i][2*j] = (cellArea * trho[i][2*j]) / d;
                phi[i][2*j+1] = (cellArea * trho[i][2*j+1]) / d;
                } else {
//                    phi[i][2*j] = trho[i][2*j];
//                    phi[i][2*j+1] = trho[i][2*j+1];
                    phi[i][2*j] = cellArea * trho[i][2*j];
                    phi[i][2*j+1] = cellArea * trho[i][2*j+1];
//                    phi[i][2*j] = 0;
//                    phi[i][2*j+1] = 0;
                }
            }
        }

There are three possibilities for the case where the denominator is 0. Here (http://www.physics.buffalo.edu/phy410-505-2004/Chapter6/ch6-lec2.pdf) they use the first version. Here are the plots of a point charge for each version:

http://htu.at/~kirill/bak/potential1.jpg
http://htu.at/~kirill/bak/potential2.jpg
http://htu.at/~kirill/bak/potential3.jpg

As you can see the second and third version are a bit darker but equal to each other. The indices in the fourier space correspond to the frequecies of the fields. The only case where the denominator is 0 is for i,j=0 which corresponds to the lowest frequencies and therefore energies.(I will remove the if from the for loop and treat this case separately) In the plotted case the following grid parameters were used
        s.width = 10;
        s.height = 6;
        g.numCellsX = 100;
        g.numCellsY = 100;

therefore cellArea < 1. My interpretation of the result is that the lowest energies were filtered out and therefore the total potential is a bit lower. But I am not sure if the filtering of the LOWEST energies can have such a visible effect...
Anyway, none of these versions makes sense mathematically! Analytically one would "make a circle" around this sigularity and then take a limes but we cant do that int the descrete case...An other option would be to use forward difference for the case i,j=0. The "problem" here would be that we would get imaginary parts and therefore a bit more code...

I am thinking about how to make a unit test that creates a lot of plots for different charge distributions and code versions... Does anyone know a good way to plot a potential with gnuplot? Currently I plot just a lot of small circles (does not look good if there only a few grid points)...
I have hardcoded the gnuplot commands and the output path for windows. This will cause assertion errors if gnuplot is not installed but the test will still pass. Is it ok if I include it in the pull request in this form?

[Kirill Streltsov]

Sorry, these are the correct links

http://htu.at/~kirill/bak/potential1.png
http://htu.at/~kirill/bak/potential2.png
http://htu.at/~kirill/bak/potential3.png

[Andreas Ipp]

I think the answer that we came up with is: It doesn't matter.

The potential is defined up to an additive constant. The important point is that the E-fields obtained from any of the three potentials equal each other (which you can test).

You could also write:

else {

   // if we reach here, i==0 and j==0
   phi[0][0] = Math.random();
   phi[0][1] = 0;  // we don't want an imaginary part

}

and still obtain a correct potential, because the fourier transform of the [0][0] entry just gives a constant contribution, which does not matter in the potential. So the simplest solution is to set it to 0.

Regarding unit tests:

It is meant that these tests can run regularly to quickly check if everything is still working (after one changed some code, for example). Ideally, the unit test should not have any output to the console or to the hard disk if everything is correct. For example, all unit tests are executed if one wants to create a jar file using "mvn package". Also, the field unit tests take a long time because many random moves are tested. It is ok to do so once during development, but once everything works, this should be reduced to a reasonable size (maybe max. 1000, ideally without output if everything works as expected) so that all unit tests can be quickly executed.

If you have code for output, you could put it for example into ui.test or ui.develop?

[Jan Kis]

I agree with Andreas about the unit tests - they should give us a fast way to partially verify that everything is working as it is supposed to. As far as the physiscs is concerned I am not of much help :)  

[Kirill Streltsov]

Ah...right, how could I forget that...:-( I will take that case out of the for loop and add a comment about it. I will not set the value to 0 but leave it untouched to have two writes less. (like they do in the example code)

Unit tests:
So we want to have the Unit tests just to test if the programm runs without errors? What about the physics? We got to test it somewhere...that is also why we need the extensive Unit tests in the Grid because some errors were only visible after 10000 or more random moves. Maybe the comparison of a point charge potential to its analytic result would be a good unit test for the PoissonSolver?

I will definitely need some analysis features as total energy, entropy, charge distribution etc. Of course it would be nice to have these parameters and plots displayed dynamicaly in the gui but I think that it will be quite time consuming to implement all these plots... I would like to put these features in the org.openpixi.pixi.analysis package and the code for the output in ui.util.

Still, it would be nice if one would have a place to play around with certain, strictly defined parts of the simulation (where one knows what happens from the beginning till the end) and test the physics behind it. I can not use the gui for that because I simply do not understand the code and do not know what it is doing. Even Ognen does not known what it is doing! (the initial conditions are not set correctly, the collisions and the timestep get overriden and maybe other things too) I have tried to run it through the debugger but right at the beginning of the MainControlApplet constructor it says source not found...
So where should I do my tests? ui.develop? I guess some of the code could later be reused in initial conditions or as demonstration material...

@Jan: are you working on the Grid? The PoissonSolver is basically finished, I need to do the charge density interpolation which is the same as CIC. I would include the PoissonSolver calculation in the grid constructors and make a pull request. But if you are changing a lot on the structure there I would not write it into the grid constructors for now.

[Jan Kis]

Unit tests:
So we want to have the Unit tests just to test if the programm runs without errors? What about the physics? We got to test it somewhere...that is also why we need the extensive Unit tests in the Grid because some errors were only visible after 10000 or more random moves.

I believe we want to have the unit tests for the following purposes.

1)  quickly verify that the changes we make do not break the essential parts of the application

2)  test quickly whether the physical calculations are as we expect them to be (of course since the unit tests have to be relatively fast and they have to have a very simple output pass/fail we can not fully test the physical behavior)

To address your problem I suggest creating a new package named eg. physics.analysis in which you can write code to produce graphs so that you can closely analyze the behavior after a large number of moves.  

I can not use the gui for that because I simply do not understand the code and do not know what it is doing. Even Ognen does not known what it is doing! (the initial conditions are not set correctly, the collisions and the timestep get overriden and maybe other things too)

This is exactly, why I would like to take look at the simulation class and improve it a bit so that it better encapsulates the functionality it should offer. Afterwards, the gui can be rewritten so that it is easier to understand. To pinpoint the problem of the gui I will ask the following question. Do you think you have problems understanding the code because you are not familiar enough with the java swing api or because the code is in general too complicated and unpredictable? 

 I have tried to run it through the debugger but right at the beginning of the MainControlApplet constructor it says source not found...

This seems to be more of an issue of the debugger. I am using Netbeans and I can debug the application without problems. If Andreas gives you freedom in what IDE to use you might want to try out also other IDEs than Eclipse and see which one suits you the best. For example, I find Netbeans much more user friendly. I understand better what it does and I think it is also easier to set it up (the maven is included in the IDE so no plugin is necessary). Another alternative is IDEA (http://www.jetbrains.com/idea/).    

@Jan: are you working on the Grid? The PoissonSolver is basically finished, I need to do the charge density interpolation which is the same as CIC. I would include the PoissonSolver calculation in the grid constructors and make a pull request. But if you are changing a lot on the structure there I would not write it into the grid constructors for now.

Currently, I am not doing any code modifications. I plan to start working on the grid on Tuesday so you can add the PoissonSolver calculation in the meantime into the grid constructor.

[Andreas Ipp]

 I have tried to run it through the debugger but right at the beginning of the MainControlApplet constructor it says source not found...

This can happen if you try to debug one of the "internal" Java classes of swing for which you don't have the source code attached. In that case, you can repeatedly press "Step return" to get out of the system classes until you are back in one of the pixi classes. Press "step over" if it is a system or swing class, and only press "step into" if it is a pixi class where you have the source code. (I think it is also possible to attach the source code of the swing classes from the JDK, but it is probably not too enlightning to debug those classes...)

I am using Netbeans and I can debug the application without problems. If Andreas gives you freedom in what IDE to use you might want to try out also other IDEs than Eclipse and see which one suits you the best. For example, I find Netbeans much more user friendly. I understand better what it does and I think it is also easier to set it up (the maven is included in the IDE so no plugin is necessary).

Of course there is freedom to use one's favorite IDE :-) Jan, if you know Netbeans better, you could add instructions to the README file for Netbeans:

https://github.com/openpixi/openpixi/blob/master/pixi/README.txt

 

Another alternative is IDEA (http://www.jetbrains.com/idea/).    

The readme could also include instructions for IDEA if someone happens to use that.

 

[Kirill Streltsov]

Yesterday we said that the PoissonSolver should not tinker with the variables in grid. So I can move the potential phi and the charge density into the PoissonSolver class but I have to output the electric fields somehow. These are two double[][] arrays. Cant really put that into a return...

I can pass them to the poisson solver instead of the grid to make it more clear what is happeneing and decouple the solver from the grid. But I also need cellWidth and cellHeigth as well as the number of cells (can get that from the length of the arrays that are passed)...So thats four variables that I would have to pass...

So whats better, one grid or four variables?

[Andreas Ipp]

Hm, given that we anticipate that the Grid may change a lot for parallelization, it is maybe better to decouple PoissonSolver from Grid as much as possible and pass the four variables.

But since PoissonSolver always works on a Grid, I think it is also fine to pass the grid, wouldn't it?

[Jan Kis]

I think you can pass the whole grid since you need to access the electric fields. These are on purpose not exposed so you can not extract from grid just Ex (double[][]). You can only access it through the accessors eg. grid.getEx(x, y).  

[Andreas Ipp]

Oh, I see, you want to write to the electric fields. I agree with Jan, better pass the grid, and use the grid.setEx(..) methods to set the electric fields from within the PoissonSolver.

Am Donnerstag, 24. Mai 2012 12:11:01 UTC+2 schrieb Jan Kis:

I think you can pass the whole grid since you need to access the electric fields. These are on purpose not exposed so you can not extract from grid just Ex (double[][]). You can only access it through the accessors eg. grid.getEx(x, y).  

[Kirill Streltsov]

Ok so I can not or should not decouple grid from the solver, but does it make sense to move phi and rho to the solver then?

If you are going to implement the GridPeriodic and GridHard etc with apropriate accessors it would be useful to have these accessors for phi and rho too (actually I think only for phi). If they are in the solver this is not possible. The advantage of having accessors is that I can use the same (short) code for the derivation of the E fields in PoissonSolverPeriodic and PoissonSolverHard. I am thinking about the following structure:

package PoissonSolver
interface PoissonSolver
class FFTPeriodic
class FFTDirichlet (hardwall)
class FFTNeumann (hardwall)
static class DeriveField

this enables an easy implementation of other solvers. Maybe we want local solvers sometime that do not rely on the FFT...all of them can rely on the deriveField method of a static class since it will be boundary sensitive because it will use the accessors of the grid.

If I pass rho and phi I would save them in PoissonSolverEmpty and make other classes extend this class. They would also have to have their own (but only slightly different) deriveField methods.

[Andreas Ipp]

If you use phi and rho only within the Solver, and never in the grid, they could stay in the solver. But it should be easy to move them later back to Grid if one needs them there.

A remark:

After looking through some of yesterday's links, I think it doesn't make sense to think in terms of GridPeriodic or GridHard. It makes more sense to think of a bare grid, to which one can attach on each side a different boundary condition, like periodic boundary, hard boundary, or transparent communication boundary to another grid. The boundary condition can be different for each direction.

But I think in PoissonSolver you have to use different fourier transforms depending on the kind of boundary, so here one is more limited.

I'm not sure if you would call local solvers in the same way as you call the PoissonSolver. So I would suggest not to make it more general than it needs to be now. (In a parallel version, one may call a local solver for a small part of the whole simulation, and I think this would work quite differently from calling a PoissonSolver once that calculates a Fourier Transform of the whole simulation area).

Instead of using an interface, DirichletPoissonSolver and NeumannPoissonSolver could be derived from PeriodicPoissonSolver so they can share the same basic method for deriving the electric field from the rho field.

But anything is fine, as long as it works and can be refactored later :-)