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

Merged latest master branch changes

parents 287c2882 bd9ec753
......@@ -42,10 +42,13 @@ tests/swift_dopair_standard.dat
tests/brute_force_perturbed.dat
tests/swift_dopair_perturbed.dat
tests/test27cells
tests/test125cells
tests/brute_force_27_standard.dat
tests/swift_dopair_27_standard.dat
tests/brute_force_27_perturbed.dat
tests/swift_dopair_27_perturbed.dat
tests/brute_force_125_standard.dat
tests/swift_dopair_125_standard.dat
tests/testGreetings
tests/testReading
tests/input.hdf5
......@@ -53,12 +56,14 @@ tests/testSingle
tests/testTimeIntegration
tests/testSPHStep
tests/testKernel
tests/testInteractions
tests/testSymmetry
tests/testMaths
tests/testParser
tests/parser_output.yml
tests/test27cells.sh
tests/test27cellsPerturbed.sh
tests/test125cells.sh
tests/testPair.sh
tests/testPairPerturbed.sh
tests/testParser.sh
......
......@@ -476,6 +476,7 @@ AC_CONFIG_FILES([tests/testPair.sh], [chmod +x tests/testPair.sh])
AC_CONFIG_FILES([tests/testPairPerturbed.sh], [chmod +x tests/testPairPerturbed.sh])
AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])
# Report general configuration.
......
......@@ -32,6 +32,7 @@
/* Local includes. */
#include "config.h"
#include "const.h"
#include "hydro.h"
#include "inline.h"
#include "part.h"
......
......@@ -20,6 +20,8 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include <float.h>
/**
* @brief Computes the hydro time-step of a given particle
*
......@@ -263,3 +265,39 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
return p->u;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
return p->force.soundspeed;
}
......@@ -283,6 +283,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
}
/**
* @brief Returns the entropy of a particle
*
......@@ -294,3 +306,15 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
return p->entropy + p->entropy_dt * dt;
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
return p->force.soundspeed;
}
......@@ -29,8 +29,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh, p->force.P_over_rho2,
p->entropy, p->entropy_dt, p->force.soundspeed, p->density.div_v,
p->density.rot_v[0], p->density.rot_v[1], p->density.rot_v[2],
p->force.balsara, p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
}
......@@ -246,3 +246,39 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
return p->u;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
return p->u * p->rho * hydro_gamma_minus_one;
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
return sqrt(p->u * hydro_gamma * hydro_gamma_minus_one);
}
......@@ -40,7 +40,7 @@ void hydro_props_init(struct hydro_props *p,
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
const float eta3 = p->eta_neighbours * p->eta_neighbours * p->eta_neighbours;
p->target_neighbours = 4.0 * M_PI * kernel_gamma3 * eta3 / 3.0;
p->target_neighbours = eta3 * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
/* Ghost stuff */
......
......@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static void kick_part(
if (p->gpart != NULL) {
a_tot[0] += p->gpart->a_grav[0];
a_tot[1] += p->gpart->a_grav[1];
a_tot[1] += p->gpart->a_grav[2];
a_tot[2] += p->gpart->a_grav[2];
}
/* Kick particles in momentum space */
......
......@@ -121,8 +121,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci,
float dxq[3 * VEC_SIZE] __attribute__((aligned(16)));
struct part *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
/* Anything to do here? */
if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
......@@ -224,7 +223,7 @@ void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) {
struct part *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
/* Anything to do here? */
if (c->ti_end_min > ti_current) return;
......@@ -326,7 +325,7 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
struct part *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
const int count_j = cj->count;
struct part *restrict parts_j = cj->parts;
......@@ -433,7 +432,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
struct part *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
const int count_j = cj->count;
struct part *restrict parts_j = cj->parts;
......@@ -627,7 +626,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
struct part *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
const int count_i = ci->count;
struct part *restrict parts_j = ci->parts;
......@@ -719,7 +718,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
/* Anything to do here? */
if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
......@@ -912,7 +911,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *piq2[VEC_SIZE], *pjq2[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
/* Anything to do here? */
if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
......@@ -1298,7 +1297,8 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
float dxq2[3 * VEC_SIZE] __attribute__((aligned(16)));
struct part *piq2[VEC_SIZE], *pjq2[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
if (c->ti_end_min > ti_current) return;
if (c->ti_end_max < ti_current) error("Cell in an impossible time-zone");
......@@ -1526,7 +1526,8 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
float dxq2[3 * VEC_SIZE] __attribute__((aligned(16)));
struct part *piq2[VEC_SIZE], *pjq2[VEC_SIZE];
#endif
TIMER_TIC
TIMER_TIC;
if (c->ti_end_min > ti_current) return;
if (c->ti_end_max < ti_current) error("Cell in an impossible time-zone");
......@@ -1721,7 +1722,7 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
struct space *s = r->e->s;
const int ti_current = r->e->ti_current;
TIMER_TIC
TIMER_TIC;
/* Should we even bother? */
if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
......@@ -1962,7 +1963,7 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
const int ti_current = r->e->ti_current;
TIMER_TIC
TIMER_TIC;
/* Should we even bother? */
if (ci->ti_end_min > ti_current) return;
......@@ -2005,7 +2006,7 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
struct space *s = r->e->s;
const int ti_current = r->e->ti_current;
TIMER_TIC
TIMER_TIC;
/* Should we even bother? */
if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
......@@ -2246,7 +2247,7 @@ void DOSUB_SELF2(struct runner *r, struct cell *ci, int gettimer) {
const int ti_current = r->e->ti_current;
TIMER_TIC
TIMER_TIC;
/* Should we even bother? */
if (ci->ti_end_min > ti_current) return;
......@@ -2278,7 +2279,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
struct space *s = r->e->s;
const int ti_current = r->e->ti_current;
TIMER_TIC
TIMER_TIC;
/* Find out in which sub-cell of ci the parts are. */
struct cell *sub = NULL;
......
......@@ -240,6 +240,70 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
}
}
void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
float r2, hi, hj, hig2, hjg2, dx[3];
struct part *pi, *pj;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
for (int j = 0; j < cj->count; ++j) {
pj = &cj->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
/* Pairwise distance */
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = ci->parts[i].x[k] - cj->parts[j].x[k];
r2 += dx[k] * dx[k];
}
/* Hit or miss? */
if (r2 < hig2 || r2 < hjg2) {
/* Interact */
runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj);
}
}
}
/* Reverse double-for loop and checks every interaction */
for (int j = 0; j < cj->count; ++j) {
pj = &cj->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
for (int i = 0; i < ci->count; ++i) {
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
/* Pairwise distance */
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = cj->parts[j].x[k] - ci->parts[i].x[k];
r2 += dx[k] * dx[k];
}
/* Hit or miss? */
if (r2 < hjg2 || r2 < hig2) {
/* Interact */
runner_iact_nonsym_force(r2, dx, hj, pi->h, pj, pi);
}
}
}
}
void self_all_density(struct runner *r, struct cell *ci) {
float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
struct part *pi, *pj;
......@@ -287,6 +351,42 @@ void self_all_density(struct runner *r, struct cell *ci) {
}
}
void self_all_force(struct runner *r, struct cell *ci) {
float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
struct part *pi, *pj;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
for (int j = i + 1; j < ci->count; ++j) {
pj = &ci->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
if (pi == pj) continue;
/* Pairwise distance */
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dxi[k] = ci->parts[i].x[k] - ci->parts[j].x[k];
r2 += dxi[k] * dxi[k];
}
/* Hit or miss? */
if (r2 < hig2 || r2 < hjg2) {
/* Interact */
runner_iact_force(r2, dxi, hi, hj, pi, pj);
}
}
}
}
void pairs_single_grav(double *dim, long long int pid,
struct gpart *restrict gparts, const struct part *parts,
int N, int periodic) {
......
......@@ -40,6 +40,8 @@ void pairs_single_density(double *dim, long long int pid,
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
void self_all_density(struct runner *r, struct cell *ci);
void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj);
void self_all_force(struct runner *r, struct cell *ci);
void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
......
......@@ -22,12 +22,13 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# List of programs and scripts to run in the test suite
TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh \
testParser.sh testSPHStep
testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh \
testParser.sh testSPHStep test125cells.sh
# List of test programs to compile
check_PROGRAMS = testGreetings testMaths testReading testSingle testKernel testSymmetry \
testTimeIntegration testSPHStep testPair test27cells testParser testInteractions
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSPHStep testPair test27cells test125cells testParser \
testKernel testInteractions testMaths testSymmetry
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
......@@ -50,6 +51,8 @@ testPair_SOURCES = testPair.c
test27cells_SOURCES = test27cells.c
test125cells_SOURCES = test125cells.c
testParser_SOURCES = testParser.c
testInteractions_SOURCES = testInteractions.c
......@@ -57,4 +60,4 @@ testInteractions_SOURCES = testInteractions.c
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \
testParserInput.yaml
test125cells.sh testParserInput.yaml
......@@ -89,6 +89,7 @@ for i in range(n_lines_to_check):
abs_diff = abs(data1[i,j] - data2[i,j])
sum = abs(data1[i,j] + data2[i,j])
if abs(data1[i,j]) + abs(data2[i,j]) < 2.5e-7: continue
if sum > 0:
rel_diff = abs(data1[i,j] - data2[i,j]) / sum
else:
......
This diff is collapsed.
#!/bin/bash
for v in {0..3}
do
for p in {0..2}
do
echo ""
rm brute_force_125_standard.dat swift_dopair_125_standard.dat
./test125cells -n 6 -r 1 -v $v -p $p -f standard
python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance_125.dat 6
done
done
exit $?
......@@ -17,11 +17,17 @@
*
******************************************************************************/
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
/* Local headers. */
#include "swift.h"
enum velocity_types {
......@@ -243,7 +249,7 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
/* And go... */
int main(int argc, char *argv[]) {
size_t runs = 0, particles = 0;
double h = 1.2348, size = 1., rho = 1.;
double h = 1.23485, size = 1., rho = 1.;
double perturbation = 0.;
char outputFileNameExtension[200] = "";
char outputFileName[200] = "";
......@@ -253,11 +259,14 @@ int main(int argc, char *argv[]) {
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
/* Get some randomness going */
srand(0);
char c;
while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:v:")) != -1) {
while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
......@@ -265,7 +274,7 @@ int main(int argc, char *argv[]) {
case 's':
sscanf(optarg, "%lf", &size);
break;
case 'p':
case 'n':
sscanf(optarg, "%zu", &particles);
break;
case 'r':
......@@ -291,9 +300,10 @@ int main(int argc, char *argv[]) {
if (h < 0 || particles == 0 || runs == 0) {
printf(
"\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
"\nGenerates a cell pair, filled with particles on a Cartesian grid."
"\nThese are then interacted using runner_dopair1_density."
"\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
"\nGenerates 27 cells, filled with particles on a Cartesian grid."
"\nThese are then interacted using runner_dopair1_density() and "
"runner_doself1_density()."
"\n\nOptions:"
"\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
"\n-m rho - Physical density in the cell"
......@@ -307,13 +317,14 @@ int main(int argc, char *argv[]) {
}
/* Help users... */
message("Adiabatic index: ga = %f", hydro_gamma);
message("Smoothing length: h = %f", h * size);
message("Kernel: %s", kernel_name);
message("Neighbour target: N = %f",
h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
message("Neighbour target: N = %f", h * h * h * kernel_norm);
message("Density target: rho = %f", rho);
message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
printf("\n");
/* Build the infrastructure */
......
#!/bin/bash
rm brute_force_27_standard.dat swift_dopair_27_standard.dat
./test27cells -p 6 -r 1 -d 0 -f standard
python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance.dat 6
for v in {0..3}
do
echo ""
rm brute_force_27_standard.dat swift_dopair_27_standard.dat
./test27cells -n 6 -r 1 -d 0 -f standard -v $v
python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance_27.dat 6
done
exit $?
#!/bin/bash
rm brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
./test27cells -p 6 -r 1 -d 0.1 -f perturbed
python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance.dat 6
for v in {0..3}
do
echo ""
rm brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
./test27cells -n 6 -r 1 -d 0.1 -f perturbed -v $v
python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tol