/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 *                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@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>

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif

/* Switch off timers. */
#ifdef TIMER
#undef TIMER
#endif

/* This object's header. */
#include "cell.h"

/* Local headers. */
#include "engine.h"
#include "error.h"
#include "multipole.h"
#include "multipole_accept.h"
#include "space.h"
#include "tools.h"

/* Global variables. */
int cell_next_tag = 0;

/** List of cell pairs for sub-cell recursion. For any sid, the entries in
 * this array contain the number of sub-cell pairs and the indices and sid
 * of the sub-cell pairs themselves. */
struct cell_split_pair cell_split_pairs[13] = {
    {1, /* (  1 ,  1 ,  1 ) */
     {{7, 0, 0}}},

    {4, /* (  1 ,  1 ,  0 ) */
     {{6, 0, 1}, {7, 1, 1}, {6, 1, 0}, {7, 0, 2}}},

    {1, /* (  1 ,  1 , -1 ) */
     {{6, 1, 2}}},

    {4, /* (  1 ,  0 ,  1 ) */
     {{5, 0, 3}, {7, 2, 3}, {5, 2, 0}, {7, 0, 6}}},

    {16, /* (  1 ,  0 ,  0 ) */
     {{4, 0, 4},
      {5, 0, 5},
      {6, 0, 7},
      {7, 0, 8},
      {4, 1, 3},
      {5, 1, 4},
      {6, 1, 6},
      {7, 1, 7},
      {4, 2, 1},
      {5, 2, 2},
      {6, 2, 4},
      {7, 2, 5},
      {4, 3, 0},
      {5, 3, 1},
      {6, 3, 3},
      {7, 3, 4}}},

    {4, /* (  1 ,  0 , -1 ) */
     {{4, 1, 5}, {6, 3, 5}, {4, 3, 2}, {6, 1, 8}}},

    {1, /* (  1 , -1 ,  1 ) */
     {{5, 2, 6}}},

    {4, /* (  1 , -1 ,  0 ) */
     {{4, 3, 6}, {5, 2, 8}, {4, 2, 7}, {5, 3, 7}}},

    {1, /* (  1 , -1 , -1 ) */
     {{4, 3, 8}}},

    {4, /* (  0 ,  1 ,  1 ) */
     {{3, 0, 9}, {7, 4, 9}, {3, 4, 0}, {7, 0, 8}}},

    {16, /* (  0 ,  1 ,  0 ) */
     {{2, 0, 10},
      {3, 0, 11},
      {6, 0, 7},
      {7, 0, 6},
      {2, 1, 9},
      {3, 1, 10},
      {6, 1, 8},
      {7, 1, 7},
      {2, 4, 1},
      {3, 4, 2},
      {6, 4, 10},
      {7, 4, 11},
      {2, 5, 0},
      {3, 5, 1},
      {6, 5, 9},
      {7, 5, 10}}},

    {4, /* (  0 ,  1 , -1 ) */
     {{2, 1, 11}, {6, 5, 11}, {2, 5, 2}, {6, 1, 6}}},

    {16, /* (  0 ,  0 ,  1 ) */
     {{1, 0, 12},
      {3, 0, 11},
      {5, 0, 5},
      {7, 0, 2},
      {1, 2, 9},
      {3, 2, 12},
      {5, 2, 8},
      {7, 2, 5},
      {1, 4, 3},
      {3, 4, 6},
      {5, 4, 12},
      {7, 4, 11},
      {1, 6, 0},
      {3, 6, 3},
      {5, 6, 9},
      {7, 6, 12}}}};

/* Cell names. */
const char *cellID_names[4] = {
    "Regular",
    "Zoom",
    "Buffer",
    "Background",
};

/* Sub-cell names. */
const char *subcellID_names[4] = {
    "Regular",
    "Neighbour",
    "Void",
};

/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */
int cell_get_tree_size(struct cell *c) {
  /* Number of cells in this subtree. */
  int count = 1;

  /* Sum up the progeny if split. */
  if (c->split)
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL) count += cell_get_tree_size(c->progeny[k]);

  /* Return the final count. */
  return count;
}

/**
 * @brief Link the cells recursively to the given #part array.
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
int cell_link_parts(struct cell *c, struct part *parts) {
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

  if (c->hydro.parts != NULL)
    error("Linking parts into a cell that was already linked");
#endif

  c->hydro.parts = parts;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->hydro.count;
}

/**
 * @brief Link the cells recursively to the given #gpart array.
 *
 * @param c The #cell.
 * @param gparts The #gpart array.
 *
 * @return The number of particles linked.
 */
