Merge remote-tracking branch 'origin/master' into scheduler_activate_root

examples/CoolingHaloWithSpin/cooling_halo.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   10.    # The end time of the simulation (in internal units).
-  dt_min:     1e-4  # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the conserved quantities statistics
@@ -32,9 +32,6 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  CoolingHalo.hdf5       # The file to read
-  shift_x:    0.                  # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    0.
-  shift_z:    0.
  
 # External potential parameters
 SoftenedIsothermalPotential:
@@ -43,12 +40,12 @@ SoftenedIsothermalPotential:
   position_z:      0.	
   vrot:            200.     # rotation speed of isothermal potential in internal units
   timestep_mult:   0.03     # controls time step
-  epsilon:         0.1    #softening for the isothermal potential
+  epsilon:         1.0      #softening for the isothermal potential
 
 # Cooling parameters
 LambdaCooling:
-  lambda_cgs:                  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)
-  cooling_tstep_mult:          1.0    # Dimensionless pre-factor for the time-step condition
+  cooling_tstep_mult:          0.1    # Dimensionless pre-factor for the time-step condition

examples/HydrostaticHalo/README

+Hydrostatic halo in equilibrium in an isothermal potential. Running
+for 10 dynamical times.
 
 To make the initial conditions we distribute gas particles randomly in
 a cube with a side length twice that of the virial radius. The density
 profile of the gas is proportional to r^(-2) where r is the distance
 from the centre of the cube.
 
-The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
+The parameter v_rot (in makeIC.py and hydrostatic.yml) sets the circular
 velocity of the halo, and by extension, the viral radius, viral mass,
 and the internal energy of the gas such that hydrostatic equilibrium
 is achieved.
@@ -12,10 +14,12 @@ is achieved.
 To run this example, make such that the code is compiled with either
 the isothermal potential or softened isothermal potential set in
 src/const.h. In the latter case, a (small) value of epsilon needs to
-be set in cooling.yml.  0.1 kpc should work well.
+be set in hydrostatic.yml.  ~1 kpc should work well.
 
 The plotting scripts produce a plot of the density, internal energy
-and radial velocity profile for each
-snapshot. test_energy_conservation.py shows the evolution of energy
+and radial velocity profile for each snapshot and divides the profile
+by the expected profile. 
+
+The script test_energy_conservation.py shows the evolution of energy
 with time. These can be used to check if the example has run properly.
 

examples/HydrostaticHalo/hydrostatic.yml

@@ -9,7 +9,7 @@ InternalUnitSystem:
 # 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).
+  time_end:   30.   # 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).
 
@@ -20,21 +20,18 @@ Statistics:
 # Parameters governing the snapshots
 Snapshots:
   basename:            Hydrostatic  # 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)
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          0.1          # Time difference between consecutive outputs (in internal units)
 
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
-  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
 
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  Hydrostatic.hdf5       # The file to read
-  shift_x:    0.                  # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    0.
-  shift_z:    0.
  
 # External potential parameters
 SoftenedIsothermalPotential:
@@ -42,6 +39,6 @@ SoftenedIsothermalPotential:
   position_y:      0.
   position_z:      0.	
   vrot:            200.     # rotation speed of isothermal potential in internal units
-  epsilon:         0.1
+  epsilon:         1.0
   timestep_mult:   0.03     # controls time step
 

examples/HydrostaticHalo/makeIC.py

@@ -116,9 +116,9 @@ print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
 print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
 print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
 
-print np.mean(coords[:,0])
-print np.mean(coords[:,1])
-print np.mean(coords[:,2])
+#print np.mean(coords[:,0])
+#print np.mean(coords[:,1])
+#print np.mean(coords[:,2])
 
 #now find the particles which are within the box
 

examples/HydrostaticHalo/run.sh

@@ -2,14 +2,23 @@
 
 # Generate the initial conditions if they are not present.
 echo "Generating initial conditions for the isothermal potential box example..."
-python makeIC.py 10000 
+python makeIC.py 100000
 
-../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
+# Run for 10 dynamical times
+../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
 
