From 61aedaecffe6fd3743daf4ff3e77f8534b161e0d Mon Sep 17 00:00:00 2001
From: d74ksy <aidan.chalk@durham.ac.uk>
Date: Mon, 2 Nov 2015 10:43:43 +0000
Subject: [PATCH] BH completed, can now output particle data, compute exact
interaction and use plot.py to compute error
---
examples/plot.py | 26 ++++-----
examples/test_bh_mpi.c | 129 +++++++++++++++++++++++++++++++++++++++--
src/qsched.c | 21 ++++---
3 files changed, 151 insertions(+), 25 deletions(-)
diff --git a/examples/plot.py b/examples/plot.py
index 941a33f..ce4c28b 100644
--- a/examples/plot.py
+++ b/examples/plot.py
@@ -37,19 +37,19 @@ import os
print "Plotting..."
# Read Quickshed accelerations
-data=loadtxt("particle_dump.dat")
+data=loadtxt("particle_dump.out")
id = data[:,0]
-accx_e=data[:,5]
-accy_e=data[:,6]
-accz_e=data[:,7]
+accx_e=data[:,4]
+accy_e=data[:,5]
+accz_e=data[:,6]
-accx_bh=data[:,8]
-accy_bh=data[:,9]
-accz_bh=data[:,10]
+accx_bh=data[:,7]
+accy_bh=data[:,8]
+accz_bh=data[:,9]
-accx_new=data[:,11]
-accy_new=data[:,12]
-accz_new=data[:,13]
+accx_new=data[:,10]
+accy_new=data[:,11]
+accz_new=data[:,12]
# Sort accelerations
rank = argsort(id)
@@ -127,7 +127,7 @@ plot(id, errx_new , 'b.')
text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" )
-ylim(-0.2, 0.2)
+ylim(-2, 2)
xlim(0,id[-1])
grid()
@@ -138,7 +138,7 @@ plot(id, erry_new , 'b.')
text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" )
-ylim(-0.2, 0.2)
+ylim(-0.02, 0.02)
xlim(0,id[-1])
grid()
@@ -150,7 +150,7 @@ plot(id, errz_bh , 'rx', label="Legacy")
#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" )
legend(loc="upper right")
-ylim(-0.2, 0.2)
+ylim(-0.02, 0.02)
xlim(0,id[-1])
grid()
diff --git a/examples/test_bh_mpi.c b/examples/test_bh_mpi.c
index 20e2c35..aff0302 100644
--- a/examples/test_bh_mpi.c
+++ b/examples/test_bh_mpi.c
@@ -35,6 +35,7 @@
/* Local includes. */
#include "quicksched.h"
+#include "res.h"
/* Some local constants. */
#define cell_pool_grow 1000
@@ -48,6 +49,7 @@
#define SANITY_CHECKS
#define NO_COM_AS_TASK
#define NO_COUNTERS
+#define EXACT
@@ -57,8 +59,13 @@ struct part {
union {
float a[3];
float a_legacy[3];
+ #ifndef EXACT
float a_exact[3];
+ #endif
};
+ #ifdef EXACT
+ float a_exact[3];
+ #endif
float mass;
int id;
}; // __attribute__((aligned(32)));
@@ -992,20 +999,103 @@ void create_tasks(struct qsched *s, struct cell *ci, struct cell *cj) {
} /* Otherwise it's a pair */
}
+inline int getindex(long long int id, struct qsched *s){
+ return id;
+}
+/**
+ * @brief Solve the particle interactions using the stupid N^2 algorithm
+ *
+ * @param N The number of particles
+ * @param parts The array of particles
+ */
+void interact_exact(int N, struct part *parts) {
+ int i, j, k;
+ float ir, w, r2, dx[3];
+ /* Loop over all particles. */
+ for (i = 0; i < N; ++i) {
+ /* Some things to store locally. */
+ double pix[3] = {parts[i].x[0], parts[i].x[1], parts[i].x[2]};
+ float mi = parts[i].mass;
+ /* Loop over every other particle. */
+ for (j = i + 1; j < N; ++j) {
+ /* Compute the pairwise distance. */
+ for (r2 = 0.0, k = 0; k < 3; k++) {
+ dx[k] = parts[j].x[k] - pix[k];
+ r2 += dx[k] * dx[k];
+ }
+ /* Apply the gravitational acceleration. */
+ ir = 1.0f / sqrtf(r2);
+ w = const_G * ir * ir * ir;
+ for (k = 0; k < 3; k++) {
+ float wdx = w * dx[k];
+ parts[j].a_exact[k] -= wdx * mi;
+ parts[i].a_exact[k] += wdx * parts[j].mass;
+ }
+ }
+ }
+#if ICHECK >= 0
+ /* for (i = 0; i < N; ++i)
+ if (parts[i].id == monitor)
+ message("[check] exact acceleration for particle %d a= %.3e %.3e %.3e\n",
+ parts[i].id, parts[i].a_exact[0], parts[i].a_exact[1],
+ parts[i].a_exact[2]);*/
+#endif
+}
-
-
-
+void output_parts(struct qsched *s)
+{
+ int i, j;
+ FILE *file;
+ file = fopen("particle_dump.out", "a");
+ for(i = 0; i < s->res_ranks[s->count_ranks]; i++)
+ {
+ /* We want to find particle resources, which are locked only. */
+ if(s->res[i].node == s->rank && s->res[i].num_lockers > 0)
+ {
+ struct res *resource = &s->res[i];
+ int parent_output = 0;
+ while(resource->parent > 0)
+ {
+ if( s->res[getindex(resource->parent, s)].num_lockers > 0)
+ {
+ parent_output = 1;
+ break;
+ }
+ resource = &s->res[getindex(resource->parent, s)];
+ }
+ /* If we didn't output any of the parents then we need to output this guy.*/
+ if(!parent_output)
+ {
+ struct part *parts = qsched_getresdata(s, s->res[i].ID);
+ if(parts == NULL)
+ {
+ message("Failed to output parts");
+ fclose(file);
+ return;
+ }
+ int count = s->res[i].size / sizeof(struct part);
+ for(j = 0; j < count; j++)
+ {
+ struct part cur = parts[j];
+ fprintf(file, "%i %e %e %e %e %e %e %e %e %e %e %e %e %e\n", cur.id, cur.mass, cur.x[0], cur.x[1], cur.x[2],
+ cur.a_exact[0], cur.a_exact[1], cur.a_exact[2], cur.a_legacy[0], cur.a_legacy[1], cur.a_legacy[2], cur.a[0], cur.a[1], cur.a[2]);
+ }
+ }
+ }
+ }
+ fclose(file);
+ message("Output parts correctly");
+}
@@ -1024,7 +1114,7 @@ void create_tasks(struct qsched *s, struct cell *ci, struct cell *cj) {
*/
void test_bh(int N, int nr_threads, int runs, char *fileName) {
- int k;
+ int i, k;
struct cell *root = NULL;
struct part *parts;
qsched_res_t parts_res;
@@ -1154,6 +1244,15 @@ if(s.rank == 0)
}
qsched_sync_resources(&s);
+#ifdef EXACT
+ if(s.rank == 0)
+ {
+ interact_exact(N, parts);
+ }
+
+ MPI_Barrier(s.comm);
+#endif
+
if(s.rank == 0)
{
create_tasks(&s, root, NULL);
@@ -1166,6 +1265,28 @@ qsched_run_MPI(&s, nr_threads, runner);
printf("Hello world from processor %s, rank = %i, count_ranks = %i\n",
processor_name, s.rank, s.count_ranks);
+if(s.rank == 0)
+{
+ file = fopen("particle_dump.out", "w");
+ fprintf(file,
+ "ID m x y z a_exact.x a_exact.y a_exact.z a_legacy.x "
+ "a_legacy.y a_legacy.z a_new.x a_new.y a_new.z\n");
+ fclose(file);
+}
+
+//Output the parts:
+for(i = 0; i < s.count_ranks; i++)
+{
+ //If it is our turn to output parts.
+ if(s.rank == i)
+ {
+ output_parts(&s);
+ }
+ MPI_Barrier(s.comm);
+}
+
+//Need to clean up everything.
+
// Finalize the MPI environment.
MPI_Finalize();
}
diff --git a/src/qsched.c b/src/qsched.c
index 25a9261..0b88d3f 100644
--- a/src/qsched.c
+++ b/src/qsched.c
@@ -858,10 +858,10 @@ void qsched_partition_compute_costs( struct qsched *s, idx_t *res_costs)
res_costs[getindex(t->locks[j],s)] += t->cost;
}
- /*for(j = 0; j < t->nr_uses; j++)
+ for(j = 0; j < t->nr_uses; j++)
{
- res_costs[getindex(t->uses[j],s)] = t->cost;
- }*/
+ res_costs[getindex(t->uses[j],s)] += 1;
+ }
}
}
@@ -884,14 +884,14 @@ void qsched_partition_build_nodelist(struct qsched *s, idx_t *nodelist, long lon
for(i = 0; i < s->res_ranks[s->count_ranks]; i++)
{
//If the resources is locked (or used).
- if(s->res[i].num_lockers > 0 /*|| s->res[i].num_users > 0*/){
+ if(s->res[i].num_lockers > 0 || s->res[i].num_users > 0){
r = &s->res[i];
//If it has a parent.
if(r->parent != -1)
{
//If the parent is locked then we don't need this.
- if(s->res[getindex(r->parent,s)].num_lockers > 0 /*|| s->res[getindex(r->parent,s)].num_users > 0*/)
+ if(s->res[getindex(r->parent,s)].num_lockers > 0 || s->res[getindex(r->parent,s)].num_users > 0)
{
noderef[node_count] = -1;
continue;
@@ -901,7 +901,7 @@ void qsched_partition_build_nodelist(struct qsched *s, idx_t *nodelist, long lon
while(r->parent != -1)
{
r = getres(r->parent, s);
- if(r->num_lockers > 0 /*|| r->num_users > 0*/)
+ if(r->num_lockers > 0 || r->num_users > 0)
{
noderef[node_count] = 0;
break;
@@ -1544,6 +1544,7 @@ void tsched_free(struct tsched *ts) {
long long int tsched_addtask ( struct tsched *ts , int type , unsigned int flags , void *data , int data_size , int cost ){
+static int i = 0;
void *temp;
struct task *t;
long long int id;
@@ -1619,7 +1620,9 @@ void *temp;
ts->count += 1;
t->node = ts->rank;
-
+ i++;
+ if(i % 1000 == 0)
+ message("Added %i send/recv tasks", i/2);
/* Return the task ID. */
return id;
@@ -2967,7 +2970,7 @@ for(i = 0; i < node_count; i++)
options[ METIS_OPTION_CONTIG ] = 0; //TODO 1
options[ METIS_OPTION_NCUTS ] = 10;
options[ METIS_OPTION_NITER ] = 10;
- options[ METIS_OPTION_UFACTOR ] = 30;
+ options[ METIS_OPTION_UFACTOR ] = 300;
options[ METIS_OPTION_SEED ] = 359872345;
idx_t one = 1;
@@ -3441,6 +3444,8 @@ printf("Rank[%i]: qsched_prepare_mpi took %lli (= %e) ticks\n", s->rank,
printf("Rank[%i]: Execution took %lli (= %e) ticks\n", s->rank,
toc - tic, (float)(toc - tic));
+
+ MPI_Barrier(s->comm);
#else
error( "QuickSched was not compiled with OpenMP support." );
#endif
--
GitLab