Added collision handler for Huayno OK algorithm

src/amuse/community/huayno/src/Makefile

@@ -26,7 +26,8 @@ RM = rm
 
 OBJS = evolve.o evolve_shared.o evolve_sf.o evolve_cc.o \
   evolve_ok.o evolve_kepler.o universal_variable_kepler.o evolve_bs.o \
-  evolve_shared_collisions.o evolve_error_control.o simple_map.o simple_hash.o
+  evolve_shared_collisions.o evolve_error_control.o simple_map.o simple_hash.o \
+  collision_check.o
 
 all: libhuayno.a
 

src/amuse/community/huayno/src/collision_check.c

@@ -0,0 +1,42 @@
+/*
+ * Algorithm to detect collisions between particles
+ */
+#include "evolve.h"
+#include "evolve_ok.h"
+#include "collision_check.h"
+// AMUSE STOPPING CONDITIONS SUPPORT
+#include <stopcond.h>
+
+
+void detect_collisions(struct sys s){
+    UINT i, j;
+    FLOAT dx[3], dr2, radius_sum;
+    struct particle *ipart, *jpart;
+    #pragma omp parallel for if((ULONG) s.n*s.n>MPWORKLIMIT && !omp_in_parallel()) default(none) \
+    private(i,j,dx,dr2,radius_sum, ipart, jpart) \
+    shared(s)
+    for (i=0; i<s.n; i++) {
+        ipart=GETPART(s,i);
+        if (ipart->mass > 0){
+          for (j=i+1; j<s.n; j++) {
+              jpart=GETPART(s,j);
+              dx[0] = ipart->pos[0] - jpart->pos[0];
+              dx[1] = ipart->pos[1] - jpart->pos[1];
+              dx[2] = ipart->pos[2] - jpart->pos[2];
+              dr2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
+              radius_sum = ipart->radius + jpart->radius;
+              if (dr2 <= radius_sum*radius_sum) {
+                  #pragma omp critical
+                  {
+                      int stopping_index = next_index_for_stopping_condition();
+                      if (stopping_index >= 0) {
+                          set_stopping_condition_info(stopping_index, COLLISION_DETECTION);
+                          set_stopping_condition_particle_index(stopping_index, 0, ipart->id);
+                          set_stopping_condition_particle_index(stopping_index, 1, jpart->id);
+                      }
+                  }
+              }
+          }
+        }
+    }
+}

src/amuse/community/huayno/src/collision_check.h

@@ -0,0 +1,6 @@
+#ifndef COLLISION_CHECK_H
+#define COLLISION_CHECK_H
+
+void detect_collisions(struct sys s);
+
+#endif // COLLISION_CHECK_H

src/amuse/community/huayno/src/evolve_ok.c

@@ -7,8 +7,13 @@
 #include <tgmath.h>
 #include <stdio.h>
 #include <stdlib.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 #include "evolve.h"
 #include "evolve_ok.h"
+// AMUSE STOPPING CONDITIONS SUPPORT
+#include <stopcond.h>
 
 struct forces zeroforces = {0, NULL};
 
@@ -150,6 +155,8 @@ static void ok_kick(int clevel,struct forces f, DOUBLE dt)
 
 void evolve_ok2(int clevel,struct sys s, struct forces f, DOUBLE stime, DOUBLE etime, DOUBLE dt, int calc_timestep)
 {
+  int is_collision_detection_enabled;
+  is_stopping_condition_enabled(COLLISION_DETECTION, &is_collision_detection_enabled);
   if (IS_ZEROFORCES(f) && clevel == 0) { f = ok_main_forces; }
   CHECK_TIMESTEP(etime,stime,dt,clevel);
   // all particles are drifted together
@@ -165,5 +172,9 @@ void evolve_ok2(int clevel,struct sys s, struct forces f, DOUBLE stime, DOUBLE e
   ok_split((FLOAT) dt, f, &slowf, &fastf);
   evolve_ok2(clevel+1,s, fastf, stime, stime+dt/2, dt/2, 0);
   ok_kick(clevel,slowf, dt);
+  if (is_collision_detection_enabled) {
+    detect_collisions(s);
+    if (set_conditions & enabled_conditions) return;
+  }
   evolve_ok2(clevel+1,s, fastf, stime+dt/2, etime, dt/2, 1);
 }

src/amuse/community/huayno/src/evolve_shared_collisions.c

@@ -10,6 +10,7 @@
 #endif
 #include "evolve.h"
 #include "integrators_shared.h"
+#include "collision_check.h"
 // AMUSE STOPPING CONDITIONS SUPPORT
 #include <stopcond.h>
 
@@ -43,37 +44,6 @@ static int update_steps_and_get_next_level(int *steps, int current_level) {
     return current_level;
 }
 
-static void detect_collisions(struct sys s) {
-    UINT i, j;
-    FLOAT dx[3], dr2, radius_sum;
-    struct particle *ipart, *jpart;
-#pragma omp parallel for if((ULONG) s.n*s.n>MPWORKLIMIT && !omp_in_parallel()) default(none) \
- private(i,j,dx,dr2,radius_sum, ipart, jpart) \
- shared(s)
-    for (i=0; i<s.n; i++) {
-        ipart=GETPART(s,i);
-        for (j=i+1; j<s.n; j++) {
-            jpart=GETPART(s,j);
-            dx[0] = ipart->pos[0] - jpart->pos[0];
-            dx[1] = ipart->pos[1] - jpart->pos[1];
-            dx[2] = ipart->pos[2] - jpart->pos[2];
-            dr2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
-            radius_sum = ipart->radius + jpart->radius;
-            if (dr2 <= radius_sum*radius_sum) {
-                #pragma omp critical
-                {
-                    int stopping_index = next_index_for_stopping_condition();
-                    if (stopping_index >= 0) {
-                        set_stopping_condition_info(stopping_index, COLLISION_DETECTION);
-                        set_stopping_condition_particle_index(stopping_index, 0, ipart->id);
-                        set_stopping_condition_particle_index(stopping_index, 1, jpart->id);
-                    }
-                }
-            }
-        }
-    }
-}
-
 static void evolve_shared_collision_detection(struct sys s, DOUBLE dt, void (*dkd_func)(int, struct sys, DOUBLE, DOUBLE, DOUBLE)) {
     FLOAT dtsys;
     int next_level, current_level = 0;
@@ -133,4 +103,3 @@ void evolve_shared8_collision_detection(struct sys s, DOUBLE dt) {
 void evolve_shared10_collision_detection(struct sys s, DOUBLE dt) {
     evolve_shared_collision_detection(s, dt, dkd10);
 }
-