Merge branch 'cosmology' into 'master'

Cosmological time integration

See merge request swift/swiftsim!509

.gitignore

@@ -48,6 +48,7 @@ tests/swift_dopair_standard.dat
 tests/brute_force_perturbed.dat
 tests/swift_dopair_perturbed.dat
 tests/test27cells
+tests/test27cells_subset
 tests/testPeriodicBC
 tests/test125cells
 tests/brute_force_27_standard.dat
@@ -121,6 +122,7 @@ theory/Multipoles/potential.pdf
 theory/Multipoles/potential_long.pdf
 theory/Multipoles/potential_short.pdf
 theory/Multipoles/force_short.pdf
+theory/Cosmology/cosmology.pdf
 
 m4/libtool.m4
 m4/ltoptions.m4

INSTALL.swift

@@ -116,6 +116,8 @@ before you can build it.
 	much like the CC one. Use this when your MPI compiler has a
 	none-standard name.
 
+ - GSL: To use cosmological time integration, a version of the GSL
+        must be available. 
 
  - libtool: The build system relies on libtool as well as the other autotools.
 

configure.ac

@@ -407,6 +407,14 @@ AC_HEADER_STDC
 # Check for the libraries we will need.
 AC_CHECK_LIB(m,sqrt,,AC_MSG_ERROR(something is wrong with the math library!))
 