-python density_profile.py 2. 200 100
+echo "Plotting density profiles"
+mkdir plots
+mkdir plots/density_profile
+python density_profile.py 2. 200 300
 
-python internal_energy_profile.py 2. 200 100
+echo "Plotting internal energy profiles"
+mkdir plots/internal_energy
+python internal_energy_profile.py 2. 200 300
 
-python velocity_profile.py 2. 200 100
+echo "Plotting radial velocity profiles"
+mkdir plots/radial_velocity_profile
+python velocity_profile.py 2. 200 300
 
-python test_energy_conservation.py
+echo "Plotting energy as a function of time"
+python test_energy_conservation.py 300

examples/HydrostaticHalo/test_energy_conservation.py

@@ -3,7 +3,7 @@ import h5py as h5
 import matplotlib.pyplot as plt
 import sys
 
-n_snaps = 5
+n_snaps = int(sys.argv[1])
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
@@ -24,7 +24,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["IsothermalPotential:vrot"])
+v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]

examples/IsothermalPotential/GravityOnly/makeIC.py → examples/IsothermalPotential/makeIC.py

@@ -141,17 +141,17 @@ ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
 ds[()] = v
 v = numpy.zeros(1)
 
-m = numpy.full((numPart, ), mass)
+m = numpy.full((numPart, ), mass, dtype='f')
 ds = grp1.create_dataset('Masses', (numPart,), 'f')
 ds[()] = m
 m = numpy.zeros(1)
 
-h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+h = numpy.full((numPart, ), 1.1255 * boxSize / L,  dtype='f')
 ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
 ds[()] = h
 h = numpy.zeros(1)
 
-u = numpy.full((numPart, ), internalEnergy)
+u = numpy.full((numPart, ), internalEnergy,  dtype='f')
 ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
 ds[()] = u
 u = numpy.zeros(1)

examples/IsothermalPotential/GravityOnly/run.sh → examples/IsothermalPotential/run.sh

@@ -7,4 +7,4 @@ then
     python makeIC.py 1000 1
 fi
 
-../../swift -g -t 2 isothermal.yml 2>&1 | tee output.log
+../swift -g -t 1 isothermal.yml 2>&1 | tee output.log

examples/Makefile.am

@@ -65,6 +65,9 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
 	     EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
 	     ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/test.pro \
 	     GreshoVortex_2D/getGlass.sh GreshoVortex_2D/gresho.yml GreshoVortex_2D/makeIC.py GreshoVortex_2D/plotSolution.py GreshoVortex_2D/run.sh \
+	     HydrostaticHalo/README HydrostaticHalo/hydrostatic.yml HydrostaticHalo/makeIC.py HydrostaticHalo/run.sh \
+	     HydrostaticHalo/density_profile.py HydrostaticHalo/velocity_profile.py HydrostaticHalo/internal_energy_profile.py HydrostaticHalo/test_energy_conservation.py \
+	     IsothermalPotential/README IsothermalPotential/run.sh IsothermalPotential/test.pro IsothermalPotential/isothermal.yml IsothermalPotential/makeIC.py \
 	     KelvinHelmholtz_2D/kelvinHelmholtz.yml KelvinHelmholtz_2D/makeIC.py KelvinHelmholtz_2D/plotSolution.py KelvinHelmholtz_2D/run.sh \
 	     MultiTypes/makeIC.py  MultiTypes/multiTypes.yml MultiTypes/run.sh \
 	     PerturbedBox_2D/makeIC.py PerturbedBox_2D/perturbedPlane.yml \

m4/ax_cc_maxopt.m4

@@ -160,6 +160,9 @@ if test "$ac_test_CFLAGS" != "set"; then
      # note that we enable "unsafe" fp optimization with other compilers, too
      AX_CHECK_COMPILE_FLAG(-ffast-math, CFLAGS="$CFLAGS -ffast-math")
 
+     # not all codes will benefit from this.
+     AX_CHECK_COMPILE_FLAG(-funroll-loops, CFLAGS="$CFLAGS -funroll-loops")
+
      AX_GCC_ARCHFLAG($acx_maxopt_portable)
      ;;
 