int cell_link_gparts(struct cell *c, struct gpart_foreign *gparts_foreign) {
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

  if (c->grav.parts_foreign != NULL)
    error("Linking gparts into a cell that was already linked");
#endif

  c->grav.parts_foreign = gparts_foreign;
  c->grav.parts_foreign_rebuild = gparts_foreign;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_gparts(c->progeny[k], &gparts_foreign[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->grav.count;
}

/**
 * @brief Link the cells recursively to the given #gpart array.
 *
 * @param c The #cell.
 * @param gparts The #gpart array.
 *
 * @return The number of particles linked.
 */
int cell_link_fof_gparts(struct cell *c,
                         struct gpart_fof_foreign *gparts_fof_foreign) {
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

  if (c->grav.parts_fof_foreign != NULL)
    error("Linking gparts into a cell that was already linked");
#endif

  c->grav.parts_fof_foreign = gparts_fof_foreign;
  c->grav.parts_fof_foreign_rebuild = gparts_fof_foreign;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset +=
            cell_link_fof_gparts(c->progeny[k], &gparts_fof_foreign[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->grav.count;
}

/**
 * @brief Link the cells recursively to the given #spart array.
 *
 * @param c The #cell.
 * @param sparts The #spart array.
 *
 * @return The number of particles linked.
 */
int cell_link_sparts(struct cell *c, struct spart *sparts) {
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

  if (c->stars.parts != NULL)
    error("Linking sparts into a cell that was already linked");
#endif

  c->stars.parts = sparts;
  c->stars.parts_rebuild = sparts;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_sparts(c->progeny[k], &sparts[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->stars.count;
}

/**
 * @brief Link the cells recursively to the given #bpart array.
 *
 * @param c The #cell.
 * @param bparts The #bpart array.
 *
 * @return The number of particles linked.
 */
int cell_link_bparts(struct cell *c, struct bpart *bparts) {

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

  if (c->black_holes.parts != NULL)
    error("Linking bparts into a cell that was already linked");
#endif

  c->black_holes.parts = bparts;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_bparts(c->progeny[k], &bparts[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->black_holes.count;
}

/**
 * @brief Link the cells recursively to the given #sink array.
 *
 * @param c The #cell.
 * @param sinks The #sink array.
 *
 * @return The number of particles linked.
 */
int cell_link_sinks(struct cell *c, struct sink *sinks) {
#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");

  if (c->sinks.parts != NULL)
    error("Linking sparts into a cell that was already linked");
#endif

  c->sinks.parts = sinks;
  c->sinks.parts_rebuild = sinks;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_sinks(c->progeny[k], &sinks[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->sinks.count;
}

/**
 * @brief Recurse down foreign cells until reaching one with hydro
 * tasks; then trigger the linking of the #part array from that
 * level.
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */
int cell_link_foreign_parts(struct cell *c, struct part *parts) {
#ifdef WITH_MPI

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");
#endif

  /* Do we have a hydro task at this level? */
  if (cell_get_recv(c, task_subtype_xv) != NULL) {

    /* Recursively attach the parts */
    const int counts = cell_link_parts(c, parts);
#ifdef SWIFT_DEBUG_CHECKS
    if (counts != c->hydro.count)
      error("Something is wrong with the foreign counts");
#endif
    return counts;
  }

  /* Go deeper to find the level where the tasks are */
  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        count += cell_link_foreign_parts(c->progeny[k], &parts[count]);
      }
    }
    return count;
  } else {
    return 0;
  }

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
}

/**
 * @brief Recurse down foreign cells until reaching one with gravity
 * tasks; then trigger the linking of the #gpart_foreign array from that
 * level.
 *
 * @param c The #cell.
 * @param gparts_foreign The #gpart_foreign array.
 *
 * @return The number of particles linked.
 */
int cell_link_foreign_gparts(struct cell *c,
                             struct gpart_foreign *gparts_foreign) {
#ifdef WITH_MPI

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");
#endif

  /* Do we have a gravity task at this level? */
  if (cell_get_recv(c, task_subtype_gpart) != NULL) {

    /* Recursively attach the gparts */
    const int counts = cell_link_gparts(c, gparts_foreign);
#ifdef SWIFT_DEBUG_CHECKS
    if (counts != c->grav.count)
      error("Something is wrong with the foreign counts");
#endif
    return counts;
  } else {
    c->grav.parts_foreign = gparts_foreign;
    c->grav.parts_foreign_rebuild = gparts_foreign;
  }

  /* Go deeper to find the level where the tasks are */
  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        count +=
            cell_link_foreign_gparts(c->progeny[k], &gparts_foreign[count]);
      }
    }
    return count;
  } else {
    return 0;
  }

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
}

/**
 * @brief Recurse down foreign cells until reaching one with gravity
 * tasks; then trigger the linking of the #gpart_fof_foreign array from that
 * level.
 *
 * @param c The #cell.
 * @param gparts_fof_foreign The #gpart_fof_foreign array.
 *
 * @return The number of particles linked.
 */
int cell_link_foreign_fof_gparts(struct cell *c,
                                 struct gpart_fof_foreign *gparts_fof_foreign) {
#ifdef WITH_MPI

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Linking foreign particles in a local cell!");
#endif

  /* Do we have a gravity task at this level? */
  if (cell_get_recv(c, task_subtype_gpart) != NULL) {

    /* Recursively attach the gparts */
    const int counts = cell_link_fof_gparts(c, gparts_fof_foreign);
#ifdef SWIFT_DEBUG_CHECKS
    if (counts != c->grav.count)
      error("Something is wrong with the foreign counts");
#endif
    return counts;
  } else {
    c->grav.parts_fof_foreign = gparts_fof_foreign;
    c->grav.parts_fof_foreign_rebuild = gparts_fof_foreign;
  }

  /* Go deeper to find the level where the tasks are */
  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        count += cell_link_foreign_fof_gparts(c->progeny[k],
                                              &gparts_fof_foreign[count]);
      }
    }
    return count;
  } else {
    return 0;
  }

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
}

/**
 * @brief Recursively nullify all the particle pointers in a cell hierarchy.
 *
 * Should only be used on foreign cells!
 *
 * This will make any task or action running on these cells likely crash.
 * Recreating the foreign links will be necessary.
 *
 * @param c The #cell to act on.
 */
void cell_unlink_foreign_particles(struct cell *c) {

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Unlinking foreign particles in a local cell!");
#endif

  c->grav.parts = NULL;
  c->hydro.parts = NULL;
  c->stars.parts = NULL;
  c->black_holes.parts = NULL;
  c->sinks.parts = NULL;

  if (c->split) {
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        cell_unlink_foreign_particles(c->progeny[k]);
      }
    }
  }
}

/**
 * @brief Recursively count the number of #part in foreign cells that
 * are in cells with hydro-related tasks.
 *
 * @param c The #cell.
 *
 * @return The number of particles linked.
 */
int cell_count_parts_for_tasks(const struct cell *c) {
#ifdef WITH_MPI

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Counting foreign particles in a local cell!");
#endif

  /* Do we have a hydro task at this level? */
  if (cell_get_recv(c, task_subtype_xv) != NULL) {
    return c->hydro.count;
  }

  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        count += cell_count_parts_for_tasks(c->progeny[k]);
      }
    }
    return count;
  } else {
    return 0;
  }

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
}

/**
 * @brief Recursively count the number of #gpart in foreign cells that
 * are in cells with gravity-related tasks.
 *
 * @param c The #cell.
 *
 * @return The number of particles linked.
 */
int cell_count_gparts_for_tasks(const struct cell *c) {
#ifdef WITH_MPI

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID == engine_rank)
    error("Counting foreign particles in a local cell!");
#endif

  /* Do we have a gravity task at this level? */
  if (cell_get_recv(c, task_subtype_gpart) != NULL) {
    return c->grav.count;
  }

  if (c->split) {
    int count = 0;
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        count += cell_count_gparts_for_tasks(c->progeny[k]);
      }
    }
    return count;
  } else {
    return 0;
  }

#else
  error("Calling linking of foregin particles in non-MPI mode.");
#endif
}

/**
 * @brief Sanitizes the smoothing length values of cells by setting large
 * outliers to more sensible values.
 *
 * Each cell with <1000 part will be processed. We limit h to be the size of
 * the cell and replace 0s with a good estimate.
 *
 * @param c The cell.
 * @param treated Has the cell already been sanitized at this level ?
 */
void cell_sanitize(struct cell *c, int treated) {
  const int count = c->hydro.count;
  const int scount = c->stars.count;
  struct part *parts = c->hydro.parts;
  struct spart *sparts = c->stars.parts;
  float h_max = 0.f;
  float h_max_active = 0.f;
  float stars_h_max = 0.f;
  float stars_h_max_active = 0.f;

  /* Treat cells will <1000 particles */
  if (count < 1000 && !treated) {
    /* Get an upper bound on h */
    const float upper_h_max = c->dmin / (1.2f * kernel_gamma);

    /* Apply it */
    for (int i = 0; i < count; ++i) {
      if (parts[i].h == 0.f || parts[i].h > upper_h_max)
        parts[i].h = upper_h_max;
    }
    for (int i = 0; i < scount; ++i) {
      if (sparts[i].h == 0.f || sparts[i].h > upper_h_max)
        sparts[i].h = upper_h_max;
    }
  }

  /* Recurse and gather the new h_max values */
  if (c->split) {
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        /* Recurse */
        cell_sanitize(c->progeny[k], (count < 1000));

        /* And collect */
        h_max = max(h_max, c->progeny[k]->hydro.h_max);
        h_max_active = max(h_max_active, c->progeny[k]->hydro.h_max_active);
        stars_h_max = max(stars_h_max, c->progeny[k]->stars.h_max);
        stars_h_max_active =
            max(stars_h_max_active, c->progeny[k]->stars.h_max_active);
      }
    }
  } else {
    /* Get the new value of h_max (note all particles are active) */
    for (int i = 0; i < count; ++i) h_max = max(h_max, parts[i].h);
    for (int i = 0; i < count; ++i)
      h_max_active = max(h_max_active, parts[i].h);
    for (int i = 0; i < scount; ++i)
      stars_h_max = max(stars_h_max, sparts[i].h);
    for (int i = 0; i < scount; ++i)
      stars_h_max_active = max(stars_h_max_active, sparts[i].h);
  }

  /* Record the change */
  c->hydro.h_max = h_max;
  c->hydro.h_max_active = h_max_active;
  c->stars.h_max = stars_h_max;
  c->stars.h_max_active = stars_h_max_active;
}

