Commit e29af86f authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'energy_conservation' into 'master'

Correct the equation for the entropy time derivative in GADGET2_SPH

This fixes #183 meaning that energy is now correctly conserved.

It also:
 - Adds a parameter file and run script for the perturbed box test case.
 - Adds entropy to the diagnostics file.
 - Replace the MPI_AllReduce in the statistics collection by an MPI_Reduce since only rank 0 writes anyway. 
 - Gives a better documentation for the `approx_exp()` function.
 - Adds an accuracy test for the approximate maths functions (currently only exp())
 - Adds a test no make sure the symmetric and non-symmetric versions of the SPH interaction routines give the same answer. That was the origin of the energy non-conservation bug.

See merge request !210
parents bd9ec753 c406bde1
......@@ -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
......
......@@ -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"
......
# 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
#!/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
......@@ -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 */
......@@ -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;
......
......@@ -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] = "";
......
......@@ -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;
}
/**
......
......@@ -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];
......
......@@ -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
......
/*******************************************************************************
* 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;
}
/*******************************************************************************
* 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;
}
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