Create ICs for cooling box in a better way.

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)

examples/CoolingBox/getGlass.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5

examples/CoolingBox/makeIC.py

 ###############################################################################
  # 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

examples/CoolingBox/run.sh

 #!/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

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);
 

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(