/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_clean_links(struct cell *c, void *data) {
  c->hydro.density = NULL;
  c->hydro.gradient = NULL;
  c->hydro.force = NULL;
  c->hydro.limiter = NULL;
  c->rt.rt_gradient = NULL;
  c->rt.rt_transport = NULL;
  c->grav.grav = NULL;
  c->grav.mm = NULL;
  c->stars.density = NULL;
  c->stars.prepare1 = NULL;
  c->stars.prepare2 = NULL;
  c->stars.feedback = NULL;
  c->sinks.swallow = NULL;
  c->sinks.density = NULL;
  c->sinks.do_sink_swallow = NULL;
  c->sinks.do_gas_swallow = NULL;
  c->black_holes.density = NULL;
  c->black_holes.swallow = NULL;
  c->black_holes.do_gas_swallow = NULL;
  c->black_holes.do_bh_swallow = NULL;
  c->black_holes.feedback = NULL;
}

/**
 * @brief Checks that the #part in a cell are at the
 * current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_part_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  /* Only check cells with content */
  if (c->hydro.count == 0) return;

  if (c->hydro.ti_old_part != ti_drift)
    error("Cell in an incorrect time-zone! c->hydro.ti_old=%lld ti_drift=%lld",
          c->hydro.ti_old_part, ti_drift);

  for (int i = 0; i < c->hydro.count; ++i)
    if (c->hydro.parts[i].ti_drift != ti_drift &&
        c->hydro.parts[i].time_bin != time_bin_inhibited)
      error("part in an incorrect time-zone! p->ti_drift=%lld ti_drift=%lld",
            c->hydro.parts[i].ti_drift, ti_drift);

  for (int i = 0; i < c->hydro.count; ++i) {
    const struct part *p = &c->hydro.parts[i];
    if (p->depth_h == c->depth) {
      if (!(p->h >= c->h_min_allowed && p->h < c->h_max_allowed) && c->split) {
        error(
            "depth_h set incorrectly! c->depth=%d p->depth_h=%d h=%e h_min=%e "
            "h_max=%e",
            c->depth, p->depth_h, p->h, c->h_min_allowed, c->h_max_allowed);
      }
    }
  }
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Checks that the #gpart in a cell are at the
 * current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_gpart_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  /* Only check cells with content */
  if (c->grav.count == 0) return;

  if (c->grav.ti_old_part != ti_drift)
    error(
        "Cell in an incorrect time-zone! c->grav.ti_old_part=%lld "
        "ti_drift=%lld",
        c->grav.ti_old_part, ti_drift);

  for (int i = 0; i < c->grav.count; ++i)
    if (c->grav.parts[i].ti_drift != ti_drift &&
        c->grav.parts[i].time_bin != time_bin_inhibited)
      error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
            c->grav.parts[i].ti_drift, ti_drift);
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Checks that the #sink in a cell are at the
 * current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_sink_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  /* Only check cells with content */
  if (c->sinks.count == 0) return;

  if (c->sinks.ti_old_part != ti_drift)
    error(
        "Cell in an incorrect time-zone! c->sinks.ti_old_part=%lld "
        "ti_drift=%lld",
        c->sinks.ti_old_part, ti_drift);

  for (int i = 0; i < c->sinks.count; ++i)
    if (c->sinks.parts[i].ti_drift != ti_drift &&
        c->sinks.parts[i].time_bin != time_bin_inhibited)
      error(
          "sink-part in an incorrect time-zone! sink->ti_drift=%lld "
          "ti_drift=%lld",
          c->sinks.parts[i].ti_drift, ti_drift);

  for (int i = 0; i < c->sinks.count; ++i) {
    const struct sink *sp = &c->sinks.parts[i];
    if (sp->depth_h == c->depth) {
      if (!(sp->h >= c->h_min_allowed && sp->h < c->h_max_allowed) &&
          c->split) {
        error(
            "depth_h set incorrectly! c->depth=%d sp->depth_h=%d h=%e h_min=%e "
            "h_max=%e",
            c->depth, sp->depth_h, sp->h, c->h_min_allowed, c->h_max_allowed);
      }
    }
  }

#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Checks that the #spart in a cell are at the
 * current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_spart_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  /* Only check cells with content */
  if (c->stars.count == 0) return;

  if (c->stars.ti_old_part != ti_drift)
    error(
        "Cell in an incorrect time-zone! c->stars.ti_old_part=%lld "
        "ti_drift=%lld",
        c->stars.ti_old_part, ti_drift);

  for (int i = 0; i < c->stars.count; ++i)
    if (c->stars.parts[i].ti_drift != ti_drift &&
        c->stars.parts[i].time_bin != time_bin_inhibited)
      error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
            c->stars.parts[i].ti_drift, ti_drift);

  for (int i = 0; i < c->stars.count; ++i) {
    const struct spart *p = &c->stars.parts[i];
    if (p->depth_h == c->depth) {
      if (!(p->h >= c->h_min_allowed && p->h < c->h_max_allowed) && c->split) {
        error(
            "depth_h set incorrectly! c->depth=%d p->depth_h=%d h=%e h_min=%e "
            "h_max=%e",
            c->depth, p->depth_h, p->h, c->h_min_allowed, c->h_max_allowed);
      }
    }
  }
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Checks that the #bpart in a cell are at the
 * current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_bpart_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  /* Only check cells with content */
  if (c->black_holes.count == 0) return;

  if (c->black_holes.ti_old_part != ti_drift)
    error(
        "Cell in an incorrect time-zone! c->black_holes.ti_old_part=%lld "
        "ti_drift=%lld",
        c->black_holes.ti_old_part, ti_drift);

  for (int i = 0; i < c->black_holes.count; ++i)
    if (c->black_holes.parts[i].ti_drift != ti_drift &&
        c->black_holes.parts[i].time_bin != time_bin_inhibited)
      error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
            c->black_holes.parts[i].ti_drift, ti_drift);

  for (int i = 0; i < c->black_holes.count; ++i) {
    const struct bpart *bp = &c->black_holes.parts[i];
    if (bp->depth_h == c->depth) {
      if (!(bp->h >= c->h_min_allowed && bp->h < c->h_max_allowed) &&
          c->split) {
        error(
            "depth_h set incorrectly! c->depth=%d p->depth_h=%d h=%e h_min=%e "
            "h_max=%e",
            c->depth, bp->depth_h, bp->h, c->h_min_allowed, c->h_max_allowed);
      }
    }
  }
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Checks that the multipole of a cell is at the current point in time
 *
 * Calls error() if the cell is not at the current time.
 *
 * @param c Cell to act upon
 * @param data The current time on the integer time-line
 */
void cell_check_multipole_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS

  const integertime_t ti_drift = *(integertime_t *)data;

  /* Only check local cells */
  if (c->nodeID != engine_rank) return;

  /* Only check cells with content */
  if (c->grav.count == 0) return;

  if (c->grav.ti_old_multipole != ti_drift)
    error(
        "Cell multipole in an incorrect time-zone! "
        "c->grav.ti_old_multipole=%lld "
        "ti_drift=%lld (depth=%d, node=%d)",
        c->grav.ti_old_multipole, ti_drift, c->depth, c->nodeID);

