1207 fix

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<int,int>::iterator iterator = mapping_from_address_j_to_index_in_update_array.find(address);
     map<int,int>::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);
+  }
 
 }
 

lib/sapporo_light/sapporo.h

@@ -22,6 +22,15 @@
 #include <map>
 
 #include <math.h>
+
+#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 <cuda_runtime.h>

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 */

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

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]);
     }

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) {