Commit ce09e653 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Split the hydro and gravity kernels into separate files

parent 9c93faae
......@@ -41,11 +41,11 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel.c tools.c part.c partition.c clocks.c
kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel.h vector.h \
runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
......
......@@ -42,7 +42,7 @@
/* Local includes. */
#include "const.h"
#include "error.h"
#include "kernel.h"
#include "kernel_hydro.h"
#include "version.h"
const char* particle_type_names[NUM_PARTICLE_TYPES] = {
......
......@@ -24,6 +24,7 @@
#include "../config.h"
/* Includes. */
#include "kernel_hydro.h"
#include "part.h"
#include "units.h"
......
......@@ -22,7 +22,7 @@
/* Includes. */
#include "const.h"
#include "kernel.h"
#include "kernel_gravity.h"
#include "vector.h"
/**
......
......@@ -22,7 +22,7 @@
/* Includes. */
#include "const.h"
#include "kernel.h"
#include "kernel_hydro.h"
#include "part.h"
#include "vector.h"
......
......@@ -22,7 +22,7 @@
/* Includes. */
#include "const.h"
#include "kernel.h"
#include "kernel_hydro.h"
#include "part.h"
#include "vector.h"
......
......@@ -21,7 +21,7 @@
/* Includes. */
#include "const.h"
#include "kernel.h"
#include "kernel_hydro.h"
#include "part.h"
#include "vector.h"
......
......@@ -21,32 +21,7 @@
#include <math.h>
#include <stdio.h>
#include "kernel.h"
/**
* @brief Test the SPH kernel function by dumping it in the interval [0,1].
*
* @param N number of intervals in [0,1].
*/
void SPH_kernel_dump(int N) {
int k;
float x, w, dw_dx;
float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
// float dw_dx4[4] __attribute__ ((aligned (16)));
for (k = 0; k <= N; k++) {
x = ((float)k) / N;
x4[3] = x4[2];
x4[2] = x4[1];
x4[1] = x4[0];
x4[0] = x;
kernel_deval(x, &w, &dw_dx);
// kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
}
}
#include "kernel_gravity.h"
/**
* @brief The Gadget-2 gravity kernel function
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_KERNEL_GRAVITY_H
#define SWIFT_KERNEL_GRAVITY_H
/* Includes. */
#include "const.h"
#include "inline.h"
#include "vector.h"
/* Gravity kernel stuff
* -----------------------------------------------------------------------------------------------
*/
/* The gravity kernel is defined as a degree 6 polynomial in the distance
r. The resulting value should be post-multiplied with r^-3, resulting
in a polynomial with terms ranging from r^-3 to r^3, which are
sufficient to model both the direct potential as well as the splines
near the origin. */
/* Coefficients for the gravity kernel. */
#define kernel_grav_degree 6
#define kernel_grav_ivals 2
#define kernel_grav_scale (2 * const_iepsilon)
static float kernel_grav_coeffs
[(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)] = {
32.0f * const_iepsilon6, -192.0f / 5.0f * const_iepsilon5,
0.0f, 32.0f / 3.0f * const_iepsilon3,
0.0f, 0.0f,
0.0f, -32.0f / 3.0f * const_iepsilon6,
192.0f / 5.0f * const_iepsilon5, -48.0f * const_iepsilon4,
64.0f / 3.0f * const_iepsilon3, 0.0f,
0.0f, -1.0f / 15.0f,
0.0f, 0.0f,
0.0f, 0.0f,
0.0f, 0.0f,
1.0f};
/**
* @brief Computes the gravity cubic spline for a given distance x.
*/
__attribute__((always_inline)) INLINE static void kernel_grav_eval(float x,
float *W) {
int ind = fmin(x * kernel_grav_scale, kernel_grav_ivals);
float *coeffs = &kernel_grav_coeffs[ind * (kernel_grav_degree + 1)];
float w = coeffs[0] * x + coeffs[1];
for (int k = 2; k <= kernel_grav_degree; k++) w = x * w + coeffs[k];
*W = w;
}
#ifdef VECTORIZE
/**
* @brief Computes the gravity cubic spline for a given distance x (Vectorized
* version).
*/
__attribute__((always_inline))
INLINE static void kernel_grav_eval_vec(vector *x, vector *w) {
vector ind, c[kernel_grav_degree + 1];
int j, k;
/* Load x and get the interval id. */
ind.m = vec_ftoi(vec_fmin(x->v * vec_set1(kernel_grav_scale),
vec_set1((float)kernel_grav_ivals)));
/* load the coefficients. */
for (k = 0; k < VEC_SIZE; k++)
for (j = 0; j < kernel_grav_degree + 1; j++)
c[j].f[k] = kernel_grav_coeffs[ind.i[k] * (kernel_grav_degree + 1) + j];
/* Init the iteration for Horner's scheme. */
w->v = (c[0].v * x->v) + c[1].v;
/* And we're off! */
for (int k = 2; k <= kernel_grav_degree; k++) w->v = (x->v * w->v) + c[k].v;
}
#endif
/* Blending function stuff
* --------------------------------------------------------------------------------------------
*/
/* Coefficients for the blending function. */
#define blender_degree 3
#define blender_ivals 3
#define blender_scale 4.0f
static float blender_coeffs[(blender_degree + 1) * (blender_ivals + 1)] = {
0.0f, 0.0f, 0.0f, 1.0f, -32.0f, 24.0f, -6.0f, 1.5f,
-32.0f, 72.0f, -54.0f, 13.5f, 0.0f, 0.0f, 0.0f, 0.0f};
/**
* @brief Computes the cubic spline blender for a given distance x.
*/
__attribute__((always_inline)) INLINE static void blender_eval(float x,
float *W) {
int ind = fmin(x * blender_scale, blender_ivals);
float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
float w = coeffs[0] * x + coeffs[1];
for (int k = 2; k <= blender_degree; k++) w = x * w + coeffs[k];
*W = w;
}
/**
* @brief Computes the cubic spline blender and its derivative for a given
* distance x.
*/
__attribute__((always_inline)) INLINE static void blender_deval(float x,
float *W,
float *dW_dx) {
int ind = fminf(x * blender_scale, blender_ivals);
float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
float w = coeffs[0] * x + coeffs[1];
float dw_dx = coeffs[0];
for (int k = 2; k <= blender_degree; k++) {
dw_dx = dw_dx * x + w;
w = x * w + coeffs[k];
}
*W = w;
*dW_dx = dw_dx;
}
#ifdef VECTORIZE
/**
* @brief Computes the cubic spline blender and its derivative for a given
* distance x (Vectorized version). Gives a sensible answer only if x<2.
*/
__attribute__((always_inline)) INLINE static void blender_eval_vec(vector *x,
vector *w) {
vector ind, c[blender_degree + 1];
int j, k;
/* Load x and get the interval id. */
ind.m = vec_ftoi(
vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
/* load the coefficients. */
for (k = 0; k < VEC_SIZE; k++)
for (j = 0; j < blender_degree + 1; j++)
c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
/* Init the iteration for Horner's scheme. */
w->v = (c[0].v * x->v) + c[1].v;
/* And we're off! */
for (int k = 2; k <= blender_degree; k++) w->v = (x->v * w->v) + c[k].v;
}
/**
* @brief Computes the cubic spline blender and its derivative for a given
* distance x (Vectorized version). Gives a sensible answer only if x<2.
*/
__attribute__((always_inline))
INLINE static void blender_deval_vec(vector *x, vector *w, vector *dw_dx) {
vector ind, c[blender_degree + 1];
int j, k;
/* Load x and get the interval id. */
ind.m = vec_ftoi(
vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
/* load the coefficients. */
for (k = 0; k < VEC_SIZE; k++)
for (j = 0; j < blender_degree + 1; j++)
c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
/* Init the iteration for Horner's scheme. */
w->v = (c[0].v * x->v) + c[1].v;
dw_dx->v = c[0].v;
/* And we're off! */
for (int k = 2; k <= blender_degree; k++) {
dw_dx->v = (dw_dx->v * x->v) + w->v;
w->v = (x->v * w->v) + c[k].v;
}
}
#endif
void gravity_kernel_dump(float r_max, int N);
#endif // SWIFT_KERNEL_GRAVITY_H
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* 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/>.
*
******************************************************************************/
#include <math.h>
#include <stdio.h>
#include "kernel_hydro.h"
/**
* @brief Test the SPH kernel function by dumping it in the interval [0,1].
*
* @param N number of intervals in [0,1].
*/
void hydro_kernel_dump(int N) {
int k;
float x, w, dw_dx;
float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
// float dw_dx4[4] __attribute__ ((aligned (16)));
for (k = 0; k <= N; k++) {
x = ((float)k) / N;
x4[3] = x4[2];
x4[2] = x4[1];
x4[1] = x4[0];
x4[0] = x;
kernel_deval(x, &w, &dw_dx);
// kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
}
}
......@@ -17,202 +17,14 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_KERNEL_H
#define SWIFT_KERNEL_H
#ifndef SWIFT_KERNEL_HYDRO_H
#define SWIFT_KERNEL_HYDRO_H
/* Includes. */
#include "const.h"
#include "inline.h"
#include "vector.h"
/**
* @file kernel.h
* @brief SPH kernel functions. Compute W(x,h) and the gradient of W(x,h),
* as well as the blending function used for gravity.
*/
/* Gravity kernel stuff
* -----------------------------------------------------------------------------------------------
*/
/* The gravity kernel is defined as a degree 6 polynomial in the distance
r. The resulting value should be post-multiplied with r^-3, resulting
in a polynomial with terms ranging from r^-3 to r^3, which are
sufficient to model both the direct potential as well as the splines
near the origin. */
/* Coefficients for the gravity kernel. */
#define kernel_grav_degree 6
#define kernel_grav_ivals 2
#define kernel_grav_scale (2 * const_iepsilon)
static float kernel_grav_coeffs
[(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)] = {
32.0f * const_iepsilon6, -192.0f / 5.0f * const_iepsilon5,
0.0f, 32.0f / 3.0f * const_iepsilon3,
0.0f, 0.0f,
0.0f, -32.0f / 3.0f * const_iepsilon6,
192.0f / 5.0f * const_iepsilon5, -48.0f * const_iepsilon4,
64.0f / 3.0f * const_iepsilon3, 0.0f,
0.0f, -1.0f / 15.0f,
0.0f, 0.0f,
0.0f, 0.0f,
0.0f, 0.0f,
1.0f};
/**
* @brief Computes the gravity cubic spline for a given distance x.
*/
__attribute__((always_inline)) INLINE static void kernel_grav_eval(float x,
float *W) {
int ind = fmin(x * kernel_grav_scale, kernel_grav_ivals);
float *coeffs = &kernel_grav_coeffs[ind * (kernel_grav_degree + 1)];
float w = coeffs[0] * x + coeffs[1];
for (int k = 2; k <= kernel_grav_degree; k++) w = x * w + coeffs[k];
*W = w;
}
#ifdef VECTORIZE
/**
* @brief Computes the gravity cubic spline for a given distance x (Vectorized
* version).
*/
__attribute__((always_inline))
INLINE static void kernel_grav_eval_vec(vector *x, vector *w) {
vector ind, c[kernel_grav_degree + 1];
int j, k;
/* Load x and get the interval id. */
ind.m = vec_ftoi(vec_fmin(x->v * vec_set1(kernel_grav_scale),
vec_set1((float)kernel_grav_ivals)));
/* load the coefficients. */
for (k = 0; k < VEC_SIZE; k++)
for (j = 0; j < kernel_grav_degree + 1; j++)
c[j].f[k] = kernel_grav_coeffs[ind.i[k] * (kernel_grav_degree + 1) + j];
/* Init the iteration for Horner's scheme. */
w->v = (c[0].v * x->v) + c[1].v;
/* And we're off! */
for (int k = 2; k <= kernel_grav_degree; k++) w->v = (x->v * w->v) + c[k].v;
}
#endif
/* Blending function stuff
* --------------------------------------------------------------------------------------------
*/
/* Coefficients for the blending function. */
#define blender_degree 3
#define blender_ivals 3
#define blender_scale 4.0f
static float blender_coeffs[(blender_degree + 1) * (blender_ivals + 1)] = {
0.0f, 0.0f, 0.0f, 1.0f, -32.0f, 24.0f, -6.0f, 1.5f,
-32.0f, 72.0f, -54.0f, 13.5f, 0.0f, 0.0f, 0.0f, 0.0f};
/**
* @brief Computes the cubic spline blender for a given distance x.
*/
__attribute__((always_inline)) INLINE static void blender_eval(float x,
float *W) {
int ind = fmin(x * blender_scale, blender_ivals);
float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
float w = coeffs[0] * x + coeffs[1];
for (int k = 2; k <= blender_degree; k++) w = x * w + coeffs[k];
*W = w;
}
/**
* @brief Computes the cubic spline blender and its derivative for a given
* distance x.
*/
__attribute__((always_inline)) INLINE static void blender_deval(float x,
float *W,
float *dW_dx) {
int ind = fminf(x * blender_scale, blender_ivals);
float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
float w = coeffs[0] * x + coeffs[1];
float dw_dx = coeffs[0];
for (int k = 2; k <= blender_degree; k++) {
dw_dx = dw_dx * x + w;
w = x * w + coeffs[k];
}
*W = w;
*dW_dx = dw_dx;
}
#ifdef VECTORIZE
/**
* @brief Computes the cubic spline blender and its derivative for a given
* distance x (Vectorized version). Gives a sensible answer only if x<2.
*/
__attribute__((always_inline)) INLINE static void blender_eval_vec(vector *x,
vector *w) {
vector ind, c[blender_degree + 1];
int j, k;
/* Load x and get the interval id. */
ind.m = vec_ftoi(
vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
/* load the coefficients. */
for (k = 0; k < VEC_SIZE; k++)
for (j = 0; j < blender_degree + 1; j++)
c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
/* Init the iteration for Horner's scheme. */
w->v = (c[0].v * x->v) + c[1].v;
/* And we're off! */
for (int k = 2; k <= blender_degree; k++) w->v = (x->v * w->v) + c[k].v;
}
/**
* @brief Computes the cubic spline blender and its derivative for a given
* distance x (Vectorized version). Gives a sensible answer only if x<2.
*/
__attribute__((always_inline))
INLINE static void blender_deval_vec(vector *x, vector *w, vector *dw_dx) {
vector ind, c[blender_degree + 1];
int j, k;
/* Load x and get the interval id. */
ind.m = vec_ftoi(
vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
/* load the coefficients. */
for (k = 0; k < VEC_SIZE; k++)
for (j = 0; j < blender_degree + 1; j++)
c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
/* Init the iteration for Horner's scheme. */
w->v = (c[0].v * x->v) + c[1].v;
dw_dx->v = c[0].v;
/* And we're off! */
for (int k = 2; k <= blender_degree; k++) {
dw_dx->v = (dw_dx->v * x->v) + w->v;
w->v = (x->v * w->v) + c[k].v;
}
}
#endif
/* --------------------------------------------------------------------------------------------------------------------
*/
#if defined(CUBIC_SPLINE_KERNEL)
/* --------------------------------------------------------------------------------------------------------------------
......@@ -611,7 +423,6 @@ __attribute__((always_inline)) INLINE static void kernel_eval(float x,
#endif // Kernel choice
/* Some cross-check functions */
void SPH_kernel_dump(int N);
void gravity_kernel_dump(float r_max, int N);
void hydro_kernel_dump(int N);
#endif // SWIFT_KERNEL_H
#endif // SWIFT_KERNEL_HYDRO_H
......@@ -25,7 +25,7 @@
/* Includes. */
#include "const.h"
#include "inline.h"
#include "kernel.h"
#include "kernel_gravity.h"
#include "part.h"
/* Some constants. */
......
......@@ -43,7 +43,7 @@
#include "cycle.h"
#include "error.h"
#include "intrinsics.h"
#include "kernel.h"
#include "kernel_hydro.h"
#include "timers.h"
/**
......
......@@ -39,7 +39,7 @@
#include "atomic.h"
#include "engine.h"
#include "error.h"
#include "kernel.h"
#include "kernel_hydro.h"
#include "lock.h"
#include "minmax.h"
#include "runner.h"
......