Merged master into shadowswift_new.

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:    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_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;
@@ -332,12 +337,6 @@ int main(int argc, char *argv[]) {
   hydro_properties.delta_neighbours = 0.0f;
 #endif
 
-  /* 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);
@@ -440,6 +439,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;
@@ -447,13 +457,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),

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
+

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;
 

src/const.h

@@ -95,6 +95,11 @@
 //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 //#define EXTERNAL_POTENTIAL_DISK_PATCH
 
+/* Cooling properties */
+//#define COOLING_CONST_DU
+#define COOLING_CONST_LAMBDA
+//#define COOLING_GRACKLE
+
 /* Are we debugging ? */
 //#define SWIFT_DEBUG_CHECKS
 

src/cooling.c

+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "cooling.h"
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * Calls cooling_init_backend for the chosen cooling function.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param cooling The cooling properties to initialize
+ */
+void cooling_init(const struct swift_params* parameter_file,
+                  const struct UnitSystem* us,
+                  const struct phys_const* phys_const,
+                  struct cooling_data* cooling) {
+
+  cooling_init_backend(parameter_file, us, phys_const, cooling);
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * Calls cooling_print_backend for the chosen cooling function.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+void cooling_print(const struct cooling_data* cooling) {
+
+  cooling_print_backend(cooling);
+}

src/cooling.h

+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_H
+#define SWIFT_COOLING_H
+
+/**
+ * @file src/cooling.h
+ * @brief Branches between the different cooling functions.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+
+/* Import the right cooling definition */
+#if defined(COOLING_CONST_DU)
+#include "./cooling/const_du/cooling.h"
+#elif defined(COOLING_CONST_LAMBDA)
+#include "./cooling/const_lambda/cooling.h"
+#elif defined(COOLING_GRACKLE)
+#include "./cooling/grackle/cooling.h"
+#else
+#error "Invalid choice of cooling function."
+#endif
+
+/* Common functions */
+void cooling_init(const struct swift_params* parameter_file,
+                  const struct UnitSystem* us,
+                  const struct phys_const* phys_const,
+                  struct cooling_data* cooling);
+
+void cooling_print(const struct cooling_data* cooling);
+
+#endif /* SWIFT_COOLING_H */

src/cooling/const_du/cooling.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *                    Richard Bower (r.g.bower@durham.ac.uk)
+ *                    Stefan Arridge  (stefan.arridge@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_COOLING_CONST_DU_H
+#define SWIFT_COOLING_CONST_DU_H
+
+/**
+ * @file src/cooling/const/cooling.h
+ * @brief Routines related to the "constant cooling" cooling function.
+ *
+ * This is the simplest possible cooling function. A constant cooling rate with
+ * a minimal energy floor is applied.
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_data {
+
+  /*! Cooling rate in internal units. du_dt = -cooling_rate */
+  float cooling_rate;
+
+  /*! Minimally allowed internal energy of the particles */
+  float min_energy;
+
+  /*! Constant multiplication factor for time-step criterion */
+  float cooling_tstep_mult;
+};
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+    const struct phys_const* const phys_const, const struct UnitSystem* us,
+    const struct cooling_data* cooling, struct part* p, float dt) {
+
+  /* Get current internal energy (dt=0) */
+  const float u_old = hydro_get_internal_energy(p, 0.f);
+
+  /* Get cooling function properties */
+  const float du_dt = -cooling->cooling_rate;
+  const float u_floor = cooling->min_energy;
+
+  /* Constant cooling with a minimal floor */
+  float u_new;
+  if (u_old - du_dt * dt > u_floor) {
+    u_new = u_old + du_dt * dt;
+  } else {
+    u_new = u_floor;
+  }
+
+  /* Update the internal energy */
+  hydro_set_internal_energy(p, u_new);
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * @param cooling The #cooling_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double cooling_timestep(
+    const struct cooling_data* cooling,
+    const struct phys_const* const phys_const, const struct part* const p) {
+
+  const float cooling_rate = cooling->cooling_rate;
+  const float internal_energy =
+      hydro_get_internal_energy(p, 0);  // dt = 0 because using current entropy
+  return cooling->cooling_tstep_mult * internal_energy / cooling_rate;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param cooling The cooling properties to initialize
+ */
+static INLINE void cooling_init_backend(
+    const struct swift_params* parameter_file, const struct UnitSystem* us,
+    const struct phys_const* phys_const, struct cooling_data* cooling) {
+
+  cooling->cooling_rate =
+      parser_get_param_double(parameter_file, "ConstCooling:cooling_rate");
+  cooling->min_energy =
+      parser_get_param_double(parameter_file, "ConstCooling:min_energy");
+  cooling->cooling_tstep_mult = parser_get_param_double(
+      parameter_file, "ConstCooling:cooling_tstep_mult");
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+static INLINE void cooling_print_backend(const struct cooling_data* cooling) {
+
+  message("Cooling function is 'Constant cooling' with rate %f and floor %f",
+          cooling->cooling_rate, cooling->min_energy);
+}
+
+#endif /* SWIFT_COOLING_CONST_DU_H */