New script to compare particle neighbour lists when running swift in debug interactions mode.

examples/EAGLE_12/check_ngbs.py

+import h5py as h
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import os.path
+
+kernel_gamma = 1.825742
+kernel_gamma2 = kernel_gamma * kernel_gamma
+kernel_gamma_dim = np.power(kernel_gamma,3)
+hydro_dimension_unit_sphere = 4. * np.pi / 3.
+kernel_norm = hydro_dimension_unit_sphere * kernel_gamma_dim
+error = False
+
+inputFile1 = ""
+inputFile2 = ""
+
+# Check list of density neighbours and check that they are correct.
+def check_density_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
+
+    error_val = False
+
+    for k in range(0,num_invalid):
+
+        # Filter neighbour lists for valid particle ids
+        filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
+        filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
+
+        # Check neighbour lists for differences
+        id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
+        
+        pid = pids[mask][k]
+
+        # Loop over discrepancies and check if they are actually neighbours
+        for pjd in id_list:
+            pi_pos = pos[np.where(pids == pid)]
+            pj_pos = pos[np.where(pids == pjd)]
+            
+            hi = h[np.where(pids == pid)]
+            
+            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]
+            
+            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:
+                error_val = True
+
+    return error_val
+
+# Check list of force neighbours and check that they are correct.
+def check_force_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
+
+    error_val = False
+
+    for k in range(0,num_invalid):
+
+        # Filter neighbour lists for valid particle ids
+        filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
+        filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
+
+        # Check neighbour lists for differences
+        id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
+        
+        pid = pids[mask][k]
+
+        # Loop over discrepancies and check if they are actually neighbours
+        for pjd in id_list:
+            pi_pos = pos[np.where(pids == pid)]
+            pj_pos = pos[np.where(pids == pjd)]
+            
+            hi = h[np.where(pids == pid)]
+            hj = h[np.where(pids == pjd)]
+            
+            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]
+            
+            r2 = dx*dx + dy*dy + dz*dz
+            
+            hig2 = hi*hi*kernel_gamma2
+            hjg2 = hj*hj*kernel_gamma2
+            
+            diff = abs(r2 - max(hig2, hjg2))
+            
+            print "Particle {} is missing {}, hig2: {}, hjg2: {}, r2: {}, |r2 - max(hig2,hjg2)|: {}".format(pid,pjd,hig2, hjg2, r2, diff)
+
+            
+            if diff < acc * max(hig2,hjg2):
+                print "Missing interaction due to precision issue will be ignored."
+            else:
+                error_val = True
+
+    return error_val
+
+# Parse command line arguments
+if len(sys.argv) < 3:
+    print "Error: pass input files as arguments"
+    sys.exit()
+else:
+    inputFile1 = sys.argv[1]
+    inputFile2 = sys.argv[2]
+    if os.path.exists(inputFile1) != 1:
+        print "\n{} does not exist!\n".format(inputFile1)
+        sys.exit()
+    if os.path.exists(inputFile2) != 1:
+        print "\n{} does not exist!\n".format(inputFile2)
+        sys.exit()
+
+# Open input files    
+file_naive = h.File(inputFile1, "r")
+file_sort = h.File(inputFile2, "r")
+
+# Read input file fields
+ids_naive = file_naive["/PartType0/ParticleIDs"][:]
+ids_sort = file_sort["/PartType0/ParticleIDs"][:]
+
+h_naive = file_naive["/PartType0/SmoothingLength"][:]
+h_sort = file_sort["/PartType0/SmoothingLength"][:]
+
+pos_naive = file_naive["/PartType0/Coordinates"][:,:]
+pos_sort = file_sort["/PartType0/Coordinates"][:,:]
+
+num_density_naive = file_naive["/PartType0/Num_ngb_density"][:]
+num_density_sort = file_sort["/PartType0/Num_ngb_density"][:]
+
+num_force_naive = file_naive["/PartType0/Num_ngb_force"][:]
+num_force_sort = file_sort["/PartType0/Num_ngb_force"][:]
+
+neighbour_ids_density_naive = file_naive["/PartType0/Ids_ngb_density"][:]
+neighbour_ids_density_sort = file_sort["/PartType0/Ids_ngb_density"][:]
+
+neighbour_ids_force_naive = file_naive["/PartType0/Ids_ngb_force"][:]
+neighbour_ids_force_sort = file_sort["/PartType0/Ids_ngb_force"][:]
+
+#wcount_naive = file_naive["/PartType0/Wcount"][:]
+#wcount_sort = file_sort["/PartType0/Wcount"][:]
+#
+#wcount_naive = wcount_naive * np.power(h_naive,3) * kernel_norm
+#wcount_sort = wcount_sort * np.power(h_sort,3) * kernel_norm
+
+# Cross check
+print "                   Min     Mean     Max "
+print "                   ---------------------"
+print "Ngbs density naiv: ", np.min(num_density_naive), np.mean(num_density_naive), np.max(num_density_naive)
+print "Ngbs density sort: ", np.min(num_density_sort), np.mean(num_density_sort), np.max(num_density_sort)
+print "Ngbs force naiv:   ", np.min(num_force_naive), np.mean(num_force_naive), np.max(num_force_naive)
+print "Ngbs force sort:   ", np.min(num_force_sort), np.mean(num_force_sort), np.max(num_force_sort)
+#print "Wcount naiv:   ", np.min(wcount_naive), np.mean(wcount_naive), np.max(wcount_naive)
+#print "Wcount sort:   ", np.min(wcount_sort), np.mean(wcount_sort), np.max(wcount_sort)
+
+# Sort
+index_naive = np.argsort(ids_naive)
+index_sort = np.argsort(ids_sort)
+
+num_density_naive = num_density_naive[index_naive]
+num_density_sort = num_density_sort[index_sort]
+num_force_naive = num_force_naive[index_naive]
+num_force_sort = num_force_sort[index_sort]
+ids_naive = ids_naive[index_naive]
+ids_sort = ids_sort[index_sort]
+neighbour_ids_density_naive = neighbour_ids_density_naive[index_naive]
+neighbour_ids_density_sort = neighbour_ids_density_sort[index_sort]
+neighbour_ids_force_naive = neighbour_ids_force_naive[index_naive]
+neighbour_ids_force_sort = neighbour_ids_force_sort[index_sort]
+#wcount_naive = wcount_naive[index_naive]
+#wcount_sort = wcount_sort[index_sort]
+h_naive = h_naive[index_naive]
+h_sort = h_sort[index_sort]
+pos_naive = pos_naive[index_naive]
+pos_sort = pos_sort[index_sort]
+
+# First check
+print "\n                   Min    Max"
+print "                   ----------"
+print "Differences for density:  ", min(num_density_naive - num_density_sort), max(num_density_naive - num_density_sort)
+print "Differences for force:    ", min(num_force_naive - num_force_sort), max(num_force_naive - num_force_sort)
+
+# Get the IDs that are different
+mask_density = num_density_naive != num_density_sort
+mask_force = num_force_naive != num_force_sort
+num_invalid_density = np.sum(mask_density)
+num_invalid_force = np.sum(mask_force)
+
+print "\nNum non-zero density: ", num_invalid_density
+print "Num non-zero force:   ", num_invalid_force
+
+print "\nParticle IDs with incorrect densities"
+print "----------------------------------------"
+print ids_naive[mask_density]
+
+# Check density neighbour lists
+error += check_density_neighbours(ids_naive, neighbour_ids_density_naive,
+        neighbour_ids_density_sort, mask_density, pos_naive, h_naive,
+        num_invalid_density, 2e-6)
+
+print "Num of density interactions", inputFile1
+print num_density_naive[mask_density]
+
+print "Num of density interactions", inputFile2
+print num_density_sort[mask_density]
+
+print "\nParticle IDs with incorrect forces"
+print "------------------------------------"
+print ids_naive[mask_force]
+
+# Check force neighbour lists
+error += check_force_neighbours(ids_naive, neighbour_ids_force_naive,
+        neighbour_ids_force_sort, mask_force, pos_naive, h_naive,
+        num_invalid_force, 2e-6)
+
+print "Num of force interactions", inputFile1
+print num_force_naive[mask_force]
+
+#print "Smoothing lengths", inputFile1
+#print h_naive[mask_force]
+
+print "Num of force interactions", inputFile2
+print num_force_sort[mask_force]
+
+#print "Smoothing lengths", inputFile2
+#print h_sort[mask_force]
+
+# Statistics of h difference
+h_relative = (h_naive - h_sort) / h_naive
+print "h statistics: {} {} (Min, 1st Percentile)".format(np.min(h_relative), np.percentile(h_relative,1))
+print "h statistics: {} {} (Mean, Median)".format(np.mean(h_relative), np.median(h_relative))
+print "h statistics: {} {} (Max, 99th Percentile)".format(np.max(h_relative), np.percentile(h_relative, 99))
+
+if error:
+    print "\n------------------"
+    print "Differences found."
+    print "------------------"
+    exit(1)
+else:
+    print "\n---------------------"
+    print "No differences found."
+    print "---------------------"
+    exit(0)