From 2ce4bba0c06ee55aaa824b5d05548c8df8841312 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 19 May 2018 15:26:28 +0200
Subject: [PATCH] Prevent any kernel value to be negative due to rounding
errors.
---
src/kernel_hydro.h | 11 ++-
tests/testKernel.c | 163 +++++++++++++++++++++++++++------------------
2 files changed, 109 insertions(+), 65 deletions(-)
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
index d3c905e50e..788a9037a6 100644
--- a/src/kernel_hydro.h
+++ b/src/kernel_hydro.h
@@ -38,6 +38,7 @@
#include "dimension.h"
#include "error.h"
#include "inline.h"
+#include "minmax.h"
#include "vector.h"
/* ------------------------------------------------------------------------- */
@@ -267,6 +268,9 @@ __attribute__((always_inline)) INLINE static void kernel_deval(
w = x * w + coeffs[k];
}
+ w = max(w, 0.f);
+ dw_dx = min(dw_dx, 0.f);
+
/* Return everything */
*W = w * kernel_constant * kernel_gamma_inv_dim;
*dW_dx = dw_dx * kernel_constant * kernel_gamma_inv_dim_plus_one;
@@ -300,6 +304,8 @@ __attribute__((always_inline)) INLINE static void kernel_eval(
/* ... and the rest of them */
for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
+ w = max(w, 0.f);
+
/* Return everything */
*W = w * kernel_constant * kernel_gamma_inv_dim;
}
@@ -331,9 +337,10 @@ __attribute__((always_inline)) INLINE static void kernel_eval_dWdx(
(float)(kernel_degree - 1) * coeffs[1];
/* ... and the rest of them */
- for (int k = 2; k < kernel_degree; k++) {
+ for (int k = 2; k < kernel_degree; k++)
dw_dx = dw_dx * x + (float)(kernel_degree - k) * coeffs[k];
- }
+
+ dw_dx = min(dw_dx, 0.f);
/* Return everything */
*dW_dx = dw_dx * kernel_constant * kernel_gamma_inv_dim_plus_one;
diff --git a/tests/testKernel.c b/tests/testKernel.c
index 0658639070..dc29d053c2 100644
--- a/tests/testKernel.c
+++ b/tests/testKernel.c
@@ -17,46 +17,72 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#include "../config.h"
+#include "align.h"
#include "kernel_hydro.h"
#include "vector.h"
+#include <fenv.h>
#include <stdlib.h>
#include <strings.h>
-#define numPoints (1 << 4)
+const int numPoints = (1 << 28);
int main() {
- const float h = 1.2348f;
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+/* Choke on FPEs */
+#ifdef HAVE_FE_ENABLE_EXCEPT
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+#endif
- float u[numPoints] = {0.f};
- float W[numPoints] = {0.f};
- float dW[numPoints] = {0.f};
+ const float h = 1.2348f;
- printf("\nSerial Output\n");
- printf("-------------\n");
+ float *u, *W, *dW;
+ if (posix_memalign((void **)&u, SWIFT_CACHE_ALIGNMENT,
+ numPoints * sizeof(float)) != 0)
+ error("Error allocating u");
+ if (posix_memalign((void **)&W, SWIFT_CACHE_ALIGNMENT,
+ numPoints * sizeof(float)) != 0)
+ error("Error allocating W");
+ if (posix_memalign((void **)&dW, SWIFT_CACHE_ALIGNMENT,
+ numPoints * sizeof(float)) != 0)
+ error("Error allocating dW");
+
+ message("Serial Output");
+ message("-------------");
const float numPoints_inv = 1. / numPoints;
- for (int i = 0; i < numPoints; ++i) {
- u[i] = i * 2.25f * numPoints_inv / h;
- }
+ for (int i = 0; i < numPoints; ++i)
+ u[i] = i * 1.2f * kernel_gamma * numPoints_inv / h;
for (int i = 0; i < numPoints; ++i) {
kernel_deval(u[i], &W[i], &dW[i]);
- printf("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i, h,
- h * kernel_gamma, u[i] * h, W[i], dW[i]);
+ if (W[i] < 0.f) error("Kernel is negative u=%e W=%e", u[i], W[i]);
+ if (dW[i] > 0.f)
+ error("Kernel derivatibe is positive u=%e dW=%e", u[i], dW[i]);
}
- printf("\nVector Output for VEC_SIZE=%d\n", VEC_SIZE);
- printf("-------------\n");
+ /* Test some additional special cases */
+ float Wtest, dWtest;
+ kernel_deval(1.930290, &Wtest, &dWtest);
+ if (Wtest < 0.f) error("Kernel is negative u=%e W=%e", 1.930290, Wtest);
+ if (dWtest > 0.f)
+ error("Kernel derivative is positive u=%e dW=%e", 1.930290, dWtest);
#ifdef WITH_VECTORIZATION
- printf("\nVector Output for kernel_deval_1_vec\n");
- printf("-------------\n");
+ message("Vector Output for VEC_SIZE=%d", VEC_SIZE);
+ message("-------------");
+
+ message("Vector Output for kernel_deval_1_vec");
+ message("-------------");
/* Test vectorised kernel that uses one vector. */
for (int i = 0; i < numPoints; i += VEC_SIZE) {
@@ -64,33 +90,35 @@ int main() {
vector vx, vx_h;
vector W_vec, dW_vec;
- for (int j = 0; j < VEC_SIZE; j++) {
- vx.f[j] = (i + j) * 2.25f / numPoints;
- }
+ for (int j = 0; j < VEC_SIZE; j++)
+ vx.f[j] = (i + j) * 1.2f * kernel_gamma / numPoints;
vx_h.v = vec_mul(vx.v, vec_set1(1.f / h));
kernel_deval_1_vec(&vx_h, &W_vec, &dW_vec);
for (int j = 0; j < VEC_SIZE; j++) {
- printf("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i + j, h,
- h * kernel_gamma, vx.f[j], W_vec.f[j], dW_vec.f[j]);
-
- if (fabsf(W_vec.f[j] - W[i + j]) > 2e-7) {
- printf("Invalid value ! scalar= %e, vector= %e\n", W[i + j],
- W_vec.f[j]);
- return 1;
- }
- if (fabsf(dW_vec.f[j] - dW[i + j]) > 2e-7) {
- printf("Invalid value ! scalar= %e, vector= %e\n", dW[i + j],
- dW_vec.f[j]);
- return 1;
- }
+
+ /* message("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i + j, h, */
+ /* h * kernel_gamma, vx.f[j], W_vec.f[j], dW_vec.f[j]); */
+
+ if (W_vec.f[j] < 0.f)
+ error("Kernel is negative u=%e W=%e", u[i + j], W_vec.f[j]);
+ if (dW_vec.f[j] > 0.f)
+ error("Kernel derivative is positive u=%e dW=%e", u[i + j],
+ dW_vec.f[j]);
+
+ if (fabsf(W_vec.f[j] - W[i + j]) > 2e-6)
+ error("Invalid Wvalue ! scalar= %e, vector= %e\n", W[i + j],
+ W_vec.f[j]);
+ if (fabsf(dW_vec.f[j] - dW[i + j]) > 2e-6)
+ error("Invalid dW value ! scalar= %e, vector= %e %e %e\n", dW[i + j],
+ dW_vec.f[j], fabsf(dW_vec.f[j] - dW[i + j]), fabsf(dW[i + j]));
}
}
- printf("\nVector Output for kernel_deval_2_vec\n");
- printf("-------------\n");
+ message("Vector Output for kernel_deval_2_vec");
+ message("-------------");
/* Test vectorised kernel that uses two vectors. */
for (int i = 0; i < numPoints; i += VEC_SIZE) {
@@ -102,8 +130,8 @@ int main() {
vector W_vec_2, dW_vec_2;
for (int j = 0; j < VEC_SIZE; j++) {
- vx.f[j] = (i + j) * 2.25f / numPoints;
- vx_2.f[j] = (i + j) * 2.25f / numPoints;
+ vx.f[j] = (i + j) * 1.2f * kernel_gamma / numPoints;
+ vx_2.f[j] = (i + j) * 1.2f * kernel_gamma / numPoints;
}
vx_h.v = vec_mul(vx.v, vec_set1(1.f / h));
@@ -113,41 +141,50 @@ int main() {
/* Check first vector results. */
for (int j = 0; j < VEC_SIZE; j++) {
- printf("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i + j, h,
- h * kernel_gamma, vx.f[j], W_vec.f[j], dW_vec.f[j]);
-
- if (fabsf(W_vec.f[j] - W[i + j]) > 2e-7) {
- printf("Invalid value ! scalar= %e, vector= %e\n", W[i + j],
- W_vec.f[j]);
- return 1;
- }
- if (fabsf(dW_vec.f[j] - dW[i + j]) > 2e-7) {
- printf("Invalid value ! scalar= %e, vector= %e\n", dW[i + j],
- dW_vec.f[j]);
- return 1;
- }
+
+ /* message("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i + j, h, */
+ /* h * kernel_gamma, vx.f[j], W_vec.f[j], dW_vec.f[j]); */
+
+ if (W_vec.f[j] < 0.f)
+ error("Kernel is negative u=%e W=%e", u[i + j], W_vec.f[j]);
+ if (dW_vec.f[j] > 0.f)
+ error("Kernel derivative is positive u=%e dW=%e", u[i + j],
+ dW_vec.f[j]);
+
+ if (fabsf(W_vec.f[j] - W[i + j]) > 2e-6)
+ error("Invalid value ! scalar= %e, vector= %e\n", W[i + j], W_vec.f[j]);
+ if (fabsf(dW_vec.f[j] - dW[i + j]) > 2e-6)
+ error("Invalid value ! scalar= %e, vector= %e\n", dW[i + j],
+ dW_vec.f[j]);
}
/* Check second vector results. */
for (int j = 0; j < VEC_SIZE; j++) {
- printf("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i + j, h,
- h * kernel_gamma, vx_2.f[j], W_vec_2.f[j], dW_vec_2.f[j]);
-
- if (fabsf(W_vec_2.f[j] - W[i + j]) > 2e-7) {
- printf("Invalid value ! scalar= %e, vector= %e\n", W[i + j],
- W_vec_2.f[j]);
- return 1;
- }
- if (fabsf(dW_vec_2.f[j] - dW[i + j]) > 2e-7) {
- printf("Invalid value ! scalar= %e, vector= %e\n", dW[i + j],
- dW_vec_2.f[j]);
- return 1;
- }
+
+ /* message("%2d: h= %f H= %f x=%f W(x,h)=%f dW(x,h)=%f\n", i + j, h, */
+ /* h * kernel_gamma, vx_2.f[j], W_vec_2.f[j], dW_vec_2.f[j]); */
+
+ if (W_vec_2.f[j] < 0.f)
+ error("Kernel is negative u=%e W=%e", u[i + j], W_vec_2.f[j]);
+ if (dW_vec_2.f[j] > 0.f)
+ error("Kernel derivative is positive u=%e dW=%e", u[i + j],
+ dW_vec_2.f[j]);
+
+ if (fabsf(W_vec_2.f[j] - W[i + j]) > 2e-6)
+ error("Invalid value ! scalar= %e, vector= %e\n", W[i + j],
+ W_vec_2.f[j]);
+ if (fabsf(dW_vec_2.f[j] - dW[i + j]) > 2e-6)
+ error("Invalid value ! scalar= %e, vector= %e\n", dW[i + j],
+ dW_vec_2.f[j]);
}
}
- printf("\nAll values are consistent\n");
+ message("All values are consistent");
#endif
+
+ free(u);
+ free(W);
+ free(dW);
return 0;
}
--
GitLab