#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Resets all the individual cell task counters to 0.
 *
 * Should only be used for debugging purposes.
 *
 * @param c The #cell to reset.
 */
void cell_reset_task_counters(struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS
  for (int t = 0; t < task_type_count; ++t) c->tasks_executed[t] = 0;
  for (int t = 0; t < task_subtype_count; ++t) c->subtasks_executed[t] = 0;
  for (int k = 0; k < 8; ++k)
    if (c->progeny[k] != NULL) cell_reset_task_counters(c->progeny[k]);
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Recursively construct all the multipoles in a cell hierarchy.
 *
 * @param c The #cell.
 * @param ti_current The current integer time.
 * @param grav_props The properties of the gravity scheme.
 */
void cell_make_multipoles(struct cell *c, integertime_t ti_current,
                          const struct gravity_props *const grav_props) {

  /* Reset everything */
  gravity_reset(c->grav.multipole);

  if (c->split) {

    /* Start by recursing */
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL)
        cell_make_multipoles(c->progeny[k], ti_current, grav_props);
    }

    /* Compute CoM of all progenies */
    double CoM[3] = {0., 0., 0.};
    double vel[3] = {0., 0., 0.};
    float max_delta_vel[3] = {0.f, 0.f, 0.f};
    float min_delta_vel[3] = {0.f, 0.f, 0.f};
    double mass = 0.;

    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        const struct gravity_tensors *m = c->progeny[k]->grav.multipole;

        mass += m->m_pole.M_000;

        CoM[0] += m->CoM[0] * m->m_pole.M_000;
        CoM[1] += m->CoM[1] * m->m_pole.M_000;
        CoM[2] += m->CoM[2] * m->m_pole.M_000;

        vel[0] += m->m_pole.vel[0] * m->m_pole.M_000;
        vel[1] += m->m_pole.vel[1] * m->m_pole.M_000;
        vel[2] += m->m_pole.vel[2] * m->m_pole.M_000;

        max_delta_vel[0] = max(m->m_pole.max_delta_vel[0], max_delta_vel[0]);
        max_delta_vel[1] = max(m->m_pole.max_delta_vel[1], max_delta_vel[1]);
        max_delta_vel[2] = max(m->m_pole.max_delta_vel[2], max_delta_vel[2]);

        min_delta_vel[0] = min(m->m_pole.min_delta_vel[0], min_delta_vel[0]);
        min_delta_vel[1] = min(m->m_pole.min_delta_vel[1], min_delta_vel[1]);
        min_delta_vel[2] = min(m->m_pole.min_delta_vel[2], min_delta_vel[2]);
      }
    }

    /* Final operation on the CoM and bulk velocity */
    const double mass_inv = 1. / mass;
    c->grav.multipole->CoM[0] = CoM[0] * mass_inv;
    c->grav.multipole->CoM[1] = CoM[1] * mass_inv;
    c->grav.multipole->CoM[2] = CoM[2] * mass_inv;
    c->grav.multipole->m_pole.vel[0] = vel[0] * mass_inv;
    c->grav.multipole->m_pole.vel[1] = vel[1] * mass_inv;
    c->grav.multipole->m_pole.vel[2] = vel[2] * mass_inv;

    /* Min max velocity along each axis */
    c->grav.multipole->m_pole.max_delta_vel[0] = max_delta_vel[0];
    c->grav.multipole->m_pole.max_delta_vel[1] = max_delta_vel[1];
    c->grav.multipole->m_pole.max_delta_vel[2] = max_delta_vel[2];
    c->grav.multipole->m_pole.min_delta_vel[0] = min_delta_vel[0];
    c->grav.multipole->m_pole.min_delta_vel[1] = min_delta_vel[1];
    c->grav.multipole->m_pole.min_delta_vel[2] = min_delta_vel[2];

    /* Now shift progeny multipoles and add them up */
    struct multipole temp;
    double r_max = 0.;
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) {
        const struct cell *cp = c->progeny[k];
        const struct multipole *m = &cp->grav.multipole->m_pole;

        /* Contribution to multipole */
        gravity_M2M(&temp, m, c->grav.multipole->CoM, cp->grav.multipole->CoM);
        gravity_multipole_add(&c->grav.multipole->m_pole, &temp);

        /* Upper limit of max CoM<->gpart distance */
        const double dx =
            c->grav.multipole->CoM[0] - cp->grav.multipole->CoM[0];
        const double dy =
            c->grav.multipole->CoM[1] - cp->grav.multipole->CoM[1];
        const double dz =
            c->grav.multipole->CoM[2] - cp->grav.multipole->CoM[2];
        const double r2 = dx * dx + dy * dy + dz * dz;
        r_max = max(r_max, cp->grav.multipole->r_max + sqrt(r2));
      }
    }
    /* Alternative upper limit of max CoM<->gpart distance */
    const double dx = c->grav.multipole->CoM[0] > c->loc[0] + c->width[0] * 0.5
                          ? c->grav.multipole->CoM[0] - c->loc[0]
                          : c->loc[0] + c->width[0] - c->grav.multipole->CoM[0];
    const double dy = c->grav.multipole->CoM[1] > c->loc[1] + c->width[1] * 0.5
                          ? c->grav.multipole->CoM[1] - c->loc[1]
                          : c->loc[1] + c->width[1] - c->grav.multipole->CoM[1];
    const double dz = c->grav.multipole->CoM[2] > c->loc[2] + c->width[2] * 0.5
                          ? c->grav.multipole->CoM[2] - c->loc[2]
                          : c->loc[2] + c->width[2] - c->grav.multipole->CoM[2];

    /* Take minimum of both limits */
    c->grav.multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));

    /* Compute the multipole power */
    gravity_multipole_compute_power(&c->grav.multipole->m_pole);

  } else {
    if (c->grav.count > 0) {

      gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count, grav_props);

      /* Compute the multipole power */
      gravity_multipole_compute_power(&c->grav.multipole->m_pole);

    } else {

      /* No gparts in that leaf cell */

      /* Set the values to something sensible */
      gravity_multipole_init(&c->grav.multipole->m_pole);
      c->grav.multipole->CoM[0] = c->loc[0] + c->width[0] * 0.5;
      c->grav.multipole->CoM[1] = c->loc[1] + c->width[1] * 0.5;
      c->grav.multipole->CoM[2] = c->loc[2] + c->width[2] * 0.5;
      c->grav.multipole->r_max = 0.;
    }
  }

  /* Also update the values at rebuild time */
  c->grav.multipole->r_max_rebuild = c->grav.multipole->r_max;
  c->grav.multipole->CoM_rebuild[0] = c->grav.multipole->CoM[0];
  c->grav.multipole->CoM_rebuild[1] = c->grav.multipole->CoM[1];
  c->grav.multipole->CoM_rebuild[2] = c->grav.multipole->CoM[2];
  c->grav.multipole->dx_max[0] = 0.f;
  c->grav.multipole->dx_max[1] = 0.f;
  c->grav.multipole->dx_max[2] = 0.f;

  c->grav.ti_old_multipole = ti_current;
}

/**
 * @brief Recursively verify that the multipoles are the sum of their progenies.
 *
 * This function does not check whether the multipoles match the particle
 * content as we may not have received the particles.
 *
 * @param c The #cell to recursively search and verify.
 */
