Dopair subset vec

Vectorises the dopair_subset function. A branched version has been implemented that only calls the vectorised version when two cells are face-on as the edge-on and corner-on cases run slower.

added 45 commits

• 52657340...5c87444c - 42 commits from branch master
• 18f70b6e - Merge branch 'master' into dopair_subset_vec
• fd3fb786 - Forgot to rename function after merge with master.
• 98dd170d - Formatting.

Compare with previous version

unmarked as a Work In Progress

Should I remove -P SPH:h_tolerance:10 from the check_interactions.sh script so it also tests the subsets? Also I couldn't read the box size from the hdf5 file, it says that boxsize doesn't exist.

/Header/BoxSize is definitely in the snapshots. Note it is an attribute and not an array so you need to read it as such with h5py.

Ah okay, I checked that it was. It just couldn't read it. I was using:

box_size = file_naive["/Header/BoxSize"]

Ah right, I'm just reading the h5py interface now and you have to read the whole header from the hdf5 file and retrieve the box size.

added 1 commit

• e5472326 - Read box size from snapshot.

Compare with previous version

That should be fixed now and ready to be looked at.

Looks good to me. The only thing I don't like is that lower tolerance. I don't see why or how we only have 200 ulps of precision for h when doing periodic BCs.

Yeah me too, this is an example of when it happens:

Particle 3674430220107 is missing 3677638524641, hig2: [  1.56662779e-07], r2: 1.56666041082e-07, |r2 - hig2|: [  3.26849658e-12]

It's when the separation is very small.

How is the distance computed in the script? With a shift or via a box-wrapping? Is it similar to the way it is computed in the code?

It is now computed with a box wrapping:

def nearest(dx):
    if(dx > 0.5 * box_size):
        return dx - box_size
    elif(dx < -0.5 * box_size):
        return dx + box_size
    else:
        return dx

...
            dx = pi_pos[0][0] - pj_pos[0][0]
            dy = pi_pos[0][1] - pj_pos[0][1]
            dz = pi_pos[0][2] - pj_pos[0][2]
    
            # Correct for BCs
            dx = nearest(dx)
            dy = nearest(dy)
            dz = nearest(dz)

            r2 = dx*dx + dy*dy + dz*dz
    
            hig2 = hi*hi*kernel_gamma2
    
            diff = abs(r2 - hig2)
    
            print "Particle {} is missing {}, hig2: {}, r2: {}, |r2 - hig2|: {}".format(pid,pjd,hig2, r2, diff)
    
            if diff < acc * hig2:
                print "Missing interaction due to precision issue will be ignored."
            else:
                return True

Also does it happen in the subset or in the normal pair? If in the subset is it because we shift by the cell where we process stuff and not necessarily by the cell where the particles are?

It only fails the test when I run without -P SPH:h_tolerance: 10 so running with the subsets on. Ah I see what you mean, how do we I get the location of the cell where the particles are?

Actually no, that discrepancy was found between the naive and serial where we don't take the location of the cell off the particles.

How so? Is the subset navie or serial not done in the frame of the cell?

No neither are.

Ah! Ok, so that may explain the lack of accuracy there. You said the vectorized version actually finds that neighbour?

Yeah, the vectorised version actually finds that neighbour which matches the naive.