src/Makefile.am

@@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
                  units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 dimension.h equation_of_state.h \
+		 dimension.h equation_of_state.h active.h \
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \

src/active.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_ACTIVE_H
+#define SWIFT_ACTIVE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "cell.h"
+#include "const.h"
+#include "engine.h"
+#include "part.h"
+
+/**
+ * @brief Does a cell contain any active particle ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_active(
+    const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->ti_end_min < e->ti_current)
+    error("cell in an impossible time-zone! c->ti_end_min=%d e->ti_current=%d",
+          c->ti_end_min, e->ti_current);
+#endif
+
+  return (c->ti_end_min == e->ti_current);
+}
+
+/**
+ * @brief Are *all* particles in a cell active ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_all_active(
+    const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->ti_end_max < e->ti_current)
+    error("cell in an impossible time-zone! c->ti_end_max=%d e->ti_current=%d",
+          c->ti_end_max, e->ti_current);
+#endif
+
+  return (c->ti_end_max == e->ti_current);
+}
+
+/**
+ * @brief Is this particle active ?
+ *
+ * @param p The #part.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int part_is_active(
+    const struct part *p, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (p->ti_end < e->ti_current)
+    error("particle in an impossible time-zone! p->ti_end=%d e->ti_current=%d",
+          p->ti_end, e->ti_current);
+#endif
+
+  return (p->ti_end == e->ti_current);
+}
+
+/**
+ * @brief Is this g-particle active ?
+ *
+ * @param gp The #gpart.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int gpart_is_active(
+    const struct gpart *gp, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (gp->ti_end < e->ti_current)
+    error(
+        "g-particle in an impossible time-zone! gp->ti_end=%d e->ti_current=%d",
+        gp->ti_end, e->ti_current);
+#endif
+
+  return (gp->ti_end == e->ti_current);
+}
+
+#endif /* SWIFT_ACTIVE_H */

src/cell.c

@@ -47,6 +47,7 @@
 #include "cell.h"
 
 /* Local headers. */
+#include "active.h"
 #include "atomic.h"
 #include "error.h"
 #include "gravity.h"
