Commit 46858619 authored by James Willis's avatar James Willis
Browse files

Removed scripts and copied to examples.

parent 50dcfedf
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)
import h5py as h
import numpy as np
import matplotlib
matplotlib.use("Agg")
from pylab import *
import os.path
kernel_gamma = 1.825742
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
inputFile1 = ""
inputFile2 = ""
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()
file_naiv = h.File(inputFile1, "r")
file_sort = h.File(inputFile2, "r")
ids_naiv = file_naiv["/PartType0/ParticleIDs"][:]
ids_sort = file_sort["/PartType0/ParticleIDs"][:]
h_naiv = file_naiv["/PartType0/SmoothingLength"][:]
h_sort = file_sort["/PartType0/SmoothingLength"][:]
pos_naiv = file_naiv["/PartType0/Coordinates"][:,:]
pos_sort = file_sort["/PartType0/Coordinates"][:,:]
num_density_naiv = file_naiv["/PartType0/Num_ngb_density"][:]
num_density_sort = file_sort["/PartType0/Num_ngb_density"][:]
num_force_naiv = file_naiv["/PartType0/Num_ngb_force"][:]
num_force_sort = file_sort["/PartType0/Num_ngb_force"][:]
#wcount_naiv = file_naiv["/PartType0/Wcount"][:]
#wcount_sort = file_sort["/PartType0/Wcount"][:]
#
#wcount_naiv = wcount_naiv * np.power(h_naiv,3) * kernel_norm
#wcount_sort = wcount_sort * np.power(h_sort,3) * kernel_norm
# Cross check
print "Ngbs density naiv: ", np.min(num_density_naiv), np.mean(num_density_naiv), np.max(num_density_naiv)
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_naiv), np.mean(num_force_naiv), np.max(num_force_naiv)
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_naiv), np.mean(wcount_naiv), np.max(wcount_naiv)
#print "Wcount sort: ", np.min(wcount_sort), np.mean(wcount_sort), np.max(wcount_sort)
# Sort
index_naiv = np.argsort(ids_naiv)
index_sort = np.argsort(ids_sort)
num_density_naiv = num_density_naiv[index_naiv]
num_density_sort = num_density_sort[index_sort]
num_force_naiv = num_force_naiv[index_naiv]
num_force_sort = num_force_sort[index_sort]
ids_naiv = ids_naiv[index_naiv]
ids_sort = ids_sort[index_sort]
#wcount_naiv = wcount_naiv[index_naiv]
#wcount_sort = wcount_sort[index_sort]
h_naiv = h_naiv[index_naiv]
h_sort = h_sort[index_sort]
# First check
print "Differences for density: ", min(num_density_naiv - num_density_sort), max(num_density_naiv - num_density_sort)
print "Differences for force: ", min(num_force_naiv - num_force_sort), max(num_force_naiv - num_force_sort)
# Get the IDs that are different
mask_density = num_density_naiv != num_density_sort
mask_force = num_force_naiv != num_force_sort
print "Num non-zero density: ", np.sum(mask_density)
print "Num non-zero force: ", np.sum(mask_force)
print "\nParticle IDs with incorrect densities"
print "----------------------------------------"
print ids_naiv[mask_density]
print "Num of density interactions", inputFile1
print num_density_naiv[mask_density]
print "Num of density interactions", inputFile2
print num_density_sort[mask_density]
print "\nParticle IDs with incorrect forces"
print "------------------------------------"
print ids_naiv[mask_force]
print "Num of force interactions", inputFile1
print num_force_naiv[mask_force]
print "Num of force interactions", inputFile2
print num_force_sort[mask_force]
# Statistics of h difference
h_relative = (h_naiv - h_sort) / h_naiv
print "h statistics:", np.min(h_relative), np.percentile(h_relative, 1)
print "h statistics:", np.mean(h_relative), np.median(h_relative)
print "h statistics:", np.percentile(h_relative, 99) , np.max(h_relative)
#print mean(h_naiv)
#exit()
# Make some plots of the wrong particles
#rcParams.update({'figure.figsize': (6,6)})
#figure()
#subplot(221)
#plot(pos_naiv[mask_force,0], pos_naiv[mask_force,1], 'bo', ms=2)
#for i in range(12):
# plot([0, 1], [i / 12., i/12.], 'k-', lw=0.5, color='0.5')
# plot([i / 12., i/12.], [0, 1], 'k-', lw=0.5, color='0.5')
#
#
#subplot(222)
#plot(pos_naiv[mask_force,0], pos_naiv[mask_force,2], 'bo', ms=2)
#for i in range(12):
# plot([0, 1], [i / 12., i/12.], 'k-', lw=0.5, color='0.5')
# plot([i / 12., i/12.], [0, 1], 'k-', lw=0.5, color='0.5')
#
#
#subplot(223)
#plot(pos_naiv[mask_force,1], pos_naiv[mask_force,2], 'bo', ms=2)
#for i in range(12):
# plot([0, 1], [i / 12., i/12.], 'k-', lw=0.5, color='0.5')
# plot([i / 12., i/12.], [0, 1], 'k-', lw=0.5, color='0.5')
#
#
#savefig("wrong_force.png")
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment