Wrote first working version of 1D moving mesh code. Got stuck in mesh steering...

Wrote first working version of 1D moving mesh code. Got stuck in mesh steering and correct slope limiting.

examples/main.c

@@ -330,6 +330,13 @@ int main(int argc, char *argv[]) {
   struct hydro_props hydro_properties;
   hydro_props_init(&hydro_properties, params);
 
+#if defined(SHADOWSWIFT)
+  /* Override the variables governing the density iteration
+     (see src/hydro/Shadowswift/hydro.h for full explanation) */
+  hydro_properties.target_neighbours = 1.0f;
+  hydro_properties.delta_neighbours = 0.0f;
+#endif
+
   /* Initialise the external potential properties */
   struct external_potential potential;
   if (with_external_gravity) potential_init(params, &us, &potential);

src/dimension.h

@@ -40,19 +40,19 @@
 
 #define hydro_dimension 3.f
 #define hydro_dimension_inv 0.3333333333f
-#define hydro_dimention_unit_sphere ((float)(4. * M_PI / 3.))
+#define hydro_dimension_unit_sphere ((float)(4. * M_PI / 3.))
 
 #elif defined(HYDRO_DIMENSION_2D)
 
 #define hydro_dimension 2.f
 #define hydro_dimension_inv 0.5f
-#define hydro_dimention_unit_sphere ((float)M_PI)
+#define hydro_dimension_unit_sphere ((float)M_PI)
 
 #elif defined(HYDRO_DIMENSION_1D)
 
 #define hydro_dimension 1.f
 #define hydro_dimension_inv 1.f
-#define hydro_dimention_unit_sphere 2.f
+#define hydro_dimension_unit_sphere 2.f
 
 #else
 
@@ -202,6 +202,34 @@ invert_dimension_by_dimension_matrix(float A[3][3]) {
 #endif
 }
 
+/**
+ * @brief Get the radius of a dimension sphere with the given volume
+ *
+ * @param volume Volume of the dimension sphere
+ * @return Radius of the dimension sphere
+ */
+__attribute__((always_inline)) INLINE static float get_radius_dimension_sphere(
+    float volume) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+  return cbrtf(volume / hydro_dimension_unit_sphere);
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+  return sqrtf(volume / hydro_dimension_unit_sphere);
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+  return volume / hydro_dimension_unit_sphere;
+
+#else
+
+  error("The dimension is not defined !");
+
+#endif
+}
+
 /* ------------------------------------------------------------------------- */
 #ifdef WITH_VECTORIZATION
 

src/engine.c

@@ -2724,7 +2724,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;

src/hydro.h

@@ -46,7 +46,7 @@
 #include "./hydro/Shadowswift/hydro.h"
 #include "./hydro/Shadowswift/hydro_iact.h"
 #define SPH_IMPLEMENTATION \
-  "Shadowfax moving mesh (Vandenbroucke & De Rijcke 2016)"
+  "Shadowfax moving mesh (Vandenbroucke and De Rijcke 2016)"
 #else
 #error "Invalid choice of SPH variant"
 #endif

src/hydro/Shadowswift/hydro.h

@@ -80,15 +80,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   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;
+
   voronoi_cell_init(&p->cell, p->x);
 }
 
@@ -115,42 +107,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
 __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);
 
   float hnew = voronoi_cell_finalize(&p->cell);
@@ -163,15 +122,15 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
      This way, we can accept an iteration by setting p->density.wcount to 1.
      To get the right correction for h, we set wcount to something else
      (say 0), and then set p->density.wcount_dh to 1/(hnew-h). */
-  if(hnew < p->h){
+  if (hnew < p->h) {
     /* Iteration succesful: we accept, but manually set h to a smaller value
        for the next time step */
     p->density.wcount = 1.0f;
-    p->h = 1.1f*hnew;
+    p->h = 1.1f * hnew;
   } else {
     /* Iteration not succesful: we force h to become 1.1*hnew */
     p->density.wcount = 0.0f;
-    p->density.wcount_dh = 1.0f/(1.1f*hnew-p->h);
+    p->density.wcount_dh = 1.0f / (1.1f * hnew - p->h);
   }
   volume = p->cell.volume;
   p->geometry.volume = volume;
@@ -306,25 +265,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
  * @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;
-}
+    struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {}
 
 /**
  * @brief Set the particle acceleration after the flux loop
@@ -341,8 +282,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __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;
+  float vcell[3];
 
   /* Update the conserved variables. We do this here and not in the kick,
      since we need the updated variables below. */
@@ -367,23 +307,49 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
            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;
+    /* We want the cell velocity to be as close as possible to the fluid
+       velocity */
+    vcell[0] = p->conserved.momentum[0] / p->conserved.mass;
+    vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
+    vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
+
+    /* To prevent stupid things like cell crossovers or generators that move
+       outside their cell, we steer the motion of the cell somewhat */
+    if (p->primitives.rho) {
+      float centroid[3], d[3];
+      float volume, csnd, R, vfac, fac, dnrm;
+      voronoi_get_centroid(&p->cell, centroid);
+      d[0] = centroid[0] - p->x[0];
+      d[1] = centroid[1] - p->x[1];
+      d[2] = centroid[2] - p->x[2];
+      dnrm = sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
+      csnd = sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+      volume = p->cell.volume;
+      R = get_radius_dimension_sphere(volume);
+      fac = 4.0f * dnrm / R;
+      if (fac > 0.9f) {
+        if (fac < 1.1f) {
+          vfac = csnd * (dnrm - 0.225f * R) / dnrm / (0.05f * R);
+        } else {
+          vfac = csnd / dnrm;
+        }
+      } else {
+        vfac = 0.0f;
+      }
+      vcell[0] += vfac * d[0];
+      vcell[1] += vfac * d[1];
+      vcell[2] += vfac * d[2];
+    }
+
+    /* We know the desired velocity; now make sure the particle is accelerated
+       accordingly during the next time step */
+    p->a_hydro[0] = (vcell[0] - p->primitives.v[0]) / p->force.dt;
+    p->a_hydro[1] = (vcell[1] - p->primitives.v[1]) / p->force.dt;
+    p->a_hydro[2] = (vcell[2] - p->primitives.v[2]) / 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;
   }
 }
 

src/hydro/Shadowswift/hydro_gradients.h

@@ -20,20 +20,14 @@
 #ifndef SWIFT_HYDRO_GRADIENTS_H
 #define SWIFT_HYDRO_GRADIENTS_H
 
-//#define SPH_GRADIENTS
-#define GIZMO_GRADIENTS
+#define SHADOWFAX_GRADIENTS
 
 #include "hydro_slope_limiters.h"
 
-#if defined(SPH_GRADIENTS)
+#if defined(SHADOWFAX_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"
+#define HYDRO_GRADIENT_IMPLEMENTATION "Shadowfax gradients (Springel 2010)"
+#include "hydro_gradients_shadowfax.h"
 
 #else
 
@@ -97,16 +91,14 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
 
   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];
+  float A[3];
+
+  if (!voronoi_get_face(&pj->cell, pi->id, A, xij_j)) {
+    /* Should never happen */
+    message("pi->ti_end: %u, pj->ti_end: %u", pi->ti_end, pj->ti_end);
+    error("pj is neighbour of pi, but pi not of pj. This should never happen!");
+    return;
+  }
 
   dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
            pi->primitives.gradients.rho[1] * xij_i[1] +

src/hydro/Shadowswift/hydro_gradients_gizmo.h

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