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

Merge branch 'testPeriodicBC' into 'master'

Test periodic bc

This tests the periodic boundary conditions which are generated in `space_getsid()` and used in `runner_dopair1_density`.

An 8x8x8 mesh is generated that is periodic and 26 pair interactions are computed for the following cells:

* All 8 corner cells
* The cell at the centre of each face
* The cell at the centre of each edge

These checks are performed for both unperturbed and perturbed cells.

See merge request !350
parents bd229c7a 7dfb9465
......@@ -35,11 +35,16 @@ examples/*/*/used_parameters.yml
examples/*/gravity_checks_*.dat
tests/testPair
tests/brute_force_periodic_BC_standard.dat
tests/swift_periodic_BC_standard.dat
tests/brute_force_periodic_BC_pertrubed.dat
tests/swift_periodic_BC_perturbed.dat
tests/brute_force_standard.dat
tests/swift_dopair_standard.dat
tests/brute_force_perturbed.dat
tests/swift_dopair_perturbed.dat
tests/test27cells
tests/testPeriodicBC
tests/test125cells
tests/brute_force_27_standard.dat
tests/swift_dopair_27_standard.dat
......@@ -64,6 +69,8 @@ tests/testMaths
tests/testThreadpool
tests/testParser
tests/parser_output.yml
tests/testPeriodicBC.sh
tests/testPeriodicBCPerturbed.sh
tests/test27cells.sh
tests/test27cellsPerturbed.sh
tests/test125cells.sh
......
......@@ -854,13 +854,23 @@ 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/test125cellsPerturbed.sh], [chmod +x tests/test125cellsPerturbed.sh])
AC_CONFIG_FILES([tests/testPeriodicBC.sh], [chmod +x tests/testPeriodicBC.sh])
AC_CONFIG_FILES([tests/testPeriodicBCPerturbed.sh], [chmod +x tests/testPeriodicBCPerturbed.sh])
AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])
# Save the compilation options
AC_DEFINE_UNQUOTED([SWIFT_CONFIG_FLAGS],["$swift_config_flags"],[Flags passed to configure])
# Make sure the latest git revision string gets included
touch src/version.c
# Generate output.
AC_OUTPUT
# Report general configuration.
AC_MSG_RESULT([
AC_MSG_RESULT([
------- Summary --------
Compiler : $CC
- vendor : $ax_cv_c_compiler_vendor
- version : $ax_cv_c_compiler_version
......@@ -891,10 +901,5 @@ AC_MSG_RESULT([
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
])
# Make sure the latest git revision string gets included
touch src/version.c
# Generate output.
AC_OUTPUT
------------------------])
......@@ -1822,7 +1822,7 @@ void *runner_main(void *data) {
runner_dopair1_branch_density(r, ci, cj);
#ifdef EXTRA_HYDRO_LOOP
else if (t->subtype == task_subtype_gradient)
runner_dopair1_gradient(r, ci, cj);
runner_dopair1_branch_gradient(r, ci, cj);
#endif
else if (t->subtype == task_subtype_force)
runner_dopair2_force(r, ci, cj);
......
......@@ -37,6 +37,7 @@
#include "gravity.h"
#include "hydro.h"
#include "part.h"
#include "periodic.h"
#include "runner.h"
/**
......@@ -181,6 +182,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
float r2, hi, hj, hig2, hjg2, dx[3];
struct part *pi, *pj;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
......@@ -197,6 +199,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = ci->parts[i].x[k] - cj->parts[j].x[k];
dx[k] = nearest(dx[k], dim[k]);
r2 += dx[k] * dx[k];
}
......@@ -224,6 +227,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = cj->parts[j].x[k] - ci->parts[i].x[k];
dx[k] = nearest(dx[k], dim[k]);
r2 += dx[k] * dx[k];
}
......@@ -241,6 +245,7 @@ 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;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
......@@ -259,6 +264,7 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = ci->parts[i].x[k] - cj->parts[j].x[k];
dx[k] = nearest(dx[k], dim[k]);
r2 += dx[k] * dx[k];
}
......@@ -288,6 +294,7 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = cj->parts[j].x[k] - ci->parts[i].x[k];
dx[k] = nearest(dx[k], dim[k]);
r2 += dx[k] * dx[k];
}
......
......@@ -25,7 +25,8 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testParser.sh testSPHStep test125cells.sh test125cellsPerturbed.sh testFFT \
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
testMatrixInversion testThreadpool testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D
testVoronoi1D testVoronoi2D testVoronoi3D \
testPeriodicBC.sh testPeriodicBCPerturbed.sh
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
......@@ -34,7 +35,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSymmetry testThreadpool benchmarkInteractions \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
......@@ -55,6 +56,8 @@ testPair_SOURCES = testPair.c
test27cells_SOURCES = test27cells.c
testPeriodicBC_SOURCES = testPeriodicBC.c
test125cells_SOURCES = test125cells.c
testParser_SOURCES = testParser.c
......@@ -91,9 +94,9 @@ testLogger_SOURCES = testLogger.c
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh \
test125cells.sh test125cellsPerturbed.sh testParserInput.yaml difffloat.py \
tolerance_125_normal.dat tolerance_125_perturbed.dat \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
testPeriodicBCPerturbed.sh test125cells.sh test125cellsPerturbed.sh testParserInput.yaml \
difffloat.py tolerance_125_normal.dat tolerance_125_perturbed.dat \
tolerance_27_normal.dat tolerance_27_perturbed.dat tolerance_27_perturbed_h.dat \
tolerance_pair_normal.dat tolerance_pair_perturbed.dat \
fft_params.yml
fft_params.yml tolerance_periodic_BC_normal.dat tolerance_periodic_BC_perturbed.dat
......@@ -42,11 +42,6 @@ if len(sys.argv) >= 4:
if len(sys.argv) >= 5:
number_to_check = int(sys.argv[4])
if len(sys.argv) == 6:
ignoreSmallRhoDh = int(sys.argv[5])
else:
ignoreSmallRhoDh = 0
# Get the particle properties being compared from the header.
with open(file1, 'r') as f:
line = f.readline()
......@@ -69,7 +64,7 @@ n_lines = shape(data1)[0]
n_columns = shape(data1)[1]
if fileTol != "":
if n_linesTol != 2:
if n_linesTol != 3:
print "Incorrect number of lines in tolerance file '%s'."%fileTol
if n_columnsTol != n_columns:
print "Incorrect number of columns in tolerance file '%s'."%fileTol
......@@ -79,10 +74,12 @@ if fileTol == "":
print "Relative difference tolerance:", rel_tol
absTol = ones(n_columns) * abs_tol
relTol = ones(n_columns) * rel_tol
limTol = zeros(n_columns)
else:
print "Tolerances read from file"
absTol = dataTol[0,:]
relTol = dataTol[1,:]
limTol = dataTol[2,:]
n_lines_to_check = 0
if number_to_check > 0:
......@@ -113,11 +110,8 @@ for i in range(n_lines_to_check):
print ""
error = True
if abs(data1[i,j]) < 4e-6 and abs(data2[i,j]) < 4e-6 : continue
if abs(data1[i,j]) + abs(data2[i,j]) < limTol[j] : continue
# Ignore pathological cases with rho_dh
if ignoreSmallRhoDh and j == 8 and abs(data1[i,j]) < 2e-4: continue
if( rel_diff > 1.1*relTol[j]):
print "Relative difference larger than tolerance (%e) for particle %d, column %s:"%(relTol[j], data1[i,0], part_props[j])
print "%10s: a = %e"%("File 1", data1[i,j])
......
......@@ -623,7 +623,7 @@ int main(int argc, char *argv[]) {
/* Do the density calculation */
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Initialise the particle cache. */
/* Initialise the particle cache. */
#ifdef WITH_VECTORIZATION
runner.ci_cache.count = 0;
cache_init(&runner.ci_cache, 512);
......
......@@ -45,7 +45,6 @@
#define DOPAIR1_NAME "runner_dopair1_density"
#endif
/* n is both particles per axis and box size:
* particles are generated on a mesh with unit spacing
*/
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (C) 2015 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/>.
*
******************************************************************************/
/* 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"
#define ACC_THRESHOLD 1e-5
#if defined(WITH_VECTORIZATION)
#define DOSELF1 runner_doself1_density_vec
#define DOPAIR1 runner_dopair1_branch_density
#define DOSELF1_NAME "runner_doself1_density_vec"
#define DOPAIR1_NAME "runner_dopair1_density_vec"
#endif
#ifndef DOSELF1
#define DOSELF1 runner_doself1_density
#define DOSELF1_NAME "runner_doself1_density"
#endif
#ifndef DOPAIR1
#define DOPAIR1 runner_dopair1_branch_density
#define DOPAIR1_NAME "runner_dopair1_density"
#endif
enum velocity_types {
velocity_zero,
velocity_random,
velocity_divergent,
velocity_rotating
};
/**
* @brief Constructs a cell and all of its particle in a valid state prior to
* a DOPAIR or DOSELF calcuation.
*
* @param n The cube root of the number of particles.
* @param offset The position of the cell offset from (0,0,0).
* @param size The cell size.
* @param h The smoothing length of the particles in units of the inter-particle
*separation.
* @param density The density of the fluid.
* @param partId The running counter of IDs.
* @param pert The perturbation to apply to the particles in the cell in units
*of the inter-particle separation.
* @param vel The type of velocity field (0, random, divergent, rotating)
*/
struct cell *make_cell(size_t n, double *offset, double size, double h,
double density, long long *partId, double pert,
enum velocity_types vel) {
const size_t count = n * n * n;
const double volume = size * size * size;
struct cell *cell = malloc(sizeof(struct cell));
bzero(cell, sizeof(struct cell));
if (posix_memalign((void **)&cell->parts, part_align,
count * sizeof(struct part)) != 0) {
error("couldn't allocate particles, no. of particles: %d", (int)count);
}
bzero(cell->parts, count * sizeof(struct part));
float h_max = 0.f;
/* Construct the parts */
struct part *part = cell->parts;
for (size_t x = 0; x < n; ++x) {
for (size_t y = 0; y < n; ++y) {
for (size_t z = 0; z < n; ++z) {
part->x[0] =
offset[0] +
size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
part->x[1] =
offset[1] +
size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
part->x[2] =
offset[2] +
size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
switch (vel) {
case velocity_zero:
part->v[0] = 0.f;
part->v[1] = 0.f;
part->v[2] = 0.f;
break;
case velocity_random:
part->v[0] = random_uniform(-0.05, 0.05);
part->v[1] = random_uniform(-0.05, 0.05);
part->v[2] = random_uniform(-0.05, 0.05);
break;
case velocity_divergent:
part->v[0] = part->x[0] - 1.5 * size;
part->v[1] = part->x[1] - 1.5 * size;
part->v[2] = part->x[2] - 1.5 * size;
break;
case velocity_rotating:
part->v[0] = part->x[1];
part->v[1] = -part->x[0];
part->v[2] = 0.f;
break;
}
part->h = size * h / (float)n;
h_max = fmax(h_max, part->h);
part->id = ++(*partId);
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#ifdef SHADOWFAX_SPH
double anchor[3] = {0., 0., 0.};
double side[3] = {1., 1., 1.};
voronoi_cell_init(&part->cell, part->x, anchor, side);
#endif
#else
part->mass = density * volume / count;
#endif
#if defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
part->entropy_one_over_gamma = 1.f;
#endif
part->time_bin = 1;
#ifdef SWIFT_DEBUG_CHECKS
part->ti_drift = 8;
part->ti_kick = 8;
#endif
++part;
}
}
}
/* Cell properties */
cell->split = 0;
cell->h_max = h_max;
cell->count = count;
cell->dx_max_part = 0.;
cell->dx_max_sort = 0.;
cell->width[0] = size;
cell->width[1] = size;
cell->width[2] = size;
cell->loc[0] = offset[0];
cell->loc[1] = offset[1];
cell->loc[2] = offset[2];
cell->ti_old_part = 8;
cell->ti_end_min = 8;
cell->ti_end_max = 8;
cell->ti_sort = 8;
shuffle_particles(cell->parts, cell->count);
cell->sorted = 0;
cell->sort = NULL;
cell->sortsize = 0;
return cell;
}
void clean_up(struct cell *ci) {
free(ci->parts);
free(ci->sort);
free(ci);
}
/**
* @brief Initializes all particles field to be ready for a density calculation
*/
void zero_particle_fields(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
hydro_init_part(&c->parts[pid], NULL);
}
}
/**
* @brief Ends the loop by adding the appropriate coefficients
*/
void end_calculation(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
hydro_end_density(&c->parts[pid]);
}
}
/**
* @brief Dump all the particles to a file
*/
void dump_particle_fields(char *fileName, struct cell *main_cell, int i, int j,
int k) {
FILE *file = fopen(fileName, "a");
/* Write header */
fprintf(file,
"# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s "
"%13s %13s %13s\n",
"ID", "pos_x", "pos_y", "pos_z", "v_x", "v_y", "v_z", "rho", "rho_dh",
"wcount", "wcount_dh", "div_v", "curl_vx", "curl_vy", "curl_vz");
fprintf(file, "# Centre cell at (i,j,k)=(%d, %d, %d) ---------------------\n",
i, j, k);
/* Write main cell */
for (int pid = 0; pid < main_cell->count; pid++) {
fprintf(file,
"%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
"%13e %13e %13e\n",
main_cell->parts[pid].id, main_cell->parts[pid].x[0],
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2],
hydro_get_density(&main_cell->parts[pid]),
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pid].density.rho_dh,
#endif
main_cell->parts[pid].density.wcount,
main_cell->parts[pid].density.wcount_dh,
#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
main_cell->parts[pid].density.div_v,
main_cell->parts[pid].density.rot_v[0],
main_cell->parts[pid].density.rot_v[1],
main_cell->parts[pid].density.rot_v[2]
#else
0., 0., 0., 0.
#endif
);
}
fclose(file);
}
/**
* @brief Compares the vectorised result against
* the serial result of the interaction.
*
* @param serial_parts Particle array that has been interacted serially
* @param vec_parts Particle array to be interacted using vectors
* @param count No. of particles that have been interacted
* @param threshold Level of accuracy needed
*
* @return Non-zero value if difference found, 0 otherwise
*/
int check_results(struct part *serial_parts, struct part *vec_parts, int count,
double threshold) {
int result = 0;
for (int i = 0; i < count; i++)
result += compare_particles(serial_parts[i], vec_parts[i], threshold);
return result;
}
/* Just a forward declaration... */
void runner_doself1_density(struct runner *r, struct cell *ci);
void runner_doself1_density_vec(struct runner *r, struct cell *ci);
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
struct cell *cj);
void test_boundary_conditions(struct cell **cells, struct runner runner,
const int loc_i, const int loc_j, const int loc_k,
const int dim, char *swiftOutputFileName,
char *bruteForceOutputFileName) {
/* Store the main cell for future use */
struct cell *main_cell = cells[loc_i * (dim * dim) + loc_j * dim + loc_k];
/* Zero the fields */
for (int j = 0; j < 512; ++j) zero_particle_fields(cells[j]);
/* Run all the pairs */
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
#ifdef WITH_VECTORIZATION
runner.ci_cache.count = 0;
cache_init(&runner.ci_cache, 512);
runner.cj_cache.count = 0;
cache_init(&runner.cj_cache, 512);
#endif
/* Now loop over all the neighbours of this cell
* and perform the pair interactions. */
for (int ii = -1; ii < 2; ii++) {
int iii = loc_i + ii;
iii = (iii + dim) % dim;
for (int jj = -1; jj < 2; jj++) {
int jjj = loc_j + jj;
jjj = (jjj + dim) % dim;
for (int kk = -1; kk < 2; kk++) {
int kkk = loc_k + kk;
kkk = (kkk + dim) % dim;
/* Get the neighbouring cell */
struct cell *cj = cells[iii * (dim * dim) + jjj * dim + kkk];
if (cj != main_cell) DOPAIR1(&runner, main_cell, cj);
}
}
}
/* And now the self-interaction */
DOSELF1(&runner, main_cell);
#endif
/* Let's get physical ! */
end_calculation(main_cell);
/* Dump particles from the main cell. */
dump_particle_fields(swiftOutputFileName, main_cell, loc_i, loc_j, loc_k);
/* Now perform a brute-force version for accuracy tests */
/* Zero the fields */
for (int i = 0; i < 512; ++i) zero_particle_fields(cells[i]);
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Now loop over all the neighbours of this cell
* and perform the pair interactions. */
for (int ii = -1; ii < 2; ii++) {
int iii = loc_i + ii;
iii = (iii + dim) % dim;
for (int jj = -1; jj < 2; jj++) {
int jjj = loc_j + jj;
jjj = (jjj + dim) % dim;
for (int kk = -1; kk < 2; kk++) {
int kkk = loc_k + kk;
kkk = (kkk + dim) % dim;
/* Get the neighbouring cell */
struct cell *cj = cells[iii * (dim * dim) + jjj * dim + kkk];
if (cj != main_cell) pairs_all_density(&runner, main_cell, cj);
}
}
}
/* And now the self-interaction */
self_all_density(&runner, main_cell);
#endif
/* Let's get physical ! */
end_calculation(main_cell);
/* Dump */
dump_particle_fields(bruteForceOutputFileName, main_cell, loc_i, loc_j,
loc_k);
}
/* And go... */
int main(int argc, char *argv[]) {
engine_pin();
size_t runs = 0, particles = 0;
double h = 1.23485, size = 1., rho = 1.;