+# Check for GSL
+have_gsl="no"
+AC_CHECK_LIB([gslcblas], [cblas_dgemm])
+AC_CHECK_LIB([gsl], [gsl_integration_qag])
+if test "x$ac_cv_lib_gslcblas_cblas_dgemm" != "x"; then
+   have_gsl="yes"
+fi
+
 # Check for pthreads.
 AX_PTHREAD([LIBS="$PTHREAD_LIBS $LIBS" CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
     CC="$PTHREAD_CC" LDFLAGS="$LDFLAGS $PTHREAD_LIBS $LIBS"],
@@ -1054,6 +1062,7 @@ AC_MSG_RESULT([
     - parallel      : $have_parallel_hdf5
    Metis enabled    : $have_metis
    FFTW3 enabled    : $have_fftw3
+   GSL enabled      : $have_gsl
    libNUMA enabled  : $have_numa
    GRACKLE enabled  : $have_grackle
    Using tcmalloc   : $have_tcmalloc

examples/CosmologicalBox/cosmo.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+#  UnitMass_in_cgs:     1   # Grams
+#  UnitLength_in_cgs:   1   # Centimeters
+#  UnitVelocity_in_cgs: 1   # Centimeters per second
+#  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.    # 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:            uniformBox # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          0.01       # 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).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./uniformBox.hdf5     # The file to read
+
+Cosmology:
+  Omega_m:        0.307
+  Omega_lambda:   0.693
+  Omega_b:        0.0455
+  h:              0.6777
+  a_begin:        0.0078125
+  a_end:          1.0
+  
+

examples/CosmologicalBox/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.
+L = int(sys.argv[1])  # Number of particles along one axis
+rho = 2.              # Density
+P = 1.                # Pressure
+gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
+fileName = "uniformBox.hdf5" 
+
+#---------------------------------------------------
+numPart = L**3
+mass = boxSize**3 * rho / numPart
+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
+grp.attrs["Dimension"] = 3
+
+#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)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+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/EAGLE_12/eagle_12.yml

@@ -13,13 +13,22 @@ TimeIntegration:
   dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
 
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.9090909     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+  
 Scheduler:
   max_top_level_cells:    15
   
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  time_first:          0.    # Time of the first output (in internal units)
+  time_first:          1.    # Time of the first output (in internal units)
   delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics

examples/main.c

@@ -126,6 +126,7 @@ int main(int argc, char *argv[]) {
    * scope.  */
   struct chemistry_data chemistry;
   struct cooling_function_data cooling_func;
+  struct cosmology cosmo;
   struct external_potential potential;
   struct gpart *gparts = NULL;
   struct gravity_props gravity_properties;
@@ -600,6 +601,13 @@ int main(int argc, char *argv[]) {
       phys_const_print(&prog_const);
     }
 
+    /* Initialise the cosmology */
+    if (with_cosmology)
+      cosmology_init(params, &us, &prog_const, &cosmo);
+    else
+      cosmology_init_no_cosmo(&cosmo);
+    if (with_cosmology) cosmology_print(&cosmo);
+
     /* Initialise the hydro properties */
     if (with_hydro) hydro_props_init(&hydro_properties, params);
     if (with_hydro) eos_init(&eos, params);
@@ -767,9 +775,9 @@ int main(int argc, char *argv[]) {
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
     engine_init(&e, &s, params, N_total[0], N_total[1], engine_policies,
-                talking, &reparttype, &us, &prog_const, &hydro_properties,
-                &gravity_properties, &potential, &cooling_func, &chemistry,
-                &sourceterms);
+                talking, &reparttype, &us, &prog_const, &cosmo,
+                &hydro_properties, &gravity_properties, &potential,
+                &cooling_func, &chemistry, &sourceterms);
     engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
                   talking, restart_file);
     if (myrank == 0) {
@@ -790,7 +798,7 @@ int main(int argc, char *argv[]) {
           "from t=%.3e until t=%.3e with %d threads and %d queues "
           "(dt_min=%.3e, "
           "dt_max=%.3e)...",
-          e.timeBegin, e.timeEnd, e.nr_threads, e.sched.nr_queues, e.dt_min,
+          e.time_begin, e.time_end, e.nr_threads, e.sched.nr_queues, e.dt_min,
           e.dt_max);
       fflush(stdout);
     }
@@ -837,9 +845,9 @@ int main(int argc, char *argv[]) {
 
   /* Legend */
   if (myrank == 0)
-    printf("# %6s %14s %14s %12s %12s %12s %16s [%s] %6s\n", "Step", "Time",
-           "Time-step", "Updates", "g-Updates", "s-Updates", "Wall-clock time",
-           clocks_getunit(), "Props");
+    printf("# %6s %14s %14s %9s %12s %12s %12s %16s [%s] %6s\n", "Step", "Time",
+           "Time-step", "Time-bins", "Updates", "g-Updates", "s-Updates",
+           "Wall-clock time", clocks_getunit(), "Props");
 
   /* File for the timers */
   if (with_verbose_timers) timers_open_file(myrank);
@@ -865,8 +873,7 @@ int main(int argc, char *argv[]) {
     if (with_verbose_timers) timers_print(e.step);
 
     /* Every so often allow the user to stop the application and dump the
-     * restart
-     * files. */
+     * restart files. */
     if (j % restart_stop_steps == 0) {
       force_stop = restart_stop_now(restart_dir, 0);
       if (myrank == 0 && force_stop)
@@ -874,8 +881,7 @@ int main(int argc, char *argv[]) {
     }
 
     /* Also if using nsteps to exit, will not have saved any restarts on exit,
-     * make
-     * sure we do that (useful in testing only). */
+     * make sure we do that (useful in testing only). */
     if (force_stop || (e.restart_onexit && e.step - 1 == nsteps))
       engine_dump_restarts(&e, 0, 1);
 

examples/parameter_example.yml

@@ -10,6 +10,36 @@ InternalUnitSystem:
 PhysicalConstants:
   G:            6.67408e-8 # (Optional) Overwrite the value of Newton's constant used internally by the code.
 
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.0078125     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+  Omega_r:        0.            # (Optional) Radiation density parameter
+  w_0:            -1.0          # (Optional) Dark-energy equation-of-state parameter at z=0.
+  w_a:            0.            # (Optional) Dark-energy equation-of-state time evolution parameter.
+
+# 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).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  h_tolerance:           1e-4     # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
+  h_max:                 10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+  max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
+  max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:          0.025    # Constant dimensionless multiplier for time integration.
+  theta:        0.7      # Opening angle (Multipole acceptance criterion)
+  epsilon:      0.1      # Softening length (in internal units).
+  a_smooth:     1.25     # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+  r_cut_max:    4.5      # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+  r_cut_min:    0.1      # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
+
 # Parameters for the task scheduling
 Scheduler:
   nr_queues:                 0         # (Optional) The number of task queues to use. Use 0  to let the system decide.
@@ -48,27 +78,6 @@ Statistics:
   energy_file_name:    energy    # (Optional) File name for energy output
   timestep_file_name:  timesteps # (Optional) File name for timing information output. Note: No underscores "_" allowed in file name
 
-# 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).
-  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-  h_tolerance:           1e-4     # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
-  h_max:                 10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
-  max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
-  max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
-
-EoS:
-  isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units).
-
-# Parameters for the self-gravity scheme
-Gravity:
-  eta:          0.025    # Constant dimensionless multiplier for time integration.
-  theta:        0.7      # Opening angle (Multipole acceptance criterion)
-  epsilon:      0.1      # Softening length (in internal units).
-  a_smooth:     1.25     # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
-  r_cut_max:    4.5      # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
-  r_cut_min:    0.1      # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
-
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  SedovBlast/sedov.hdf5 # The file to read
@@ -108,6 +117,11 @@ DomainDecomposition:
   minfrac:          0.9     # (Optional) Fractional of all particles that should be updated in previous step when
                             # using CPU time trigger
 
+# Parameters related to the equation of state ------------------------------------------
+
+EoS:
+  isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units).
+
 # Parameters related to external potentials --------------------------------------------
 
 # Point mass external potentials

examples/plot_scaling_results.py

@@ -148,8 +148,8 @@ def parse_files():
     
     # Loop over all files for a given series and load the times
     for j in range(0,len(file_list)):
-      times = np.loadtxt(file_list[j],usecols=(6,))
-      updates = np.loadtxt(file_list[j],usecols=(3,))
+      times = np.loadtxt(file_list[j],usecols=(8,))
+      updates = np.loadtxt(file_list[j],usecols=(5,))
       totalTime[i].append(np.sum(times))
       
     sumTotal.append(np.sum(totalTime[i]))

examples/plot_scaling_results_breakdown.py

@@ -149,8 +149,8 @@ def parse_files():
    
     # Loop over all files for a given series and load the times
     for j in range(0,len(file_list)):
-      times = np.loadtxt(file_list[j],usecols=(6,), skiprows=11)
-      updates = np.loadtxt(file_list[j],usecols=(3,), skiprows=11)
+      times = np.loadtxt(file_list[j],usecols=(8,))
+      updates = np.loadtxt(file_list[j],usecols=(5,))
       totalTime[i].append(np.sum(times))
       
     sumTotal.append(np.sum(totalTime[i]))

src/Makefile.am

@@ -47,7 +47,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
     dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
     gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
-    chemistry.h chemistry_io.h chemistry_struct.h restart.h
+    chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -59,7 +59,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
     part_type.c xmf.c gravity_properties.c gravity.c \
     collectgroup.c hydro_space.c equation_of_state.c \
-    chemistry.c restart.c
+    chemistry.c cosmology.c restart.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \

src/active.h

@@ -42,9 +42,9 @@ __attribute__((always_inline)) INLINE static int cell_are_part_drifted(
   if (c->ti_old_part > e->ti_current)
     error(
         "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
-        "and e->ti_current=%lld (t=%e)",
-        c->ti_old_part, c->ti_old_part * e->timeBase, e->ti_current,
-        e->ti_current * e->timeBase);
+        "and e->ti_current=%lld (t=%e, a=%e)",
+        c->ti_old_part, c->ti_old_part * e->time_base, e->ti_current,
+        e->ti_current * e->time_base, e->cosmology->a);
 #endif
 
   return (c->ti_old_part == e->ti_current);
@@ -66,8 +66,8 @@ __attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
     error(
         "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
         "and e->ti_current=%lld (t=%e)",
-        c->ti_old_gpart, c->ti_old_gpart * e->timeBase, e->ti_current,
-        e->ti_current * e->timeBase);
+        c->ti_old_gpart, c->ti_old_gpart * e->time_base, e->ti_current,
+        e->ti_current * e->time_base);
 #endif
 
   return (c->ti_old_gpart == e->ti_current);
@@ -89,9 +89,9 @@ __attribute__((always_inline)) INLINE static int cell_is_active_hydro(
   if (c->ti_hydro_end_min < e->ti_current)
     error(
         "cell in an impossible time-zone! c->ti_end_min=%lld (t=%e) and "
-        "e->ti_current=%lld (t=%e)",
-        c->ti_hydro_end_min, c->ti_hydro_end_min * e->timeBase, e->ti_current,
-        e->ti_current * e->timeBase);
+        "e->ti_current=%lld (t=%e, a=%e)",
+        c->ti_hydro_end_min, c->ti_hydro_end_min * e->time_base, e->ti_current,
+        e->ti_current * e->time_base, e->cosmology->a);
 #endif
 
   return (c->ti_hydro_end_min == e->ti_current);
@@ -132,9 +132,9 @@ __attribute__((always_inline)) INLINE static int cell_is_active_gravity(
   if (c->ti_gravity_end_min < e->ti_current)
     error(
         "cell in an impossible time-zone! c->ti_end_min=%lld (t=%e) and "
-        "e->ti_current=%lld (t=%e)",
-        c->ti_gravity_end_min, c->ti_gravity_end_min * e->timeBase,
-        e->ti_current, e->ti_current * e->timeBase);
+        "e->ti_current=%lld (t=%e, a=%e)",
+        c->ti_gravity_end_min, c->ti_gravity_end_min * e->time_base,
+        e->ti_current, e->ti_current * e->time_base, e->cosmology->a);
 #endif
 
   return (c->ti_gravity_end_min == e->ti_current);
@@ -265,9 +265,9 @@ __attribute__((always_inline)) INLINE static int cell_is_starting_hydro(
   if (c->ti_hydro_beg_max > e->ti_current)
     error(
         "cell in an impossible time-zone! c->ti_beg_max=%lld (t=%e) and "
-        "e->ti_current=%lld (t=%e)",
-        c->ti_hydro_beg_max, c->ti_hydro_beg_max * e->timeBase, e->ti_current,
-        e->ti_current * e->timeBase);
+        "e->ti_current=%lld (t=%e, a=%e)",
+        c->ti_hydro_beg_max, c->ti_hydro_beg_max * e->time_base, e->ti_current,
+        e->ti_current * e->time_base, e->cosmology->a);
 #endif
 
   return (c->ti_hydro_beg_max == e->ti_current);
@@ -287,9 +287,9 @@ __attribute__((always_inline)) INLINE static int cell_is_starting_gravity(
   if (c->ti_gravity_beg_max > e->ti_current)
     error(
         "cell in an impossible time-zone! c->ti_beg_max=%lld (t=%e) and "
-        "e->ti_current=%lld (t=%e)",
-        c->ti_gravity_beg_max, c->ti_gravity_beg_max * e->timeBase,
-        e->ti_current, e->ti_current * e->timeBase);
+        "e->ti_current=%lld (t=%e, a=%e)",
+        c->ti_gravity_beg_max, c->ti_gravity_beg_max * e->time_base,
+        e->ti_current, e->ti_current * e->time_base, e->cosmology->a);
 #endif
 
   return (c->ti_gravity_beg_max == e->ti_current);

src/adiabatic_index.h

@@ -431,4 +431,73 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
 #endif
 }
 
+/**
+ * @brief Return the argument to the power three adiabatic index minus two.
+ *
+ * Computes \f$x^{3\gamma - 2}\f$.
+ *
+ * @param x Argument
+ */
+__attribute__((always_inline)) INLINE static float pow_three_gamma_minus_two(
+    float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return x * x * x; /* x^(3) */
+
+#elif defined(HYDRO_GAMMA_7_5)
+
+  return powf(x, 2.2f); /* x^(11/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return x * x; /* x^(2) */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x * x * x * x; /* x^(4) */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Return the argument to the power three adiabatic index minus five over
+ * two.
+ *
+ * Computes \f$x^{(3\gamma - 5)/2}\f$.
+ *
+ * @param x Argument
+ */
+__attribute__((always_inline)) INLINE static float
+pow_three_gamma_minus_five_over_two(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return 1.f; /* x^(0) */
+
+#elif defined(HYDRO_GAMMA_7_5)
+
+  return powf(x, -0.4f); /* x^(-2/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return 1.f / sqrtf(x); /* x^(-1/2) */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return sqrtf(x); /* x^(1/2) */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
 #endif /* SWIFT_ADIABATIC_INDEX_H */

src/cell.c

@@ -2344,14 +2344,11 @@ int cell_has_tasks(struct cell *c) {
 void cell_drift_part(struct cell *c, const struct engine *e, int force) {
 
   const float hydro_h_max = e->hydro_properties->h_max;
-  const double timeBase = e->timeBase;
   const integertime_t ti_old_part = c->ti_old_part;
   const integertime_t ti_current = e->ti_current;
   struct part *const parts = c->parts;
   struct xpart *const xparts = c->xparts;
 
-  /* Drift from the last time the cell was drifted to the current time */
-  const double dt = (ti_current - ti_old_part) * timeBase;
   float dx_max = 0.f, dx2_max = 0.f;
   float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
   float cell_h_max = 0.f;
@@ -2371,7 +2368,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
   if (c->split && (force || c->do_sub_drift)) {
 
     /* Loop over the progeny and collect their data. */
-    for (int k = 0; k < 8; k++)
+    for (int k = 0; k < 8; k++) {
       if (c->progeny[k] != NULL) {
         struct cell *cp = c->progeny[k];
 
@@ -2383,6 +2380,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
         dx_max_sort = max(dx_max_sort, cp->dx_max_sort);
         cell_h_max = max(cell_h_max, cp->h_max);
       }
+    }
 
     /* Store the values */
     c->h_max = cell_h_max;
@@ -2394,6 +2392,24 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
 
   } else if (!c->split && force && ti_current > ti_old_part) {
 
+    /* Drift from the last time the cell was drifted to the current time */