@@ -867,14 +868,14 @@ void cell_clean(struct cell *c) {
  * @brief Checks whether a given cell needs drifting or not.
  *
  * @param c the #cell.
- * @param ti_current The current time on the integer time-line.
+ * @param e The #engine (holding current time information).
  *
  * @return 1 If the cell needs drifting, 0 otherwise.
  */
-int cell_is_drift_needed(struct cell *c, int ti_current) {
+int cell_is_drift_needed(struct cell *c, const struct engine *e) {
 
   /* Do we have at least one active particle in the cell ?*/
-  if (c->ti_end_min == ti_current) return 1;
+  if (cell_is_active(c, e)) return 1;
 
   /* Loop over the pair tasks that involve this cell */
   for (struct link *l = c->density; l != NULL; l = l->next) {
@@ -882,9 +883,9 @@ int cell_is_drift_needed(struct cell *c, int ti_current) {
     if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair)
       continue;
 
-    /* Does the other cell in the pair have an active particle ? */
-    if ((l->t->ci == c && l->t->cj->ti_end_min == ti_current) ||
-        (l->t->cj == c && l->t->ci->ti_end_min == ti_current))
+    /* Is the other cell in the pair active ? */
+    if ((l->t->ci == c && cell_is_active(l->t->cj, e)) ||
+        (l->t->cj == c && cell_is_active(l->t->ci, e)))
       return 1;
   }
 

src/cell.h

@@ -37,6 +37,7 @@
 #include "task.h"
 
 /* Avoid cyclic inclusions */
+struct engine;
 struct space;
 struct scheduler;
 
@@ -290,7 +291,7 @@ int cell_are_neighbours(const struct cell *restrict ci,
                         const struct cell *restrict cj);
 void cell_check_multipole(struct cell *c, void *data);
 void cell_clean(struct cell *c);
-int cell_is_drift_needed(struct cell *c, int ti_current);
+int cell_is_drift_needed(struct cell *c, const struct engine *e);
 int cell_unskip_tasks(struct cell *c, struct scheduler *s);
 void cell_set_super(struct cell *c, struct cell *super);
 

src/cooling/const_lambda/cooling.h

@@ -24,6 +24,7 @@
 #define SWIFT_COOLING_CONST_LAMBDA_H
 
 /* Some standard headers. */
+#include <float.h>
 #include <math.h>
 
 /* Local includes. */
@@ -84,7 +85,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   /* Calculate du_dt */
   const float du_dt = cooling_rate(phys_const, us, cooling, p);
 
-  /* Intergrate cooling equation, but enforce energy floor */
+  /* Integrate cooling equation, but enforce energy floor */
   float u_new;
   if (u_old + du_dt * dt > u_floor) {
     u_new = u_old + du_dt * dt;
@@ -92,6 +93,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
     u_new = u_floor;
   }
 
+  /* Don't allow particle to cool too much in one timestep */
+  if (u_new < 0.5f * u_old) u_new = 0.5f * u_old;
+
   /* Update the internal energy */
   hydro_set_internal_energy(p, u_new);
 
@@ -112,13 +116,16 @@ __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 du_dt */
-  const float du_dt = cooling_rate(phys_const, us, cooling, p);
-
   /* Get current internal energy (dt=0) */
   const float u = hydro_get_internal_energy(p, 0.f);
+  const float du_dt = cooling_rate(phys_const, us, cooling, p);
 
-  return u / fabsf(du_dt);
+  /* If we are close to (or below) the energy floor, we ignore cooling timestep
+   */
+  if (u < 1.01f * cooling->min_energy)
+    return FLT_MAX;
+  else
+    return cooling->cooling_tstep_mult * u / fabsf(du_dt);
 }
 
 /**

src/debug.c

@@ -193,6 +193,59 @@ int checkSpacehmax(struct space *s) {
   return 0;
 }
 
+/**
+ * @brief Check if the h_max and dx_max values of a cell's hierarchy are
+ * consistent with the particles. Report verbosely if not.
+ *
+ * @param c the top cell of the hierarchy.
+ * @param depth the recursion depth for use in messages. Set to 0 initially.
+ * @result 1 or 0
+ */
+int checkCellhdxmax(const struct cell *c, int *depth) {
+
+  *depth = *depth + 1;
+
+  float h_max = 0.0f;
+  float dx_max = 0.0f;
+  if (!c->split) {
+    const size_t nr_parts = c->count;
+    struct part *parts = c->parts;
+    for (size_t k = 0; k < nr_parts; k++) {
+      h_max = (h_max > parts[k].h) ? h_max : parts[k].h;
+    }
+  } else {
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) {
+        struct cell *cp = c->progeny[k];
+        checkCellhdxmax(cp, depth);
+        dx_max = max(dx_max, cp->dx_max);
+        h_max = max(h_max, cp->h_max);
+      }
+  }
+
+  /* Check. */
+  int result = 1;
+  if (c->h_max != h_max) {
+    message("%d Inconsistent h_max: cell %f != parts %f", *depth, c->h_max, h_max);
+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+    result = 0;
+  }
+  if (c->dx_max != dx_max) {
+    message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+    result = 0;
+  }
+
+  /* Check rebuild criterion. */
+  if (h_max > c->dmin) {
+    message("%d Inconsistent c->dmin: %f > %f", *depth, h_max, c->dmin);
+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+    result = 0;
+  }
+
+  return result;
+}
+
 #ifdef HAVE_METIS
 
 /**

src/debug.h

@@ -32,6 +32,7 @@ void printParticle_single(const struct part *p, const struct xpart *xp);
 void printgParticle_single(struct gpart *gp);
 
 int checkSpacehmax(struct space *s);
+int checkCellhdxmax(const struct cell *c, int *depth);
 
 #ifdef HAVE_METIS
 #include "metis.h"