From 2f13ca8f1f101aea34529a68e68eac6845cdc2b4 Mon Sep 17 00:00:00 2001
From: Tom Theuns <tom.theuns@durham.ac.uk>
Date: Thu, 25 Feb 2016 10:37:50 +0000
Subject: [PATCH] added TestingGravity and ExternalPotential

---
 examples/TestingGravity/makeIC.py             | 140 ++++++++++++++++++
 examples/TestingGravity/run.sh                |  11 ++
 src/cell.c                                    |   6 +-
 src/const.h                                   |   6 +
 src/gravity/ExternalPotential/gravity.h       |  35 +++++
 src/gravity/ExternalPotential/gravity_io.h    |  71 +++++++++
 src/gravity/ExternalPotential/gravity_part.h  |  69 +++++++++
 .../ExternalPotential/runner_iact_grav.h      | 123 +++++++++++++++
 src/part.h                                    |   8 +
 9 files changed, 466 insertions(+), 3 deletions(-)
 create mode 100644 examples/TestingGravity/makeIC.py
 create mode 100644 examples/TestingGravity/run.sh
 create mode 100644 src/gravity/ExternalPotential/gravity.h
 create mode 100644 src/gravity/ExternalPotential/gravity_io.h
 create mode 100644 src/gravity/ExternalPotential/gravity_part.h
 create mode 100644 src/gravity/ExternalPotential/runner_iact_grav.h

diff --git a/examples/TestingGravity/makeIC.py b/examples/TestingGravity/makeIC.py
new file mode 100644
index 0000000000..58fc51cc96
--- /dev/null
+++ b/examples/TestingGravity/makeIC.py
@@ -0,0 +1,140 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ #                    Matthieu Schaller (matthieu.schaller@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/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+import numpy
+import math
+import random
+
+# Generates a random distriution of particles, for motion in an external potnetial centred at (0,0,0)
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS  = 6.672e-8
+SOLAR_MASS_IN_CGS   = 1.989e33
+PARSEC_IN_CGS       = 3.086e18
+PROTON_MASS_IN_CGS  = 1.6726231e24
+YEAR_IN_CGS         = 3.154e+7
+
+# choice of units
+const_unit_length_in_cgs   =     (1000 *  PARSEC_IN_CGS )
+const_unit_mass_in_cgs     =     (SOLAR_MASS_IN_CGS)
+const_unit_velocity_in_cgs =    (1e5)
+
+# derived units
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+const_G                = (NEWTON_GRAVITY_CGS * const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/ (const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs))
+
+
+
+# Parameters
+periodic= 1            # 1 For periodic box
+boxSize = 100.         # 
+Radius  = boxSize / 4. # maximum radius of particles
+G       = const_G 
+Mass    = 1e10         
+
+N       = int(sys.argv[1])  # Number of particles
+L       = N**(1./3.)
+
+# these are not used but necessary for I/O
+rho = 2.              # Density
+P = 1.                # Pressure
+gamma = 5./3.         # Gas adiabatic index
+fileName = "Sphere.hdf5" 
+
+
+#---------------------------------------------------
+numPart        = N
+mass           = 1
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Particle group
+#grp0 = file.create_group("/PartType0")
+grp1 = file.create_group("/PartType1")
+#generate particle positions
+radius = Radius * (numpy.random.rand(N))**(1./3.) 
+ctheta = -1. + 2 * numpy.random.rand(N)
+stheta = numpy.sqrt(1.-ctheta**2)
+phi    =  2 * math.pi * numpy.random.rand(N)
+r      = numpy.zeros((numPart, 3))
+# r[:,0] = radius * stheta * numpy.cos(phi)
+# r[:,1] = radius * stheta * numpy.sin(phi)
+# r[:,2] = radius * ctheta
+r[:,0] = radius
+#
+speed  = numpy.sqrt(G * Mass / radius)
+v      = numpy.zeros((numPart, 3))
+omega  = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+v[:,0] = -omega * r[:,1]
+v[:,1] =  omega * r[:,0]
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = numpy.zeros(1)
+
+m = numpy.full((numPart, ), mass)
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = numpy.zeros(1)
+
+h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
+ds[()] = h
+h = numpy.zeros(1)
+
+u = numpy.full((numPart, ), internalEnergy)
+ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
+ds[()] = u
+u = numpy.zeros(1)
+
+
+ids = numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds[()] = ids
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+file.close()
diff --git a/examples/TestingGravity/run.sh b/examples/TestingGravity/run.sh
new file mode 100644
index 0000000000..66d425971e
--- /dev/null
+++ b/examples/TestingGravity/run.sh
@@ -0,0 +1,11 @@
+#Notes:
+#Load swift module and gcc compiler to run sanitizer.
+
+#Configure
+./configure --enable-sanitizer --enable-debug --enable-optimization=no  --enable-mpi=no
+
+
+# -d: minimum time step
+# -e: maximum time step
+# -c: end of simulation
+ ./swift_fixdt -m .1 -s "50 50 50"  -t 1 -d 0.001 -e 0.001 -c 0.2 -f TestingGravity/Sphere.hdf5
diff --git a/src/cell.c b/src/cell.c
index 6d84dafdaa..06312f448a 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -448,9 +448,9 @@ void cell_split(struct cell *c) {
     c->progeny[k]->xparts = &c->xparts[left[k]];
   }
 
