Skip to content
Snippets Groups Projects
Commit ca916c92 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added the calculation of the orthogonal distance to the diagonal linking the... 

Added the calculation of the orthogonal distance to the diagonal linking the cells for each particles. Should help understand the bias.
parent 3e84351a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/plot_sorted.py
+ 102
− 51
View file @ ca916c92
......@@ -33,6 +33,9 @@ rcParams.update(params)
rc('font', family='serif')
import sys
import os
from scipy import stats
dist_cutoff_ratio=1.2
print "Plotting..."
......@@ -68,19 +71,28 @@ axis = [
#names = ["side", "edge", "corner"]
for orientation in range( 26 ):
# for jjj in range(2):
# if jjj == 0:
# orientation = 0
# if jjj == 1:
# orientation = 8
#for orientation in range( 26 ):
for jjj in range(3):
if jjj == 0:
orientation = 0
if jjj == 1:
orientation = 1
if jjj == 2:
orientation = 4
# Read Quickshed accelerations
data=loadtxt( "interaction_dump_%d.dat"%orientation )
id = data[:,0]
pos = data[:,2]
pos -= mean(pos)
x = data[:,3]
y = data[:,4]
z = data[:,5]
dist = data[:,12]
accx_u=data[:,6]
accy_u=data[:,7]
accz_u=data[:,8]
......@@ -112,56 +124,95 @@ for orientation in range( 26 ):
stdz_s = std(errz_s[abs(errz_s) < 0.1])
# sample_pos = pos[dist<0.2]
# sample_x = e_errx_s[dist<0.2]
# sample_y = e_erry_s[dist<0.2]
# sample_z = e_errz_s[dist<0.2]
# numBins = 100
# binEdges = linspace(-1.5-1.5/numBins, 1.5+1.5/numBins, numBins+2)
# bins = linspace(-1.5, 1.5, numBins+1)
# corr_x, a, b = stats.binned_statistic(sample_pos, sample_x, statistic='mean', bins=binEdges)
# corr_y, a, b = stats.binned_statistic(sample_pos, sample_y, statistic='mean', bins=binEdges)
# corr_z, a, b = stats.binned_statistic(sample_pos, sample_z, statistic='mean', bins=binEdges)
# a,b, sample_bin = stats.binned_statistic(pos, pos, statistic='mean', bins=binEdges)
# sample_bin -= 2
# # for j in range(size(pos)):
# # e_errx_s /= corr_x[sample_bin[j]]
# # e_erry_s /= corr_y[sample_bin[j]]
# # e_errz_s /= corr_z[sample_bin[j]]
# Plot -------------------------------------------------------
figure(frameon=True)
subplot(311, title="Acceleration along X")
#plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
plot(pos, e_errx_s , 'ro')
# for j in range(size(pos)):
# if ( pos[j-1] != pos[j] ):
# text(pos[j], e_errx_s[j]-0.005, "%d"%id[j], fontsize=10)
# else:
# text(pos[j], e_errx_s[j]-0.015, "%d"%id[j], fontsize=10)
text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
scatter(pos[dist<1.], e_errx_s[dist<1.] ,c=dist[dist<1.], marker='o', s=1, linewidth=0, cmap='jet')
plot([-dist_cutoff_ratio/2., -dist_cutoff_ratio/2.], [-0.01, 0.01], 'k--')
plot([dist_cutoff_ratio/2., dist_cutoff_ratio/2.], [-0.01, 0.01], 'k--')
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
ylim(-0.05, 0.05)
ylim(-0.03, 0.03)
grid()
subplot(312, title="Acceleration along Y")
#plot(id[abs(erry_s) > 0.001], e_erry_s , 'ro')
plot(pos, e_erry_s , 'ro')
# for j in range(size(pos)):
# if ( pos[j-1] != pos[j] ):
# text(pos[j], e_errx_s[j]-0.005, "%d"%id[j], fontsize=10)
# else:
# text(pos[j], e_errx_s[j]-0.015, "%d"%id[j], fontsize=10)
text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
scatter(pos[dist<1.], e_erry_s[dist<1.] , c=dist[dist<1.], marker='o', s=1, linewidth=0, cmap='jet')
plot([-dist_cutoff_ratio/2., -dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
plot([dist_cutoff_ratio/2., dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
ylim(-0.05, 0.05)
ylim(-0.03, 0.03)
grid()
subplot(313, title="Acceleration along Z")
#plot(id[abs(errz_s) > 0.001], e_errz_s , 'ro', label="Sorted")
plot(pos, e_errz_s , 'ro', label="Sorted")
# for j in range(size(pos)):
# if ( pos[j-1] != pos[j] ):
# text(pos[j], e_errx_s[j]-0.005, "%d"%id[j], fontsize=10)
# else:
# text(pos[j], e_errx_s[j]-0.015, "%d"%id[j], fontsize=10)
text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
legend(loc="upper right")
scatter(pos[dist<1.], e_errz_s[dist<1.] , c=dist[dist<1.], label="Sorted", marker='o', s=1, linewidth=0, cmap='jet')
plot([-dist_cutoff_ratio/2., -dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
plot([dist_cutoff_ratio/2., dist_cutoff_ratio/2.], [-0.03, 0.03], 'k--')
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
#legend(loc="upper right")
xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
ylim(-0.05, 0.05)
#xlim(0, size(id)-1)
ylim(-0.03, 0.03)
grid()
savefig("1quadrupole_accelerations_relative_%d.png"%orientation)
savefig("quadrupole_accelerations_relative_%d.png"%orientation)
close()
figure(frameon=True)
subplot(311, title="Acceleration along X")
scatter(dist, e_errx_s )
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01)
grid()
subplot(312, title="Acceleration along Y")
scatter(dist, e_erry_s )
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01)
grid()
subplot(313, title="Acceleration along Z")
scatter(dist, e_errz_s , label="Sorted")
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
#legend(loc="upper right")
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01)
grid()
savefig("error_distance_%d.png"%orientation)
close()
# # Plot -------------------------------------------------------
# figure(frameon=True)
......@@ -207,22 +258,22 @@ for orientation in range( 26 ):
# figure(frameon=True)
# subplot(311, title="Acceleration along X")
# hist(e_errx_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Sorted")
# legend(loc="upper right")
# xlim(-0.12, 0.12)
figure(frameon=True)
subplot(311, title="Acceleration along X")
hist(e_errx_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Sorted")
legend(loc="upper right")
xlim(-0.12, 0.12)
# subplot(312, title="Acceleration along Y")
# hist(e_erry_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
# xlim(-0.12, 0.12)
subplot(312, title="Acceleration along Y")
hist(e_erry_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
xlim(-0.12, 0.12)
# subplot(313, title="Acceleration along Z")
# hist(e_errz_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
# xlim(-0.12, 0.12)
subplot(313, title="Acceleration along Z")
hist(e_errz_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
xlim(-0.12, 0.12)
# savefig("histogram_%d.png"%orientation)
# close()
savefig("histogram_%d.png"%orientation)
close()
......
examples/test_bh_sorted.c
+ 26
− 12
View file @ ca916c92
......@@ -190,6 +190,7 @@ static inline void multipole_iact(struct multipole *d, const float *x, float* a)
w = const_G * ir * ir * ir * d->mass;
for (k = 0; k < 3; k++) a[k] -= w * dx[k];
/* Compute some helpful terms */
float w1 = const_G * ir * ir * ir * ir * ir;
float w2 = -2.5f * const_G * ir * ir * ir * ir * ir * ir * ir;
......@@ -2207,18 +2208,30 @@ void test_direct_neighbour(int N_parts, int orientation) {
}
/* Position along the axis */
double *position;
double *position, *ortho_dist;
if ((position = (double *)malloc(sizeof(double) * N_parts * 2)) == NULL)
error("Failed to allocate position buffer.");
if ((ortho_dist = (double *)malloc(sizeof(double) * N_parts * 2)) == NULL)
error("Failed to allocate orthogonal distance buffer.");
float w = 1.0f / sqrtf(shift[0]*shift[0] + shift[1]*shift[1] + shift[2]*shift[2]);
shift[0] *= w;
shift[1] *= w;
shift[2] *= w;
for (k = 0; k < 2 * N_parts; ++k)
position[k] = parts[k].x[0] * shift[0] + parts[k].x[1] * shift[1] +
parts[k].x[2] * shift[2];
for (k = 0; k < 2 * N_parts; ++k) {
float dx = parts[k].x[0] * shift[0];
float dy = parts[k].x[1] * shift[1];
float dz = parts[k].x[2] * shift[2];
position[k] = dx + dy + dz;
dx = (parts[k].x[1] - 0.5f)*shift[2] - (parts[k].x[2] - 0.5f)*shift[1];
dy = (parts[k].x[2] - 0.5f)*shift[0] - (parts[k].x[0] - 0.5f)*shift[2];
dz = (parts[k].x[0] - 0.5f)*shift[1] - (parts[k].x[1] - 0.5f)*shift[0];
ortho_dist[k] = sqrtf(dx*dx + dy*dy + dz*dz);
}
/* Now, output everything */
char fileName[100];
......@@ -2226,12 +2239,13 @@ void test_direct_neighbour(int N_parts, int orientation) {
message("Writing file '%s'", fileName);
FILE *file = fopen(fileName, "w");
fprintf(file,
"# ID m r x y z a_u.x a_u.y a_u.z a_s.x a_s.y a_s.z\n");
"# ID m r x y z a_u.x a_u.y a_u.z a_s.x a_s.y a_s.z ortho\n");
for (k = 0; k < 2 * N_parts; k++) {
fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,
fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,
parts[k].mass, position[k], parts[k].x[0], parts[k].x[1],
parts[k].x[2], parts[k].a_exact[0], parts[k].a_exact[1],
parts[k].a_exact[2], parts[k].a[0], parts[k].a[1], parts[k].a[2]);
parts[k].a_exact[2], parts[k].a[0], parts[k].a[1], parts[k].a[2],
ortho_dist[k]);
}
fclose(file);
......@@ -2490,11 +2504,11 @@ int main(int argc, char *argv[]) {
N_parts);
/* Run the test */
for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
/* k++; */
/* test_direct_neighbour(N_parts, 0); */
/* test_direct_neighbour(N_parts, 8); */
//for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
test_direct_neighbour(N_parts, 0);
test_direct_neighbour(N_parts, 1);
test_direct_neighbour(N_parts, 4);
k++;
/* Dump arguments */
message("Interacting one cell with %d particles with its 27 neighbours"
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment