diff --git a/.gitignore b/.gitignore
index 907d6b4685189f98f4d0895f23788cf6da382317..dba88cee0989e00baef51d9e75a42b1141d15714 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,6 +57,8 @@ tests/testTimeIntegration
tests/testSPHStep
tests/testKernel
tests/testInteractions
+tests/testSymmetry
+tests/testMaths
tests/testParser
tests/parser_output.yml
tests/test27cells.sh
diff --git a/examples/PerturbedBox/makeIC.py b/examples/PerturbedBox/makeIC.py
index ee1d845fc2149892909a54bf588046b0b1691b03..cc7fffe14d4f361153a07101ddcec20a3c979b4a 100644
--- a/examples/PerturbedBox/makeIC.py
+++ b/examples/PerturbedBox/makeIC.py
@@ -19,6 +19,7 @@
##############################################################################
import h5py
+import sys
import random
from numpy import *
@@ -26,13 +27,13 @@ from numpy import *
# at a constant density and pressure in a cubic box
# Parameters
-periodic= 1 # 1 For periodic box
+periodic= 1 # 1 For periodic box
boxSize = 1.
-L = 50 # Number of particles along one axis
-rho = 1. # Density
-P = 1. # Pressure
-gamma = 5./3. # Gas adiabatic index
-pert = 0.01 # Perturbation scale (in units of the interparticle separation)
+L = int(sys.argv[1]) # Number of particles along one axis
+rho = 1. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+pert = 0.01 # Perturbation scale (in units of the interparticle separation)
fileName = "perturbedBox.hdf5"
diff --git a/examples/PerturbedBox/perturbedBox.yml b/examples/PerturbedBox/perturbedBox.yml
new file mode 100644
index 0000000000000000000000000000000000000000..3a445e27e74fb83fb8deb56fde6003c22f7dedf1
--- /dev/null
+++ b/examples/PerturbedBox/perturbedBox.yml
@@ -0,0 +1,35 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1. # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: perturbedBox # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-1 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./perturbedBox.hdf5 # The file to read
diff --git a/examples/PerturbedBox/run.sh b/examples/PerturbedBox/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c15011ca526efb7250e93dc8cb002e49647b82b6
--- /dev/null
+++ b/examples/PerturbedBox/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e perturbedBox.hdf5 ]
+then
+ echo "Generating initial conditions for the perturbed box example..."
+ python makeIC.py 50
+fi
+
+../swift -s -t 16 perturbedBox.yml
diff --git a/src/approx_math.h b/src/approx_math.h
index cbca602b3fafcc5044b0939b2207b8f9d50a7446..ad07adeb4f3b1b54ca5f33d80eabb6a004d2a3aa 100644
--- a/src/approx_math.h
+++ b/src/approx_math.h
@@ -24,12 +24,16 @@
/**
* @brief Approximate version of expf(x) using a 4th order Taylor expansion
*
- * The absolute error is of order 10^-6 for -0.2 < x < 0.2.
+ * The absolute error is smaller than 3 * 10^-6 for -0.2 < x < 0.2.
+ * The absolute error is smaller than 2 * 10^-7 for -0.1 < x < 0.1.
+
+ * The relative error is smaller than 1 * 10^-6 for -0.2 < x < 0.2.
+ * The relative error is smaller than 4 * 10^-8 for -0.1 < x < 0.1.
*
* @param x The number to take the exponential of.
*/
__attribute__((always_inline)) INLINE static float approx_expf(float x) {
- return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.0f + 1.f / 24.0f * x)));
+ return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.f + 1.f / 24.f * x)));
}
#endif /* SWIFT_APPROX_MATH_H */
diff --git a/src/cell.h b/src/cell.h
index 9f35f4386cb161517b1b0c57468f1c0317c0c6a7..b9e060050827e6e0bbd62c87024bbf066d3807f8 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -152,7 +152,7 @@ struct cell {
double mom[3], ang_mom[3];
/* Mass, potential, internal and kinetic energy of particles in this cell. */
- double mass, e_pot, e_int, e_kin;
+ double mass, e_pot, e_int, e_kin, entropy;
/* Number of particles updated in this cell. */
int updated, g_updated;
diff --git a/src/engine.c b/src/engine.c
index 8e471466fad5bac7c17603a42a59fba50d8c024e..6609f7f0d53453c6649a8089d2b656f555e00f75 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2098,7 +2098,7 @@ void engine_collect_drift(struct cell *c) {
if (c->drift != NULL) return;
/* Counters for the different quantities. */
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
/* Only do something is the cell is non-empty */
@@ -2120,6 +2120,7 @@ void engine_collect_drift(struct cell *c) {
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
+ entropy += cp->entropy;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
@@ -2135,6 +2136,7 @@ void engine_collect_drift(struct cell *c) {
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
+ c->entropy = entropy;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
@@ -2151,7 +2153,7 @@ void engine_print_stats(struct engine *e) {
const struct space *s = e->s;
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
/* Collect the cell data. */
@@ -2167,6 +2169,7 @@ void engine_print_stats(struct engine *e) {
e_kin += c->e_kin;
e_int += c->e_int;
e_pot += c->e_pot;
+ entropy += c->entropy;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
mom[2] += c->mom[2];
@@ -2178,8 +2181,8 @@ void engine_print_stats(struct engine *e) {
/* Aggregate the data from the different nodes. */
#ifdef WITH_MPI
{
- double in[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
- double out[10];
+ double in[11] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
+ double out[11];
out[0] = e_kin;
out[1] = e_int;
out[2] = e_pot;
@@ -2190,7 +2193,8 @@ void engine_print_stats(struct engine *e) {
out[7] = ang_mom[1];
out[8] = ang_mom[2];
out[9] = mass;
- if (MPI_Allreduce(out, in, 10, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
+ out[10] = entropy;
+ if (MPI_Reduce(out, in, 11, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate stats.");
e_kin = out[0];
@@ -2203,6 +2207,7 @@ void engine_print_stats(struct engine *e) {
ang_mom[1] = out[7];
ang_mom[2] = out[8];
mass = out[9];
+ entropy = out[10];
}
#endif
@@ -2210,10 +2215,11 @@ void engine_print_stats(struct engine *e) {
/* Print info */
if (e->nodeID == 0) {
- fprintf(e->file_stats,
- " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
- e->time, mass, e_tot, e_kin, e_int, e_pot, mom[0], mom[1], mom[2],
- ang_mom[0], ang_mom[1], ang_mom[2]);
+ fprintf(
+ e->file_stats,
+ " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
+ e->time, mass, e_tot, e_kin, e_int, e_pot, entropy, mom[0], mom[1],
+ mom[2], ang_mom[0], ang_mom[1], ang_mom[2]);
fflush(e->file_stats);
}
}
@@ -2973,10 +2979,11 @@ void engine_init(struct engine *e, struct space *s,
engine_default_energy_file_name);
sprintf(energyfileName + strlen(energyfileName), ".txt");
e->file_stats = fopen(energyfileName, "w");
- fprintf(e->file_stats,
- "# %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
- "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "p_x", "p_y",
- "p_z", "ang_x", "ang_y", "ang_z");
+ fprintf(
+ e->file_stats,
+ "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
+ "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "Entropy", "p_x",
+ "p_y", "p_z", "ang_x", "ang_y", "ang_z");
fflush(e->file_stats);
char timestepsfileName[200] = "";
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 7af9a22c7d12dad36f40e5eea1f69715ffd35c35..323f066db6a047db7ae3401a3525afa941bfc047 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -470,7 +470,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Change in entropy */
pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
- pj->entropy_dt -= 0.5f * mi * visc_term * dvdr;
+ pj->entropy_dt += 0.5f * mi * visc_term * dvdr;
}
/**
diff --git a/src/runner.c b/src/runner.c
index e48d2ccaba69315f4295a0dbacd36e9c7d26a0c4..63243ab66d24a2efc045e19e866e24b7a0f4868c 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -593,7 +593,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gparts = c->gparts;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0};
double ang_mom[3] = {0.0, 0.0, 0.0};
@@ -670,6 +670,9 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
e_pot += 0.;
e_int += m * hydro_get_internal_energy(p, half_dt);
+
+ /* Collect entropy */
+ entropy += m * hydro_get_entropy(p, half_dt);
}
/* Now, get the maximal particle motion from its square */
@@ -694,6 +697,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
+ entropy += cp->entropy;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
@@ -710,6 +714,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
+ c->entropy = entropy;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 90220e972065a19878aa330af3d03c6988e4e2cf..4d6cb65ca6971d30fd6e39f52ce582c59cecb8b8 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -21,20 +21,26 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) -DTIMER
AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# List of programs and scripts to run in the test suite
-TESTS = testGreetings testReading.sh testSingle testPair.sh testPairPerturbed.sh \
- test27cells.sh test27cellsPerturbed.sh testParser.sh testKernel testSPHStep \
- test125cells.sh
+TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
+ testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh \
+ testParser.sh testSPHStep test125cells.sh
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSPHStep testPair test27cells test125cells testParser \
- testKernel testInteractions
+ testKernel testInteractions testMaths testSymmetry
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
+testMaths_SOURCES = testMaths.c
+
testReading_SOURCES = testReading.c
+testKernel_SOURCES = testKernel.c
+
+testSymmetry_SOURCES = testSymmetry.c
+
testTimeIntegration_SOURCES = testTimeIntegration.c
testSPHStep_SOURCES = testSPHStep.c
@@ -49,8 +55,6 @@ test125cells_SOURCES = test125cells.c
testParser_SOURCES = testParser.c
-testKernel_SOURCES = testKernel.c
-
testInteractions_SOURCES = testInteractions.c
# Files necessary for distribution
diff --git a/tests/testMaths.c b/tests/testMaths.c
new file mode 100644
index 0000000000000000000000000000000000000000..96c75313db18edc653aa3a48cb6ac34913297806
--- /dev/null
+++ b/tests/testMaths.c
@@ -0,0 +1,67 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 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 "../config.h"
+
+#include "approx_math.h"
+#include "vector.h"
+
+#include <math.h>
+#include <stdio.h>
+
+int main() {
+
+ const int numPoints = 60000;
+
+ for (int i = 0; i < numPoints; ++i) {
+
+ const float x = 0.6f * (i / (float)numPoints) - 0.3f;
+ const float exp_correct = expf(x);
+ const float exp_approx = approx_expf(x);
+
+ const float abs = fabs(exp_correct - exp_approx);
+ const float rel =
+ 0.5f * fabs(exp_correct - exp_approx) / fabs(exp_correct + exp_approx);
+
+ printf("%2d: x= %f exp(x)= %e approx_exp(x)=%e abs=%e rel=%e\n", i, x,
+ exp_correct, exp_approx, abs, rel);
+
+ if (abs > 3e-6 && fabsf(x) <= 0.2) {
+ printf("Absolute difference too large !\n");
+ return 1;
+ }
+ if (abs > 1.2e-7 && fabsf(x) <= 0.1) {
+ printf("Absolute difference too large !\n");
+ return 1;
+ }
+
+ if (rel > 1e-6 && fabsf(x) <= 0.2) {
+ printf("Relative difference too large !\n");
+ return 1;
+ }
+ if (rel > 4e-8 && fabsf(x) <= 0.1) {
+ printf("Relative difference too large !\n");
+ return 1;
+ }
+ }
+
+ printf("\nAll values are consistent\n");
+
+ return 0;
+}
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
new file mode 100644
index 0000000000000000000000000000000000000000..eb3fab6becca08e9ef87e7c60cc8c04bd2a0290c
--- /dev/null
+++ b/tests/testSymmetry.c
@@ -0,0 +1,111 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 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 "../config.h"
+
+#include <fenv.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "swift.h"
+
+int main(int argc, char *argv[]) {
+
+ /* Choke if need be */
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+
+ /* Create two random particles (don't do this at home !) */
+ struct part pi, pj;
+ for (size_t i = 0; i < sizeof(struct part) / sizeof(float); ++i) {
+ *(((float *)&pi) + i) = (float)random_uniform(0., 2.);
+ *(((float *)&pj) + i) = (float)random_uniform(0., 2.);
+ }
+
+ /* Make the particle smoothing length and position reasonable */
+ for (size_t i = 0; i < 3; ++i) pi.x[0] = random_uniform(-1., 1.);
+ for (size_t i = 0; i < 3; ++i) pj.x[0] = random_uniform(-1., 1.);
+ pi.h = 2.f;
+ pj.h = 2.f;
+ pi.id = 1;
+ pj.id = 2;
+
+ /* Make an xpart companion */
+ struct xpart xpi, xpj;
+ bzero(&xpi, sizeof(struct xpart));
+ bzero(&xpj, sizeof(struct xpart));
+
+ /* Make some copies */
+ struct part pi2, pj2;
+ memcpy(&pi2, &pi, sizeof(struct part));
+ memcpy(&pj2, &pj, sizeof(struct part));
+
+ int i_ok = memcmp(&pi, &pi2, sizeof(struct part));
+ int j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+
+ if (i_ok != 0) error("Particles 'pi' do not match after copy");
+ if (j_ok != 0) error("Particles 'pj' do not match after copy");
+
+ /* Compute distance vector */
+ float dx[3];
+ dx[0] = pi.x[0] - pj.x[0];
+ dx[1] = pi.x[1] - pj.x[1];
+ dx[2] = pi.x[2] - pj.x[2];
+ float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ /* --- Test the density loop --- */
+
+ /* Call the symmetric version */
+ runner_iact_density(r2, dx, pi.h, pj.h, &pi, &pj);
+
+ /* Call the non-symmetric version */
+ runner_iact_nonsym_density(r2, dx, pi2.h, pj2.h, &pi2, &pj2);
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
+ runner_iact_nonsym_density(r2, dx, pj2.h, pi2.h, &pj2, &pi2);
+
+ /* Check that the particles are the same */
+ i_ok = memcmp(&pi, &pi2, sizeof(struct part));
+ j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+
+ if (i_ok) error("Particles 'pi' do not match after density");
+ if (j_ok) error("Particles 'pj' do not match after density");
+
+ /* --- Test the force loop --- */
+
+ /* Call the symmetric version */
+ runner_iact_force(r2, dx, pi.h, pj.h, &pi, &pj);
+
+ /* Call the non-symmetric version */
+ runner_iact_nonsym_force(r2, dx, pi2.h, pj2.h, &pi2, &pj2);
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
+ runner_iact_nonsym_force(r2, dx, pj2.h, pi2.h, &pj2, &pi2);
+
+ /* Check that the particles are the same */
+ i_ok = memcmp(&pi, &pi2, sizeof(struct part));
+ j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+
+ if (i_ok) error("Particles 'pi' do not match after force");
+ if (j_ok) error("Particles 'pj' do not match after force");
+
+ return 0;
+}