void cell_check_foreign_multipole(const struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS

  if (c->split) {
    double M_000 = 0.;
    long long num_gpart = 0;

    for (int k = 0; k < 8; k++) {
      const struct cell *cp = c->progeny[k];

      if (cp != NULL) {
        /* Check the mass */
        M_000 += cp->grav.multipole->m_pole.M_000;

        /* Check the number of particles */
        num_gpart += cp->grav.multipole->m_pole.num_gpart;

        /* Now recurse */
        cell_check_foreign_multipole(cp);
      }
    }

    if (num_gpart != c->grav.multipole->m_pole.num_gpart)
      error("Sum of particles in progenies does not match");

    if (fabs(M_000 / c->grav.multipole->m_pole.M_000 - 1.) > 1e-2)
      error("Mass in progenies does not match!");
  }

#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Computes the multi-pole brutally and compare to the
 * recursively computed one.
 *
 * @param c Cell to act upon
 * @param grav_props The properties of the gravity scheme.
 */
void cell_check_multipole(struct cell *c,
                          const struct gravity_props *const grav_props) {

#ifdef SWIFT_DEBUG_CHECKS
  struct gravity_tensors ma;
  const double tolerance = 1e-3; /* Relative */

  /* If the cell is a void, exit immediately. We don't want
   * to double count particles in the zoom region which are also in the void
   * cell multipoles. This is because the void cell multipoles are populated
   * bottom up from the zoom cells. The void cell tree itself is tested
   * elsewhere (NOTE: This issue and these cell types only appear when running
   * with a zoom region). */
  if (c->subtype == cell_subtype_void) return;

  /* First recurse */
  if (c->split)
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL)
        cell_check_multipole(c->progeny[k], grav_props);

  if (c->grav.count > 0) {
    /* Brute-force calculation */
    gravity_P2M(&ma, c->grav.parts, c->grav.count, grav_props);
    gravity_multipole_compute_power(&ma.m_pole);

    /* Now  compare the multipole expansion */
    if (!gravity_multipole_equal(&ma, c->grav.multipole, tolerance)) {
      message("Multipoles are not equal at depth=%d! tol=%f", c->depth,
              tolerance);
      message("Correct answer:");
      gravity_multipole_print(&ma.m_pole);
      message("Recursive multipole:");
      gravity_multipole_print(&c->grav.multipole->m_pole);
      error("Aborting");
    }

    /* Check that the upper limit of r_max is good enough */
    if (!(1.1 * c->grav.multipole->r_max >= ma.r_max)) {
      error("Upper-limit r_max=%e too small. Should be >=%e.",
            c->grav.multipole->r_max, ma.r_max);
    } else if (c->grav.multipole->r_max * c->grav.multipole->r_max >
               3. * c->width[0] * c->width[0]) {
      error("r_max=%e larger than cell diagonal %e.", c->grav.multipole->r_max,
            sqrt(3. * c->width[0] * c->width[0]));
    }

    if (fabsf(c->grav.multipole->dx_max[0]) > 0.f ||
        fabsf(c->grav.multipole->dx_max[1]) > 0.f ||
        fabsf(c->grav.multipole->dx_max[2]) > 0.f) {
      error("dx_max should be zero at rebuild: %e %e %e",
            c->grav.multipole->dx_max[0], c->grav.multipole->dx_max[1],
            c->grav.multipole->dx_max[2]);
    }
  }
#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

/**
 * @brief Frees up the memory allocated for this #cell.
 *
 * @param c The #cell.
 */
void cell_clean(struct cell *c) {
  /* Hydro */
  cell_free_hydro_sorts(c);

  /* Stars */
  cell_free_stars_sorts(c);

  /* Grid */
  cell_free_grid(c);

  /* Recurse */
  for (int k = 0; k < 8; k++)
    if (c->progeny[k]) cell_clean(c->progeny[k]);
}

/**
 * @brief Clear the drift flags on the given cell.
 */
void cell_clear_drift_flags(struct cell *c, void *data) {
  cell_clear_flag(c, cell_flag_do_hydro_drift | cell_flag_do_hydro_sub_drift |
                         cell_flag_do_grav_drift | cell_flag_do_grav_sub_drift |
                         cell_flag_do_bh_drift | cell_flag_do_bh_sub_drift |
                         cell_flag_do_stars_drift |
                         cell_flag_do_stars_sub_drift |
                         cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
}

/**
 * @brief Clear the limiter flags on the given cell.
 */
void cell_clear_limiter_flags(struct cell *c, void *data) {
  cell_clear_flag(c,
                  cell_flag_do_hydro_limiter | cell_flag_do_hydro_sub_limiter);
}

/**
 * @brief Set the super-cell pointers for all cells in a hierarchy.
 *
 * @param c The top-level #cell to play with.
 * @param super Pointer to the deepest cell with tasks in this part of the
 * tree.
 * @param with_hydro Are we running with hydrodynamics on?
 * @param with_grav Are we running with gravity on?
 */
void cell_set_super(struct cell *c, struct cell *super, const int with_hydro,
                    const int with_grav) {
  /* Are we in a cell which is either the hydro or gravity super? */
  if (super == NULL && ((with_hydro && c->hydro.super != NULL) ||
                        (with_grav && c->grav.super != NULL)))
    super = c;

  /* Set the super-cell */
  c->super = super;

  /* Recurse */
  if (c->split)
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL)
        cell_set_super(c->progeny[k], super, with_hydro, with_grav);
}

/**
 * @brief Set the super-cell pointers for all cells in a hierarchy.
 *
 * @param c The top-level #cell to play with.
 * @param super_hydro Pointer to the deepest cell with tasks in this part of
 * the tree.
 */
void cell_set_super_hydro(struct cell *c, struct cell *super_hydro) {
  /* Are we in a cell with some kind of self/pair task ? */
  if (super_hydro == NULL && c->hydro.density != NULL) super_hydro = c;

  /* Set the super-cell */
  c->hydro.super = super_hydro;

  /* Recurse */
  if (c->split)
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL)
        cell_set_super_hydro(c->progeny[k], super_hydro);
}

/**
 * @brief Set the super-cell pointers for all cells in a hierarchy.
 *
 * @param c The top-level #cell to play with.
 * @param super_gravity Pointer to the deepest cell with tasks in this part of
 * the tree.
 */
void cell_set_super_gravity(struct cell *c, struct cell *super_gravity) {
  /* Are we in a cell with some kind of self/pair task ? */
  if (super_gravity == NULL && (c->grav.grav != NULL || c->grav.mm != NULL)) {
#ifdef SWIFT_DEBUG_CHECKS
    /* Make sure in zoom land we don't get any confusing empty top level cells
     * with tasks (this breaks hierarchical task creation) */
    if (cell_is_empty_grav(c))
      error(
          "Setting super_gravity to non-void cell at depth %d with no gparts "
          "(%s/%s) "
          "c->grav.grav=%p c->grav.mm=%p, c=%p",
          c->depth, cellID_names[c->type], subcellID_names[c->subtype],
          (void *)c->grav.grav, (void *)c->grav.mm, (void *)c);
#endif

    super_gravity = c;
  }

  /* Set the super-cell */
  c->grav.super = super_gravity;

#ifdef SWIFT_DEBUG_CHECKS
  if (super_gravity != NULL && super_gravity->subtype == cell_subtype_void &&
      c->type == cell_type_zoom)
    error("Zoom cell has a void cell super-gravity pointer!");
#endif

  /* Recurse */
  if (c->split)
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL)
        cell_set_super_gravity(c->progeny[k], super_gravity);
}

/**
 * @brief Mapper function to set the super pointer of the cells.
 *
 * @param map_data The top-level cells.
 * @param num_elements The number of top-level cells.
 * @param extra_data Unused parameter.
 */
void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data) {
  const struct engine *e = (const struct engine *)extra_data;

  const int with_hydro = (e->policy & engine_policy_hydro);
  const int with_grav = (e->policy & engine_policy_self_gravity) ||
                        (e->policy & engine_policy_external_gravity);

  for (int ind = 0; ind < num_elements; ind++) {
    struct cell *c = &((struct cell *)map_data)[ind];

    /* All top-level cells get an MPI tag. */
#ifdef WITH_MPI
    cell_ensure_tagged(c);
#endif

    /* Super-pointer for hydro */
    if (with_hydro) cell_set_super_hydro(c, NULL);

    /* Super-pointer for gravity */
    if (with_grav) cell_set_super_gravity(c, NULL);

    /* Super-pointer for common operations */
    cell_set_super(c, NULL, with_hydro, with_grav);
  }
}

/**
 * @brief Does this cell or any of its children have any task ?
 *
 * We use the timestep-related tasks to probe this as these always
 * exist in a cell hierarchy that has any kind of task.
 *
 * @param c The #cell to probe.
 */
int cell_has_tasks(struct cell *c) {
#ifdef WITH_MPI
  return (c->timestep_collect != NULL || c->mpi.recv != NULL);
#else
  return (c->timestep_collect != NULL);
#endif
}

/**
 * @brief Resets all the sorting properties for the stars in a given cell
 * hierarchy.
 *
 * The clear_unused_flags argument can be used to additionally clean up all
 * the flags demanding a sort for the given cell. This should be used with
 * caution as it will prevent the sort tasks from doing anything on that cell
 * until these flags are reset.
 *
 * @param c The #cell to clean.
 * @param clear_unused_flags Do we also clean the flags demanding a sort?
 */
void cell_clear_stars_sort_flags(struct cell *c, const int clear_unused_flags) {

  /* Clear the flags that have not been reset by the sort task? */
  if (clear_unused_flags) {
    c->stars.requires_sorts = 0;
    c->stars.do_sort = 0;
    cell_clear_flag(c, cell_flag_do_stars_sub_sort);
  }

  /* Indicate that the cell is not sorted and cancel the pointer sorting
   * arrays.
   */
  c->stars.sorted = 0;
  cell_free_stars_sorts(c);

  /* Recurse if possible */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL)
        cell_clear_stars_sort_flags(c->progeny[k], clear_unused_flags);
  }
}

/**
 * @brief Resets all the sorting properties for the hydro in a given cell
 * hierarchy.
 *
 * The clear_unused_flags argument can be used to additionally clean up all
 * the flags demanding a sort for the given cell. This should be used with
 * caution as it will prevent the sort tasks from doing anything on that cell
 * until these flags are reset.
 *
 * @param c The #cell to clean.
 * @param clear_unused_flags Do we also clean the flags demanding a sort?
 */
void cell_clear_hydro_sort_flags(struct cell *c, const int clear_unused_flags) {

  /* Clear the flags that have not been reset by the sort task? */
  if (clear_unused_flags) {
    c->hydro.do_sort = 0;
    c->hydro.requires_sorts = 0;
    cell_clear_flag(c, cell_flag_do_hydro_sub_sort);
  }

  /* Indicate that the cell is not sorted and cancel the pointer sorting
   * arrays.
   */
  c->hydro.sorted = 0;
  cell_free_hydro_sorts(c);

  /* Recurse if possible */
  if (c->split) {
    for (int k = 0; k < 8; k++)
      if (c->progeny[k] != NULL)
        cell_clear_hydro_sort_flags(c->progeny[k], clear_unused_flags);
  }
}

/**
 * @brief Recursively checks that all particles in a cell have a time-step
 */
void cell_check_timesteps(const struct cell *c, const integertime_t ti_current,
                          const timebin_t max_bin) {
#ifdef SWIFT_DEBUG_CHECKS

  if (c->hydro.ti_end_min == 0 && c->grav.ti_end_min == 0 &&
      c->stars.ti_end_min == 0 && c->black_holes.ti_end_min == 0 &&
      c->sinks.ti_end_min == 0 && c->rt.ti_rt_end_min == 0 && c->nr_tasks > 0)
    error("Cell without assigned time-step");

  if (c->split) {
    for (int k = 0; k < 8; ++k)
      if (c->progeny[k] != NULL)
        cell_check_timesteps(c->progeny[k], ti_current, max_bin);
  } else {
    if (c->nodeID == engine_rank)
      for (int i = 0; i < c->hydro.count; ++i)
        if (c->hydro.parts[i].time_bin == 0)
          error("Particle without assigned time-bin");
  }

  /* Other checks not relevent when starting-up */
  if (ti_current == 0) return;

  integertime_t ti_end_min = max_nr_timesteps;
  integertime_t ti_beg_max = 0;

  int count = 0;

  for (int i = 0; i < c->hydro.count; ++i) {

    const struct part *p = &c->hydro.parts[i];
    if (p->time_bin == time_bin_inhibited) continue;
    if (p->time_bin == time_bin_not_created) continue;

    ++count;

    integertime_t ti_end, ti_beg;

    if (p->time_bin <= max_bin) {
      integertime_t time_step = get_integer_timestep(p->time_bin);
      ti_end = get_integer_time_end(ti_current, p->time_bin) + time_step;
      ti_beg = get_integer_time_begin(ti_current + 1, p->time_bin);
    } else {
      ti_end = get_integer_time_end(ti_current, p->time_bin);
      ti_beg = get_integer_time_begin(ti_current + 1, p->time_bin);
    }

    ti_end_min = min(ti_end, ti_end_min);
    ti_beg_max = max(ti_beg, ti_beg_max);
  }

  /* Only check cells that have at least one non-inhibited particle */
  if (count > 0) {

    if (count != c->hydro.count) {

      /* Note that we use a < as the particle with the smallest time-bin
         might have been swallowed. This means we will run this cell with
         0 active particles but that's not wrong */
      if (ti_end_min < c->hydro.ti_end_min)
        error(
            "Non-matching ti_end_min. Cell=%lld true=%lld ti_current=%lld "
            "depth=%d",
            c->hydro.ti_end_min, ti_end_min, ti_current, c->depth);

    } else /* Normal case: nothing was swallowed/converted */ {
      if (ti_end_min != c->hydro.ti_end_min)
        error(
            "Non-matching ti_end_min. Cell=%lld true=%lld ti_current=%lld "
            "depth=%d",
            c->hydro.ti_end_min, ti_end_min, ti_current, c->depth);
    }

    if (ti_beg_max != c->hydro.ti_beg_max)
      error(
          "Non-matching ti_beg_max. Cell=%lld true=%lld ti_current=%lld "
          "depth=%d",
          c->hydro.ti_beg_max, ti_beg_max, ti_current, c->depth);
  }

#else
  error("Calling debugging code without debugging flag activated.");
#endif
}

void cell_check_spart_pos(const struct cell *c,
                          const struct spart *global_sparts) {
#ifdef SWIFT_DEBUG_CHECKS

  /* Recurse */
  if (c->split) {
    for (int k = 0; k < 8; ++k)
      if (c->progeny[k] != NULL)
        cell_check_spart_pos(c->progeny[k], global_sparts);
  }

  struct spart *sparts = c->stars.parts;
  const int count = c->stars.count;
  for (int i = 0; i < count; ++i) {
    const struct spart *sp = &sparts[i];
    if ((sp->x[0] < c->loc[0] / space_stretch) ||
        (sp->x[1] < c->loc[1] / space_stretch) ||
        (sp->x[2] < c->loc[2] / space_stretch) ||
        (sp->x[0] >= (c->loc[0] + c->width[0]) * space_stretch) ||
        (sp->x[1] >= (c->loc[1] + c->width[1]) * space_stretch) ||
        (sp->x[2] >= (c->loc[2] + c->width[2]) * space_stretch))
      error("spart not in its cell!");

    if (sp->time_bin != time_bin_not_created &&
        sp->time_bin != time_bin_inhibited) {
      const struct gpart *gp = sp->gpart;
      if (gp == NULL && sp->time_bin != time_bin_not_created)
        error("Unlinked spart!");

      if (&global_sparts[-gp->id_or_neg_offset] != sp)
        error("Incorrectly linked spart!");
    }
  }

#else
  error("Calling a degugging function outside debugging mode.");
#endif
}

/**
 * @brief Checks that a cell and all its progenies have cleared their sort
 * flags.
 *
 * Should only be used for debugging purposes.
 *
 * @param c The #cell to check.
 */
void cell_check_sort_flags(const struct cell *c) {

#ifdef SWIFT_DEBUG_CHECKS
  const int do_hydro_sub_sort = cell_get_flag(c, cell_flag_do_hydro_sub_sort);
  const int do_stars_sub_sort = cell_get_flag(c, cell_flag_do_stars_sub_sort);

  if (do_hydro_sub_sort)
    error(
        "cell %lld has a hydro sub_sort flag set. Node=%d depth=%d maxdepth=%d",
        c->cellID, c->nodeID, c->depth, c->maxdepth);

  if (do_stars_sub_sort)
    error(
        "cell %lld has a stars sub_sort flag set. Node=%d depth=%d maxdepth=%d",
        c->cellID, c->nodeID, c->depth, c->maxdepth);

  if (c->split) {
    for (int k = 0; k < 8; ++k) {
      if (c->progeny[k] != NULL) cell_check_sort_flags(c->progeny[k]);
    }
  }
#endif
}

/**
 * @brief Can we use the MM interactions fo a given pair of cells?
 *
 * The two cells have to be different!
 *
 * @param ci The first #cell.
 * @param cj The second #cell.
 * @param e The #engine.
 * @param s The #space.
 * @param use_rebuild_data Are we considering the data at the last tree-build
 * (1) or current data (0)?
 * @param is_tree_walk Are we calling this in the tree walk (1) or for the
 * top-level task construction (0)?
 */
int cell_can_use_pair_mm(const struct cell *restrict ci,
                         const struct cell *restrict cj, const struct engine *e,
                         const struct space *s, const int use_rebuild_data,
                         const int is_tree_walk, const int periodic,
                         const int use_mesh) {

  const struct gravity_props *props = e->gravity_properties;
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};

  /* Check for trivial cases */
  if (is_tree_walk && ci->grav.count <= 1) return 0;
  if (is_tree_walk && cj->grav.count <= 1) return 0;

  /* Recover the multipole information */
  const struct gravity_tensors *restrict multi_i = ci->grav.multipole;
  const struct gravity_tensors *restrict multi_j = cj->grav.multipole;

  const double r2 =
      cell_mpole_CoM_dist2(multi_i, multi_j, use_rebuild_data, periodic, dim);

  return gravity_M2L_accept_symmetric(props, multi_i, multi_j, r2,
                                      use_rebuild_data, use_mesh);
}

/**
 * @brief Test whether two cells can use PM interactions.
 *
 * This will test if particles in the two cells are far enough apart to use
 * the mesh for their interaction. If so we won't need an expensive pair
 * task or multipole-multipole interaction.
 *
 * @param e The #engine.
 * @param ci The first #cell.
 * @param cj The second #cell.
 *
 * @return 1 if the mesh can be used, 0 otherwise.
 */
int cell_can_use_mesh(struct engine *e, const struct cell *ci,
                      const struct cell *cj) {

  struct space *s = e->s;
  const double max_distance = e->mesh->r_cut_max;
  const double max_distance2 = max_distance * max_distance;

  /* If not periodic then we cannot use the mesh */
  if (!s->periodic) {
    return 0;
  }

  /* Minimal distance between any pair of particles */
  const double min_radius2 = cell_min_dist2(ci, cj, s->periodic, s->dim);

  /* Are we beyond the distance where the truncated forces are 0 ?*/
  return (min_radius2 > max_distance2);
}

/**
 * @brief Do we need a rebuild based on whether two cells can still use PM
 * interactions despite motion since the last rebuild.
 *
 * This will test if particles in the two cells are far enough apart to use
 * the mesh with and without including the maximal displacement of a gpart since
 * the last rebuild. At the rebuild we made tasks based on the particle which
 * are wholly inside cells.
 *
 * @param e The #engine.
 * @param ci The first #cell.
 * @param cj The second #cell.
 * @return 1 if we could use the mesh and now can't, 0 otherwise.
 */
int cell_cant_use_mesh_anymore(struct engine *e, const struct cell *ci,
                               const struct cell *cj) {

  struct space *s = e->s;
  const double max_distance = e->mesh->r_cut_max;
  const double max_distance2 = max_distance * max_distance;

  /* If not periodic then we cannot use the mesh */
  if (!s->periodic) {
    return 0;
  }

  /* Could we use the mesh at rebuild time? */
  const int could_use_mesh_at_rebuild = cell_can_use_mesh(e, ci, cj);

  /* If we could not use the mesh at rebuild time then no one cares */
  if (!could_use_mesh_at_rebuild) {
    return 0;
  }

  /*Minimal distance between any pair of particles including max
   * displacement since rebuild */
  const double min_radius2 =
      cell_min_dist2_with_max_dx(ci, cj, s->periodic, s->dim);

  /* Are we beyond the distance where the truncated forces are 0 ?*/
  const int can_use_mesh_now = (min_radius2 > max_distance2);

  /* Report related information if we could use the mesh at rebuild but can no
   * longer do so */
  if (could_use_mesh_at_rebuild && !can_use_mesh_now) {
    message(
        "Forcing rebuild due to particle motion. Cell pair: "
        "min_radius2=%e ci->grav.multipole->r_max=%e "
        "ci->grav.multipole->r_max_rebuild=%e ci->grav.multipole->dx_max=%e %e "
        "%e "
        "cj->grav.multipole->r_max=%e cj->grav.multipole->r_max_rebuild=%e "
        "cj->grav.multipole->dx_max=%e %e %e, ci->grav.count=%d "
        "cj->grav.count=%d",
        min_radius2, ci->grav.multipole->r_max,
        ci->grav.multipole->r_max_rebuild, ci->grav.multipole->dx_max[0],
        ci->grav.multipole->dx_max[1], ci->grav.multipole->dx_max[2],
        cj->grav.multipole->r_max, cj->grav.multipole->r_max_rebuild,
        cj->grav.multipole->dx_max[0], cj->grav.multipole->dx_max[1],
        cj->grav.multipole->dx_max[2], ci->grav.count, cj->grav.count);
  }

  return (could_use_mesh_at_rebuild && !can_use_mesh_now);
}

