Commit 7ce42890 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Moved the functions that dump kernels for cross-check to a new file...

Moved the functions that dump kernels for cross-check to a new file src/kernel.c and the factorisation routine to src/tools.c
parent 4338674e
......@@ -240,70 +240,6 @@ void pairs_single_grav(double *dim, long long int pid,
pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
}
/**
* @brief Test the kernel function by dumping it in the interval [0,1].
*
* @param N number of intervals in [0,1].
*/
void 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]);
}
}
float gadget(float r) {
float fac, h_inv, u, r2 = r * r;
if (r >= const_epsilon)
fac = 1.0f / (r2 * r);
else {
h_inv = 1. / const_epsilon;
u = r * h_inv;
if (u < 0.5)
fac = const_iepsilon3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
else
fac = const_iepsilon3 *
(21.333333333333 - 48.0 * u + 38.4 * u * u -
10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
return const_G * fac;
}
void gravity_dump(float r_max, int N) {
int k;
float x, w;
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 = 1; k <= N; k++) {
x = (r_max * k) / N;
x4[3] = x4[2];
x4[2] = x4[1];
x4[1] = x4[0];
x4[0] = x;
kernel_grav_eval(x, &w);
w *= const_G / (x * x * x);
// blender_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
printf(" %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n", x, w * x, w4[0],
w4[1], w4[2], w4[3], gadget(x) * x);
}
}
/**
* @brief Test the density function by dumping it for two random parts.
......@@ -354,22 +290,6 @@ void density_dump(int N) {
}
}
/**
* Factorize a given integer, attempts to keep larger pair of factors.
*/
void factor(int value, int *f1, int *f2) {
int j;
int i;
j = (int)sqrt(value);
for (i = j; i > 0; i--) {
if ((value % i) == 0) {
*f1 = i;
*f2 = value / i;
break;
}
}
}
/**
* @brief Main routine that loads a few particles and generates some output.
......
......@@ -35,12 +35,13 @@ endif
# List required headers
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
common_io.h single_io.h multipole.h map.h
common_io.h single_io.h multipole.h map.h tools.h
# Common source files
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
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel.c tools.c
# Include files for distribution, not installation.
noinst_HEADERS = atomic.h cycle.h error.h inline.h kernel.h vector.h \
......
......@@ -23,14 +23,12 @@
#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 kernel_dump(int N) {
void SPH_kernel_dump(int N) {
int k;
float x, w, dw_dx;
......@@ -50,7 +48,6 @@ void kernel_dump(int N) {
}
}
/**
* @brief The Gadget-2 gravity kernel function
*
......@@ -73,14 +70,13 @@ float gadget(float r) {
return const_G * fac;
}
/**
* @brief Test the gravity kernel function by dumping it in the interval [0,1].
*
* @param r_max The radius up to which the kernel is dumped.
* @param N number of intervals in [0,1].
*/
void gravity_dump(float r_max, int N) {
void gravity_kernel_dump(float r_max, int N) {
int k;
float x, w;
......
......@@ -609,4 +609,8 @@ __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);
#endif // SWIFT_KERNEL_H
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2015 Peter W. Draper (p.w.draper@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>
/**
* Factorize a given integer, attempts to keep larger pair of factors.
*/
void factor(int value, int *f1, int *f2) {
int j;
int i;
j = (int)sqrt(value);
for (i = j; i > 0; i--) {
if ((value % i) == 0) {
*f1 = i;
*f2 = value / i;
break;
}
}
}
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2015 Peter W. Draper (p.w.draper@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/>.
*
******************************************************************************/
void factor(int value, int *f1, int *f2);
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment