Commit f2b5d359 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Merge branch 'master' into gizmo

parents f50fbe0f b6cdc5e5
......@@ -668,9 +668,6 @@ case "$with_hydro" in
minimal)
AC_DEFINE([MINIMAL_SPH], [1], [Minimal SPH])
;;
debug)
AC_DEFINE([DEBUG_INTERACTIONS_SPH], [1], [Debug SPH])
;;
hopkins)
AC_DEFINE([HOPKINS_PE_SPH], [1], [Pressure-Entropy SPH])
;;
......@@ -689,6 +686,29 @@ case "$with_hydro" in
;;
esac
# Check if debugging interactions is switched on.
AC_ARG_ENABLE([debug-interactions],
[AS_HELP_STRING([--enable-debug-interactions],
[Activate interaction debugging, logging a maximum of @<:@N@:>@ neighbours. Defaults to 256 if no value set.]
)],
[enable_debug_interactions="$enableval"],
[enable_debug_interactions="no"]
)
if test "$enable_debug_interactions" != "no"; then
if test "$with_hydro" = "gadget2"; then
AC_DEFINE([DEBUG_INTERACTIONS_SPH],1,[Enable interaction debugging])
if test "$enable_debug_interactions" == "yes"; then
AC_DEFINE([MAX_NUM_OF_NEIGHBOURS],256,[The maximum number of particle neighbours to be logged])
[enable_debug_interactions="yes (Logging up to 256 neighbours)"]
else
AC_DEFINE_UNQUOTED([MAX_NUM_OF_NEIGHBOURS], [$enableval] ,[The maximum number of particle neighbours to be logged])
[enable_debug_interactions="yes (Logging up to $enableval neighbours)"]
fi
else
[enable_debug_interactions="no (only available for gadget2 hydro scheme)"]
fi
fi
# SPH Kernel function
AC_ARG_WITH([kernel],
[AS_HELP_STRING([--with-kernel=<kernel>],
......@@ -945,11 +965,12 @@ AC_MSG_RESULT([
Multipole order : $with_multipole_order
No gravity below ID : $no_gravity_below_id
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
Threadpool debugging : $enable_threadpool_debugging
Debugging checks : $enable_debugging_checks
Naive interactions : $enable_naive_interactions
Gravity checks : $gravity_force_checks
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
Threadpool debugging : $enable_threadpool_debugging
Debugging checks : $enable_debugging_checks
Interaction debugging : $enable_debug_interactions
Naive interactions : $enable_naive_interactions
Gravity checks : $gravity_force_checks
------------------------])
#!/bin/bash -l
#
# Runs the SedovBlast_3D and EAGLE_12 examples using naive, serial and vectorised particle interactions. Then compares the output between versions using check_ngbs.py
echo
echo "# Running SedovBlast_3D and EAGLE_12 with naive interactions and neighbour logging, 16 thread"
echo
cd ../
./autogen.sh
# Naive interactions run
./configure --disable-mpi --enable-debug-interactions=1024 --disable-vec --enable-naive-interactions
make clean; make -j 6
cd examples/SedovBlast_3D/
./getGlass.sh
python makeIC.py
../swift -s -t 16 -n 5 sedov.yml -P SPH:h_tolerance:10
mv sedov_0000.hdf5 sedov_naive.hdf5
cd ../EAGLE_12/
# Link to ICs
ln -s /gpfs/data/Swift/web-storage/ICs/EAGLE_ICs_12.hdf5 EAGLE_ICs_12.hdf5
../swift -s -t 16 -n 5 eagle_12.yml -P SPH:h_tolerance:10
mv eagle_0000.hdf5 eagle_12_naive.hdf5
cd ../../
echo
echo "# Running SedovBlast_3D and EAGLE_12 with serial interactions and neighbour logging, 16 thread"
echo
# Serial interactions run
./configure --disable-mpi --enable-debug-interactions=1024 --disable-vec
make clean; make -j 6
cd examples/SedovBlast_3D/
../swift -s -t 16 -n 5 sedov.yml -P SPH:h_tolerance:10
mv sedov_0000.hdf5 sedov_serial.hdf5
cd ../EAGLE_12/
../swift -s -t 16 -n 5 eagle_12.yml -P SPH:h_tolerance:10
mv eagle_0000.hdf5 eagle_12_serial.hdf5
cd ../../
echo
echo "# Running SedovBlast_3D and EAGLE_12 with vectorised interactions and neighbour logging, 16 thread"
echo
# Vectorised interactions run
./configure --disable-mpi --enable-debug-interactions=1024
make clean; make -j 6
cd examples/SedovBlast_3D/
../swift -s -t 16 -n 5 sedov.yml -P SPH:h_tolerance:10
mv sedov_0000.hdf5 sedov_vec.hdf5
# Compare outputs
if python ../check_ngbs.py sedov_naive.hdf5 sedov_serial.hdf5
then
echo "SedovBlast_3D comparison between naive and serial passed"
else
echo "SedovBlast_3D comparison between naive and serial failed"
exit 1
fi
if python ../check_ngbs.py sedov_naive.hdf5 sedov_vec.hdf5
then
echo "SedovBlast_3D comparison between naive and vectorised passed"
else
echo "SedovBlast_3D comparison between naive and vectorised failed"
exit 1
fi
if python ../check_ngbs.py sedov_serial.hdf5 sedov_vec.hdf5
then
echo "SedovBlast_3D comparison between serial and vectorised passed"
else
echo "SedovBlast_3D comparison between serial and vectorised failed"
exit 1
fi
cd ../EAGLE_12/
../swift -s -t 16 -n 5 eagle_12.yml -P SPH:h_tolerance:10
mv eagle_0000.hdf5 eagle_12_vec.hdf5
# Compare outputs
if python ../check_ngbs.py eagle_12_naive.hdf5 eagle_12_serial.hdf5
then
echo "EAGLE_12 comparison between naive and serial passed"
else
echo "EAGLE_12 comparison between naive and serial failed"
exit 1
fi
if python ../check_ngbs.py eagle_12_naive.hdf5 eagle_12_vec.hdf5
then
echo "EAGLE_12 comparison between naive and vectorised passed"
else
echo "EAGLE_12 comparison between naive and vectorised failed"
exit 1
fi
if python ../check_ngbs.py eagle_12_serial.hdf5 eagle_12_vec.hdf5
then
echo "EAGLE_12 comparison between serial and vectorised passed"
exit 0
else
echo "EAGLE_12 comparison between serial and vectorised failed"
exit 1
fi
import h5py as h
import numpy as np
import matplotlib
matplotlib.use("Agg")
from pylab import *
import os.path
kernel_gamma = 1.825742
kernel_gamma2 = kernel_gamma * kernel_gamma
kernel_gamma_dim = np.power(kernel_gamma,3)
hydro_dimension_unit_sphere = 4. * np.pi / 3.
kernel_norm = hydro_dimension_unit_sphere * kernel_gamma_dim
error = False
inputFile1 = ""
inputFile2 = ""
# Check list of density neighbours and check that they are correct.
def check_density_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
error_val = False
for k in range(0,num_invalid):
# Filter neighbour lists for valid particle ids
filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
# Check neighbour lists for differences
id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
pid = pids[mask][k]
# Loop over discrepancies and check if they are actually neighbours
for pjd in id_list:
pi_pos = pos[np.where(pids == pid)]
pj_pos = pos[np.where(pids == pjd)]
hi = h[np.where(pids == pid)]
dx = pi_pos[0][0] - pj_pos[0][0]
dy = pi_pos[0][1] - pj_pos[0][1]
dz = pi_pos[0][2] - pj_pos[0][2]
r2 = dx*dx + dy*dy + dz*dz
hig2 = hi*hi*kernel_gamma2
diff = abs(r2 - hig2)
print "Particle {} is missing {}, hig2: {}, r2: {}, |r2 - hig2|: {}".format(pid,pjd,hig2, r2, diff)
if diff < acc * hig2:
print "Missing interaction due to precision issue will be ignored."
else:
error_val = True
return error_val
# Check list of force neighbours and check that they are correct.
def check_force_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
error_val = False
for k in range(0,num_invalid):
# Filter neighbour lists for valid particle ids
filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
# Check neighbour lists for differences
id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
pid = pids[mask][k]
# Loop over discrepancies and check if they are actually neighbours
for pjd in id_list:
pi_pos = pos[np.where(pids == pid)]
pj_pos = pos[np.where(pids == pjd)]
hi = h[np.where(pids == pid)]
hj = h[np.where(pids == pjd)]
dx = pi_pos[0][0] - pj_pos[0][0]
dy = pi_pos[0][1] - pj_pos[0][1]
dz = pi_pos[0][2] - pj_pos[0][2]
r2 = dx*dx + dy*dy + dz*dz
hig2 = hi*hi*kernel_gamma2
hjg2 = hj*hj*kernel_gamma2
diff = abs(r2 - max(hig2, hjg2))
print "Particle {} is missing {}, hig2: {}, hjg2: {}, r2: {}, |r2 - max(hig2,hjg2)|: {}".format(pid,pjd,hig2, hjg2, r2, diff)
if diff < acc * max(hig2,hjg2):
print "Missing interaction due to precision issue will be ignored."
else:
error_val = True
return error_val
# Parse command line arguments
if len(sys.argv) < 3:
print "Error: pass input files as arguments"
sys.exit()
else:
inputFile1 = sys.argv[1]
inputFile2 = sys.argv[2]
if os.path.exists(inputFile1) != 1:
print "\n{} does not exist!\n".format(inputFile1)
sys.exit()
if os.path.exists(inputFile2) != 1:
print "\n{} does not exist!\n".format(inputFile2)
sys.exit()
# Open input files
file_naive = h.File(inputFile1, "r")
file_sort = h.File(inputFile2, "r")
# Read input file fields
ids_naive = file_naive["/PartType0/ParticleIDs"][:]
ids_sort = file_sort["/PartType0/ParticleIDs"][:]
h_naive = file_naive["/PartType0/SmoothingLength"][:]
h_sort = file_sort["/PartType0/SmoothingLength"][:]
pos_naive = file_naive["/PartType0/Coordinates"][:,:]
pos_sort = file_sort["/PartType0/Coordinates"][:,:]
num_density_naive = file_naive["/PartType0/Num_ngb_density"][:]
num_density_sort = file_sort["/PartType0/Num_ngb_density"][:]
num_force_naive = file_naive["/PartType0/Num_ngb_force"][:]
num_force_sort = file_sort["/PartType0/Num_ngb_force"][:]
neighbour_ids_density_naive = file_naive["/PartType0/Ids_ngb_density"][:]
neighbour_ids_density_sort = file_sort["/PartType0/Ids_ngb_density"][:]
neighbour_ids_force_naive = file_naive["/PartType0/Ids_ngb_force"][:]
neighbour_ids_force_sort = file_sort["/PartType0/Ids_ngb_force"][:]
#wcount_naive = file_naive["/PartType0/Wcount"][:]
#wcount_sort = file_sort["/PartType0/Wcount"][:]
#
#wcount_naive = wcount_naive * np.power(h_naive,3) * kernel_norm
#wcount_sort = wcount_sort * np.power(h_sort,3) * kernel_norm
# Cross check
max_density_ngbs_naive = np.max(num_density_naive)
max_density_ngbs_sort = np.max(num_density_sort)
max_force_ngbs_naive = np.max(num_force_naive)
max_force_ngbs_sort = np.max(num_force_sort)
print " Min Mean Max "
print " ---------------------"
print "Ngbs density naiv: ", np.min(num_density_naive), np.mean(num_density_naive), max_density_ngbs_naive
print "Ngbs density sort: ", np.min(num_density_sort), np.mean(num_density_sort), max_density_ngbs_sort
print "Ngbs force naiv: ", np.min(num_force_naive), np.mean(num_force_naive), max_force_ngbs_naive
print "Ngbs force sort: ", np.min(num_force_sort), np.mean(num_force_sort), max_force_ngbs_sort
#print "Wcount naiv: ", np.min(wcount_naive), np.mean(wcount_naive), np.max(wcount_naive)
#print "Wcount sort: ", np.min(wcount_sort), np.mean(wcount_sort), np.max(wcount_sort)
# Sort
index_naive = np.argsort(ids_naive)
index_sort = np.argsort(ids_sort)
num_density_naive = num_density_naive[index_naive]
num_density_sort = num_density_sort[index_sort]
num_force_naive = num_force_naive[index_naive]
num_force_sort = num_force_sort[index_sort]
ids_naive = ids_naive[index_naive]
ids_sort = ids_sort[index_sort]
neighbour_ids_density_naive = neighbour_ids_density_naive[index_naive]
neighbour_ids_density_sort = neighbour_ids_density_sort[index_sort]
neighbour_ids_force_naive = neighbour_ids_force_naive[index_naive]
neighbour_ids_force_sort = neighbour_ids_force_sort[index_sort]
#wcount_naive = wcount_naive[index_naive]
#wcount_sort = wcount_sort[index_sort]
h_naive = h_naive[index_naive]
h_sort = h_sort[index_sort]
pos_naive = pos_naive[index_naive]
pos_sort = pos_sort[index_sort]
neighbour_length_naive = len(neighbour_ids_density_naive[0])
neighbour_length_sort = len(neighbour_ids_density_sort[0])
# Check that input files are logging the same number of neighbours
if neighbour_length_naive != neighbour_length_sort:
print "Input files have logged different numbers of neighbour lengths!"
print "{} has logged: {} neighbours".format(inputFile1, neighbour_length_naive)
print "{} has logged: {} neighbours".format(inputFile2, neighbour_length_sort)
exit(1)
if (max_density_ngbs_naive > neighbour_length_naive or max_force_ngbs_naive > neighbour_length_naive or
max_density_ngbs_sort > neighbour_length_sort or max_force_ngbs_sort > neighbour_length_sort):
print "The number of neighbours has exceeded the number of neighbours logged."
print "Modify NUM_OF_NEIGHBOURS in hydro_part.h to log more neighbours."
print "The highest neighbour count is: ", max(max_density_ngbs_naive,max_force_ngbs_naive, max_density_ngbs_sort,max_force_ngbs_sort)
exit(1)
# First check
print "\n Min Max"
print " ----------"
print "Differences for density: ", min(num_density_naive - num_density_sort), max(num_density_naive - num_density_sort)
print "Differences for force: ", min(num_force_naive - num_force_sort), max(num_force_naive - num_force_sort)
# Get the IDs that are different
mask_density = num_density_naive != num_density_sort
mask_force = num_force_naive != num_force_sort
num_invalid_density = np.sum(mask_density)
num_invalid_force = np.sum(mask_force)
print "\nNum non-zero density: ", num_invalid_density
print "Num non-zero force: ", num_invalid_force
print "\nParticle IDs with incorrect densities"
print "----------------------------------------"
print ids_naive[mask_density]
# Check density neighbour lists
error += check_density_neighbours(ids_naive, neighbour_ids_density_naive,
neighbour_ids_density_sort, mask_density, pos_naive, h_naive,
num_invalid_density, 2e-6)
print "Num of density interactions", inputFile1
print num_density_naive[mask_density]
print "Num of density interactions", inputFile2
print num_density_sort[mask_density]
print "\nParticle IDs with incorrect forces"
print "------------------------------------"
print ids_naive[mask_force]
# Check force neighbour lists
error += check_force_neighbours(ids_naive, neighbour_ids_force_naive,
neighbour_ids_force_sort, mask_force, pos_naive, h_naive,
num_invalid_force, 2e-6)
print "Num of force interactions", inputFile1
print num_force_naive[mask_force]
#print "Smoothing lengths", inputFile1
#print h_naive[mask_force]
print "Num of force interactions", inputFile2
print num_force_sort[mask_force]
#print "Smoothing lengths", inputFile2
#print h_sort[mask_force]
# Statistics of h difference
h_relative = (h_naive - h_sort) / h_naive
print "h statistics: {} {} (Min, 1st Percentile)".format(np.min(h_relative), np.percentile(h_relative,1))
print "h statistics: {} {} (Mean, Median)".format(np.mean(h_relative), np.median(h_relative))
print "h statistics: {} {} (Max, 99th Percentile)".format(np.max(h_relative), np.percentile(h_relative, 99))
if error:
print "\n------------------"
print "Differences found."
print "------------------"
exit(1)
else:
print "\n---------------------"
print "No differences found."
print "---------------------"
exit(0)
......@@ -40,9 +40,7 @@
#include "space.h"
/* Import the right hydro definition */
#if defined(DEBUG_INTERACTIONS_SPH)
#include "./hydro/DebugInteractions/hydro_debug.h"
#elif defined(MINIMAL_SPH)
#if defined(MINIMAL_SPH)
#include "./hydro/Minimal/hydro_debug.h"
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_debug.h"
......
......@@ -27,11 +27,7 @@
#include "kernel_hydro.h"
/* Import the right functions */
#if defined(DEBUG_INTERACTIONS_SPH)
#include "./hydro/DebugInteractions/hydro.h"
#include "./hydro/DebugInteractions/hydro_iact.h"
#define SPH_IMPLEMENTATION "Debug SELF/PAIR"
#elif defined(MINIMAL_SPH)
#if defined(MINIMAL_SPH)
#include "./hydro/Minimal/hydro.h"
#include "./hydro/Minimal/hydro_iact.h"
#define SPH_IMPLEMENTATION "Minimal version of SPH (e.g. Price 2010)"
......
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2017 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/>.
*
******************************************************************************/
#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_H
#define SWIFT_DEBUG_INTERACTIONS_HYDRO_H
/**
* @file DebugInteractions/hydro.h
* @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
*/
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "part.h"
#include <float.h>
/**
* @brief Returns the internal energy of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the density of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_density(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the mass of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_mass(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the time derivative of internal energy of a particle
*
* We assume a constant density.
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
const struct part *restrict p) {
return 0.f;
}
/**
* @brief Returns the time derivative of internal energy of a particle
*
* We assume a constant density.
*
* @param p The particle of interest.
* @param du_dt The new time derivative of the internal energy.
*/
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
struct part *restrict p, float du_dt) {}
/**
* @brief Computes the hydro time-step of a given particle
*
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const struct part *restrict p, const struct xpart *restrict xp,
const struct hydro_props *restrict hydro_properties) {
return FLT_MAX;
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* @param p The particle to act upon.
* @param dt Physical time step of the particle during the next step.