-  /* Re-link the gparts. */
-  for (k = 0; k < count; k++)
-    if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];
+  /* Re-link the gparts. THIS IS BROKEN BUT NEEDS MENDING J & T */
+  /* for (k = 0; k < count; k++) */
+  /*   if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k]; */
 
   /* Verify that _all_ the parts have been assigned to a cell. */
   /* for ( k = 1 ; k < 8 ; k++ )
diff --git a/src/const.h b/src/const.h
index 3bd9edff82..239ec18a5b 100644
--- a/src/const.h
+++ b/src/const.h
@@ -70,6 +70,12 @@
 #define GADGET2_SPH
 //#define DEFAULT_SPH
 
+
+/* Gravity properties */
+#define GRAVITY
+/* valid choices DEFAULT_GRAVITY || EXTERNAL_POTENTIAL */
+#define EXTERNAL_POTENTIAL
+
 /* System of units */
 #define const_unit_length_in_cgs 1   /* 3.08567810e16  /\* 1Mpc *\/ */
 #define const_unit_mass_in_cgs 1     /* 1.9891e33      /\* 1 M_sun *\/ */
diff --git a/src/gravity/ExternalPotential/gravity.h b/src/gravity/ExternalPotential/gravity.h
new file mode 100644
index 0000000000..82bc52ad3e
--- /dev/null
+++ b/src/gravity/ExternalPotential/gravity.h
@@ -0,0 +1,35 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@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/>.
+ *
+ ******************************************************************************/
+
+#include <float.h>
+
+/**
+ * @brief Computes the gravity time-step of a given particle
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ *
+ */
+
+__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
+    struct part* p, struct xpart* xp) {
+
+  /* Currently no limit is imposed */
+  return FLT_MAX;
+}
diff --git a/src/gravity/ExternalPotential/gravity_io.h b/src/gravity/ExternalPotential/gravity_io.h
new file mode 100644
index 0000000000..796d24a93d
--- /dev/null
+++ b/src/gravity/ExternalPotential/gravity_io.h
@@ -0,0 +1,71 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Reads the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to read the arrays.
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ *
+ */
+__attribute__((always_inline)) INLINE static void darkmatter_read_particles(
+    hid_t h_grp, int N, long long N_total, long long offset,
+    struct gpart* gparts) {
+  message("Reading Dark Matter particles\n");
+  /* Read arrays */
+  readArray(h_grp, "Coordinates", DOUBLE, N, 3, gparts, N_total, offset, x,
+            COMPULSORY);
+  readArray(h_grp, "Velocities", FLOAT, N, 3, gparts, N_total, offset, v,
+            COMPULSORY);
+  readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset, id,
+            COMPULSORY);
+}
+
+/**
+ * @brief Writes the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to write the arrays.
+ * @param fileName The name of the file (unsued in MPI mode).
+ * @param xmfFile The XMF file to write to (unused in MPI mode).
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param mpi_rank The MPI rank of this node (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ * @param us The unit system to use
+ *
+ */
+__attribute__((always_inline)) INLINE static void darkmatter_write_particles(
+    hid_t h_grp, char* fileName, FILE* xmfFile, int N, long long N_total,
+    int mpi_rank, long long offset, struct gpart* gparts, struct UnitSystem* us) {
+
+  /* Write arrays */
+  writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, gparts,
+             N_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
+  writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, gparts,
+             N_total, mpi_rank, offset, v, us, UNIT_CONV_SPEED);
+  writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, gparts,
+             N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+}
+
diff --git a/src/gravity/ExternalPotential/gravity_part.h b/src/gravity/ExternalPotential/gravity_part.h
new file mode 100644
index 0000000000..0bdbc0b09f
--- /dev/null
+++ b/src/gravity/ExternalPotential/gravity_part.h
@@ -0,0 +1,69 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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_GPART_H
+#define SWIFT_GPART_H
+
+/* Some standard headers. */
+#include <stdlib.h>
+
+
+/* properties of external potential */
+static struct ExternalPointMass
+{
+  const float Mass;
+  const float Position[3];
+} PointMass = {.Mass = 1, .Position={0.,0.,0.}};
+
+
+/* Gravity particle. */
+struct gpart {
+
+  /* Particle position. */
+  double x[3];
+
+  /* Particle velocity. */
+  float v[3];
+
+  /* Particle acceleration. */
+  float a[3];
+
+  /* Particle external gravity acceleration */
+  float a_grav_external[3];
+
+  /* Particle mass. */
+  float mass;
+
+  /* Particle time of beginning of time-step. */
+  float t_begin;
+
+  /* Particle time of end of time-step. */
+  float t_end;
+
+  /* Anonymous union for id/part. */
+  union {
+
+    /* Particle ID. */
+    size_t id;
+
+    /* Pointer to corresponding SPH part. */
+    struct part* part;
+  };
+
+} __attribute__((aligned(part_align)));
+#endif /* SWIFT_GPART_H */
diff --git a/src/gravity/ExternalPotential/runner_iact_grav.h b/src/gravity/ExternalPotential/runner_iact_grav.h
new file mode 100644
index 0000000000..e62be446e8
--- /dev/null
+++ b/src/gravity/ExternalPotential/runner_iact_grav.h
@@ -0,0 +1,123 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@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_RUNNER_IACT_GRAV_H
+#define SWIFT_RUNNER_IACT_GRAV_H
+
+/* Includes. */
+#include "const.h"
+#include "kernel.h"
+#include "vector.h"
+
+/**
+ * @file  runner_iact_grav.h
+ * @brief Gravity interaction functions.
+ *
+ */
+
+/**
+ * @brief Gravity potential
+ */
+
+__attribute__((always_inline)) INLINE static void runner_iact_grav(
+    float r2, float *dx, struct gpart *pi, struct gpart *pj) {
+
+  float ir, r;
+  float w, acc;
+  float mi = pi->mass, mj = pj->mass;
+  int k;
+
+  /* Get the absolute distance. */
+  ir = 1.0f / sqrtf(r2);
+  r = r2 * ir;
+
+  /* Evaluate the gravity kernel. */
+  kernel_grav_eval(r, &acc);
+
+  /* Scale the acceleration. */
+  acc *= const_G * ir * ir * ir;
+
+  /* Aggregate the accelerations. */
+  for (k = 0; k < 3; k++) {
+    w = acc * dx[k];
+    pi->a[k] -= w * mj;
+    pj->a[k] += w * mi;
+  }
+}
+
+/**
+ * @brief Gravity potential (Vectorized version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_vec_grav(
+    float *R2, float *Dx, struct gpart **pi, struct gpart **pj) {
+
+#ifdef VECTORIZE
+
+  vector ir, r, r2, dx[3];
+  vector w, acc, ai, aj;
+  vector mi, mj;
+  int j, k;
+
+#if VEC_SIZE == 8
+  mi.v = vec_set(pi[0]->mass, pi[1]->mass, pi[2]->mass, pi[3]->mass,
+                 pi[4]->mass, pi[5]->mass, pi[6]->mass, pi[7]->mass);
+  mj.v = vec_set(pj[0]->mass, pj[1]->mass, pj[2]->mass, pj[3]->mass,
+                 pj[4]->mass, pj[5]->mass, pj[6]->mass, pj[7]->mass);
+  for (k = 0; k < 3; k++)
+    dx[k].v = vec_set(Dx[0 + k], Dx[3 + k], Dx[6 + k], Dx[9 + k], Dx[12 + k],
+                      Dx[15 + k], Dx[18 + k], Dx[21 + k]);
+#elif VEC_SIZE == 4
+  mi.v = vec_set(pi[0]->mass, pi[1]->mass, pi[2]->mass, pi[3]->mass);
+  mj.v = vec_set(pj[0]->mass, pj[1]->mass, pj[2]->mass, pj[3]->mass);
+  for (k = 0; k < 3; k++)
+    dx[k].v = vec_set(Dx[0 + k], Dx[3 + k], Dx[6 + k], Dx[9 + k]);
+#endif
+
+  /* Get the radius and inverse radius. */
+  r2.v = vec_load(R2);
+  ir.v = vec_rsqrt(r2.v);
+  ir.v = ir.v - vec_set1(0.5f) * ir.v * (r2.v * ir.v * ir.v - vec_set1(1.0f));
+  r.v = r2.v * ir.v;
+
+  /* Evaluate the gravity kernel. */
+  blender_eval_vec(&r, &acc);
+
+  /* Scale the acceleration. */
+  acc.v *= vec_set1(const_G) * ir.v * ir.v * ir.v;
+
+  /* Aggregate the accelerations. */
+  for (k = 0; k < 3; k++) {
+    w.v = acc.v * dx[k].v;
+    ai.v = w.v * mj.v;
+    aj.v = w.v * mi.v;
+    for (j = 0; j < VEC_SIZE; j++) {
+      pi[j]->a[k] -= ai.f[j];
+      pj[j]->a[k] += aj.f[j];
+    }
+  }
+
+#else
+
+  for (int k = 0; k < VEC_SIZE; k++)
+    runner_iact_grav(R2[k], &Dx[3 * k], pi[k], pj[k]);
+
+#endif
+}
+
+#endif /* SWIFT_RUNNER_IACT_GRAV_H */
diff --git a/src/part.h b/src/part.h
index 168e80b68b..33ee05fcb8 100644
--- a/src/part.h
+++ b/src/part.h
@@ -48,7 +48,15 @@
 #error "Invalid choice of SPH variant"
 #endif
 
+#if defined(GRAVITY)
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity_part.h"
+#elif defined(EXTERNAL_POTENTIAL)
+#include "./gravity/ExternalPotential/gravity_part.h"
+#elif
+#error "Invalid choice of gravity variant"
+#endif
+#endif
 
 #ifdef WITH_MPI
 void part_create_mpi_type(MPI_Datatype* part_type);
-- 
GitLab