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

Prevent any kernel value to be negative due to rounding errors.

parent 383ec6d9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/kernel_hydro.h
+ 9
− 2
View file @ 2ce4bba0
......@@ -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;
......
tests/testKernel.c
+ 100
− 63
View file @ 2ce4bba0
......@@ -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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment