Skip to content
Snippets Groups Projects
Commit d2f1f017 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'shadowswift_new' into 'master'

Voronoi tests work now.

I think I added all missing files used by `SHADOWFAX_SPH` and `GIZMO_SPH` to `Makefile.am`.

I have also introduced a tolerance for geometric tests in the 2D Voronoi algorithm, which fixes round off issues that caused `testVoronoi2D` to fail if compiled with optimizations enabled.

This fixes #282.

See merge request !325
parents b8f5850f 0bed4834
No related branches found
No related tags found
1 merge request!325Voronoi tests work now.
...@@ -77,8 +77,31 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h ...@@ -77,8 +77,31 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \ hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \ hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \ hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \ hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h \
hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \ hydro/Gizmo/hydro_io.h hydro/Gizmo/hydro_debug.h \
hydro/Gizmo/hydro_part.h \
hydro/Gizmo/hydro_gradients_gizmo.h \
hydro/Gizmo/hydro_gradients.h \
hydro/Gizmo/hydro_gradients_sph.h \
hydro/Gizmo/hydro_slope_limiters_cell.h \
hydro/Gizmo/hydro_slope_limiters_face.h \
hydro/Gizmo/hydro_slope_limiters.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
hydro/Shadowswift/hydro_io.h \
hydro/Shadowswift/hydro_part.h \
hydro/Shadowswift/hydro_slope_limiters_cell.h \
hydro/Shadowswift/hydro_slope_limiters_face.h \
hydro/Shadowswift/hydro_slope_limiters.h \
hydro/Shadowswift/voronoi1d_algorithm.h \
hydro/Shadowswift/voronoi1d_cell.h \
hydro/Shadowswift/voronoi2d_algorithm.h \
hydro/Shadowswift/voronoi2d_cell.h \
hydro/Shadowswift/voronoi3d_algorithm.h \
hydro/Shadowswift/voronoi3d_cell.h \
hydro/Shadowswift/voronoi_algorithm.h \
hydro/Shadowswift/voronoi_cell.h \
riemann/riemann_hllc.h riemann/riemann_trrs.h \ riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \ riemann/riemann_exact.h riemann/riemann_vacuum.h \
stars.h stars_io.h \ stars.h stars_io.h \
......
...@@ -43,6 +43,8 @@ ...@@ -43,6 +43,8 @@
#define VORONOI2D_BOX_TOP 18446744073709551604llu #define VORONOI2D_BOX_TOP 18446744073709551604llu
#define VORONOI2D_BOX_BOTTOM 18446744073709551605llu #define VORONOI2D_BOX_BOTTOM 18446744073709551605llu
#define VORONOI2D_TOLERANCE 1.e-6f
/** /**
* @brief Initialize a 2D Voronoi cell. * @brief Initialize a 2D Voronoi cell.
* *
...@@ -152,7 +154,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact( ...@@ -152,7 +154,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
/* start testing a random vertex: the first one */ /* start testing a random vertex: the first one */
test = cell->vertices[0][0] * half_dx[0] + cell->vertices[0][1] * half_dx[1] - test = cell->vertices[0][0] * half_dx[0] + cell->vertices[0][1] * half_dx[1] -
r2; r2;
if (test < 0.) { if (test < -VORONOI2D_TOLERANCE) {
/* vertex is below midline */ /* vertex is below midline */
#ifdef VORONOI_VERBOSE #ifdef VORONOI_VERBOSE
message("First vertex is below midline (%g %g --> %g)!", message("First vertex is below midline (%g %g --> %g)!",
...@@ -215,7 +217,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact( ...@@ -215,7 +217,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
i = 1; i = 1;
test = cell->vertices[i][0] * half_dx[0] + test = cell->vertices[i][0] * half_dx[0] +
cell->vertices[i][1] * half_dx[1] - r2; cell->vertices[i][1] * half_dx[1] - r2;
while (i < cell->nvert && test >= 0.) { while (i < cell->nvert && test > -VORONOI2D_TOLERANCE) {
/* make sure we always store the most recent test result */ /* make sure we always store the most recent test result */
a1 = test; a1 = test;
++i; ++i;
...@@ -272,7 +274,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact( ...@@ -272,7 +274,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
r2; r2;
/* this loop can never deadlock, as we know there is at least 1 vertex below /* this loop can never deadlock, as we know there is at least 1 vertex below
the midline */ the midline */
while (test >= 0.) { while (test > -VORONOI2D_TOLERANCE) {
/* make sure we always store the most recent test result */ /* make sure we always store the most recent test result */
a2 = test; a2 = test;
i += increment; i += increment;
......
...@@ -206,7 +206,6 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent( ...@@ -206,7 +206,6 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
float fa, fb, fc, fs; float fa, fb, fc, fs;
float tmp, tmp2; float tmp, tmp2;
int mflag; int mflag;
int i;
a = lower_limit; a = lower_limit;
b = upper_limit; b = upper_limit;
...@@ -243,7 +242,6 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent( ...@@ -243,7 +242,6 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
c = a; c = a;
fc = fa; fc = fa;
mflag = 1; mflag = 1;
i = 0;
while (!(fb == 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b))) { while (!(fb == 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b))) {
if ((fa != fc) && (fb != fc)) /* Inverse quadratic interpolation */ if ((fa != fc) && (fb != fc)) /* Inverse quadratic interpolation */
...@@ -286,7 +284,6 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent( ...@@ -286,7 +284,6 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
fa = fb; fa = fb;
fb = tmp; fb = tmp;
} }
i++;
} }
return b; return b;
} }
......
...@@ -197,12 +197,11 @@ int main() { ...@@ -197,12 +197,11 @@ int main() {
/* Check the neighbour relations for an arbitrary cell: cell 44 /* Check the neighbour relations for an arbitrary cell: cell 44
We plotted the grid and manually found the correct neighbours and their We plotted the grid and manually found the correct neighbours and their
order. */ order. */
assert(cells[44].nvert == 5); assert(cells[44].nvert == 4);
assert(cells[44].ngbs[0] == 34); assert(cells[44].ngbs[0] == 34);
assert(cells[44].ngbs[1] == 35); assert(cells[44].ngbs[1] == 45);
assert(cells[44].ngbs[2] == 45); assert(cells[44].ngbs[2] == 54);
assert(cells[44].ngbs[3] == 54); assert(cells[44].ngbs[3] == 43);
assert(cells[44].ngbs[4] == 43);
message("Done."); message("Done.");
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment