periodic boundary conditions

[glao]

Hello , i am trying to understand how periodic boundary condition works.Of course, i understand what it is and how it works,just i am trying to implement it in c++ and i have some problems.( for now i want to make it work in 1d but i want it for 2d also).

Lets say for example i have this line :

 -2    -1    0    1     2         The size(length) of the line is 4.
 |-----|-----|----|-----|

Ok , now for example i want to create a function which take as inputs the position of the 'particle' ,the step that it does and the size of the line.

I create for example this table:

initial position |  step  |  final position
--------------------------------------------
  0              |     2  |         2  <size        we are ok here , the particle will go to '2'
  1              |     2  |         3  <size  

  here the particle must go to the end of line (2) and then   begin from start (-2) and go finally to -2.


 2              |     4  |         6  >size        here it must go finally to position 1.

I can see that if the final position is for example >size/2 (the second example above) or <-size/2 ,
i must do a calculation (i can't figure what!) in order to start from the beginning?

Any help?

Also, i know that this works with the modulo operator also , but right now i can't figure.I must understand this first.

Thank you!

[ravenous]

To make a simple periodic-boundary-condition-like access to an array, you can do

size_t size = 5;
double array[ size ];

// ... Fill the array with data

std::cout << "enter index:" << std::endl;
int index;
std::cin >> index;

double value = array[ ( size + ( index % size ) ) % size ];
std::cout << "array[ " << index << "] = " << value << std::endl;

This is pretty ugly, so it's probably best to encapsulate it in a function:

size_t GetPeriodicIndex( int index, unsigned size )
{
    return static_cast< unsigned >( size + ( index % size ) ) % size;
}

Using this method, you can do a python-like array access of using -1 to get the last element of the array. Since you have to do a bunch of extra calculations at each access though, it will be slower than just accessing the elements directly.

[glao]

Thanks for your answer.
I just need sth close to my approach in order to understand.

[glao]

Well, i managed to do sth :

//function to compute periodic boundary conditions for 1d size
double boundary1d(int position,int step , int size){
//size : size of the line
//position : initial position of particle
//step : step that the particle does

/*


     -2   -1  0   1   2    -> size=5
        |---|---|---|---|



    position | step | final position
    -------------------------------
    2           1       -2  (2+1=3 is > size/2 ,so 3-5=-2)
    2           2       -1  (2+2=4 is > size/2 , so 4-5=-1)
    2           7       -1  (2+7=9 is > size/2 , so 9-5=4 ,4 is >size/2 , so 4-5=-1)
    -1          1        0  (-1+1=0 is < size/2 , so leave it)
    -1          3        2  (-1+3=2 is < size/2 , so leave it)
*/
int final_pos=position+step;//final position

if (final_pos < size/2) final_pos=position+step;
if (final_pos > size/2) {

do {
    final_pos=(position+step)-size;
}while (final_pos < size/2);

}


return final_pos;
}

The problem is that when the (position+step) is >size/2 the program doesn't give the right results or the loop doesn't stop.For example :

boundary1d(2,2,5) , it should give '-1' , but instead the loop doesn't end.
boundary1d(2,7,5) , it should give '-1' , but it gives '4'.It seems it doesn't do the next step (to subtract size again).

Any ideas? Thanks!

[glao]

Any ideas? Thanks!

[sfuo]

I don't know what this "sth" is and I don't know what assumptions you are making as far as the size coming in goes, but here is a method you could use.

int boundary( int pos, int step, int size )
{
    return (pos + static_cast<int>(size/2) + step)%size - static_cast<int>(size/2);
}

The assumptions that are being made here are that size will be odd and that the range of values will always be equal in both directions off 0. And that step is going to be positive.

[glao]

Hello ,

I managed to do it finally!

double boundary1d(int position,int step , int size){

int final_pos=position+step;//final position

//for step <0
if (step<0) {
    if (abs(final_pos)> static_cast<int>(size/2) )

    do {

    final_pos+=size;

}while (final_pos > static_cast<int>(size/2));

}

if (final_pos > static_cast<int>(size/2)) {

do {

    final_pos-=size;

}while (final_pos > static_cast<int>(size/2));

}

return final_pos;
}

I had wrong the while statement before.

return (pos + static_cast<int>(size/2) + step)%size - static_cast<int>(size/2);

Also, i used the case where step<0 (you gave me the idea)

I have 2 questions though about your way (with the modulo)

1) I can't understand what's the idea of using the modulo , i mean how did you think of using modulo in this case?
2) I didn't manage to use it with step <0.

( the static_cast<int> is necessary here? Because size is int)

Thanks

[sfuo]

Modulos gives the remainder of the two numbers if they were used in division.

Places where I use the mod operator are where I see wrapping. So if I had a problem where I wanted to know what angle a line was making from the center of a cricle, I would probably want it to be from [0-360) and not some larger or smaller angle.

The reason why I used static_cast<int> is because I was just trying to show that it will give an int and not pop out a x.5 (and for some reason you are returning double on your function when you are dealing with ints).

Here are two methods that you can use. One is just a cleaned up version of what you posted (an if with the same condition of the do{}while() can be replaced with just a while()) and a method that uses % for both positive and negative steps.

Cleaned up while() method:

int boundaryWHILE( int pos, int step, int size )
{
    if( step == 0 )
        return pos;

    int fpos = pos + step;

    while( fpos < -(size/2) )
        fpos += size;

    while( fpos > (size/2) )
        fpos -= size;

    return fpos;
}

mod method:

int boundaryMOD( int pos, int step, int size )
{
    if( step == 0 )
        return pos;

    step %= size;

    if( step < 0 )
        step += size;

    return (pos + size/2 + step)%size - size/2;
}

The modulus method can probably be cleaned up a bit but I hope you get the general idea.

[ravenous]

Sorry, missed the ongoing discussion. You don't really need any branching logic or loops for this. It's a one-liner (albeit a slightly confusing one):

int boundary1d( int position, int step , int size )
{
    return ( size + ( ( position + step + (size/2) ) % size ) ) % size - (size/2);
}

The position + step part generates the new position (in terms of the problem reference frame, not array indices). The problem reference frame is off-set by size/2 from array indices, so position + step - (size/2) moves us into an array index frame. Say

i_array = position + step - (size/2)

then doing i_array % size gives us a number in the range [-size - 1,size - 1], since i_array could, and will often, be negative. We can't have a negative array index, so we add size to this number, giving a number in the range [0, 2*size - 1]. Lets say that

i_array_2 = size + ( i_array % size )

So, a final i_array_2 % size gives us a value in the range [0, size - 1]. The final thing is to move back into the problem reference frame by subtracting the size/2 that we added at the start, giving the final

( size + ( ( position + step + (size/2) ) % size ) ) % size - (size/2)

I hope that makes some sense.

I also changed the return-type of your function to int, since that makes more sense than a double.

[glao]

Ok, thank you both!
Now ,its better!
Thanks!

(I am trying to do it now in 2d dimensions)