From c40f411dc9e83d2834dde29c8bc06c19579f3a2f Mon Sep 17 00:00:00 2001
From: James Willis <james.s.willis@durham.ac.uk>
Date: Thu, 19 Apr 2018 20:26:39 +0800
Subject: [PATCH] Added a naive implementation of a 3D FOF.
---
src/Makefile.am | 4 +-
src/fof.c | 112 ++++++++++++++++++++++++++++++++++++++++++++++++
src/fof.h | 39 +++++++++++++++++
3 files changed, 153 insertions(+), 2 deletions(-)
create mode 100644 src/fof.c
create mode 100644 src/fof.h
diff --git a/src/Makefile.am b/src/Makefile.am
index 1bd80b6693..36f2921025 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -47,7 +47,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
- chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h
+ chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h fof.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -59,7 +59,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c \
collectgroup.c hydro_space.c equation_of_state.c \
- chemistry.c cosmology.c restart.c
+ chemistry.c cosmology.c restart.c fof.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/fof.c b/src/fof.c
new file mode 100644
index 0000000000..913f81fa59
--- /dev/null
+++ b/src/fof.c
@@ -0,0 +1,112 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "fof.h"
+
+/* Local headers. */
+//#include "active.h"
+
+void fof_search(struct engine *e) {
+
+ const size_t nr_gparts = e->s->nr_gparts;
+ struct gpart *gparts = e->s->gparts;
+ const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+ const double l_x = 0.2 * (dim[0] / pow(nr_gparts, 1./3.));
+ const double l_x2 = l_x * l_x;
+ int *pid;
+ int *num_in_groups;
+ int num_groups = nr_gparts;
+
+ message("Searching %ld gravity particles for links with l_x2: %lf", nr_gparts, l_x2);
+
+ /* Allocate and populate array of particle group IDs. */
+ if (posix_memalign((void **)&pid, 32,
+ nr_gparts * sizeof(int)) != 0)
+ error("Failed to allocate list of particle group IDs for FOF search.");
+
+ for(size_t i=0; i<nr_gparts; i++) pid[i] = i;
+
+ if (posix_memalign((void **)&num_in_groups, 32,
+ nr_gparts * sizeof(int)) != 0)
+ error("Failed to allocate list of number in groups for FOF search.");
+
+ for(size_t i=0; i<nr_gparts; i++) num_in_groups[i] = 1;
+
+ /* Loop over particles and find which particles belong in the same group. */
+ for(size_t i=0; i<nr_gparts; i++) {
+
+ //message("Searching for particle: %ld groups.", i);
+
+ struct gpart *pi = &gparts[i];
+ const double pix = pi->x[0];
+ const double piy = pi->x[1];
+ const double piz = pi->x[2];
+
+ for(size_t j=0; j<nr_gparts; j++) {
+
+ /* Skip yourself. */
+ if(i == j) continue;
+
+ struct gpart *pj = &gparts[j];
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute pairwise distance, remembering to account for boundary conditions. */
+ float dx[3], r2 = 0.0f;
+ dx[0] = pix - pjx;
+ dx[1] = piy - pjy;
+ dx[2] = piz - pjz;
+
+ for (int k = 0; k < 3; k++) {
+ dx[k] = nearest(dx[k], dim[k]);
+ r2 += dx[k] * dx[k];
+ }
+
+ /* Hit or miss? */
+ if (r2 < l_x2 && pid[j] < pid[i]) {
+ pid[i] = pid[j];
+ num_in_groups[i]--;
+ num_in_groups[j]++;
+ num_groups--;
+ }
+ }
+ }
+
+ int num_parts_in_groups = 0;
+ int max_group_size = 0, max_group_id = 0;
+ for(size_t i=0; i<nr_gparts; i++) {
+
+ if(num_in_groups[i] > 1) num_parts_in_groups += num_in_groups[i];
+ if( num_in_groups[i] > max_group_size) {
+ max_group_size = num_in_groups[i];
+ max_group_id = i;
+ }
+ }
+
+ message("No. of groups: %d. No. of particles in groups: %d. No. of particles not in groups: %d.", num_groups, num_parts_in_groups, nr_gparts - num_parts_in_groups);
+ message("Biggest group size: %d with ID: %d", max_group_size, max_group_id);
+
+ free(pid);
+ free(num_in_groups);
+}
diff --git a/src/fof.h b/src/fof.h
new file mode 100644
index 0000000000..1d3d0aa231
--- /dev/null
+++ b/src/fof.h
@@ -0,0 +1,39 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_FOF_H
+#define SWIFT_FOF_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers */
+//#include "active.h"
+//#include "cell.h"
+#include "engine.h"
+//#include "hydro.h"
+//#include "part.h"
+//#include "runner.h"
+//#include "timers.h"
+//#include "vector.h"
+
+/* Function prototypes. */
+void fof_search(struct engine *e);
+
+#endif /* SWIFT_FOF_H */
--
GitLab