diff --git a/src/const.h b/src/const.h
index daa573db25c0b554d5b0f62598682111dc9aa640..9693aebb4f15d969ed80d2f501cc2a7b6b21b0d3 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
diff --git a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
index 9977c36d94d9647da071b867de5a8d101457de6d..9ca40a604da3dc12bbb48ac033cd078f0561d8ab 100644
--- a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
+++ b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
@@ -86,15 +86,14 @@ __attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- float A[3], midpoint[3];
+ float A, midpoint[3];
- if (!voronoi_get_face(&pi->cell, pj->id, A, midpoint)) {
+ A = voronoi_get_face(&pi->cell, pj->id, midpoint);
+ if (!A) {
/* particle is not a cell neighbour: do nothing */
return;
}
- float Anorm = sqrtf(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
-
float c[3];
/* midpoint is relative w.r.t. pi->x, as is dx */
/* c is supposed to be the vector pointing from the midpoint of pi and pj to
@@ -108,15 +107,15 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r = sqrtf(r2);
hydro_gradients_single_quantity(pi->primitives.rho, pj->primitives.rho, c, dx,
- r, Anorm, pi->primitives.gradients.rho);
+ r, A, pi->primitives.gradients.rho);
hydro_gradients_single_quantity(pi->primitives.v[0], pj->primitives.v[0], c,
- dx, r, Anorm, pi->primitives.gradients.v[0]);
+ dx, r, A, pi->primitives.gradients.v[0]);
hydro_gradients_single_quantity(pi->primitives.v[1], pj->primitives.v[1], c,
- dx, r, Anorm, pi->primitives.gradients.v[1]);
+ dx, r, A, pi->primitives.gradients.v[1]);
hydro_gradients_single_quantity(pi->primitives.v[2], pj->primitives.v[2], c,
- dx, r, Anorm, pi->primitives.gradients.v[2]);
+ dx, r, A, pi->primitives.gradients.v[2]);
hydro_gradients_single_quantity(pi->primitives.P, pj->primitives.P, c, dx, r,
- Anorm, pi->primitives.gradients.P);
+ A, pi->primitives.gradients.P);
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -125,19 +124,15 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
mindx[1] = -dx[1];
mindx[2] = -dx[2];
hydro_gradients_single_quantity(pj->primitives.rho, pi->primitives.rho, c,
- mindx, r, Anorm,
- pj->primitives.gradients.rho);
+ mindx, r, A, pj->primitives.gradients.rho);
hydro_gradients_single_quantity(pj->primitives.v[0], pi->primitives.v[0], c,
- mindx, r, Anorm,
- pj->primitives.gradients.v[0]);
+ mindx, r, A, pj->primitives.gradients.v[0]);
hydro_gradients_single_quantity(pj->primitives.v[1], pi->primitives.v[1], c,
- mindx, r, Anorm,
- pj->primitives.gradients.v[1]);
+ mindx, r, A, pj->primitives.gradients.v[1]);
hydro_gradients_single_quantity(pj->primitives.v[2], pi->primitives.v[2], c,
- mindx, r, Anorm,
- pj->primitives.gradients.v[2]);
+ mindx, r, A, pj->primitives.gradients.v[2]);
hydro_gradients_single_quantity(pj->primitives.P, pi->primitives.P, c, mindx,
- r, Anorm, pj->primitives.gradients.P);
+ r, A, pj->primitives.gradients.P);
hydro_slope_limit_cell_collect(pj, pi, r);
}
@@ -156,15 +151,14 @@ __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 A[3], midpoint[3];
+ float A, midpoint[3];
- if (!voronoi_get_face(&pi->cell, pj->id, A, midpoint)) {
+ A = voronoi_get_face(&pi->cell, pj->id, midpoint);
+ if (!A) {
/* particle is not a cell neighbour: do nothing */
return;
}
- float Anorm = sqrtf(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
-
float c[3];
/* midpoint is relative w.r.t. pi->x, as is dx */
/* c is supposed to be the vector pointing from the midpoint of pi and pj to
@@ -178,15 +172,15 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
float r = sqrtf(r2);
hydro_gradients_single_quantity(pi->primitives.rho, pj->primitives.rho, c, dx,
- r, Anorm, pi->primitives.gradients.rho);
+ r, A, pi->primitives.gradients.rho);
hydro_gradients_single_quantity(pi->primitives.v[0], pj->primitives.v[0], c,
- dx, r, Anorm, pi->primitives.gradients.v[0]);
+ dx, r, A, pi->primitives.gradients.v[0]);
hydro_gradients_single_quantity(pi->primitives.v[1], pj->primitives.v[1], c,
- dx, r, Anorm, pi->primitives.gradients.v[1]);
+ dx, r, A, pi->primitives.gradients.v[1]);
hydro_gradients_single_quantity(pi->primitives.v[2], pj->primitives.v[2], c,
- dx, r, Anorm, pi->primitives.gradients.v[2]);
+ dx, r, A, pi->primitives.gradients.v[2]);
hydro_gradients_single_quantity(pi->primitives.P, pj->primitives.P, c, dx, r,
- Anorm, pi->primitives.gradients.P);
+ A, pi->primitives.gradients.P);
hydro_slope_limit_cell_collect(pi, pj, r);
}
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 1c1ec9e1ad1f833d3bf73986026d7d0c96d8ff5f..8997155a0cb30848a74b220811970421b6d3396a 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -134,8 +134,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r = sqrtf(r2);
int k;
- float A[3];
- float Anorm;
+ float A;
float xij_i[3];
float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
@@ -143,7 +142,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float dti, dtj, mindt;
float n_unit[3];
- if (!voronoi_get_face(&pi->cell, pj->id, A, xij_i)) {
+ A = voronoi_get_face(&pi->cell, pj->id, xij_i);
+ if (!A) {
/* this neighbour does not share a face with the cell, return */
return;
}
@@ -190,11 +190,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dti = mindt;
dtj = mindt;
- Anorm = A[0] * A[0] + A[1] * A[1] + A[2] * A[2];
-
/* compute the normal vector of the interface */
- Anorm = sqrtf(Anorm);
- for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
+ for (k = 0; k < 3; ++k) {
+ n_unit[k] = -dx[k] / r;
+ }
/* Compute interface velocity */
float fac = (vi[0] - vj[0]) * (xij_i[0] + 0.5f * dx[0]) +
@@ -261,19 +260,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* 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];
+ pi->conserved.flux.mass -= dti * A * totflux[0];
+ pi->conserved.flux.momentum[0] -= dti * A * totflux[1];
+ pi->conserved.flux.momentum[1] -= dti * A * totflux[2];
+ pi->conserved.flux.momentum[2] -= dti * A * totflux[3];
+ pi->conserved.flux.energy -= dti * A * 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;
+ pi->conserved.flux.energy += dti * A * totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += dti * A * totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += dti * A * totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= dti * A * totflux[0] * ekin;
/* here is how it works:
Mode will only be 1 if both particles are ACTIVE and they are in the same
@@ -288,19 +287,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> 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];
+ pj->conserved.flux.mass += dtj * A * totflux[0];
+ pj->conserved.flux.momentum[0] += dtj * A * totflux[1];
+ pj->conserved.flux.momentum[1] += dtj * A * totflux[2];
+ pj->conserved.flux.momentum[2] += dtj * A * totflux[3];
+ pj->conserved.flux.energy += dtj * A * 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;
+ pj->conserved.flux.energy -= dtj * A * totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= dtj * A * totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= dtj * A * totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += dtj * A * totflux[0] * ekin;
}
}
diff --git a/src/hydro/Shadowswift/voronoi1d_algorithm.h b/src/hydro/Shadowswift/voronoi1d_algorithm.h
index 222c97b617512b4fcaa1bfd935674bd5b7d8eca5..e8a6ccd61cdf90f4887028fa042fde93ecd90d5e 100644
--- a/src/hydro/Shadowswift/voronoi1d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi1d_algorithm.h
@@ -128,43 +128,37 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
* deal with it appropriately.
*
* For this specific case, we simply check if the neighbour is the left or
- * right neighbour and set the surface area to the unit vector pointing either
- * to the right or to the left. The midpoint is set to the relative position
- * vector of the appropriate face.
+ * right neighbour and set the surface area to 1. The midpoint is set to the
+ * relative position vector of the appropriate face.
*
* @param cell 1D Voronoi cell.
* @param ngb ID of a particle that is possibly a neighbour of this cell.
- * @param A Array to store the oriented surface area in.
* @param midpoint Array to store the relative position of the face in.
- * @return 0 if the given neighbour is not a neighbour, 1 otherwise.
+ * @return 0 if the given neighbour is not a neighbour, surface area 1.0f
+ * otherwise.
*/
-__attribute__((always_inline)) INLINE int voronoi_get_face(
- struct voronoi_cell *cell, unsigned long long ngb, float *A,
- float *midpoint) {
+__attribute__((always_inline)) INLINE float voronoi_get_face(
+ struct voronoi_cell *cell, unsigned long long ngb, float *midpoint) {
if (ngb != cell->idL && ngb != cell->idR) {
/* this is perfectly possible: we interact with all particles within the
smoothing length, and they do not need to be all neighbours.
If this happens, we return 0, so that the flux method can return */
- return 0;
+ return 0.0f;
}
if (ngb == cell->idL) {
/* Left face */
- A[0] = -1.0f;
midpoint[0] = cell->xL;
} else {
/* Right face */
- A[0] = 1.0f;
midpoint[0] = cell->xR;
}
- /* The other components of A and midpoint are just zero */
- A[1] = 0.0f;
- A[2] = 0.0f;
+ /* The other components of midpoint are just zero */
midpoint[1] = 0.0f;
midpoint[2] = 0.0f;
- return 1;
+ return 1.0f;
}
/**
diff --git a/src/hydro/Shadowswift/voronoi3d_algorithm.h b/src/hydro/Shadowswift/voronoi3d_algorithm.h
index 165798bad1584173cb0df42894a1ca01091e2c0a..253db7e54fd7ea832924f93d4369eb32bc759fcf 100644
--- a/src/hydro/Shadowswift/voronoi3d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi3d_algorithm.h
@@ -21,6 +21,9 @@
#define SWIFT_VORONOI3D_ALGORITHM_H
#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "error.h"
#include "inline.h"
#include "voronoi3d_cell.h"
@@ -1219,6 +1222,13 @@ __attribute__((always_inline)) INLINE void voronoi_calculate_cell(
cell->centroid[0] += cell->x[0];
cell->centroid[1] += cell->x[1];
cell->centroid[2] += cell->x[2];
+
+ /* Reset the edges: we still need them for the face calculation */
+ for (i = 0; i < VORONOI3D_MAXNUMEDGE; ++i) {
+ if (cell->edges[i] < 0) {
+ cell->edges[i] = -1 - cell->edges[i];
+ }
+ }
}
/**
@@ -1389,24 +1399,36 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
}
max_radius = sqrtf(max_radius);
- return max_radius;
+ return 2.0f * max_radius;
}
/**
- * @brief Get the oriented surface area and midpoint of the face between a
- * 3D Voronoi cell and the given neighbour
+ * @brief Get the surface area and midpoint of the face between a 3D Voronoi
+ * cell and the given neighbour
*
* @param cell 3D Voronoi cell.
* @param ngb ID of a particle that is possibly a neighbour of this cell->
- * @param A Array to store the oriented surface area in.
* @param midpoint Array to store the relative position of the face in.
- * @return 0 if the given neighbour is not a neighbour, 1 otherwise.
+ * @return 0 if the given neighbour is not a neighbour, the surface area of
+ * the face otherwise.
*/
-__attribute__((always_inline)) INLINE int voronoi_get_face(
- struct voronoi_cell *cell, unsigned long long ngb, float *A,
- float *midpoint) {
+__attribute__((always_inline)) INLINE float voronoi_get_face(
+ struct voronoi_cell *cell, unsigned long long ngb, float *midpoint) {
- return 0;
+ int i = 0;
+ while (i < cell->nface && cell->ngbs[i] != ngb) {
+ i++;
+ }
+ if (i == cell->nface) {
+ /* Ngb not found */
+ return 0.0f;
+ }
+
+ midpoint[0] = cell->face_midpoints[i][0];
+ midpoint[1] = cell->face_midpoints[i][1];
+ midpoint[2] = cell->face_midpoints[i][2];
+
+ return cell->face_areas[i];
}
/**
diff --git a/tests/testVoronoi1D.c b/tests/testVoronoi1D.c
index 1b96dc1d1d98f8a7d7ce3044822631c1c4515496..dfc0d77b35d83976cfeeb4aed68f9c69d5ed395d 100644
--- a/tests/testVoronoi1D.c
+++ b/tests/testVoronoi1D.c
@@ -55,18 +55,18 @@ int main() {
}
/* Check face method */
- float A[3];
+ float A;
float midpoint[3];
- voronoi_get_face(&cell, 1, A, midpoint);
- if (A[0] != -1.0f) {
- error("Wrong surface normal returned for left neighbour!");
+ A = voronoi_get_face(&cell, 1, midpoint);
+ if (A != 1.0f) {
+ error("Wrong surface area returned for left neighbour!");
}
if (midpoint[0] != -0.25f) {
error("Wrong midpoint returned for left neighbour!");
}
- voronoi_get_face(&cell, 2, A, midpoint);
- if (A[0] != 1.0f) {
- error("Wrong surface normal returned for right neighbour!");
+ A = voronoi_get_face(&cell, 2, midpoint);
+ if (A != 1.0f) {
+ error("Wrong surface area returned for right neighbour!");
}
if (midpoint[0] != 0.25f) {
error("Wrong midpoint returned for right neighbour!");
diff --git a/tests/testVoronoi3D.c b/tests/testVoronoi3D.c
index 644a9508af2c8b4e1c2630367f972e8af5f40ee4..87c6def86a2a90ee94e6300b89660764667bc272 100644
--- a/tests/testVoronoi3D.c
+++ b/tests/testVoronoi3D.c
@@ -64,12 +64,54 @@ void test_calculate_cell() {
voronoi_initialize(&cell);
/* Calculate the volume and centroid of the large cube. */
voronoi_calculate_cell(&cell);
+ /* Calculate the faces. */
+ voronoi_calculate_faces(&cell);
/* Update these values if you ever change to another large cube! */
assert(cell.volume == 27.0f);
assert(cell.centroid[0] = 0.5f);
assert(cell.centroid[1] = 0.5f);
assert(cell.centroid[2] = 0.5f);
+
+ /* Check cell neighbours. */
+ assert(cell.nface == 6);
+ assert(cell.ngbs[0] == VORONOI3D_BOX_FRONT);
+ assert(cell.ngbs[1] == VORONOI3D_BOX_LEFT);
+ assert(cell.ngbs[2] == VORONOI3D_BOX_BOTTOM);
+ assert(cell.ngbs[3] == VORONOI3D_BOX_TOP);
+ assert(cell.ngbs[4] == VORONOI3D_BOX_BACK);
+ assert(cell.ngbs[5] == VORONOI3D_BOX_RIGHT);
+
+ /* Check cell faces */
+ assert(cell.face_areas[0] == 9.0f);
+ assert(cell.face_midpoints[0][0] == 0.5f);
+ assert(cell.face_midpoints[0][1] == -1.0f);
+ assert(cell.face_midpoints[0][2] == 0.5f);
+
+ assert(cell.face_areas[1] == 9.0f);
+ assert(cell.face_midpoints[1][0] == -1.0f);
+ assert(cell.face_midpoints[1][1] == 0.5f);
+ assert(cell.face_midpoints[1][2] == 0.5f);
+
+ assert(cell.face_areas[2] == 9.0f);
+ assert(cell.face_midpoints[2][0] == 0.5f);
+ assert(cell.face_midpoints[2][1] == 0.5f);
+ assert(cell.face_midpoints[2][2] == -1.0f);
+
+ assert(cell.face_areas[3] == 9.0f);
+ assert(cell.face_midpoints[3][0] == 0.5f);
+ assert(cell.face_midpoints[3][1] == 0.5f);
+ assert(cell.face_midpoints[3][2] == 2.0f);
+
+ assert(cell.face_areas[4] == 9.0f);
+ assert(cell.face_midpoints[4][0] == 0.5f);
+ assert(cell.face_midpoints[4][1] == 2.0f);
+ assert(cell.face_midpoints[4][2] == 0.5f);
+
+ assert(cell.face_areas[5] == 9.0f);
+ assert(cell.face_midpoints[5][0] == 2.0f);
+ assert(cell.face_midpoints[5][1] == 0.5f);
+ assert(cell.face_midpoints[5][2] == 0.5f);
}
int main() {
@@ -116,14 +158,22 @@ int main() {
error("Wrong centroid: %g %g %g!", cell.centroid[0], cell.centroid[1],
cell.centroid[2]);
}
- /* Check neighbour order */
- // TODO
-
- /* Check face method */
- float A[3];
- float midpoint[3];
- voronoi_get_face(&cell, 1, A, midpoint);
- // TODO
+
+ /* Check faces. */
+ float A, midpoint[3];
+ A = voronoi_get_face(&cell, 1, midpoint);
+ if (A) {
+ if (A != 0.25f) {
+ error("Wrong surface area: %g!", A);
+ }
+ if (fabs(midpoint[0] - 0.25f) > 1.e-5 || fabs(midpoint[1] - 0.5f) > 1.e-5 ||
+ fabs(midpoint[2] - 0.5f) > 1.e-5) {
+ error("Wrong face midpoint: %g %g %g!", midpoint[0], midpoint[1],
+ midpoint[2]);
+ }
+ } else {
+ error("Neighbour not found!");
+ }
return 0;
}