diff --git a/lib/sapporo_light/sapporo.cpp b/lib/sapporo_light/sapporo.cpp index 1afd0c1c4e..4c1b32963d 100644 --- a/lib/sapporo_light/sapporo.cpp +++ b/lib/sapporo_light/sapporo.cpp @@ -79,6 +79,9 @@ int sapporo::set_j_particle(int cluster_id, exit(-1); } + // Skip test particles (massless particles) as they don't contribute to gravity + if (mass <= SAPPORO_TEST_PARTICLE_MASS) return 0; + DS Dmass = (DS){mass, INT_AS_FLOAT(id)}; map::iterator iterator = mapping_from_address_j_to_index_in_update_array.find(address); map::iterator end = mapping_from_address_j_to_index_in_update_array.end(); @@ -134,6 +137,7 @@ void sapporo::calc_firsthalf(int cluster_id, } if (address_j.size() > 0) { + // Since we skip adding test particles, all particles in address_j are massive send_j_particles_to_device(device_id); } @@ -147,7 +151,9 @@ void sapporo::calc_firsthalf(int cluster_id, acc_j.clear(); jrk_j.clear(); } - evaluate_gravity(ni, nj); + if (nj > 0) { + evaluate_gravity(ni, nj); + } } diff --git a/lib/sapporo_light/sapporo.h b/lib/sapporo_light/sapporo.h index 6969d16dfa..f6218f7326 100644 --- a/lib/sapporo_light/sapporo.h +++ b/lib/sapporo_light/sapporo.h @@ -22,6 +22,15 @@ #include #include + +#ifndef SAPPORO_TEST_PARTICLE_MASS +# if defined(_TINY_) +# define SAPPORO_TEST_PARTICLE_MASS ((float)(_TINY_)) +# else +# define SAPPORO_TEST_PARTICLE_MASS 1e-30f +# endif +#endif + using namespace std; #include diff --git a/lib/sapporo_light/sapporo_defs.h b/lib/sapporo_light/sapporo_defs.h index 0ae0876fe4..217721f313 100644 --- a/lib/sapporo_light/sapporo_defs.h +++ b/lib/sapporo_light/sapporo_defs.h @@ -3,7 +3,7 @@ #define MAXCUDADEVICES 4 #define NBODIES_MAX 524288 -#define NBLOCKS 16 /* number of block which can be run simultaneously */ +#define NBLOCKS 256 /* number of block which can be run simultaneously */ #ifdef NGB #define NTHREADS 256 /* max number of threads which can run per block */ diff --git a/lib/sapporo_light/send_fetch_data.cpp b/lib/sapporo_light/send_fetch_data.cpp index e47ae146d4..c8118884fd 100644 --- a/lib/sapporo_light/send_fetch_data.cpp +++ b/lib/sapporo_light/send_fetch_data.cpp @@ -111,6 +111,7 @@ int sapporo::fetch_ngb_list_from_device(int ignore) { double sapporo::evaluate_gravity(int ni, int nj) { + #ifdef NGB bool ngb = true; #else diff --git a/src/amuse_ph4/src/idata.cc b/src/amuse_ph4/src/idata.cc index 4f97d4f376..6a48ff172e 100644 --- a/src/amuse_ph4/src/idata.cc +++ b/src/amuse_ph4/src/idata.cc @@ -204,33 +204,30 @@ void idata::get_partial_acc_and_jerk() ldnn[i] = _INFINITY_; for (int k = 0; k < 3; k++) lacc[i][k] = ljerk[i][k] = 0; for (int j = j_start; j < j_end; j++) { - // Skip particles with zero mass. - if (jdat->mass[j] > _TINY_){ - r2 = xv = 0; - for (int k = 0; k < 3; k++) { - dx[k] = jdat->pred_pos[j][k] - ipos[i][k]; - dv[k] = jdat->pred_vel[j][k] - ivel[i][k]; - r2 += dx[k]*dx[k]; - xv += dx[k]*dv[k]; - } - r2i = 1/(r2+eps2+_TINY_); - ri = sqrt(r2i); - mri = jdat->mass[j]*ri; - mr3i = mri*r2i; - a3 = -3*xv*r2i; - // PRC(jdat->mpi_rank); PRC(ri); PRL(mri); - if (r2 > _TINY_) { - lpot[i] -= mri; - if (r2 < ldnn[i]) { - ldnn[i] = r2; - lnn[i] = j; - } - } - for (int k = 0; k < 3; k++) { - lacc[i][k] += mr3i*dx[k]; - ljerk[i][k] += mr3i*(dv[k]+a3*dx[k]); + r2 = xv = 0; + for (int k = 0; k < 3; k++) { + dx[k] = jdat->pred_pos[j][k] - ipos[i][k]; + dv[k] = jdat->pred_vel[j][k] - ivel[i][k]; + r2 += dx[k]*dx[k]; + xv += dx[k]*dv[k]; + } + r2i = 1/(r2+eps2+_TINY_); + ri = sqrt(r2i); + mri = jdat->mass[j]*ri; + mr3i = mri*r2i; + a3 = -3*xv*r2i; + // PRC(jdat->mpi_rank); PRC(ri); PRL(mri); + if (r2 > _TINY_) { + lpot[i] -= mri; + if (r2 < ldnn[i]) { + ldnn[i] = r2; + lnn[i] = j; } } + for (int k = 0; k < 3; k++) { + lacc[i][k] += mr3i*dx[k]; + ljerk[i][k] += mr3i*(dv[k]+a3*dx[k]); + } } ldnn[i] = sqrt(ldnn[i]); } diff --git a/src/amuse_ph4/src/jdata.cc b/src/amuse_ph4/src/jdata.cc index 9acef22d21..cf227ef39b 100644 --- a/src/amuse_ph4/src/jdata.cc +++ b/src/amuse_ph4/src/jdata.cc @@ -139,6 +139,11 @@ int jdata::add_particle(real pmass, real pradius, real dt) // default = -1 { const char *in_function = "jdata::add_particle"; + + if (pmass <= _TINY_) { + return -1; + } + if (DEBUG > 2 && mpi_rank == 0) PRL(in_function); if (nj >= njbuf) {