Patch vecorised pair interactions

Implements:

• Small optimisations to runner_dopair1_density_vec.

Added 1 commit:

• 854c66b9 - Formatting.

mentioned in issue #354 (closed)

Have you tested it ? I'd bet a pint that the test-suite fails.

I will not take that bet. I will run the test-suite.

testPeriodicBC should show you what's wrong.

Added 1 commit:

• ee6316b8 - Sign change bug.

Try again. Glad my test worked.

Indeed, good test of the test. That's the bug I had spotted.

How about the same bit of code for the DOPAIR2 case ?

While we are at it: Any reason why we shift both cells to the frame of ci->loc ? Wouldn't cj->loc be a better choice given that's where we compute the shift for in the case of periodic BC?

What do you mean about DOPAIR2? The reference frame is arbitrary, but how is it better in the frame of cj?

2 separate things, sorry for the confusion.

• Just thinking that the same will have to be done for the cache construction in DOPAIR2.
• For the other aspect I am not 100% sure yet. Let me try to detail this.

The idea of changing frame is that we want to be somewhere where the particle coordinates are close to 0, leaving as many digits as possible for the differences in particle distance within the PAIR. If the two cells are actual neighbours then either ci or cj's position is fine as they are basically the same. Now, if we want to apply periodic BC, we start by moving the particles from one side of the box to the other. If I am not mistaken, we do this by moving the particles in ci to the frame of cj using the shift vector. At this stage the positions are commensurate with the position of cj. So removing cj->loc from them would bring all positions in ci and cj to something close to 0, leaving all digits for precise position differences. If we instead remove ci->loc, it is still correct, but we are burning some digits unnecessarily I think. Does that make sense ? Basically, all the positions in the cache should ideally be in the range [0, 2*c->width]. And I am not sure this is the case now with periodic BCs.

Shoudn't the range be: [0,2*sqrt(2)*c->width] or [-sqrt(2)*c->width, sqrt(2)*c->width]?

why sqrt(2) ? I guess you are talking about the diagonal case so sqrt(3) should be the limit.

Oh yeah in 3D it's sqrt(3).

Okay, so I have tested this and when we use ci as the location we get cases that are outside of: [-sqrt(3)*c->width, sqrt(3)*c->width], but when we use cj its fine. So I will change the location to that of cj.

Excellent! Out of interest do you know which cases ? An by how much ?

ci->loc[7.000000,7.000000,0.000000], cj->loc[0.000000,0.000000,0.000000]
shift[8.000000,8.000000,-0.000000]
Error: cj->loc[7.000000,7.000000,0.000000] Particle 0 x pos is not within [-ci->width, ci->width]. x=-7.083333, ci->width[0]=1.000000

That's the first case, but I throw an error so I'm not sure how many more there are.

Ok. So that's exactly the case we had drown on the piece of paper.

Yeah, exactly

I get an error in test125cells now though:

ci->loc[0.000000,0.000000,1.000000], cj->loc[0.000000,1.000000,0.000000]
shift[0.000000,0.000000,0.000000]
Error: cj->loc[0.000000,1.000000,0.000000] Particle 0 z pos is not within [-ci->width, ci->width]. z=1.916667, ci->width[2]=1.000000

but there are no periodic conditions in effect.

Can't explain this...

Does that mean testPeriodicBC is not testing that configuration ?

Correct, but it's not supposed to because it is not a periodic boundary condition.

I think it's fine, because the particles won't be shifted in the z dimension. I think this condition: [0,2*sqrt(3)*c->width] should apply in this case.

So the test should be [-sqrt(3)*c->width, 2*sqrt(3)*c->width] ?

I think so, for all possible cases

Cool

That passes when running test125cells.

Added 4 commits:

• e7e6c6c6 - Updated tolerances when cj->loc is removed from particle positions.
• 26622f44 - Shift the particle positions in the cache by cj instead of ci to prevent loss of…
• 88b2163a - Compute h_max correctly for cells.
• 2be912fe - Defined a dimension variable and used it in loops and to allocate memory.

Added 1 commit:

• 0fd2b62e - Formatting.

Added 1 commit:

• 966bed0c - Missed sqrt(3) term in error message.

Of course EAGLE_25 throws an error.

Error: cj->loc[11.148849,16.597804,12.040922] Particle 0 z pos is not within [-sqrt(3)*ci->width, 2*sqrt(3)*ci->width]. z=-0.002415, ci->width[2]=0.001379

after the 67th step. The values are closeish, z=-0.002415 sqrt(3)*ci->width = -0.002388498.

Mmmmh... Not sure what to say...

Is it the test that's not accurate enough because of a mix of precision ?

I think it should be between: [-2*sqrt(3)ci->width, 2sqrt(3)*ci->width], shouldn't it? That updated condition still fails:

Error: cj->loc[11.148849,16.597804,12.040922] Particle 0 z pos is not within [-2*sqrt(3)*ci->width, 2*sqrt(3)*ci->width]. z=-0.004795, ci->width[2]=0.001379

The values are closer, z=-0.004795 2*sqrt(3)*ci->width = -0.004776996.

I'm wondering if I need to add a small factor to the comparison because the values are so close

I guess 2*sqrt(3) can be replaced by 4.

Now, if the condition has to be [-4c->widh, 4c->width], then why is this case not appearing in the testPeriodicBC ?

[-4c->widh, 4c->width] still doesn't work

z=-0.005564, ci->width[2]=0.001379. So -0.005564 vs -0.005516

Wait, after how many steps is that ? Remember that we allow particles to move a bit out of their cells so we should take this into account when writing the condition.

279

Ok, so particles are allowed to move. That factor needs taking into account. It's space_stretch.

@jwillis actually, it can be more than that. Check cell_need_rebuild_for_pair() to see what the maximal motion of a particle outside of its starting cell can be.

I have used 4. * ci->width[0] * (1. + 2.*space_maxreldx) and EAGLE_25 ran fine. And the reason it's not triggered in testPeriodicBC is because the particles never move outside their cells.

Sure but didn't we conclude that [-c->width, 2c->width] was the range ? The movement should just multiply in not double only the negative range.

I doubled the negative bound, because the shift could be negative also. Is that not right?

That's what I am confused about. In that case it should be part of testPeriodicBC, no ?

Okay, so:

• [-2*sqrt(3) * ci->width[0] * (1. + 2.*space_maxreldx), 2*sqrt(3) * ci->width[0] * (1. + 2.*space_maxreldx)] fails EAGLE_25,
• [-4 * ci->width[0] * (1. + 2.*space_maxreldx), 4 * ci->width[0] * (1. + 2.*space_maxreldx)] passes EAGLE_25,
• [-2 * ci->width[0] * (1. + 2.*space_maxreldx), 4 * ci->width[0] * (1. + 2.*space_maxreldx)] fails EAGLE_25.

I will try and figure out why the 3rd bullet point does not fail in testPeriodicBC