diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index 4bf621dfbf870d712e3b994aa1cfc13644ff24a3..90247203793765933e9a90b9d01650ab1069c7ba 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -35,7 +35,7 @@ InitialConditions:
# Dimensionless pre-factor for the time-step condition
LambdaCooling:
- lambda_cgs: 0. #1.0e-22 # Cooling rate (in cgs units)
+ lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
diff --git a/examples/CoolingBox/getGlass.sh b/examples/CoolingBox/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46
--- /dev/null
+++ b/examples/CoolingBox/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index 5de012a17af4eef71e56548602e7956faef529f5..f863e174b1fcd404ae178fe324c7a165598b4af0 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -1,6 +1,6 @@
###############################################################################
# This file is part of SWIFT.
- # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durhama.ac.uk)
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
@@ -22,13 +22,11 @@ import h5py
import sys
from numpy import *
-# Generates a swift IC file containing a cartesian distribution of particles
-# at a constant density and pressure in a cubic box
+# Generates a SWIFT IC file with a constant density and pressure
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1 # 1 kiloparsec
-L = int(sys.argv[1]) # Number of particles along one axis
rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
P = 4.5e6 # Pressure in code units (at 10^5K)
gamma = 5./3. # Gas adiabatic index
@@ -36,12 +34,17 @@ eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
#---------------------------------------------------
-numPart = L**3
-mass = boxSize**3 * rho / numPart
-print mass
-internalEnergy = P / ((gamma - 1.)*rho)
-#--------------------------------------------------
+# Read id, position and h from glass
+glass = h5py.File("glassCube_32.hdf5", "r")
+ids = glass["/PartType0/ParticleIDs"][:]
+pos = glass["/PartType0/Coordinates"][:,:] * boxSize
+h = glass["/PartType0/SmoothingLength"][:] * boxSize
+
+# Compute basic properties
+numPart = size(pos) / 3
+mass = boxSize**3 * rho / numPart
+internalEnergy = P / ((gamma - 1.) * rho)
#File
file = h5py.File(fileName, 'w')
@@ -57,11 +60,11 @@ grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
-#Runtime parameters
+# Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
-#Units
+# Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 3.0857e21
grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33
@@ -75,35 +78,26 @@ grp = file.create_group("/PartType0")
v = zeros((numPart, 3))
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
-v = zeros(1)
m = full((numPart, 1), mass)
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
-m = zeros(1)
-h = full((numPart, 1), eta * boxSize / L)
-ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+h = reshape(h, (numPart, 1))
+ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
ds[()] = h
-h = zeros(1)
u = full((numPart, 1), internalEnergy)
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
-u = zeros(1)
-
-ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ids = reshape(ids, (numPart, 1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids + 1
-x = ids % L;
-y = ((ids - x) / L) % L;
-z = (ids - x - L * y) / L**2;
-coords = zeros((numPart, 3))
-coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
-coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
-coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds[()] = ids
+
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
-ds[()] = coords
+ds[()] = pos
file.close()
+
+print numPart
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 7a8e804d5a858c695b44319f161a3eb69a8945e8..fae7fd4ac5a760a3dc3b85dbddf541586ece79b8 100755
--- a/examples/CoolingBox/run.sh
+++ b/examples/CoolingBox/run.sh
@@ -1,14 +1,18 @@
#!/bin/bash
# Generate the initial conditions if they are not present.
-echo "Generating initial conditions for the cooling box example..."
-
-python makeIC.py 10
+if [ ! -e glassCube_32.hdf5 ]
+then
+ echo "Fetching initial glass file for the cooling box example..."
+ ./getGlass.sh
+fi
+if [ ! -e coolingBox.hdf5 ]
+then
+ echo "Generating initial conditions for the cooling box example..."
+ python makeIC.py
+fi
+# Run SWIFT
../swift -s -C -t 1 coolingBox.yml
-#-C 2>&1 | tee output.log
-
-python energy_plot.py 0
-
-#python test_energy_conservation.py 0
+# Check energy conservation and cooling rate
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 50a21a68b65539785369b29f35c0ae72664eecce..eced1910e80f58daeeeeda61f263e0ba09811614 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -97,7 +97,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
/* Store the radiated energy */
- xp->cooling_data.radiated_energy += hydro_get_mass(p) * cooling_du_dt * dt;
+ xp->cooling_data.radiated_energy += -hydro_get_mass(p) * cooling_du_dt * dt;
}
/**
@@ -113,7 +113,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
const struct phys_const* restrict phys_const,
const struct UnitSystem* restrict us, const struct part* restrict p) {
- /* Get current internal energy (dt=0) */
+ /* Get current internal energy */
const float u = hydro_get_internal_energy(p);
const float du_dt = cooling_rate(phys_const, us, cooling, p);
diff --git a/src/statistics.c b/src/statistics.c
index 297d88c1f25c1b5b42be8edcd1282fd437964894..b40f7afb6ed621f6a9ae3e3539b0c7c7ba2cbd17 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -163,7 +163,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
/* Collect energies. */
stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
stats.E_int += m * hydro_get_internal_energy(p);
- stats.E_rad += cooling_get_radiated_energy(xp);
+ stats.E_rad += cooling_get_radiated_energy(xp) / 2.;
stats.E_pot_self += 0.f;
if (gp != NULL)
stats.E_pot_ext += m * external_gravity_get_potential_energy(