Merge branch 'master' into drift-fixes

AUTHORS

@@ -8,3 +8,5 @@ John A. Regan 		john.a.regan@durham.ac.uk
 Angus Lepper		angus.lepper@ed.ac.uk
 Tom Theuns 		tom.theuns@durham.ac.uk
 Richard G. Bower	r.g.bower@durham.ac.uk
+Stefan Arridge		stefan.arridge@durham.ac.uk
+Massimiliano Culpo	massimiliano.culpo@googlemail.com

README

@@ -19,6 +19,7 @@ Usage: swift [OPTION]... PARAMFILE
 Valid options are:
   -a          Pin runners using processor affinity
   -c          Run with cosmological time integration
+  -C          Run with cooling
   -d          Dry run. Read the parameter file, allocate memory but does not read 
               the particles from ICs and exit before the start of time integration.
               Allows user to check validy of parameter and IC files as well as memory limits.
@@ -30,8 +31,9 @@ Valid options are:
   -n    {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. 
   -s          Run with SPH
   -t    {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-  -v     [12] Increase the level of verbosity 1: MPI-rank 0 writes 
-              2: All MPI-ranks write
+  -v     [12] Increase the level of verbosity
+  	      1: MPI-rank 0 writes
+	      2: All MPI-ranks write
   -y    {int} Time-step frequency at which task graphs are dumped
   -h          Print this help message and exit
 

doc/Doxyfile.in

@@ -760,8 +760,10 @@ WARN_LOGFILE           =
 # Note: If this tag is empty the current directory is searched.
 
 INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
-INPUT		       += @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
-INPUT		       += @top_srcdir@/src/riemann 
+INPUT		       += @top_srcdir@/src/hydro/Minimal
+INPUT		       += @top_srcdir@/src/gravity/Default
+INPUT		       += @top_srcdir@/src/riemann
+INPUT		       += @top_srcdir@/src/cooling/const_lambda
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

examples/CoolingBox/coolingBox.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     2.0e33   # Solar masses
+  UnitLength_in_cgs:   3.01e21   # Kilparsecs
+  UnitVelocity_in_cgs: 1.0e5   # Time unit is cooling time
+  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.0    # 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-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            coolingBox # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          1.0e-1       # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # 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:  ./coolingBox.hdf5     # The file to read
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+  lambda:                      0.0    # 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)
+  cooling_tstep_mult:          1.0    # Dimensionless pre-factor for the time-step condition
+  

examples/CoolingBox/energy_plot.py

+import numpy as np
+import matplotlib.pyplot as plt
+import h5py as h5
+import sys
+
+stats_filename = "./energy.txt"
+snap_filename = "coolingBox_000.hdf5"
+#plot_dir = "./"
+
+#some constants in cgs units
+k_b = 1.38E-16 #boltzmann
+m_p = 1.67e-24 #proton mass
+#initial conditions set in makeIC.py
+rho = 3.2e3
+P = 4.5e6
+n_H_cgs = 0.0001
+gamma = 5./3.
+T_init = 1.0e5
+
+#Read the units parameters from the snapshot
+f = h5.File(snap_filename,'r')
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+parameters = f["Parameters"]
+cooling_lambda = float(parameters.attrs["LambdaCooling:lambda"])
+min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"])
+mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"])
+X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
+
+#get number of particles
+header = f["Header"]
+n_particles = header.attrs["NumPart_ThisFile"][0]
+#read energy and time arrays
+array = np.genfromtxt(stats_filename,skip_header = 1)
+time = array[:,0]
+total_energy = array[:,2]
+total_mass = array[:,1]
+
+time = time[1:]
+total_energy = total_energy[1:]
+total_mass = total_mass[1:]
+
+#conversions to cgs
+rho_cgs = rho * unit_mass / (unit_length)**3
+time_cgs = time * unit_time
+u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2 
+
+#find the energy floor
+print min_T
+u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
+#find analytic solution
+analytic_time = np.linspace(time_cgs[0],time_cgs[-1],1000)
+print time_cgs[1]
+print analytic_time[1]
+du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
+u_analytic = du_dt_cgs*(analytic_time - analytic_time[0]) + u_init_cgs
+cooling_time = u_init_cgs/(-du_dt_cgs)
+#rescale energy to initial energy
+total_energy /= total_energy[0]
+u_analytic /= u_init_cgs
+u_floor_cgs /= u_init_cgs
+# plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
+# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
+#analytic_solution = np.zeros(n_snaps-1)
+for i in range(u_analytic.size):
+    if u_analytic[i]<u_floor_cgs:
+        u_analytic[i] = u_floor_cgs
+plt.plot(time_cgs,total_energy,'k',label = "Numerical solution")
+plt.plot(analytic_time,u_analytic,'--r',lw = 2.0,label = "Analytic Solution")
+plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
+plt.plot((time_cgs[0],time_cgs[0]),(0,1),'m',label = "First output")
+plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
+plt.xlabel("Time (seconds)")
+plt.ylabel("Energy/Initial energy")
+plt.ylim(0.999,1.001)
+plt.xlim(0,min(10.0*cooling_time,time_cgs[-1]))
+plt.legend(loc = "upper right")    
+if (int(sys.argv[1])==0):
+    plt.show()
+else:
+    plt.savefig(full_plot_filename,format = "png")
+    plt.close()

examples/CoolingBox/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ #                    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/>.
+ # 
+ ##############################################################################
+
+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
+
+# 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 (0.01 hydrogen atoms per cm^3)
+P = 4.5e6          # Pressure in code units (at 10^5K)
+gamma = 5./3.         # Gas adiabatic index
+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)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+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
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 3.08e21 
+grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33 
+grp.attrs["Unit time in cgs (U_t)"] = 3.08e16 
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+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')
+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))
+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 = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+
+file.close()

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
+
+../swift -s -t 1 coolingBox.yml -C 2>&1 | tee output.log
+
+python energy_plot.py 0