/**
 * @brief Recursively check if gravity mesh pairs can still be skipped.
 *
 * This function recursively checks pairs of cells (or a single cell for
 * self-interactions) to determine if they could use mesh at rebuild but
 * can no longer do so due to particle motion. If such a pair is found,
 * it sets e->forcerebuild and returns early.
 *
 * @param ci The first #cell.
 * @param cj The second #cell (can be NULL for self-interactions).
 * @param e The #engine.
 * @return 1 if forcerebuild was set, 0 otherwise.
 */
static int cell_check_grav_mesh_pairs_recursive(struct cell *ci,
                                                struct cell *cj,
                                                struct engine *e) {

  struct space *s = e->s;

  /* Self-interaction? */
  if (cj == NULL) {

    /* Should this self task be split?
     * Matches scheduler_splittask_gravity logic */
    if (cell_can_split_self_gravity_task(ci)) {

      /* Check particle count threshold - matches scheduler_splittask_gravity */
      if (ci->grav.count < space_subsize_self_grav) {
        return 0;
      }

      /* Recurse on self interactions for each progeny */
      for (int k = 0; k < 8; k++) {
        if (ci->progeny[k] != NULL) {
          if (cell_check_grav_mesh_pairs_recursive(ci->progeny[k], NULL, e)) {
            return 1;
          }
        }
      }

      /* Now handle pair interactions between progeny */
      for (int j = 0; j < 8; j++) {
        if (ci->progeny[j] == NULL) continue;
        struct cell *cpj = ci->progeny[j];
        for (int k = j + 1; k < 8; k++) {
          if (ci->progeny[k] == NULL) continue;
          struct cell *cpk = ci->progeny[k];

          /* Can we use the mesh for this pair? */
          if (cell_can_use_mesh(e, cpj, cpk)) {
            /* Check if we can no longer use the mesh */
            if (cell_cant_use_mesh_anymore(e, cpj, cpk)) {
              atomic_inc(&e->forcerebuild);
              return 1;
            }
            /* Mesh still valid, continue */
            continue;
          }

          /* Otherwise recurse as a pair interaction */
          if (cell_check_grav_mesh_pairs_recursive(cpj, cpk, e)) {
            return 1;
          }
        }
      }
    }
    /* else: We have a real task that doesn't split further */

  } else {
    /* Pair interaction */

    /* Should this pair be split?
     * Matches scheduler_splittask_gravity logic */
    if (cell_can_split_pair_gravity_task(ci, cj)) {

      /* Check particle count threshold - matches scheduler_splittask_gravity */
      const long long gcount_i = ci->grav.count;
      const long long gcount_j = cj->grav.count;
      if (gcount_i * gcount_j < ((long long)space_subsize_pair_grav)) {
        return 0;
      }

      /* Recurse on all progeny pairs */
      for (int i = 0; i < 8; i++) {
        if (ci->progeny[i] == NULL) continue;
        struct cell *cpi = ci->progeny[i];
        for (int j = 0; j < 8; j++) {
          if (cj->progeny[j] == NULL) continue;
          struct cell *cpj = cj->progeny[j];

          /* Can we use the mesh for this pair? */
          if (cell_can_use_mesh(e, cpi, cpj)) {
            /* Check if we can no longer use the mesh */
            if (cell_cant_use_mesh_anymore(e, cpi, cpj)) {
              atomic_inc(&e->forcerebuild);
              return 1;
            }
            /* Mesh still valid, continue */
            continue;
          }

          /* Can we use M-M for this pair? */
          if (cell_can_use_pair_mm(cpi, cpj, e, s, /*use_rebuild_data=*/1,
                                   /*is_tree_walk=*/1, s->periodic,
                                   /*use_mesh=*/s->periodic)) {
            /* This would be handled by a M-M task, skip */
            continue;
          }

          /* We would create real tasks, so recurse */
          if (cell_check_grav_mesh_pairs_recursive(cpi, cpj, e)) {
            return 1;
          }
        }
      }
    }
    /* else: We have a real task that doesn't split further */
  }

  return 0;
}

/**
 * @brief Check if gravity mesh pairs from a top-level cell can still be
 * skipped.
 *
 * This function finds all top-level pairs involving the given cell (using
 * the same logic as engine_make_self_gravity_tasks_mapper) and recursively
 * checks if pairs that were skipped due to the mesh at rebuild can still
 * be skipped. If any pair can no longer use the mesh, a rebuild is triggered.
 *
 * @param c The top-level #cell.
 * @param e The #engine.
 */
void cell_check_grav_mesh_pairs(struct cell *c, struct engine *e) {

  struct space *s = e->s;
  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
  struct cell *cells = s->cells_top;

  /* Skip cells without gravity particles */
  if (c->grav.count == 0) return;

  /* Early exit if not doing self-gravity */
  if (!(e->policy & engine_policy_self_gravity)) {
    return;
  }

  /* Early exit if not using mesh */
  if (!s->periodic) {
    return;
  }

  /* Early exit if forcerebuild already set */
  if (e->forcerebuild) {
    return;
  }

  /* Compute maximal distance where we can expect a direct interaction */
  const float distance = gravity_M2L_min_accept_distance(
      e->gravity_properties, sqrtf(3) * cells[0].width[0], s->max_softening,
      s->min_a_grav, s->max_mpole_power, /*periodic=*/1);

  /* Convert the maximal search distance to a number of cells */
  const int delta = max((int)(sqrt(3) * distance / cells[0].width[0]) + 1, 2);
  int delta_m = delta;
  int delta_p = delta;

  /* Special case where every cell is in range of every other one */
  if (delta >= cdim[0] / 2) {
    if (cdim[0] % 2 == 0) {
      delta_m = cdim[0] / 2;
      delta_p = cdim[0] / 2 - 1;
    } else {
      delta_m = cdim[0] / 2;
      delta_p = cdim[0] / 2;
    }
  }

  /* Get the cell index and integer indices in the top-level grid */
  const int cid = c - cells;
  const int i = cid / (cdim[1] * cdim[2]);
  const int j = (cid / cdim[2]) % cdim[1];
  const int k = cid % cdim[2];

  /* Check self-interaction. If it returns true (forcerebuild set), we don't
   * need to check any pairs. */
  if (cell_check_grav_mesh_pairs_recursive(c, NULL, e)) {
    message("Self interaction triggers a rebuild due to mesh pair failure.");
    return;
  }

  /* Loop over every other cell within (Manhattan) range delta */
  for (int ii = i - delta_m; ii <= i + delta_p; ii++) {
    for (int jj = j - delta_m; jj <= j + delta_p; jj++) {
      for (int kk = k - delta_m; kk <= k + delta_p; kk++) {

        /* Apply periodic BC */
        const int iii = (ii + cdim[0]) % cdim[0];
        const int jjj = (jj + cdim[1]) % cdim[1];
        const int kkk = (kk + cdim[2]) % cdim[2];

        /* Get the second cell */
        const int cjd = cell_getid(cdim, iii, jjj, kkk);
        struct cell *cj = &cells[cjd];

        /* Avoid duplicates, empty cells */
        if (cid >= cjd || cj->grav.count == 0) continue;

        /* Check this pair recursively */
        if (cell_check_grav_mesh_pairs_recursive(c, cj, e)) {
          message(
              "Pair interaction triggers a rebuild due to mesh pair failure");
          return;
        }
      }
    }
  }
}
