diff --git a/.gitignore b/.gitignore
index 5564b9e1e00c735e6f27935f8443decd7c61647b..e883fac300bb92a5fade0c778b513e709cbd2854 100644
--- a/.gitignore
+++ b/.gitignore
@@ -77,6 +77,7 @@ tests/testRiemannExact
tests/testRiemannTRRS
tests/testRiemannHLLC
tests/testMatrixInversion
+tests/testVoronoi1D
theory/latex/swift.pdf
theory/kernel/kernels.pdf
diff --git a/src/const.h b/src/const.h
index 4d510ceffad23dbf0a50d4a871af95b494e1fdd0..e9d64e384c1e0db4a7fa556bd67658706695d1c1 100644
--- a/src/const.h
+++ b/src/const.h
@@ -37,9 +37,9 @@
#define const_max_u_change 0.1f
/* Dimensionality of the problem */
-#define HYDRO_DIMENSION_3D
+//#define HYDRO_DIMENSION_3D
//#define HYDRO_DIMENSION_2D
-//#define HYDRO_DIMENSION_1D
+#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
@@ -62,7 +62,8 @@
//#define MINIMAL_SPH
//#define GADGET2_SPH
//#define DEFAULT_SPH
-#define GIZMO_SPH
+//#define GIZMO_SPH
+#define SHADOWSWIFT
/* Self gravity stuff. */
#define const_gravity_multipole_order 2
diff --git a/src/debug.c b/src/debug.c
index ca77745ec9f905d85403445ff3d8ef2de16034e7..cbe1eecf5299a652aa50d2fb859d0dac9ea1a8ae 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -45,6 +45,8 @@
#include "./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_debug.h"
+#elif defined(SHADOWSWIFT)
+#include "./hydro/Shadowswift/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/engine.c b/src/engine.c
index dcdd956250045d6274ebb1ee2cbd5c43be216e83..003e6b3d42e40ffef5dc3250116dd0c99440f0ed 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2720,7 +2720,7 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_sub_pair;
mask |= 1 << task_type_ghost;
mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask
- does not harm */
+ does not harm */
submask |= 1 << task_subtype_density;
submask |= 1 << task_subtype_gradient;
diff --git a/src/hydro.h b/src/hydro.h
index 4a2b0bd029d494e2091b9081d22b7949cec5648c..9ca1bfd6e71eb372f763f939b93db454f06ad12c 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -42,6 +42,11 @@
#include "./hydro/Gizmo/hydro.h"
#include "./hydro/Gizmo/hydro_iact.h"
#define SPH_IMPLEMENTATION "GIZMO (Hopkins 2015)"
+#elif defined(SHADOWSWIFT)
+#include "./hydro/Shadowswift/hydro.h"
+#include "./hydro/Shadowswift/hydro_iact.h"
+#define SPH_IMPLEMENTATION \
+ "Shadowfax moving mesh (Vandenbroucke & De Rijcke 2016)"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a05ae68495158ff6a4a4e726eb23acf24add39c
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro.h
@@ -0,0 +1,428 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * 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>
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "hydro_gradients.h"
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ * @param hydro_properties Pointer to the hydro parameters.
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ const struct part* restrict p, const struct xpart* restrict xp,
+ const struct hydro_props* restrict hydro_properties) {
+
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * In this case, we copy the particle velocities into the corresponding
+ * primitive variable field. We do this because the particle velocities in GIZMO
+ * can be independent of the actual fluid velocity. The latter is stored as a
+ * primitive variable and integrated using the linear momentum, a conserved
+ * variable.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part* p, struct xpart* xp) {
+
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+
+ p->primitives.v[0] = p->v[0];
+ p->primitives.v[1] = p->v[1];
+ p->primitives.v[2] = p->v[2];
+}
+
+/**
+ * @brief Prepares a particle for the volume calculation.
+ *
+ * Simply makes sure all necessary variables are initialized to zero.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part* p) {
+
+ p->density.wcount = 0.0f;
+ p->density.wcount_dh = 0.0f;
+ p->rho = 0.0f;
+ p->geometry.volume = 0.0f;
+ p->geometry.matrix_E[0][0] = 0.0f;
+ p->geometry.matrix_E[0][1] = 0.0f;
+ p->geometry.matrix_E[0][2] = 0.0f;
+ p->geometry.matrix_E[1][0] = 0.0f;
+ p->geometry.matrix_E[1][1] = 0.0f;
+ p->geometry.matrix_E[1][2] = 0.0f;
+ p->geometry.matrix_E[2][0] = 0.0f;
+ p->geometry.matrix_E[2][1] = 0.0f;
+ p->geometry.matrix_E[2][2] = 0.0f;
+}
+
+/**
+ * @brief Finishes the volume calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and adds the self-contribution term. Calculates the volume and uses it to
+ * update the primitive variables (based on the conserved variables). The latter
+ * should only be done for active particles. This is okay, since this method is
+ * only called for active particles.
+ *
+ * Multiplies the components of the matrix E with the appropriate constants and
+ * inverts it. Initializes the variables used during the gradient loop. This
+ * cannot be done in hydro_prepare_force, since that method is called for all
+ * particles, and not just the active ones. If we would initialize the
+ * variables there, gradients for passive particles would be zero, while we
+ * actually use the old gradients in the flux calculation between active and
+ * passive particles.
+ *
+ * @param p The particle to act upon.
+ * @param The current physical time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part* restrict p, float time) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float ih = 1.0f / h;
+
+ /* Final operation on the density. */
+ p->density.wcount += kernel_root;
+ p->density.wcount *= kernel_norm;
+
+ p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
+
+ const float ihdim = pow_dimension(ih);
+
+ p->rho += p->mass * kernel_root;
+ p->rho *= ihdim;
+
+ float volume;
+ float m, momentum[3], energy;
+
+ /* Final operation on the geometry. */
+ /* we multiply with the smoothing kernel normalization ih3 and calculate the
+ * volume */
+ volume = ihdim * (p->geometry.volume + kernel_root);
+ p->geometry.volume = volume = 1. / volume;
+ /* we multiply with the smoothing kernel normalization */
+ p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
+ p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
+ p->geometry.matrix_E[0][2] = ihdim * p->geometry.matrix_E[0][2];
+ p->geometry.matrix_E[1][0] = ihdim * p->geometry.matrix_E[1][0];
+ p->geometry.matrix_E[1][1] = ihdim * p->geometry.matrix_E[1][1];
+ p->geometry.matrix_E[1][2] = ihdim * p->geometry.matrix_E[1][2];
+ p->geometry.matrix_E[2][0] = ihdim * p->geometry.matrix_E[2][0];
+ p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
+ p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
+
+ invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
+
+ hydro_gradients_init(p);
+
+ /* compute primitive variables */
+ /* eqns (3)-(5) */
+ m = p->conserved.mass;
+ if (m) {
+ momentum[0] = p->conserved.momentum[0];
+ momentum[1] = p->conserved.momentum[1];
+ momentum[2] = p->conserved.momentum[2];
+ p->primitives.rho = m / volume;
+ p->primitives.v[0] = momentum[0] / m;
+ p->primitives.v[1] = momentum[1] / m;
+ p->primitives.v[2] = momentum[2] / m;
+ energy = p->conserved.energy;
+ p->primitives.P = hydro_gamma_minus_one * energy / volume;
+ }
+}
+
+/**
+ * @brief Prepare a particle for the gradient calculation.
+ *
+ * The name of this method is confusing, as this method is really called after
+ * the density loop and before the gradient loop.
+ *
+ * We use it to set the physical timestep for the particle and to copy the
+ * actual velocities, which we need to boost our interfaces during the flux
+ * calculation. We also initialize the variables used for the time step
+ * calculation.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param ti_current Current integer time.
+ * @param timeBase Conversion factor between integer time and physical time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part* restrict p, struct xpart* restrict xp, int ti_current,
+ double timeBase) {
+
+ /* Set the physical time step */
+ p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
+
+ /* Initialize time step criterion variables */
+ p->timestepvars.vmax = 0.0f;
+
+ /* Set the actual velocity of the particle */
+ p->force.v_full[0] = xp->v_full[0];
+ p->force.v_full[1] = xp->v_full[1];
+ p->force.v_full[2] = xp->v_full[2];
+}
+
+/**
+ * @brief Finishes the gradient calculation.
+ *
+ * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
+ * in which case no gradients are used.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+ struct part* p) {
+
+ hydro_gradients_finalize(p);
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is actually not necessary for GIZMO, since we just set the accelerations
+ * after the flux calculation.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part* p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->force.h_dt = 0.0f;
+}
+
+/**
+ * @brief Converts the hydrodynamic variables from the initial condition file to
+ * conserved variables that can be used during the integration
+ *
+ * Requires the volume to be known.
+ *
+ * The initial condition file contains a mixture of primitive and conserved
+ * variables. Mass is a conserved variable, and we just copy the particle
+ * mass into the corresponding conserved quantity. We need the volume to
+ * also derive a density, which is then used to convert the internal energy
+ * to a pressure. However, we do not actually use these variables anymore.
+ * We do need to initialize the linear momentum, based on the mass and the
+ * velocity of the particle.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p) {
+
+ float volume;
+ float m;
+ float momentum[3];
+ volume = p->geometry.volume;
+
+ p->conserved.mass = m = p->mass;
+ p->primitives.rho = m / volume;
+
+ /* P actually contains internal energy at this point */
+ p->primitives.P *= hydro_gamma_minus_one * p->primitives.rho;
+
+ p->conserved.momentum[0] = momentum[0] = m * p->primitives.v[0];
+ p->conserved.momentum[1] = momentum[1] = m * p->primitives.v[1];
+ p->conserved.momentum[2] = momentum[2] = m * p->primitives.v[2];
+ p->conserved.energy = p->primitives.P / hydro_gamma_minus_one * volume;
+}
+
+/**
+ * @brief Extra operations to be done during the drift
+ *
+ * Not used for GIZMO.
+ *
+ * @param p Particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param t0 Integer start time of the drift interval.
+ * @param t1 Integer end time of the drift interval.
+ * @param timeBase Conversion factor between integer and physical time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
+
+ // return;
+ float dt = (t1 - t0) * timeBase;
+ float h_inv = 1.0f / p->h;
+ float w = -hydro_dimension * p->force.h_dt * h_inv * dt;
+
+ if (fabsf(w) < 0.2f) {
+ p->primitives.rho *= approx_expf(w);
+ } else {
+ p->primitives.rho *= expf(w);
+ }
+ p->primitives.v[0] += p->a_hydro[0] * dt;
+ p->primitives.v[1] += p->a_hydro[1] * dt;
+ p->primitives.v[2] += p->a_hydro[2] * dt;
+ float u = p->conserved.energy + p->du_dt * dt;
+ p->primitives.P =
+ hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+}
+
+/**
+ * @brief Set the particle acceleration after the flux loop
+ *
+ * We use the new conserved variables to calculate the new velocity of the
+ * particle, and use that to derive the change of the velocity over the particle
+ * time step.
+ *
+ * If the particle time step is zero, we set the accelerations to zero. This
+ * should only happen at the start of the simulation.
+ *
+ * @param p Particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part* p) {
+
+ /* Add normalization to h_dt. */
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+
+ /* Update the conserved variables. We do this here and not in the kick,
+ since we need the updated variables below. */
+ p->conserved.mass += p->conserved.flux.mass;
+ p->conserved.momentum[0] += p->conserved.flux.momentum[0];
+ p->conserved.momentum[1] += p->conserved.flux.momentum[1];
+ p->conserved.momentum[2] += p->conserved.flux.momentum[2];
+ p->conserved.energy += p->conserved.flux.energy;
+
+ /* reset fluxes */
+ /* we can only do this here, since we need to keep the fluxes for inactive
+ particles */
+ p->conserved.flux.mass = 0.0f;
+ p->conserved.flux.momentum[0] = 0.0f;
+ p->conserved.flux.momentum[1] = 0.0f;
+ p->conserved.flux.momentum[2] = 0.0f;
+ p->conserved.flux.energy = 0.0f;
+
+ /* Set the hydro acceleration, based on the new momentum and mass */
+ /* NOTE: the momentum and mass are only correct for active particles, since
+ only active particles have received flux contributions from all their
+ neighbours. Since this method is only called for active particles,
+ this is indeed the case. */
+ if (p->force.dt) {
+ p->a_hydro[0] =
+ (p->conserved.momentum[0] / p->conserved.mass - p->primitives.v[0]) /
+ p->force.dt;
+ p->a_hydro[1] =
+ (p->conserved.momentum[1] / p->conserved.mass - p->primitives.v[1]) /
+ p->force.dt;
+ p->a_hydro[2] =
+ (p->conserved.momentum[2] / p->conserved.mass - p->primitives.v[2]) /
+ p->force.dt;
+
+ p->du_dt = p->conserved.flux.energy / p->force.dt;
+ } else {
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ p->du_dt = 0.0f;
+ }
+}
+
+/**
+ * @brief Extra operations done during the kick
+ *
+ * Not used for GIZMO.
+ *
+ * @param p Particle to act upon.
+ * @param xp Extended particle data to act upon.
+ * @param dt Physical time step.
+ * @param half_dt Half the physical time step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part* p, struct xpart* xp, float dt, float half_dt) {}
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param dt Time since the last kick.
+ * @return Internal energy of the particle.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part* restrict p, float dt) {
+
+ return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest.
+ * @param dt Time since the last kick.
+ * @return Entropy of the particle.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part* restrict p, float dt) {
+
+ return p->primitives.P / pow_gamma(p->primitives.rho);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest.
+ * @param dt Time since the last kick.
+ * @param Sound speed of the particle.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part* restrict p, float dt) {
+
+ return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ * @param Pressure of the particle.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part* restrict p, float dt) {
+
+ return p->primitives.P;
+}
diff --git a/src/hydro/Shadowswift/hydro_debug.h b/src/hydro/Shadowswift/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..db672bab16662961a99558b563b9320d4e75ecc4
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_debug.h
@@ -0,0 +1,83 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ printf(
+ "x=[%.16e,%.16e,%.16e], "
+ "v=[%.3e,%.3e,%.3e], "
+ "a=[%.3e,%.3e,%.3e], "
+ "h=%.3e, "
+ "ti_begin=%d, "
+ "ti_end=%d, "
+ "primitives={"
+ "v=[%.3e,%.3e,%.3e], "
+ "rho=%.3e, "
+ "P=%.3e, "
+ "gradients={"
+ "rho=[%.3e,%.3e,%.3e], "
+ "v=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]], "
+ "P=[%.3e,%.3e,%.3e]}, "
+ "limiter={"
+ "rho=[%.3e,%.3e], "
+ "v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
+ "P=[%.3e,%.3e], "
+ "maxr=%.3e}}, "
+ "conserved={"
+ "momentum=[%.3e,%.3e,%.3e], "
+ "mass=%.3e, "
+ "energy=%.3e}, "
+ "geometry={"
+ "volume=%.3e, "
+ "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
+ "timestepvars={"
+ "vmax=%.3e}, "
+ "density={"
+ "div_v=%.3e, "
+ "wcount_dh=%.3e, "
+ "curl_v=[%.3e,%.3e,%.3e], "
+ "wcount=%.3e}, "
+ "mass=%.3e\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
+ p->a_hydro[1], p->a_hydro[2], p->h, p->ti_begin, p->ti_end,
+ p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
+ p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
+ p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
+ p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
+ p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
+ p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
+ p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
+ p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
+ p->primitives.gradients.P[1], p->primitives.gradients.P[2],
+ p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
+ p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
+ p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
+ p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
+ p->primitives.limiter.P[0], p->primitives.limiter.P[1],
+ p->primitives.limiter.maxr, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
+ p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+ p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
+ p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
+ p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
+ p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
+ p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
+ p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
+ p->density.wcount, p->mass);
+}
diff --git a/src/hydro/Shadowswift/hydro_gradients.h b/src/hydro/Shadowswift/hydro_gradients.h
new file mode 100644
index 0000000000000000000000000000000000000000..90c09a0621d7feccc835c612512e699664fd7376
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_gradients.h
@@ -0,0 +1,211 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_HYDRO_GRADIENTS_H
+#define SWIFT_HYDRO_GRADIENTS_H
+
+//#define SPH_GRADIENTS
+#define GIZMO_GRADIENTS
+
+#include "hydro_slope_limiters.h"
+
+#if defined(SPH_GRADIENTS)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "SPH gradients (Price 2012)"
+#include "hydro_gradients_sph.h"
+
+#elif defined(GIZMO_GRADIENTS)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "GIZMO gradients (Hopkins 2015)"
+#include "hydro_gradients_gizmo.h"
+
+#else
+
+/* No gradients. Perfectly acceptable, but we have to provide empty functions */
+#define HYDRO_GRADIENT_IMPLEMENTATION "No gradients (first order scheme)"
+
+/**
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part* p) {}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float* dx, float hi, float hj, struct part* pi, struct part* pj) {
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, float* dx, float hi, float hj,
+ struct part* pi, struct part* pj) {}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part* p) {}
+
+#endif
+
+/**
+ * @brief Gradients reconstruction. Is the same for all gradient types (although
+ * gradients_none does nothing, since all gradients are zero -- are they?).
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
+ struct part* pi, struct part* pj, float hi, float hj, float* dx, float r,
+ float* xij_i, float* Wi, float* Wj, float mindt) {
+
+ float dWi[5], dWj[5];
+ float xij_j[3];
+ int k;
+ float xfac;
+
+ /* perform gradient reconstruction in space and time */
+ /* space */
+ /* Compute interface position (relative to pj, since we don't need the actual
+ * position) */
+ /* eqn. (8) */
+ xfac = hj / (hi + hj);
+ for (k = 0; k < 3; k++) xij_j[k] = xfac * dx[k];
+
+ dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
+ pi->primitives.gradients.rho[1] * xij_i[1] +
+ pi->primitives.gradients.rho[2] * xij_i[2];
+ dWi[1] = pi->primitives.gradients.v[0][0] * xij_i[0] +
+ pi->primitives.gradients.v[0][1] * xij_i[1] +
+ pi->primitives.gradients.v[0][2] * xij_i[2];
+ dWi[2] = pi->primitives.gradients.v[1][0] * xij_i[0] +
+ pi->primitives.gradients.v[1][1] * xij_i[1] +
+ pi->primitives.gradients.v[1][2] * xij_i[2];
+ dWi[3] = pi->primitives.gradients.v[2][0] * xij_i[0] +
+ pi->primitives.gradients.v[2][1] * xij_i[1] +
+ pi->primitives.gradients.v[2][2] * xij_i[2];
+ dWi[4] = pi->primitives.gradients.P[0] * xij_i[0] +
+ pi->primitives.gradients.P[1] * xij_i[1] +
+ pi->primitives.gradients.P[2] * xij_i[2];
+
+ dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
+ pj->primitives.gradients.rho[1] * xij_j[1] +
+ pj->primitives.gradients.rho[2] * xij_j[2];
+ dWj[1] = pj->primitives.gradients.v[0][0] * xij_j[0] +
+ pj->primitives.gradients.v[0][1] * xij_j[1] +
+ pj->primitives.gradients.v[0][2] * xij_j[2];
+ dWj[2] = pj->primitives.gradients.v[1][0] * xij_j[0] +
+ pj->primitives.gradients.v[1][1] * xij_j[1] +
+ pj->primitives.gradients.v[1][2] * xij_j[2];
+ dWj[3] = pj->primitives.gradients.v[2][0] * xij_j[0] +
+ pj->primitives.gradients.v[2][1] * xij_j[1] +
+ pj->primitives.gradients.v[2][2] * xij_j[2];
+ dWj[4] = pj->primitives.gradients.P[0] * xij_j[0] +
+ pj->primitives.gradients.P[1] * xij_j[1] +
+ pj->primitives.gradients.P[2] * xij_j[2];
+
+ hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
+
+ /* time */
+ dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+ Wi[2] * pi->primitives.gradients.rho[1] +
+ Wi[3] * pi->primitives.gradients.rho[2] +
+ Wi[0] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+ dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
+ Wi[2] * pi->primitives.gradients.v[0][1] +
+ Wi[3] * pi->primitives.gradients.v[0][2] +
+ pi->primitives.gradients.P[0] / Wi[0]);
+ dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
+ Wi[2] * pi->primitives.gradients.v[1][1] +
+ Wi[3] * pi->primitives.gradients.v[1][2] +
+ pi->primitives.gradients.P[1] / Wi[0]);
+ dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
+ Wi[2] * pi->primitives.gradients.v[2][1] +
+ Wi[3] * pi->primitives.gradients.v[2][2] +
+ pi->primitives.gradients.P[2] / Wi[0]);
+ dWi[4] -=
+ 0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
+ Wi[2] * pi->primitives.gradients.P[1] +
+ Wi[3] * pi->primitives.gradients.P[2] +
+ hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+
+ dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+ Wj[2] * pj->primitives.gradients.rho[1] +
+ Wj[3] * pj->primitives.gradients.rho[2] +
+ Wj[0] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+ dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
+ Wj[2] * pj->primitives.gradients.v[0][1] +
+ Wj[3] * pj->primitives.gradients.v[0][2] +
+ pj->primitives.gradients.P[0] / Wj[0]);
+ dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
+ Wj[2] * pj->primitives.gradients.v[1][1] +
+ Wj[3] * pj->primitives.gradients.v[1][2] +
+ pj->primitives.gradients.P[1] / Wj[0]);
+ dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
+ Wj[2] * pj->primitives.gradients.v[2][1] +
+ Wj[3] * pj->primitives.gradients.v[2][2] +
+ pj->primitives.gradients.P[2] / Wj[0]);
+ dWj[4] -=
+ 0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
+ Wj[2] * pj->primitives.gradients.P[1] +
+ Wj[3] * pj->primitives.gradients.P[2] +
+ hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+
+ Wi[0] += dWi[0];
+ Wi[1] += dWi[1];
+ Wi[2] += dWi[2];
+ Wi[3] += dWi[3];
+ Wi[4] += dWi[4];
+
+ Wj[0] += dWj[0];
+ Wj[1] += dWj[1];
+ Wj[2] += dWj[2];
+ Wj[3] += dWj[3];
+ Wj[4] += dWj[4];
+}
+
+#endif // SWIFT_HYDRO_GRADIENTS_H
diff --git a/src/hydro/Shadowswift/hydro_gradients_gizmo.h b/src/hydro/Shadowswift/hydro_gradients_gizmo.h
new file mode 100644
index 0000000000000000000000000000000000000000..aa6e4406b94e7a5cafcd0ca556162476003477de
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_gradients_gizmo.h
@@ -0,0 +1,341 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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 Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {
+
+ p->primitives.gradients.rho[0] = 0.0f;
+ p->primitives.gradients.rho[1] = 0.0f;
+ p->primitives.gradients.rho[2] = 0.0f;
+
+ p->primitives.gradients.v[0][0] = 0.0f;
+ p->primitives.gradients.v[0][1] = 0.0f;
+ p->primitives.gradients.v[0][2] = 0.0f;
+
+ p->primitives.gradients.v[1][0] = 0.0f;
+ p->primitives.gradients.v[1][1] = 0.0f;
+ p->primitives.gradients.v[1][2] = 0.0f;
+
+ p->primitives.gradients.v[2][0] = 0.0f;
+ p->primitives.gradients.v[2][1] = 0.0f;
+ p->primitives.gradients.v[2][2] = 0.0f;
+
+ p->primitives.gradients.P[0] = 0.0f;
+ p->primitives.gradients.P[1] = 0.0f;
+ p->primitives.gradients.P[2] = 0.0f;
+
+ hydro_slope_limit_cell_init(p);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float r = sqrtf(r2);
+ float xi, xj;
+ float hi_inv, hj_inv;
+ float wi, wj, wi_dx, wj_dx;
+ int k, l;
+ float Bi[3][3];
+ float Bj[3][3];
+ float Wi[5], Wj[5];
+
+ /* Initialize local variables */
+ for (k = 0; k < 3; k++) {
+ for (l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ Bj[k][l] = pj->geometry.matrix_E[k][l];
+ }
+ }
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ /* Compute kernel of pi. */
+ hi_inv = 1.0 / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+
+ /* Compute kernel of pj. */
+ hj_inv = 1.0 / hj;
+ xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* Compute gradients for pj */
+ /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
+ * want
+ * it to be */
+ pj->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ pj->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ pj->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ hydro_slope_limit_cell_collect(pj, pi, r);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
+
+ float r = sqrtf(r2);
+ float xi;
+ float hi_inv;
+ float wi, wi_dx;
+ int k, l;
+ float Bi[3][3];
+ float Wi[5], Wj[5];
+
+ /* Initialize local variables */
+ for (k = 0; k < 3; k++) {
+ for (l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ }
+ }
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ /* Compute kernel of pi. */
+ hi_inv = 1.0 / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {
+
+ float h, ih;
+
+ /* add kernel normalization to gradients */
+ h = p->h;
+ ih = 1.0f / h;
+ const float ihdim = pow_dimension(ih);
+
+ p->primitives.gradients.rho[0] *= ihdim;
+ p->primitives.gradients.rho[1] *= ihdim;
+ p->primitives.gradients.rho[2] *= ihdim;
+
+ p->primitives.gradients.v[0][0] *= ihdim;
+ p->primitives.gradients.v[0][1] *= ihdim;
+ p->primitives.gradients.v[0][2] *= ihdim;
+ p->primitives.gradients.v[1][0] *= ihdim;
+ p->primitives.gradients.v[1][1] *= ihdim;
+ p->primitives.gradients.v[1][2] *= ihdim;
+ p->primitives.gradients.v[2][0] *= ihdim;
+ p->primitives.gradients.v[2][1] *= ihdim;
+ p->primitives.gradients.v[2][2] *= ihdim;
+
+ p->primitives.gradients.P[0] *= ihdim;
+ p->primitives.gradients.P[1] *= ihdim;
+ p->primitives.gradients.P[2] *= ihdim;
+
+ hydro_slope_limit_cell(p);
+}
diff --git a/src/hydro/Shadowswift/hydro_gradients_sph.h b/src/hydro/Shadowswift/hydro_gradients_sph.h
new file mode 100644
index 0000000000000000000000000000000000000000..f635faecea549f7da280ade9b944021a5e4aeb4c
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_gradients_sph.h
@@ -0,0 +1,248 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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 Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {
+
+ p->primitives.gradients.rho[0] = 0.0f;
+ p->primitives.gradients.rho[1] = 0.0f;
+ p->primitives.gradients.rho[2] = 0.0f;
+
+ p->primitives.gradients.v[0][0] = 0.0f;
+ p->primitives.gradients.v[0][1] = 0.0f;
+ p->primitives.gradients.v[0][2] = 0.0f;
+
+ p->primitives.gradients.v[1][0] = 0.0f;
+ p->primitives.gradients.v[1][1] = 0.0f;
+ p->primitives.gradients.v[1][2] = 0.0f;
+ p->primitives.gradients.v[2][0] = 0.0f;
+ p->primitives.gradients.v[2][1] = 0.0f;
+ p->primitives.gradients.v[2][2] = 0.0f;
+
+ p->primitives.gradients.P[0] = 0.0f;
+ p->primitives.gradients.P[1] = 0.0f;
+ p->primitives.gradients.P[2] = 0.0f;
+
+ hydro_slope_limit_cell_init(p);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float wi, wi_dx, xi, hi_inv;
+ float wj, wj_dx, xj, hj_inv;
+ float r = sqrtf(r2);
+
+ hi_inv = 1.0f / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* very basic gradient estimate */
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+
+ hj_inv = 1.0f / hj;
+ xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* signs are the same as before, since we swap i and j twice */
+ pj->primitives.gradients.rho[0] -=
+ wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pj->primitives.gradients.rho[1] -=
+ wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pj->primitives.gradients.rho[2] -=
+ wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pj->primitives.gradients.v[0][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pj->primitives.gradients.v[0][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pj->primitives.gradients.v[0][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+ pj->primitives.gradients.v[1][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[1][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[1][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[2][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pj->primitives.gradients.v[2][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pj->primitives.gradients.v[2][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pj->primitives.gradients.P[0] -=
+ wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pj->primitives.gradients.P[1] -=
+ wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pj->primitives.gradients.P[2] -=
+ wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ hydro_slope_limit_cell_collect(pj, pi, r);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
+
+ float wi, wi_dx, xi, hi_inv;
+ float r = sqrtf(r2);
+
+ hi_inv = 1.0f / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* very basic gradient estimate */
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {
+
+ const float h = p->h;
+ const float ih = 1.0f / h;
+ const float ihdimp1 = pow_dimension_plus_one(ih);
+
+ float volume = p->geometry.volume;
+
+ /* finalize gradients by multiplying with volume */
+ p->primitives.gradients.rho[0] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[1] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.P[0] *= ihdimp1 * volume;
+ p->primitives.gradients.P[1] *= ihdimp1 * volume;
+ p->primitives.gradients.P[2] *= ihdimp1 * volume;
+
+ hydro_slope_limit_cell(p);
+}
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..b68afc77a889385dd8c909084ded66dfe06f38db
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -0,0 +1,486 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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 "adiabatic_index.h"
+#include "hydro_gradients.h"
+#include "riemann.h"
+
+/**
+ * @brief Calculate the volume interaction between particle i and particle j
+ *
+ * The volume is in essence the same as the weighted number of neighbours in a
+ * classical SPH density calculation.
+ *
+ * We also calculate the components of the matrix E, which is used for second
+ * order accurate gradient calculations and for the calculation of the interface
+ * surface areas.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float r = sqrtf(r2);
+ float xi, xj;
+ float h_inv;
+ float wi, wj, wi_dx, wj_dx;
+ int k, l;
+
+ /* Compute density of pi. */
+ h_inv = 1.0 / hi;
+ xi = r * h_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= xi * wi_dx;
+
+ /* Density. Needed for h_dt. */
+ pi->rho += pj->mass * wi;
+
+ /* these are eqns. (1) and (2) in the summary */
+ pi->geometry.volume += wi;
+ for (k = 0; k < 3; k++)
+ for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+ /* Compute density of pj. */
+ h_inv = 1.0 / hj;
+ xj = r * h_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= xj * wj_dx;
+
+ /* Density. Needed for h_dt. */
+ pj->rho += pi->mass * wi;
+
+ /* these are eqns. (1) and (2) in the summary */
+ pj->geometry.volume += wj;
+ for (k = 0; k < 3; k++)
+ for (l = 0; l < 3; l++) pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+}
+
+/**
+ * @brief Calculate the volume interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * The volume is in essence the same as the weighted number of neighbours in a
+ * classical SPH density calculation.
+ *
+ * We also calculate the components of the matrix E, which is used for second
+ * order accurate gradient calculations and for the calculation of the interface
+ * surface areas.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float r;
+ float xi;
+ float h_inv;
+ float wi, wi_dx;
+ int k, l;
+
+ /* Get r and r inverse. */
+ r = sqrtf(r2);
+
+ h_inv = 1.0 / hi;
+ xi = r * h_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= xi * wi_dx;
+
+ pi->rho += pj->mass * wi;
+
+ /* these are eqns. (1) and (2) in the summary */
+ pi->geometry.volume += wi;
+ for (k = 0; k < 3; k++)
+ for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
+}
+
+/**
+ * @brief Common part of the flux calculation between particle i and j
+ *
+ * Since the only difference between the symmetric and non-symmetric version
+ * of the flux calculation is in the update of the conserved variables at the
+ * very end (which is not done for particle j if mode is 0 and particle j is
+ * active), both runner_iact_force and runner_iact_nonsym_force call this
+ * method, with an appropriate mode.
+ *
+ * This method calculates the surface area of the interface between particle i
+ * and particle j, as well as the interface position and velocity. These are
+ * then used to reconstruct and predict the primitive variables, which are then
+ * fed to a Riemann solver that calculates a flux. This flux is used to update
+ * the conserved variables of particle i or both particles.
+ *
+ * This method also calculates the maximal velocity used to calculate the time
+ * step.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
+ int mode) {
+
+ float r = sqrtf(r2);
+ float xi, xj;
+ float hi_inv, hi_inv_dim;
+ float hj_inv, hj_inv_dim;
+ float wi, wj, wi_dx, wj_dx;
+ int k, l;
+ float A[3];
+ float Anorm;
+ float Bi[3][3];
+ float Bj[3][3];
+ float Vi, Vj;
+ float xij_i[3], xfac, xijdotdx;
+ float vmax, dvdotdx;
+ float vi[3], vj[3], vij[3];
+ float Wi[5], Wj[5];
+ float dti, dtj, mindt;
+ float n_unit[3];
+
+ /* Initialize local variables */
+ for (k = 0; k < 3; k++) {
+ for (l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ Bj[k][l] = pj->geometry.matrix_E[k][l];
+ }
+ vi[k] = pi->force.v_full[k]; /* particle velocities */
+ vj[k] = pj->force.v_full[k];
+ }
+ Vi = pi->geometry.volume;
+ Vj = pj->geometry.volume;
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ dti = pi->force.dt;
+ dtj = pj->force.dt;
+
+ /* calculate the maximal signal velocity */
+ if (Wi[0] && Wj[0]) {
+ vmax =
+ sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+ } else {
+ vmax = 0.0f;
+ }
+ dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+ (Wi[3] - Wj[3]) * dx[2];
+ if (dvdotdx > 0.) {
+ vmax -= dvdotdx / r;
+ }
+ pi->timestepvars.vmax = fmaxf(pi->timestepvars.vmax, vmax);
+ if (mode == 1) {
+ pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax);
+ }
+
+ /* The flux will be exchanged using the smallest time step of the two
+ * particles */
+ mindt = fminf(dti, dtj);
+ dti = mindt;
+ dtj = mindt;
+
+ /* Compute kernel of pi. */
+ hi_inv = 1.0 / hi;
+ hi_inv_dim = pow_dimension(hi_inv);
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Compute kernel of pj. */
+ hj_inv = 1.0 / hj;
+ hj_inv_dim = pow_dimension(hj_inv);
+ xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* Compute h_dt */
+ float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+ float ri = 1.0f / r;
+ float hidp1 = pow_dimension_plus_one(hi_inv);
+ float hjdp1 = pow_dimension_plus_one(hj_inv);
+ float wi_dr = hidp1 * wi_dx;
+ float wj_dr = hjdp1 * wj_dx;
+ dvdr *= ri;
+ pi->force.h_dt -= pj->mass * dvdr / pj->rho * wi_dr;
+ if (mode == 1) {
+ pj->force.h_dt -= pi->mass * dvdr / pi->rho * wj_dr;
+ }
+
+ /* Compute area */
+ /* eqn. (7) */
+ Anorm = 0.0f;
+ for (k = 0; k < 3; k++) {
+ /* we add a minus sign since dx is pi->x - pj->x */
+ A[k] = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
+ hi_inv_dim -
+ Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
+ hj_inv_dim;
+ Anorm += A[k] * A[k];
+ }
+
+ if (!Anorm) {
+ /* if the interface has no area, nothing happens and we return */
+ /* continuing results in dividing by zero and NaN's... */
+ return;
+ }
+
+ /* compute the normal vector of the interface */
+ Anorm = sqrtf(Anorm);
+ for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
+
+ /* Compute interface position (relative to pi, since we don't need the actual
+ * position) */
+ /* eqn. (8) */
+ xfac = hi / (hi + hj);
+ for (k = 0; k < 3; k++) xij_i[k] = -xfac * dx[k];
+
+ /* Compute interface velocity */
+ /* eqn. (9) */
+ xijdotdx = xij_i[0] * dx[0] + xij_i[1] * dx[1] + xij_i[2] * dx[2];
+ for (k = 0; k < 3; k++) vij[k] = vi[k] + (vi[k] - vj[k]) * xijdotdx / r2;
+
+ /* complete calculation of position of interface */
+ /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
+ correction terms for periodicity. If we do the interpolation,
+ we have to use xij w.r.t. the actual particle.
+ => we need a separate xij for pi and pj... */
+ /* tldr: we do not need the code below, but we do need the same code as above
+ but then
+ with i and j swapped */
+ // for ( k = 0 ; k < 3 ; k++ )
+ // xij[k] += pi->x[k];
+
+ /* Boost the primitive variables to the frame of reference of the interface */
+ /* Note that velocities are indices 1-3 in W */
+ Wi[1] -= vij[0];
+ Wi[2] -= vij[1];
+ Wi[3] -= vij[2];
+ Wj[1] -= vij[0];
+ Wj[2] -= vij[1];
+ Wj[3] -= vij[2];
+
+ hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt);
+
+ /* we don't need to rotate, we can use the unit vector in the Riemann problem
+ * itself (see GIZMO) */
+
+ if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
+ printf("mindt: %g\n", mindt);
+ printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
+ pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
+#ifdef USE_GRADIENTS
+ printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
+#endif
+ printf("gradWL[0]: %g %g %g\n", pi->primitives.gradients.rho[0],
+ pi->primitives.gradients.rho[1], pi->primitives.gradients.rho[2]);
+ printf("gradWL[1]: %g %g %g\n", pi->primitives.gradients.v[0][0],
+ pi->primitives.gradients.v[0][1], pi->primitives.gradients.v[0][2]);
+ printf("gradWL[2]: %g %g %g\n", pi->primitives.gradients.v[1][0],
+ pi->primitives.gradients.v[1][1], pi->primitives.gradients.v[1][2]);
+ printf("gradWL[3]: %g %g %g\n", pi->primitives.gradients.v[2][0],
+ pi->primitives.gradients.v[2][1], pi->primitives.gradients.v[2][2]);
+ printf("gradWL[4]: %g %g %g\n", pi->primitives.gradients.P[0],
+ pi->primitives.gradients.P[1], pi->primitives.gradients.P[2]);
+ printf("WL': %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
+ printf("WR: %g %g %g %g %g\n", pj->primitives.rho, pj->primitives.v[0],
+ pj->primitives.v[1], pj->primitives.v[2], pj->primitives.P);
+#ifdef USE_GRADIENTS
+ printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
+#endif
+ printf("gradWR[0]: %g %g %g\n", pj->primitives.gradients.rho[0],
+ pj->primitives.gradients.rho[1], pj->primitives.gradients.rho[2]);
+ printf("gradWR[1]: %g %g %g\n", pj->primitives.gradients.v[0][0],
+ pj->primitives.gradients.v[0][1], pj->primitives.gradients.v[0][2]);
+ printf("gradWR[2]: %g %g %g\n", pj->primitives.gradients.v[1][0],
+ pj->primitives.gradients.v[1][1], pj->primitives.gradients.v[1][2]);
+ printf("gradWR[3]: %g %g %g\n", pj->primitives.gradients.v[2][0],
+ pj->primitives.gradients.v[2][1], pj->primitives.gradients.v[2][2]);
+ printf("gradWR[4]: %g %g %g\n", pj->primitives.gradients.P[0],
+ pj->primitives.gradients.P[1], pj->primitives.gradients.P[2]);
+ printf("WR': %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
+ error("Negative density or pressure!\n");
+ }
+
+ float totflux[5];
+ riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
+
+ /* Update conserved variables */
+ /* eqn. (16) */
+ pi->conserved.flux.mass -= dti * Anorm * totflux[0];
+ pi->conserved.flux.momentum[0] -= dti * Anorm * totflux[1];
+ pi->conserved.flux.momentum[1] -= dti * Anorm * totflux[2];
+ pi->conserved.flux.momentum[2] -= dti * Anorm * totflux[3];
+ pi->conserved.flux.energy -= dti * Anorm * totflux[4];
+
+ float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
+ pi->primitives.v[1] * pi->primitives.v[1] +
+ pi->primitives.v[2] * pi->primitives.v[2]);
+ pi->conserved.flux.energy += dti * Anorm * totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += dti * Anorm * totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += dti * Anorm * totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= dti * Anorm * totflux[0] * ekin;
+
+ /* here is how it works:
+ Mode will only be 1 if both particles are ACTIVE and they are in the same
+ cell. In this case, this method IS the flux calculation for particle j, and
+ we HAVE TO UPDATE it.
+ Mode 0 can mean several things: it can mean that particle j is INACTIVE, in
+ which case we NEED TO UPDATE it, since otherwise the flux is lost from the
+ system and the conserved variable is not conserved.
+ It can also mean that particle j sits in another cell and is ACTIVE. In
+ this case, the flux exchange for particle j is done TWICE and we SHOULD NOT
+ UPDATE particle j.
+ ==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
+ */
+ if (mode == 1 || pj->ti_end > pi->ti_end) {
+ pj->conserved.flux.mass += dtj * Anorm * totflux[0];
+ pj->conserved.flux.momentum[0] += dtj * Anorm * totflux[1];
+ pj->conserved.flux.momentum[1] += dtj * Anorm * totflux[2];
+ pj->conserved.flux.momentum[2] += dtj * Anorm * totflux[3];
+ pj->conserved.flux.energy += dtj * Anorm * totflux[4];
+
+ ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
+ pj->primitives.v[1] * pj->primitives.v[1] +
+ pj->primitives.v[2] * pj->primitives.v[2]);
+ pj->conserved.flux.energy -= dtj * Anorm * totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += dtj * Anorm * totflux[0] * ekin;
+ }
+}
+
+/**
+ * @brief Flux calculation between particle i and particle j
+ *
+ * This method calls runner_iact_fluxes_common with mode 1.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
+}
+
+/**
+ * @brief Flux calculation between particle i and particle j: non-symmetric
+ * version
+ *
+ * This method calls runner_iact_fluxes_common with mode 0.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
+}
+
+//// EMPTY VECTORIZED VERSIONS (gradients methods are missing...)
+
+__attribute__((always_inline)) INLINE static void runner_iact_vec_density(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
+ struct part **pi, struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void runner_iact_vec_force(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {}
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..aaae075fe8b8614333a95032424cda746bc52313
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -0,0 +1,133 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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 "adiabatic_index.h"
+#include "hydro_gradients.h"
+#include "hydro_slope_limiters.h"
+#include "io_properties.h"
+#include "riemann.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+void hydro_read_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 8;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, parts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] =
+ io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, primitives.P);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_DENSITY, parts, primitives.rho);
+}
+
+float convert_u(struct engine* e, struct part* p) {
+ return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+}
+
+float convert_A(struct engine* e, struct part* p) {
+ return p->primitives.P / pow_gamma(p->primitives.rho);
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+void hydro_write_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 13;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
+ parts, x);
+ list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
+ primitives.v);
+ list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
+ conserved.mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ parts, primitives.P, convert_u);
+ list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_output_field("Acceleration", FLOAT, 3,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
+ primitives.rho);
+ list[8] = io_make_output_field("Volume", FLOAT, 1, UNIT_CONV_VOLUME, parts,
+ geometry.volume);
+ list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
+ parts, primitives.gradients.rho);
+ list[10] = io_make_output_field_convert_part(
+ "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
+ list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+ parts, primitives.P);
+ list[12] = io_make_output_field("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
+ parts, conserved.energy);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+void writeSPHflavour(hid_t h_grpsph) {
+ /* Gradient information */
+ writeAttribute_s(h_grpsph, "Gradient reconstruction model",
+ HYDRO_GRADIENT_IMPLEMENTATION);
+
+ /* Slope limiter information */
+ writeAttribute_s(h_grpsph, "Cell wide slope limiter model",
+ HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
+ writeAttribute_s(h_grpsph, "Piecewise slope limiter model",
+ HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
+
+ /* Riemann solver information */
+ writeAttribute_s(h_grpsph, "Riemann solver type",
+ RIEMANN_SOLVER_IMPLEMENTATION);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..ee0a1fc13f1db85fa596018cf2e0c3eb9ca737ad
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -0,0 +1,196 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Extra particle data not needed during the computation. */
+struct xpart {
+
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /* Velocity at the last full step. */
+ float v_full[3];
+
+ /* Old density. */
+ float omega;
+
+} __attribute__((aligned(xpart_align)));
+
+/* Data of a single particle. */
+struct part {
+
+ /* Particle position. */
+ double x[3];
+
+ /* Particle predicted velocity. */
+ float v[3];
+
+ /* Particle acceleration. */
+ float a_hydro[3];
+
+ /* Particle cutoff radius. */
+ float h;
+
+ /* Particle time of beginning of time-step. */
+ int ti_begin;
+
+ /* Particle time of end of time-step. */
+ int ti_end;
+
+ /* The primitive hydrodynamical variables. */
+ struct {
+
+ /* Fluid velocity. */
+ float v[3];
+
+ /* Density. */
+ float rho;
+
+ /* Pressure. */
+ float P;
+
+ /* Gradients of the primitive variables. */
+ struct {
+
+ /* Density gradients. */
+ float rho[3];
+
+ /* Fluid velocity gradients. */
+ float v[3][3];
+
+ /* Pressure gradients. */
+ float P[3];
+
+ } gradients;
+
+ /* Quantities needed by the slope limiter. */
+ struct {
+
+ /* Extreme values of the density among the neighbours. */
+ float rho[2];
+
+ /* Extreme values of the fluid velocity among the neighbours. */
+ float v[3][2];
+
+ /* Extreme values of the pressure among the neighbours. */
+ float P[2];
+
+ /* Maximal distance to all neighbouring faces. */
+ float maxr;
+
+ } limiter;
+
+ } primitives;
+
+ /* The conserved hydrodynamical variables. */
+ struct {
+
+ /* Fluid momentum. */
+ float momentum[3];
+
+ /* Fluid mass (this field already exists outside of this struct as well). */
+ float mass;
+
+ /* Fluid thermal energy (not per unit mass!). */
+ float energy;
+
+ /* Fluxes. */
+ struct {
+
+ /* Mass flux. */
+ float mass;
+
+ /* Momentum flux. */
+ float momentum[3];
+
+ /* Energy flux. */
+ float energy;
+
+ } flux;
+
+ } conserved;
+
+ /* Geometrical quantities used for hydro. */
+ struct {
+
+ /* Volume of the particle. */
+ float volume;
+
+ /* Geometrical shear matrix used to calculate second order accurate
+ gradients */
+ float matrix_E[3][3];
+
+ } geometry;
+
+ /* Variables used for timestep calculation (currently not used). */
+ struct {
+
+ /* Maximum fluid velocity among all neighbours. */
+ float vmax;
+
+ } timestepvars;
+
+ /* Quantities used during the volume (=density) loop. */
+ struct {
+
+ /* Particle velocity divergence. */
+ float div_v;
+
+ /* Derivative of particle number density. */
+ float wcount_dh;
+
+ /* Particle velocity curl. */
+ float curl_v[3];
+
+ /* Particle number density. */
+ float wcount;
+
+ } density;
+
+ /* Quantities used during the force loop. */
+ struct {
+
+ /* Needed to drift the primitive variables. */
+ float h_dt;
+
+ /* Physical time step of the particle. */
+ float dt;
+
+ /* Actual velocity of the particle. */
+ float v_full[3];
+
+ } force;
+
+ /* Particle mass (this field is also part of the conserved quantities...). */
+ float mass;
+
+ /* Particle ID. */
+ unsigned long long id;
+
+ /* Associated gravitas. */
+ struct gpart *gpart;
+
+ /* Variables needed for the code to compile (should be removed/replaced). */
+ float rho;
+
+ /* Old internal energy flux */
+ float du_dt;
+
+} __attribute__((aligned(part_align)));
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters.h b/src/hydro/Shadowswift/hydro_slope_limiters.h
new file mode 100644
index 0000000000000000000000000000000000000000..e1f841411a65fc5fce78ecab07e272c9c6df475b
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_slope_limiters.h
@@ -0,0 +1,97 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_HYDRO_SLOPE_LIMITERS_H
+#define SWIFT_HYDRO_SLOPE_LIMITERS_H
+
+#define PER_FACE_LIMITER
+#define CELL_WIDE_LIMITER
+
+#include "dimension.h"
+#include "kernel_hydro.h"
+
+#ifdef PER_FACE_LIMITER
+
+#define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION \
+ "GIZMO piecewise slope limiter (Hopkins 2015)"
+#include "hydro_slope_limiters_face.h"
+
+#else
+
+#define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION "No piecewise slope limiter"
+
+/**
+ * @brief Slope limit the slopes at the interface between two particles
+ *
+ * @param Wi Hydrodynamic variables of particle i.
+ * @param Wj Hydrodynamic variables of particle j.
+ * @param dWi Difference between the hydrodynamic variables of particle i at the
+ * position of particle i and at the interface position.
+ * @param dWj Difference between the hydrodynamic variables of particle j at the
+ * position of particle j and at the interface position.
+ * @param xij_i Relative position vector of the interface w.r.t. particle i.
+ * @param xij_j Relative position vector of the interface w.r.t. partilce j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+ float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+ float r) {}
+
+#endif
+
+#ifdef CELL_WIDE_LIMITER
+
+#define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION \
+ "Cell wide slope limiter (Springel 2010)"
+#include "hydro_slope_limiters_cell.h"
+
+#else
+
+#define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION "No cell wide slope limiter"
+
+/**
+ * @brief Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+ struct part *p) {}
+
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_slope_limit_cell_collect(struct part *pi, struct part *pj, float r) {}
+
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+ struct part *p) {}
+
+#endif
+
+#endif // SWIFT_HYDRO_SLOPE_LIMITERS_H
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
new file mode 100644
index 0000000000000000000000000000000000000000..aa99b43721f669f47a7888a5da0b1933ca1ebd62
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
@@ -0,0 +1,173 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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 Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+ struct part* p) {
+
+ p->primitives.limiter.rho[0] = FLT_MAX;
+ p->primitives.limiter.rho[1] = -FLT_MAX;
+ p->primitives.limiter.v[0][0] = FLT_MAX;
+ p->primitives.limiter.v[0][1] = -FLT_MAX;
+ p->primitives.limiter.v[1][0] = FLT_MAX;
+ p->primitives.limiter.v[1][1] = -FLT_MAX;
+ p->primitives.limiter.v[2][0] = FLT_MAX;
+ p->primitives.limiter.v[2][1] = -FLT_MAX;
+ p->primitives.limiter.P[0] = FLT_MAX;
+ p->primitives.limiter.P[1] = -FLT_MAX;
+
+ p->primitives.limiter.maxr = -FLT_MAX;
+}
+
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
+
+ /* basic slope limiter: collect the maximal and the minimal value for the
+ * primitive variables among the ngbs */
+ pi->primitives.limiter.rho[0] =
+ fmin(pj->primitives.rho, pi->primitives.limiter.rho[0]);
+ pi->primitives.limiter.rho[1] =
+ fmax(pj->primitives.rho, pi->primitives.limiter.rho[1]);
+
+ pi->primitives.limiter.v[0][0] =
+ fmin(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
+ pi->primitives.limiter.v[0][1] =
+ fmax(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
+ pi->primitives.limiter.v[1][0] =
+ fmin(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
+ pi->primitives.limiter.v[1][1] =
+ fmax(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
+ pi->primitives.limiter.v[2][0] =
+ fmin(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
+ pi->primitives.limiter.v[2][1] =
+ fmax(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
+
+ pi->primitives.limiter.P[0] =
+ fmin(pj->primitives.P, pi->primitives.limiter.P[0]);
+ pi->primitives.limiter.P[1] =
+ fmax(pj->primitives.P, pi->primitives.limiter.P[1]);
+
+ pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
+}
+
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+ struct part* p) {
+
+ float gradrho[3], gradv[3][3], gradP[3];
+ float gradtrue, gradmax, gradmin, alpha;
+
+ gradrho[0] = p->primitives.gradients.rho[0];
+ gradrho[1] = p->primitives.gradients.rho[1];
+ gradrho[2] = p->primitives.gradients.rho[2];
+
+ gradv[0][0] = p->primitives.gradients.v[0][0];
+ gradv[0][1] = p->primitives.gradients.v[0][1];
+ gradv[0][2] = p->primitives.gradients.v[0][2];
+
+ gradv[1][0] = p->primitives.gradients.v[1][0];
+ gradv[1][1] = p->primitives.gradients.v[1][1];
+ gradv[1][2] = p->primitives.gradients.v[1][2];
+
+ gradv[2][0] = p->primitives.gradients.v[2][0];
+ gradv[2][1] = p->primitives.gradients.v[2][1];
+ gradv[2][2] = p->primitives.gradients.v[2][2];
+
+ gradP[0] = p->primitives.gradients.P[0];
+ gradP[1] = p->primitives.gradients.P[1];
+ gradP[2] = p->primitives.gradients.P[2];
+
+ gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
+ gradrho[2] * gradrho[2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
+ gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
+ alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ p->primitives.gradients.rho[0] *= alpha;
+ p->primitives.gradients.rho[1] *= alpha;
+ p->primitives.gradients.rho[2] *= alpha;
+ }
+
+ gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
+ gradv[0][2] * gradv[0][2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
+ gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
+ alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ p->primitives.gradients.v[0][0] *= alpha;
+ p->primitives.gradients.v[0][1] *= alpha;
+ p->primitives.gradients.v[0][2] *= alpha;
+ }
+
+ gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
+ gradv[1][2] * gradv[1][2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
+ gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
+ alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ p->primitives.gradients.v[1][0] *= alpha;
+ p->primitives.gradients.v[1][1] *= alpha;
+ p->primitives.gradients.v[1][2] *= alpha;
+ }
+
+ gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
+ gradv[2][2] * gradv[2][2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
+ gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
+ alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ p->primitives.gradients.v[2][0] *= alpha;
+ p->primitives.gradients.v[2][1] *= alpha;
+ p->primitives.gradients.v[2][2] *= alpha;
+ }
+
+ gradtrue =
+ sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ gradmax = p->primitives.limiter.P[1] - p->primitives.P;
+ gradmin = p->primitives.P - p->primitives.limiter.P[0];
+ alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ p->primitives.gradients.P[0] *= alpha;
+ p->primitives.gradients.P[1] *= alpha;
+ p->primitives.gradients.P[2] *= alpha;
+ }
+}
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters_face.h b/src/hydro/Shadowswift/hydro_slope_limiters_face.h
new file mode 100644
index 0000000000000000000000000000000000000000..7ae5dd2eb073d9aae8ab6f2efffdf8df15b4bb4a
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_slope_limiters_face.h
@@ -0,0 +1,121 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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 Slope limit a single quantity at the interface
+ *
+ * @param phi_i Value of the quantity at the particle position.
+ * @param phi_j Value of the quantity at the neighbouring particle position.
+ * @param phi_mid0 Extrapolated value of the quantity at the interface position.
+ * @param xij_norm Distance between the particle position and the interface
+ * position.
+ * @param r Distance between the particle and its neighbour.
+ * @return The slope limited difference between the quantity at the particle
+ * position and the quantity at the interface position.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
+ float xij_norm, float r) {
+
+ float delta1, delta2, phimin, phimax, phibar, phiplus, phiminus, phi_mid;
+ const float psi1 = 0.5f;
+ const float psi2 = 0.25f;
+
+ if (phi_i == phi_j) {
+ return 0.0f;
+ }
+
+ delta1 = psi1 * fabs(phi_i - phi_j);
+ delta2 = psi2 * fabs(phi_i - phi_j);
+
+ phimin = fmin(phi_i, phi_j);
+ phimax = fmax(phi_i, phi_j);
+
+ phibar = phi_i + xij_norm / r * (phi_j - phi_i);
+
+ /* if sign(phimax+delta1) == sign(phimax) */
+ if ((phimax + delta1) * phimax > 0.0f) {
+ phiplus = phimax + delta1;
+ } else {
+ phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+ }
+
+ /* if sign(phimin-delta1) == sign(phimin) */
+ if ((phimin - delta1) * phimin > 0.0f) {
+ phiminus = phimin - delta1;
+ } else {
+ phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+ }
+
+ if (phi_i < phi_j) {
+ phi_mid = fmax(phiminus, fmin(phibar + delta2, phi_mid0));
+ } else {
+ phi_mid = fmin(phiplus, fmax(phibar - delta2, phi_mid0));
+ }
+
+ return phi_mid - phi_i;
+}
+
+/**
+ * @brief Slope limit the slopes at the interface between two particles
+ *
+ * @param Wi Hydrodynamic variables of particle i.
+ * @param Wj Hydrodynamic variables of particle j.
+ * @param dWi Difference between the hydrodynamic variables of particle i at the
+ * position of particle i and at the interface position.
+ * @param dWj Difference between the hydrodynamic variables of particle j at the
+ * position of particle j and at the interface position.
+ * @param xij_i Relative position vector of the interface w.r.t. particle i.
+ * @param xij_j Relative position vector of the interface w.r.t. partilce j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+ float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+ float r) {
+
+ float xij_i_norm, xij_j_norm;
+
+ xij_i_norm =
+ sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] + xij_i[2] * xij_i[2]);
+
+ xij_j_norm =
+ sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] + xij_j[2] * xij_j[2]);
+
+ dWi[0] = hydro_slope_limit_face_quantity(Wi[0], Wj[0], Wi[0] + dWi[0],
+ xij_i_norm, r);
+ dWi[1] = hydro_slope_limit_face_quantity(Wi[1], Wj[1], Wi[1] + dWi[1],
+ xij_i_norm, r);
+ dWi[2] = hydro_slope_limit_face_quantity(Wi[2], Wj[2], Wi[2] + dWi[2],
+ xij_i_norm, r);
+ dWi[3] = hydro_slope_limit_face_quantity(Wi[3], Wj[3], Wi[3] + dWi[3],
+ xij_i_norm, r);
+ dWi[4] = hydro_slope_limit_face_quantity(Wi[4], Wj[4], Wi[4] + dWi[4],
+ xij_i_norm, r);
+
+ dWj[0] = hydro_slope_limit_face_quantity(Wj[0], Wi[0], Wj[0] + dWj[0],
+ xij_j_norm, r);
+ dWj[1] = hydro_slope_limit_face_quantity(Wj[1], Wi[1], Wj[1] + dWj[1],
+ xij_j_norm, r);
+ dWj[2] = hydro_slope_limit_face_quantity(Wj[2], Wi[2], Wj[2] + dWj[2],
+ xij_j_norm, r);
+ dWj[3] = hydro_slope_limit_face_quantity(Wj[3], Wi[3], Wj[3] + dWj[3],
+ xij_j_norm, r);
+ dWj[4] = hydro_slope_limit_face_quantity(Wj[4], Wi[4], Wj[4] + dWj[4],
+ xij_j_norm, r);
+}
diff --git a/src/hydro/Shadowswift/voronoi1d_algorithm.h b/src/hydro/Shadowswift/voronoi1d_algorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..92f07918ffabefdfb1e45ee0024e275cd7b3e0f2
--- /dev/null
+++ b/src/hydro/Shadowswift/voronoi1d_algorithm.h
@@ -0,0 +1,78 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * 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_VORONOI1D_ALGORITHM_H
+#define SWIFT_VORONOI1D_ALGORITHM_H
+
+#include <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include "error.h"
+#include "inline.h"
+#include "voronoi1d_cell.h"
+
+__attribute__((always_inline)) INLINE void voronoi_cell_init(
+ struct voronoi_cell *cell, double *x, float max_radius) {
+ cell->x = x[0];
+ cell->xL = -DBL_MAX;
+ cell->xR = DBL_MAX;
+ cell->idL = 0;
+ cell->idR = 0;
+ cell->volume = 0.0f;
+ cell->centroid = 0.0f;
+ cell->max_radius = max_radius;
+}
+
+__attribute__((always_inline)) INLINE void voronoi_cell_interact(
+ struct voronoi_cell *cell, double *x, unsigned long long id) {
+
+ /* Check for stupidity */
+ if (x[0] == cell->x) {
+ error("Cannot interact a Voronoi cell generator with itself!");
+ }
+
+ if (x[0] < cell->x) {
+ /* New left neighbour? */
+ if (x[0] > cell->xL) {
+ cell->xL = x[0];
+ cell->idL = id;
+ }
+ } else {
+ /* New right neighbour? */
+ if (x[0] < cell->xR) {
+ cell->xR = x[0];
+ cell->idR = id;
+ }
+ }
+}
+
+__attribute__((always_inline)) INLINE void voronoi_cell_finalize(
+ struct voronoi_cell *cell) {
+
+ double x, xL, xR;
+ x = cell->x;
+ cell->max_radius = fmax(cell->xL, cell->xR);
+ cell->xL = xL = 0.5f * (cell->xL + x);
+ cell->xR = xR = 0.5f * (x + cell->xR);
+
+ cell->volume = xR - xL;
+ cell->centroid = 0.5f * (xL + xR);
+}
+
+#endif // SWIFT_VORONOI1D_ALGORITHM_H
diff --git a/src/hydro/Shadowswift/voronoi1d_cell.h b/src/hydro/Shadowswift/voronoi1d_cell.h
new file mode 100644
index 0000000000000000000000000000000000000000..0b1d62e8e6209a8e00dff39d99bf8f04a45b1966
--- /dev/null
+++ b/src/hydro/Shadowswift/voronoi1d_cell.h
@@ -0,0 +1,53 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * 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_VORONOI1D_CELL_H
+#define SWIFT_VORONOI1D_CELL_H
+
+#include "voronoi1d_cell.h"
+
+struct voronoi_cell {
+
+ /* The position of the generator of the cell. */
+ double x;
+
+ /* The position of the left neighbour of the cell. */
+ double xL;
+
+ /* The position of the right neighbour of the cell. */
+ double xR;
+
+ /* The particle ID of the left neighbour. */
+ unsigned long long idL;
+
+ /* The particle ID of the right neighbour. */
+ unsigned long long idR;
+
+ /* The "volume" of the 1D cell. */
+ float volume;
+
+ /* The centroid of the cell. */
+ float centroid;
+
+ /* The maximal search radius for the cell. If we interacted the cell with all
+ generators within this radius, the cell is complete. */
+ float max_radius;
+};
+
+#endif // SWIFT_VORONOI1D_CELL_H
diff --git a/src/hydro_io.h b/src/hydro_io.h
index f0a619b90b774574c434007b1c01a0e55e75e464..0cc167a481db75653b07890fa99c3afaa15a37a1 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -30,6 +30,8 @@
#include "./hydro/Default/hydro_io.h"
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_io.h"
+#elif defined(SHADOWSWIFT)
+#include "./hydro/Shadowswift/hydro_io.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/part.h b/src/part.h
index ea895e6e0295d6a8b63309c7bd6855daa2cf7d64..0c969aad6b5bc8f703dfa0198e6eab17ca624d73 100644
--- a/src/part.h
+++ b/src/part.h
@@ -48,6 +48,9 @@
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_part.h"
#define EXTRA_HYDRO_LOOP
+#elif defined(SHADOWSWIFT)
+#include "./hydro/Shadowswift/hydro_part.h"
+#define EXTRA_HYDRO_LOOP
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/tests/Makefile.am b/tests/Makefile.am
index ac99573af98da665be254defd11beac91920e865..f930ac2e13d5e32eed24afcd214a1592d0932732 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -25,14 +25,14 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh \
testParser.sh testSPHStep test125cells.sh testKernelGrav testFFT \
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
- testMatrixInversion
+ testMatrixInversion testVoronoi1D
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSPHStep testPair test27cells test125cells testParser \
testKernel testKernelGrav testFFT testInteractions testMaths testSymmetry \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
- testRiemannHLLC testMatrixInversion
+ testRiemannHLLC testMatrixInversion testVoronoi1D
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
@@ -75,6 +75,8 @@ testRiemannHLLC_SOURCES = testRiemannHLLC.c
testMatrixInversion_SOURCES = testMatrixInversion.c
+testVoronoi1D_SOURCES = testVoronoi1D.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 4619d1fce34e67d4a1f62af59792390310dcfccb..3d09c1c67d8c8f2240c04c6765e2787ca67b9701 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -99,7 +99,7 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(GIZMO_SPH)
+#elif defined(GIZMO_SPH) || defined(SHADOWSWIFT)
part->primitives.P = pressure;
#else
error("Need to define pressure here !");
@@ -196,7 +196,7 @@ void reset_particles(struct cell *c, enum velocity_field vel,
set_velocity(p, vel, size);
set_energy_state(p, press, size, density);
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
p->geometry.volume = p->mass / density;
p->primitives.rho = density;
p->primitives.v[0] = p->v[0];
@@ -269,7 +269,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
hydro_first_init_part(part, xpart);
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
part->geometry.volume = part->mass / density;
part->primitives.rho = density;
part->primitives.v[0] = part->v[0];
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 3d1d56d9245f9a9aee9e705754c46bb22ba0ef0a..dd7cc19962ff84c42473bfe4a34543e9f1c27b21 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -185,7 +185,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2], main_cell->parts[pid].rho,
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
0.f,
#else
main_cell->parts[pid].rho_dh,
@@ -227,7 +227,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
cj->parts[pjd].v[2], cj->parts[pjd].rho,
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
0.f,
#else
main_cell->parts[pjd].rho_dh,
diff --git a/tests/testPair.c b/tests/testPair.c
index 8afc2250434de7fdfeaa19a36a213dd9ea79d861..00d454f1cc80eeb858e7785d6a9664f5169275f7 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -132,7 +132,7 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
ci->parts[pid].id, ci->parts[pid].x[0], ci->parts[pid].x[1],
ci->parts[pid].x[2], ci->parts[pid].v[0], ci->parts[pid].v[1],
ci->parts[pid].v[2], ci->parts[pid].rho,
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
0.f,
#else
cj->parts[pid].rho_dh,
@@ -156,7 +156,7 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
cj->parts[pjd].v[2], cj->parts[pjd].rho,
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
0.f,
#else
cj->parts[pjd].rho_dh,
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index 6469d314fb8b1438cc2c9737669c1a13a97bd803..5db424991823a5752800d1b6d05f52592b5ae456 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -46,7 +46,7 @@ int main(int argc, char *argv[]) {
pi.id = 1;
pj.id = 2;
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
/* Give the primitive variables sensible values, since the Riemann solver does
not like negative densities and pressures */
pi.primitives.rho = random_uniform(0.1f, 1.0f);
@@ -150,7 +150,7 @@ int main(int argc, char *argv[]) {
runner_iact_nonsym_force(r2, dx, pj2.h, pi2.h, &pj2, &pi2);
/* Check that the particles are the same */
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
i_ok = 0;
j_ok = 0;
for (size_t i = 0; i < sizeof(struct part) / sizeof(float); ++i) {
diff --git a/tests/testVoronoi1D.c b/tests/testVoronoi1D.c
new file mode 100644
index 0000000000000000000000000000000000000000..4d556b93fb246d772b93bfa60893960478ab6221
--- /dev/null
+++ b/tests/testVoronoi1D.c
@@ -0,0 +1,52 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * 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 "hydro/Shadowswift/voronoi1d_algorithm.h"
+
+int main() {
+
+ /* Create a Voronoi cell */
+ double x[1] = {0.5f};
+ struct voronoi_cell cell;
+ voronoi_cell_init(&cell, x, 1.0f);
+
+ /* Interact with a left and right neighbour */
+ double xL[1] = {0.0f};
+ double xR[1] = {1.0f};
+ voronoi_cell_interact(&cell, xL, 1);
+ voronoi_cell_interact(&cell, xR, 2);
+
+ /* Finalize cell and check results */
+ voronoi_cell_finalize(&cell);
+
+ if (cell.volume != 0.5f) {
+ error("Wrong volume: %g!", cell.volume);
+ }
+ if (cell.centroid != 0.5f) {
+ error("Wrong centroid: %g!", cell.centroid);
+ }
+ if (cell.idL != 1) {
+ error("Wrong left neighbour: %llu!", cell.idL);
+ }
+ if (cell.idR != 2) {
+ error("Wrong right neighbour: %llu!", cell.idR);
+ }
+
+ return 0;
+}