examples/CosmoVolume/cosmoVolume.yml

@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
-# Parameters for the task scheduling
-Scheduler:
-  cell_sub_size:    6000     # Value used for the original scaling tests
-  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).

examples/EAGLE_12/eagle_12.yml

@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters for the task scheduling
-#Scheduler:
-  cell_sub_size:    64000000     # Value used for the original scaling tests
-#  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).

examples/EAGLE_25/eagle_25.yml

@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters for the task scheduling
-Scheduler:
-  cell_sub_size:    6000     # Value used for the original scaling tests
-  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).

examples/EAGLE_50/eagle_50.yml

@@ -6,11 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters for the task scheduling
-Scheduler:
-  cell_sub_size:    6000     # Value used for the original scaling tests
-  cell_split_size:  300      # Value used for the original scaling tests
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).

examples/Makefile.am

@@ -68,6 +68,7 @@ swift_fixdt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(MPI_LIBS) $(EXTRA_LIBS
 EXTRA_DIST = BigCosmoVolume/makeIC.py \
 	     BigPerturbedBox/makeIC_fcc.py \
 	     CosmoVolume/cosmoVolume.yml CosmoVolume/getIC.sh CosmoVolume/run.sh \
+	     CoolingBox/coolingBox.yml CoolingBox/energy_plot.py CoolingBox/makeIC.py CoolingBox/run.sh \
 	     EAGLE_12/eagle_12.yml EAGLE_12/getIC.sh EAGLE_12/README EAGLE_12/run.sh \
 	     EAGLE_25/eagle_25.yml EAGLE_25/getIC.sh EAGLE_25/README EAGLE_25/run.sh \
 	     EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \

examples/main.c

@@ -58,6 +58,7 @@ void print_help_message() {
   printf("Valid options are:\n");
   printf("  %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
   printf("  %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
+  printf("  %2s %8s %s\n", "-C", "", "Run with cooling");
   printf(
       "  %2s %8s %s\n", "-d", "",
       "Dry run. Read the parameter file, allocate memory but does not read ");
@@ -83,8 +84,8 @@ void print_help_message() {
   printf("  %2s %8s %s\n", "-t", "{int}",
          "The number of threads to use on each MPI rank. Defaults to 1 if not "
          "specified.");
-  printf("  %2s %8s %s\n", "-v", "[12]",
-         "Increase the level of verbosity 1: MPI-rank 0 writes ");
+  printf("  %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity");
+  printf("  %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
   printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
   printf("  %2s %8s %s\n", "-y", "{int}",
          "Time-step frequency at which task graphs are dumped");
@@ -146,6 +147,7 @@ int main(int argc, char *argv[]) {
   int nsteps = -2;
   int with_cosmology = 0;
   int with_external_gravity = 0;
+  int with_cooling = 0;
   int with_self_gravity = 0;
   int with_hydro = 0;
   int with_fp_exceptions = 0;
@@ -157,13 +159,16 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acdDef:gGhn:st:v:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acCdDef:gGhn:st:v:y:")) != -1) switch (c) {
       case 'a':
         with_aff = 1;
         break;
       case 'c':
         with_cosmology = 1;
         break;
+      case 'C':
+        with_cooling = 1;
+        break;
       case 'd':
         dry_run = 1;
         break;
@@ -325,12 +330,6 @@ int main(int argc, char *argv[]) {
   struct hydro_props hydro_properties;
   hydro_props_init(&hydro_properties, params);
 
-  /* Initialise the external potential properties */
-  struct external_potential potential;
-  if (with_external_gravity)
-    potential_init(params, &prog_const, &us, &potential);
-  if (with_external_gravity && myrank == 0) potential_print(&potential);
-
   /* Read particles and space information from (GADGET) ICs */
   char ICfileName[200] = "";
   parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -433,6 +432,17 @@ int main(int argc, char *argv[]) {
     message("nr of cells at depth %i is %i.", data[0], data[1]);
   }
 
+  /* Initialise the external potential properties */
+  struct external_potential potential;
+  if (with_external_gravity)
+    potential_init(params, &prog_const, &us, &potential);
+  if (with_external_gravity && myrank == 0) potential_print(&potential);
+
+  /* Initialise the cooling function properties */
+  struct cooling_data cooling;
+  if (with_cooling) cooling_init(params, &us, &prog_const, &cooling);
+  if (with_cooling && myrank == 0) cooling_print(&cooling);
+
   /* Construct the engine policy */
   int engine_policies = ENGINE_POLICY | engine_policy_steal;
   if (with_drift_all) engine_policies |= engine_policy_drift_all;
@@ -440,13 +450,14 @@ int main(int argc, char *argv[]) {
   if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
   if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
   if (with_cosmology) engine_policies |= engine_policy_cosmology;
+  if (with_cooling) engine_policies |= engine_policy_cooling;
 
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
   struct engine e;
   engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
               engine_policies, talking, &us, &prog_const, &hydro_properties,
-              &potential);
+              &potential, &cooling);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -519,7 +530,7 @@ int main(int argc, char *argv[]) {
 
       /* Make sure output file is empty, only on one rank. */
       char dumpfile[30];
-      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j);
+      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j + 1);
       FILE *file_thread;
       if (myrank == 0) {
         file_thread = fopen(dumpfile, "w");
@@ -571,7 +582,7 @@ int main(int argc, char *argv[]) {
 
 #else
       char dumpfile[30];
-      snprintf(dumpfile, 30, "thread_info-step%d.dat", j);
+      snprintf(dumpfile, 30, "thread_info-step%d.dat", j + 1);
       FILE *file_thread;
       file_thread = fopen(dumpfile, "w");
       /* Add some information to help with the plots */

examples/parameter_example.yml

@@ -9,8 +9,8 @@ InternalUnitSystem:
 # Parameters for the task scheduling
 Scheduler:
   nr_queues:        0        # (Optional) The number of task queues to use. Use 0  to let the system decide.
-  cell_max_size:    8000000  # (Optional) Maximal number of interactions per task (this is the default value).
-  cell_sub_size:    8000000  # (Optional) Maximal number of interactions per sub-task  (this is the default value).
+  cell_max_size:    8000000  # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
+  cell_sub_size:    64000000 # (Optional) Maximal number of interactions per sub-task  (this is the default value).
   cell_split_size:  400      # (Optional) Maximal number of particles per cell (this is the default value).
 
 # Parameters governing the time integration
@@ -63,7 +63,7 @@ DomainDecomposition:
   initial_grid_z:    10
   repartition_type:   b     # (Optional) The re-decomposition strategy ("n", "b", "v", "e" or "x").
  
-# Parameters related to external potentials
+# Parameters related to external potentials --------------------------------------------
   
 # Point mass external potentials
 PointMass:
@@ -73,6 +73,7 @@ PointMass:
   mass:            1e10     # mass of external point mass (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
 
+# Isothermal potential parameters
 IsothermalPotential:
   position_x:      100.     # Location of centre of isothermal potential (internal units)
   position_y:      100.
@@ -80,9 +81,27 @@ IsothermalPotential:
   vrot:            200.     # Rotation speed of isothermal potential (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
 
+# Disk-patch potential parameters
 Disk-PatchPotential:
   surface_density: 10.      # Surface density of the disk (internal units)
   scale_height:    100.     # Scale height of the disk (internal units)
   z_disk:          200.     # Disk height (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
   growth_time:     5.       # (Optional) Time for the disk to grow to its final size (multiple of the dynamical time)
+
+# Parameters related to cooling function  ----------------------------------------------
+
+# Constant du/dt cooling function
+ConstCooling:
+  cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
+  min_energy:   1.          # Minimal internal energy per unit mass (internal units)
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition
+
+# Constant lambda cooling function
+LambdaCooling:
+  lambda:                      2.0   # 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)
+  cooling_tstep_mult:          1.0   # Dimensionless pre-factor for the time-step condition
+

examples/plot_scaling_results.py

@@ -48,7 +48,7 @@ threadList = []
 hexcols = ['#332288', '#88CCEE', '#44AA99', '#117733', '#999933', '#DDCC77',
            '#CC6677', '#882255', '#AA4499', '#661100', '#6699CC', '#AA4466',
            '#4477AA']
-linestyle = (hexcols[0],hexcols[1],hexcols[3],hexcols[5],hexcols[6])
+linestyle = (hexcols[0],hexcols[1],hexcols[3],hexcols[5],hexcols[6],hexcols[8])
 #cmdLine = './swift_fixdt -s -t 16 cosmoVolume.yml'
 #platform = 'KNL'
 
@@ -68,6 +68,9 @@ elif len(sys.argv) == 5:
 elif len(sys.argv) == 6:
   inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])
   numOfSeries = 5
+elif len(sys.argv) == 7:
+  inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5],sys.argv[6])
+  numOfSeries = 6
 
 # Get the names of the branch, Git revision, hydro scheme and hydro kernel
 def parse_header(inputFile):
@@ -76,7 +79,8 @@ def parse_header(inputFile):
     for line in f:
       if 'Branch:' in line:
         s = line.split()
-        branch.append(s[2])
+        line = s[2:]
+        branch.append(" ".join(line))
       elif 'Revision:' in line:
         s = line.split() 
         revision.append(s[2])
@@ -129,8 +133,8 @@ def parse_files():
     parse_header(file_list[0])
     
     version.append(branch[i] + " " + revision[i] + "\n" + hydro_scheme[i] + 
-                   "\n" + hydro_kernel[i] + r", $N_{ngb}$=" + hydro_neighbours[i] + 
-                   r", $\eta$=" + hydro_eta[i])                  
+                   "\n" + hydro_kernel[i] + r", $N_{ngb}=%d$"%float(hydro_neighbours[i]) + 
+                   r", $\eta=%.3f$"%float(hydro_eta[i]))
     times.append([])
     totalTime.append([])
     speedUp.append([])
@@ -176,7 +180,7 @@ def print_results(times,totalTime,parallelEff,version):
 
 def plot_results(times,totalTime,speedUp,parallelEff):
   
-  fig, axarr = plt.subplots(2, 2, figsize=(10,10), frameon=False)
+  fig, axarr = plt.subplots(2, 2, figsize=(10,10), frameon=True)
   speedUpPlot = axarr[0, 0]
   parallelEffPlot = axarr[0, 1]
   totalTimePlot = axarr[1, 0]
@@ -218,7 +222,7 @@ def plot_results(times,totalTime,speedUp,parallelEff):
   totalTimePlot.set_xlim([0.9, 10**(np.floor(np.log10(threadList[i][-1]))+0.5)])
   totalTimePlot.set_ylim([10**np.floor(np.log10(np.min(totalTime)*0.6)), 1.2*10**np.floor(np.log10(np.max(totalTime) * 1.5)+1)])
   
-  totalTimePlot.legend(bbox_to_anchor=(1.14, 1), loc=2, borderaxespad=0.,prop={'size':12})
+  totalTimePlot.legend(bbox_to_anchor=(1.14, 0.97), loc=2, borderaxespad=0.,prop={'size':12}, frameon=False)
   emptyPlot.axis('off')
   
   for i, txt in enumerate(threadList[0]):

src/Makefile.am

@@ -43,7 +43,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
     common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
     physical_constants.h physical_constants_cgs.h potentials.h version.h \
-    hydro_properties.h threadpool.h
+    hydro_properties.h threadpool.h cooling.h
+
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -51,7 +52,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     units.c common_io.c single_io.c multipole.c version.c map.c \
     kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
     physical_constants.c potentials.c hydro_properties.c \
-    runner_doiact_fft.c threadpool.c
+    runner_doiact_fft.c threadpool.c cooling.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -71,7 +72,9 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
 		 hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \
                  hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \
 	         riemann.h riemann/riemann_hllc.h riemann/riemann_trrs.h \
-		 riemann/riemann_exact.h riemann/riemann_vacuum.h 
+		 riemann/riemann_exact.h riemann/riemann_vacuum.h \
+	         cooling/const_du/cooling.h cooling/const_lambda/cooling.h
+
 
 # Sources and flags for regular library
 libswiftsim_la_SOURCES = $(AM_SOURCES)

src/cell.h

@@ -163,6 +163,9 @@ struct cell {
   /* Task for external gravity */
   struct task *grav_external;
 
+  /* Task for cooling */
+  struct task *cooling;
+
   /* Number of tasks that are associated with this cell. */
   int nr_tasks;