/*******************************************************************************
 * 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)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

/* Config parameters. */
#include <config.h>

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

/* Local headers. */
#include "active.h"
#include "adaptive_softening.h"
#include "black_holes.h"
#include "cell.h"
#include "engine.h"
#include "feedback.h"
#include "mhd.h"
#include "rt.h"
#include "sink.h"
#include "space_getsid.h"
#include "star_formation.h"
#include "stars.h"
#include "timers.h"
#include "timestep_limiter.h"
#include "tracers.h"

/* Import the density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_hydro.h"
#undef FUNCTION
#undef FUNCTION_TASK_LOOP

/* Import the stars density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_stars.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

/* Import the black hole density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_black_holes.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

/* Import the sink density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_sinks.h"
#undef FUNCTION_TASK_LOOP
#undef FUNCTION

/**
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param offset First particle in the cell to treat (for split tasks).
 * @param ntasks Interval between successive particles that are treated.
 * @param timer Are we timing this ?
 */
void runner_do_stars_ghost(struct runner *r, struct cell *c, const int offset,
                           const int ntasks, const int timer) {

  struct spart *restrict sparts = c->stars.parts;
  const struct engine *e = r->e;
  const struct unit_system *us = e->internal_units;
  const struct phys_const *phys_const = e->physical_constants;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
  const int with_rt = (e->policy & engine_policy_rt);
  const struct cosmology *cosmo = e->cosmology;
  const struct feedback_props *feedback_props = e->feedback_props;
  const struct rt_props *rt_props = e->rt_props;
  const float stars_h_max = e->hydro_properties->h_max;
  const float stars_h_min = e->hydro_properties->h_min;
  const float eps = e->stars_properties->h_tolerance;
  const float stars_eta_dim =
      pow_dimension(e->stars_properties->eta_neighbours);
  const int max_smoothing_iter = e->stars_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;

  /* Running value of the maximal smoothing length */
  float h_max = c->stars.h_max;
  float h_max_active = 0.f;

  TIMER_TIC;

#ifdef SWIFT_DEBUG_CHECKS
  if (c->nodeID != e->nodeID)
    error("Running the star ghost on a foreign node!");
#endif

  /* Anything to do here? */
  if (c->stars.count == 0) return;
  if (!cell_is_active_stars(c, e)) return;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        runner_do_stars_ghost(r, c->progeny[k], offset, ntasks, /*timer=*/0);

        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->stars.h_max);
        h_max_active = max(h_max_active, c->progeny[k]->stars.h_max_active);
      }
    }
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
    if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
      error("Can't allocate memory for sid.");
    if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
      error("Can't allocate memory for h_0.");
    if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
      error("Can't allocate memory for left.");
    if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
      error("Can't allocate memory for right.");
    for (int k = offset; k < c->stars.count; k += ntasks)
      if (spart_is_active(&sparts[k], e) &&
          (feedback_is_active(&sparts[k], e) || with_rt)) {
        sid[scount] = k;
        h_0[scount] = sparts[k].h;
        left[scount] = 0.f;
        right[scount] = stars_h_max;
        ++scount;
      }

    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {

      ghost_stats_account_for_stars(&c->ghost_statistics, num_reruns, scount,
                                    sparts, sid);

      /* Reset the redo-count. */
      redo = 0;

      /* Loop over the remaining active parts in this cell. */
      for (int i = 0; i < scount; i++) {

        /* Get a direct pointer on the part. */
        struct spart *sp = &sparts[sid[i]];

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!spart_is_active(sp, e))
          error("Ghost applied to inactive particle");
        if (!feedback_is_active(sp, e) && !with_rt)
          error("Ghost applied to particle inactive for feedback and RT");
#endif

        /* Get some useful values */
        const float h_init = h_0[i];
        const float h_old = sp->h;
        const float h_old_dim = pow_dimension(h_old);
        const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);

        float h_new;
        int has_no_neighbours = 0;

        if (sp->density.wcount < 1.e-5 * kernel_root) { /* No neighbours case */

          ghost_stats_no_ngb_star_iteration(&c->ghost_statistics, num_reruns);

          /* Flag that there were no neighbours */
          has_no_neighbours = 1;

          /* Double h and try again */
          h_new = 2.f * h_old;

        } else {

          /* Finish the density calculation */
          stars_end_density(sp, cosmo);

          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = sp->density.wcount * h_old_dim;
          const float n_target = stars_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              sp->density.wcount_dh * h_old_dim +
              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;

          /* Improve the bisection bounds */
          if (n_sum < n_target)
            left[i] = max(left[i], h_old);
          else if (n_sum > n_target)
            right[i] = min(right[i], h_old);

#ifdef SWIFT_DEBUG_CHECKS
          /* Check the validity of the left and right bounds */
          if (left[i] > right[i])
            error("Invalid left (%e) and right (%e)", left[i], right[i]);
#endif

          /* Skip if h is already h_max and we don't have enough neighbours
           */
          /* Same if we are below h_min */
          if (((sp->h >= stars_h_max) && (f < 0.f)) ||
              ((sp->h <= stars_h_min) && (f > 0.f))) {

            stars_reset_feedback(sp);

            /* Only do feedback if stars have a reasonable birth time */
            if (feedback_is_active(sp, e)) {

              const integertime_t ti_step = get_integer_timestep(sp->time_bin);
              const integertime_t ti_begin =
                  get_integer_time_begin(e->ti_current - 1, sp->time_bin);

              /* Get particle time-step */
              double dt_star;
              if (with_cosmology) {
                dt_star = cosmology_get_delta_time(e->cosmology, ti_begin,
                                                   ti_begin + ti_step);
              } else {
                dt_star = get_timestep(sp->time_bin, e->time_base);
              }

              /* Calculate age of the star at current time */
              const double star_age_end_of_step =
                  stars_compute_age(sp, e->cosmology, e->time, with_cosmology);

              /* Has this star been around for a while ? */
              if (star_age_end_of_step > 0.) {

                /* Get the length of the enrichment time-step */
                const double dt_enrichment = feedback_get_enrichment_timestep(
                    sp, with_cosmology, cosmo, e->time, dt_star);
                const double star_age_beg_of_step =
                    star_age_end_of_step - dt_enrichment;

                /* Compute the stellar evolution  */
                feedback_prepare_feedback(sp, feedback_props, cosmo, us,
                                          phys_const, star_age_beg_of_step,
                                          dt_enrichment, e->time, ti_begin,
                                          with_cosmology);
              } else {

                /* Reset the feedback fields of the star particle */
                feedback_reset_feedback(sp, feedback_props);
              }
            } else {

              feedback_reset_feedback(sp, feedback_props);
            }

            if (with_rt) {

              rt_reset_spart(sp);

              /* Get particle time-step */
              double dt_star;
              if (with_cosmology) {

                /* get star's age and time step for stellar emission rates */
                const integertime_t ti_begin =
                    get_integer_time_begin(e->ti_current - 1, sp->time_bin);
                const integertime_t ti_step =
                    get_integer_timestep(sp->time_bin);
                dt_star = cosmology_get_delta_time(e->cosmology, ti_begin,
                                                   ti_begin + ti_step);

              } else {
                dt_star = get_timestep(sp->time_bin, e->time_base);
              }

              /* Calculate age of the star at current time */
              const double star_age_end_of_step =
                  stars_compute_age(sp, e->cosmology, e->time, with_cosmology);

              rt_compute_stellar_emission_rate(sp, e->time,
                                               star_age_end_of_step, dt_star,
                                               rt_props, phys_const, us);
            }

            if (feedback_is_active(sp, e) || with_rt) {

              /* Check if h_max has increased */
              h_max = max(h_max, sp->h);
              h_max_active = max(h_max_active, sp->h);
            }

            /* Ok, we are done with this particle */
            continue;
          }

          /* Normal case: Use Newton-Raphson to get a better value of h */

          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);

          /* Be verbose about the particles that struggle to converge */
          if (num_reruns > max_smoothing_iter - 10) {

            message(
                "Smoothing length convergence problem: iter=%d p->id=%lld "
                "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
                "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
                num_reruns, sp->id, h_init, h_old, h_new, f, f_prime, n_sum,
                n_target, left[i], right[i]);
          }

          /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
          h_new = min(h_new, 2.f * h_old);
          h_new = max(h_new, 0.5f * h_old);

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
        }

        /* Check whether the particle has an inappropriate smoothing length
         */
        if (fabsf(h_new - h_old) > eps * h_old) {

          /* Ok, correct then */

          /* Case where we have been oscillating around the solution */
          if ((h_new == left[i] && h_old == right[i]) ||
              (h_old == left[i] && h_new == right[i])) {

            /* Bisect the remaining interval */
            sp->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            sp->h = h_new;
          }

          /* If below the absolute maximum, try again */
          if (sp->h < stars_h_max && sp->h > stars_h_min) {

            /* Flag for another round of fun */
            sid[redo] = sid[i];
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
            redo += 1;

            /* Re-initialise everything */
            stars_init_spart(sp);
            feedback_init_spart(sp);
            rt_init_spart(sp);

            /* Off we go ! */
            continue;

          } else if (sp->h <= stars_h_min) {

            /* Ok, this particle is a lost cause... */
            sp->h = stars_h_min;

          } else if (sp->h >= stars_h_max) {

            /* Ok, this particle is a lost cause... */
            sp->h = stars_h_max;

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
              ghost_stats_no_ngb_star_converged(&c->ghost_statistics);
              stars_spart_has_no_neighbours(sp, cosmo);
              rt_spart_has_no_neighbours(sp);
            }

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
          }
        }

        /* We now have a particle whose smoothing length has converged */

        /* Set the correct depth */
        cell_set_spart_h_depth(sp, c);

        /* Check if h_max has increased */
        h_max = max(h_max, sp->h);
        h_max_active = max(h_max_active, sp->h);

        ghost_stats_converged_star(&c->ghost_statistics, sp);

        stars_reset_feedback(sp);

        /* Only do feedback if stars have a reasonable birth time */
        if (feedback_is_active(sp, e)) {

          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, sp->time_bin);

          /* Get particle time-step */
          double dt_star;
          if (with_cosmology) {
            dt_star = cosmology_get_delta_time(e->cosmology, ti_begin,
                                               ti_begin + ti_step);
          } else {
            dt_star = get_timestep(sp->time_bin, e->time_base);
          }

          /* Calculate age of the star at current time */
          const double star_age_end_of_step =
              stars_compute_age(sp, e->cosmology, e->time, with_cosmology);

          /* Has this star been around for a while ? */
          if (star_age_end_of_step > 0.) {

            /* Get the length of the enrichment time-step */
            const double dt_enrichment = feedback_get_enrichment_timestep(
                sp, with_cosmology, cosmo, e->time, dt_star);
            const double star_age_beg_of_step =
                star_age_end_of_step - dt_enrichment;

            /* Compute the stellar evolution  */
            feedback_prepare_feedback(sp, feedback_props, cosmo, us, phys_const,
                                      star_age_beg_of_step, dt_enrichment,
                                      e->time, ti_begin, with_cosmology);
          } else {

            /* Reset the feedback fields of the star particle */
            feedback_reset_feedback(sp, feedback_props);
          }
        } else {

          /* Reset the feedback fields of the star particle */
          feedback_reset_feedback(sp, feedback_props);
        }

        if (with_rt) {

          rt_reset_spart(sp);

          /* Get particle time-step */
          double dt_star;
          if (with_cosmology) {

            /* get star's age and time step for stellar emission rates */
            const integertime_t ti_begin =
                get_integer_time_begin(e->ti_current - 1, sp->time_bin);
            const integertime_t ti_step = get_integer_timestep(sp->time_bin);
            dt_star = cosmology_get_delta_time(e->cosmology, ti_begin,
                                               ti_begin + ti_step);

          } else {
            dt_star = get_timestep(sp->time_bin, e->time_base);
          }

          /* Calculate age of the star at current time */
          const double star_age_end_of_step =
              stars_compute_age(sp, e->cosmology, e->time, with_cosmology);

          rt_compute_stellar_emission_rate(sp, e->time, star_age_end_of_step,
                                           dt_star, rt_props, phys_const, us);
        }
      }

      /* We now need to treat the particles whose smoothing length had not
       * converged again */

      /* Re-set the counter for the next loop (potentially). */
      scount = redo;
      if (scount > 0) {

        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {

          /* Run through this cell's density interactions. */
          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {

#ifdef SWIFT_DEBUG_CHECKS
            if (l->t->ti_run < r->e->ti_current)
              error("Density task should have been run.");
#endif

            /* Self interaction? */
            if (l->t->type == task_type_self) {

              /* if we have split the self tasks, then we only want to call
               *  doself_subset once, for the split task with flag 0 */
#if STARS_SELF_NTASK > 1
              const int run_on_all_selfs = 0;
#else
              const int run_on_all_selfs = 1;
#endif
              if (run_on_all_selfs || l->t->flags == 0) {

                runner_dosub_self_subset_stars_density(r, finger, sparts, sid,
                                                       scount, 1);
              }
            }

            /* Otherwise, pair interaction? */
            else if (l->t->type == task_type_pair) {

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dosub_pair_subset_stars_density(r, finger, sparts, sid,
                                                       scount, l->t->cj, 1);
              else
                runner_dosub_pair_subset_stars_density(r, finger, sparts, sid,
                                                       scount, l->t->ci, 1);
            } else {
#ifdef SWIFT_DEBUG_CHECKS
              error("Invalid sub-type!");
#endif
            }
          }
        }
      }
    }

    if (scount) {
      warning(
          "Smoothing length failed to converge for the following star "
          "particles:");
      for (int i = 0; i < scount; i++) {
        struct spart *sp = &sparts[sid[i]];
        warning("ID: %lld, h: %g, wcount: %g", sp->id, sp->h,
                sp->density.wcount);
      }

      error("Smoothing length failed to converge on %i particles.", scount);
    }

    /* Be clean */
    free(left);
    free(right);
    free(sid);
    free(h_0);

#ifdef SWIFT_DEBUG_CHECKS
    /* Update h_max */
    atomic_max_f(&c->stars.h_max, h_max);
    atomic_max_f(&c->stars.h_max_active, h_max_active);

    for (int i = offset; i < c->stars.count; i += ntasks) {
      const struct spart *sp = &c->stars.parts[i];
      const float h = c->stars.parts[i].h;
      if (spart_is_inhibited(sp, e)) continue;

      if (h > c->stars.h_max)
        error("Particle has h larger than h_max (id=%lld)", sp->id);

      if (spart_is_active(sp, e) && (feedback_is_active(sp, e) || with_rt) &&
          (h > c->stars.h_max_active))
        error("Active particle has h larger than h_max_active (id=%lld)",
              sp->id);
    }
#endif
  }

  /* Update h_max */
  atomic_max_f(&c->stars.h_max, h_max);
  atomic_max_f(&c->stars.h_max_active, h_max_active);

  /* The ghost may not always be at the top level.
   * Therefore we need to update h_max between the super- and top-levels
   * Note: We can check the 0th entry as they are either all NULL or none */
  if (c->stars.density_ghost[0]) {
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
      atomic_max_f(&tmp->stars.h_max, h_max);
      atomic_max_f(&tmp->stars.h_max_active, h_max_active);
    }
  }

  if (timer) TIMER_TOC(timer_do_stars_ghost);
}

/**
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
                                         int timer) {

  struct bpart *restrict bparts = c->black_holes.parts;
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
  const float black_holes_h_max = e->hydro_properties->h_max;
  const float black_holes_h_min = e->hydro_properties->h_min;
  const float eps = e->black_holes_properties->h_tolerance;
  const float black_holes_eta_dim =
      pow_dimension(e->black_holes_properties->eta_neighbours);
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
  int redo = 0, bcount = 0;

  /* Running value of the maximal smoothing length */
  float h_max = c->black_holes.h_max;
  float h_max_active = 0.f;

  TIMER_TIC;

  /* Anything to do here? */
  if (c->black_holes.count == 0) return;
  if (!cell_is_active_black_holes(c, e)) return;

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

        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->black_holes.h_max);
        h_max_active =
            max(h_max_active, c->progeny[k]->black_holes.h_max_active);
      }
    }
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
    if ((sid = (int *)malloc(sizeof(int) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for sid.");
    if ((h_0 = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for h_0.");
    if ((left = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for left.");
    if ((right = (float *)malloc(sizeof(float) * c->black_holes.count)) == NULL)
      error("Can't allocate memory for right.");
    for (int k = 0; k < c->black_holes.count; k++)
      if (bpart_is_active(&bparts[k], e)) {
        sid[bcount] = k;
        h_0[bcount] = bparts[k].h;
        left[bcount] = 0.f;
        right[bcount] = black_holes_h_max;
        ++bcount;
      }

    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; bcount > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {

      ghost_stats_account_for_black_holes(&c->ghost_statistics, num_reruns,
                                          bcount, bparts, sid);

      /* Reset the redo-count. */
      redo = 0;

      /* Loop over the remaining active parts in this cell. */
      for (int i = 0; i < bcount; i++) {

        /* Get a direct pointer on the part. */
        struct bpart *bp = &bparts[sid[i]];

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!bpart_is_active(bp, e))
          error("Ghost applied to inactive particle");
#endif

        /* Get some useful values */
        const float h_init = h_0[i];
        const float h_old = bp->h;
        const float h_old_dim = pow_dimension(h_old);
        const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);

        float h_new;
        int has_no_neighbours = 0;

        if (bp->density.wcount < 1.e-5 * kernel_root) { /* No neighbours case */

          ghost_stats_no_ngb_black_hole_iteration(&c->ghost_statistics,
                                                  num_reruns);

          /* Flag that there were no neighbours */
          has_no_neighbours = 1;

          /* Double h and try again */
          h_new = 2.f * h_old;

        } else {

          /* Finish the density calculation */
          black_holes_end_density(bp, cosmo);

          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = bp->density.wcount * h_old_dim;
          const float n_target = black_holes_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              bp->density.wcount_dh * h_old_dim +
              hydro_dimension * bp->density.wcount * h_old_dim_minus_one;

          /* Improve the bisection bounds */
          if (n_sum < n_target)
            left[i] = max(left[i], h_old);
          else if (n_sum > n_target)
            right[i] = min(right[i], h_old);

#ifdef SWIFT_DEBUG_CHECKS
          /* Check the validity of the left and right bounds */
          if (left[i] > right[i])
            error("Invalid left (%e) and right (%e)", left[i], right[i]);
#endif

          /* Skip if h is already h_max and we don't have enough neighbours
           */
          /* Same if we are below h_min */
          if (((bp->h >= black_holes_h_max) && (f < 0.f)) ||
              ((bp->h <= black_holes_h_min) && (f > 0.f))) {

            black_holes_reset_feedback(bp);

            /* Check if h_max has increased */
            h_max = max(h_max, bp->h);
            h_max_active = max(h_max_active, bp->h);

            /* Ok, we are done with this particle */
            continue;
          }

          /* Normal case: Use Newton-Raphson to get a better value of h */

          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);

          /* Be verbose about the particles that struggle to converge */
          if (num_reruns > max_smoothing_iter - 10) {

            message(
                "Smoothing length convergence problem: iter=%d p->id=%lld "
                "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
                "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
                num_reruns, bp->id, h_init, h_old, h_new, f, f_prime, n_sum,
                n_target, left[i], right[i]);
          }

          /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
          h_new = min(h_new, 2.f * h_old);
          h_new = max(h_new, 0.5f * h_old);

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
        }

        /* Check whether the particle has an inappropriate smoothing length
         */
        if (fabsf(h_new - h_old) > eps * h_old) {

          /* Ok, correct then */

          /* Case where we have been oscillating around the solution */
          if ((h_new == left[i] && h_old == right[i]) ||
              (h_old == left[i] && h_new == right[i])) {

            /* Bisect the remaining interval */
            bp->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            bp->h = h_new;
          }

          /* If below the absolute maximum, try again */
          if (bp->h < black_holes_h_max && bp->h > black_holes_h_min) {

            /* Flag for another round of fun */
            sid[redo] = sid[i];
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
            redo += 1;

            /* Re-initialise everything */
            black_holes_init_bpart(bp);

            /* Off we go ! */
            continue;

          } else if (bp->h <= black_holes_h_min) {

            /* Ok, this particle is a lost cause... */
            bp->h = black_holes_h_min;

          } else if (bp->h >= black_holes_h_max) {

            /* Ok, this particle is a lost cause... */
            bp->h = black_holes_h_max;

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
              ghost_stats_no_ngb_black_hole_converged(&c->ghost_statistics);
              black_holes_bpart_has_no_neighbours(bp, cosmo);
            }

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
          }
        }

        /* We now have a particle whose smoothing length has converged */

        ghost_stats_converged_black_hole(&c->ghost_statistics, bp);

        black_holes_reset_feedback(bp);

        /* Set the correct depth */
        cell_set_bpart_h_depth(bp, c);

        /* Check if h_max has increased */
        h_max = max(h_max, bp->h);
        h_max_active = max(h_max_active, bp->h);
      }

      /* We now need to treat the particles whose smoothing length had not
       * converged again */

      /* Re-set the counter for the next loop (potentially). */
      bcount = redo;
      if (bcount > 0) {

        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {

          /* Run through this cell's density interactions. */
          for (struct link *l = finger->black_holes.density; l != NULL;
               l = l->next) {

#ifdef SWIFT_DEBUG_CHECKS
            if (l->t->ti_run < r->e->ti_current)
              error("Density task should have been run.");
#endif

            /* Self-interaction? */
            if (l->t->type == task_type_self) {
              runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
                                             NULL, 1);
            }

            /* Otherwise, pair interaction? */
            else if (l->t->type == task_type_pair) {

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
                                               l->t->cj, 1);
              else
                runner_dosub_subset_bh_density(r, finger, bparts, sid, bcount,
                                               l->t->ci, 1);
            } else {
#ifdef SWIFT_DEBUG_CHECKS
              error("Invalid sub-type!");
#endif
            }
          }
        }
      }
    }

    if (bcount) {
      warning(
          "Smoothing length failed to converge for the following BH "
          "particles:");
      for (int i = 0; i < bcount; i++) {
        struct bpart *bp = &bparts[sid[i]];
        warning("ID: %lld, h: %g, wcount: %g", bp->id, bp->h,
                bp->density.wcount);
      }

      error("Smoothing length failed to converge on %i particles.", bcount);
    }

    /* Be clean */
    free(left);
    free(right);
    free(sid);
    free(h_0);
  }

  /* Update h_max */
  c->black_holes.h_max = h_max;
  c->black_holes.h_max_active = h_max_active;

#ifdef SWIFT_DEBUG_CHECKS
  for (int i = 0; i < c->black_holes.count; ++i) {
    const struct bpart *bp = &c->black_holes.parts[i];
    const float h = c->black_holes.parts[i].h;
    if (bpart_is_inhibited(bp, e)) continue;

    if (h > c->black_holes.h_max)
      error("Particle has h larger than h_max (id=%lld)", bp->id);
    if (bpart_is_active(bp, e) && h > c->black_holes.h_max_active)
      error("Active particle has h larger than h_max_active (id=%lld)", bp->id);
  }
#endif

  /* The ghost may not always be at the top level.
   * Therefore we need to update h_max between the super- and top-levels */
  if (c->black_holes.density_ghost) {
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
      atomic_max_f(&tmp->black_holes.h_max, h_max);
      atomic_max_f(&tmp->black_holes.h_max_active, h_max_active);
    }
  }

  if (timer) TIMER_TOC(timer_do_black_holes_ghost);
}

/**
 * @brief Intermediate task after the BHs have done their swallowing step.
 * This is used to update the BH quantities if necessary.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
                                         int timer) {

  struct bpart *restrict bparts = c->black_holes.parts;
  const int count = c->black_holes.count;
  const struct engine *e = r->e;
  const int with_cosmology = e->policy & engine_policy_cosmology;

  TIMER_TIC;

  /* Anything to do here? */
  if (c->black_holes.count == 0) return;
  if (!cell_is_active_black_holes(c, e)) return;

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

    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {

      /* Get a direct pointer on the part. */
      struct bpart *bp = &bparts[i];

      if (bpart_is_active(bp, e)) {

        /* Get particle time-step */
        double dt;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(bp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, bp->time_bin);

          dt = cosmology_get_delta_time(e->cosmology, ti_begin,
                                        ti_begin + ti_step);
        } else {
          dt = get_timestep(bp->time_bin, e->time_base);
        }

        /* Compute the final operations for repositioning of this BH */
        black_holes_end_reposition(bp, e->black_holes_properties,
                                   e->physical_constants, e->cosmology, dt,
                                   e->ti_current);

        /* Compute variables required for the feedback loop */
        black_holes_prepare_feedback(bp, e->black_holes_properties,
                                     e->physical_constants, e->cosmology,
                                     e->cooling_func, e->entropy_floor, e->time,
                                     with_cosmology, dt, e->ti_current);
      }
    }
  }

  if (timer) TIMER_TOC(timer_do_black_holes_ghost);
}

/**
 * @brief Intermediate task after the gradient loop that does final operations
 * on the gradient quantities and optionally slope limits the gradients
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {

#ifdef EXTRA_HYDRO_LOOP

  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
  const int count = c->hydro.count;
  const struct engine *e = r->e;
  const integertime_t ti_current = e->ti_current;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
  const double time_base = e->time_base;
  const struct cosmology *cosmo = e->cosmology;
  const struct hydro_props *hydro_props = e->hydro_properties;
  const struct pressure_floor_props *pressure_floor = e->pressure_floor_props;

  TIMER_TIC;

  /* Anything to do here? */
  if (!cell_is_active_hydro(c, e)) return;

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

    /* Loop over the parts in this cell. */
    for (int i = 0; i < count; i++) {

      /* Get a direct pointer on the part. */
      struct part *restrict p = &parts[i];
      struct xpart *restrict xp = &xparts[i];

      if (part_is_active(p, e)) {

        /* Finish the gradient calculation */
        hydro_end_gradient(p);
        chemistry_end_gradient(p, e->chemistry);
        mhd_end_gradient(p);

        /* As of here, particle force variables will be set. */

        /* Calculate the time-step for passing to hydro_prepare_force.
         * This is the physical time between the start and end of the time-step
         * without any scale-factor powers. */
        double dt_alpha, dt_therm;

        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(ti_current - 1, p->time_bin);

          dt_alpha =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
          dt_therm = cosmology_get_therm_kick_factor(cosmo, ti_begin,
                                                     ti_begin + ti_step);
        } else {
          dt_alpha = get_timestep(p->time_bin, time_base);
          dt_therm = get_timestep(p->time_bin, time_base);
        }

        /* Compute variables required for the force loop */
        hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, dt_alpha,
                            dt_therm);
        chemistry_prepare_force(p, xp, cosmo, dt_alpha, dt_therm, e->chemistry);
        mhd_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
        timestep_limiter_prepare_force(p, xp);
        rt_prepare_force(p);
        rt_timestep_prepare_force(p);

        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
        mhd_reset_acceleration(p);
      }
    }
  }

  if (timer) TIMER_TOC(timer_do_extra_ghost);

#else
  error("SWIFT was not compiled with the extra hydro loop activated.");
#endif
}

/**
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param offset First particle in the cell to treat (for split tasks).
 * @param ntasks Interval between successive particles that are treated.
 * @param timer Are we timing this ?
 */
void runner_do_ghost(struct runner *r, struct cell *c, const int offset,
                     const int ntasks, const int timer) {

  struct part *restrict parts = c->hydro.parts;
  struct xpart *restrict xparts = c->hydro.xparts;
  const struct engine *e = r->e;
  const struct space *s = e->s;
  const struct hydro_space *hs = &s->hs;
  const struct cosmology *cosmo = e->cosmology;
  const struct chemistry_global_data *chemistry = e->chemistry;
  const struct star_formation *star_formation = e->star_formation;
  const struct hydro_props *hydro_props = e->hydro_properties;
  const struct pressure_floor_props *pressure_floor = e->pressure_floor_props;

  const int with_cosmology = (e->policy & engine_policy_cosmology);
  const int with_rt = (e->policy & engine_policy_rt);

  const float hydro_h_max = e->hydro_properties->h_max;
  const float hydro_h_min = e->hydro_properties->h_min;
  const float eps = e->hydro_properties->h_tolerance;
  const float hydro_eta_dim =
      pow_dimension(e->hydro_properties->eta_neighbours);
  const int use_mass_weighted_num_ngb =
      e->hydro_properties->use_mass_weighted_num_ngb;
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
  int redo = 0, count = 0;

  /* Running value of the maximal smoothing length */
  float h_max = c->hydro.h_max;
  float h_max_active = 0.f;

  TIMER_TIC;

  /* Anything to do here? */
  if (c->hydro.count == 0) return;
  if (!cell_is_active_hydro(c, e)) return;

  /* Recurse? */
  if (c->split) {
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL) {
        runner_do_ghost(r, c->progeny[k], offset, ntasks, /*timer=*/0);

        /* Update h_max */
        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);
      }
    }
  } else {

    /* Init the list of active particles that have to be updated and their
     * current smoothing lengths. */
    int *pid = NULL;
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
    if ((pid = (int *)malloc(sizeof(int) * c->hydro.count)) == NULL)
      error("Can't allocate memory for pid.");
    if ((h_0 = (float *)malloc(sizeof(float) * c->hydro.count)) == NULL)
      error("Can't allocate memory for h_0.");
    if ((left = (float *)malloc(sizeof(float) * c->hydro.count)) == NULL)
      error("Can't allocate memory for left.");
    if ((right = (float *)malloc(sizeof(float) * c->hydro.count)) == NULL)
      error("Can't allocate memory for right.");
    for (int k = offset; k < c->hydro.count; k += ntasks)
      if (part_is_active(&parts[k], e)) {
        pid[count] = k;
        h_0[count] = parts[k].h;
        left[count] = 0.f;
        right[count] = hydro_h_max;
        ++count;
      }

    /* While there are particles that need to be updated... */
    for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
         num_reruns++) {

      ghost_stats_account_for_hydro(&c->ghost_statistics, num_reruns, count,
                                    parts, pid);

      /* Reset the redo-count. */
      redo = 0;

      /* Loop over the remaining active parts in this cell. */
      for (int i = 0; i < count; i++) {

        /* Get a direct pointer on the part. */
        struct part *p = &parts[pid[i]];
        struct xpart *xp = &xparts[pid[i]];

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif

        /* Get some useful values */
        const float h_init = h_0[i];
        const float h_old = p->h;
        const float h_old_dim = pow_dimension(h_old);
        const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);

        float h_new;
        int has_no_neighbours = 0;

        if (p->density.wcount < 1.e-5 * kernel_root) { /* No neighbours case */

          ghost_stats_no_ngb_hydro_iteration(&c->ghost_statistics, num_reruns);

          /* Flag that there were no neighbours */
          has_no_neighbours = 1;

          /* Double h and try again */
          h_new = 2.f * h_old;

        } else {

          /* Finish the density calculation */
          hydro_end_density(p, cosmo);
          adaptive_softening_end_density(p, e->gravity_properties);
          mhd_end_density(p, cosmo);
          chemistry_end_density(p, chemistry, cosmo);
          star_formation_end_density(p, xp, star_formation, cosmo);

          /* Are we using the alternative definition of the
             number of neighbours? */
          if (use_mass_weighted_num_ngb) {
#if defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH) || defined(SHADOWFAX_SPH)
            error(
                "Can't use alternative neighbour definition with this scheme!");
#else
            const float inv_mass = 1.f / hydro_get_mass(p);
            p->density.wcount = p->rho * inv_mass;
            p->density.wcount_dh = p->density.rho_dh * inv_mass;
#endif
          }

          /* Compute one step of the Newton-Raphson scheme */
          const float n_sum = p->density.wcount * h_old_dim;
          const float n_target = hydro_eta_dim;
          const float f = n_sum - n_target;
          const float f_prime =
              p->density.wcount_dh * h_old_dim +
              hydro_dimension * p->density.wcount * h_old_dim_minus_one;

          /* Improve the bisection bounds */
          if (n_sum < n_target)
            left[i] = max(left[i], h_old);
          else if (n_sum > n_target)
            right[i] = min(right[i], h_old);

#ifdef SWIFT_DEBUG_CHECKS
          /* Check the validity of the left and right bounds */
          if (left[i] > right[i])
            error("Invalid left (%e) and right (%e)", left[i], right[i]);
#endif

          /* Skip if h is already h_max and we don't have enough neighbours */
          /* Same if we are below h_min */
          if (((p->h >= hydro_h_max) && (f < 0.f)) ||
              ((p->h <= hydro_h_min) && (f > 0.f))) {

            /* We have a particle whose smoothing length is already set (wants
             * to be larger but has already hit the maximum OR wants to be
             * smaller but has already reached the minimum). So, just tidy up
             * as if the smoothing length had converged correctly  */

#ifdef EXTRA_HYDRO_LOOP

            /* As of here, particle gradient variables will be set. */
            /* The force variables are set in the extra ghost. */

            /* Compute variables required for the gradient loop */
            hydro_prepare_gradient(p, xp, cosmo, hydro_props, pressure_floor);
            mhd_prepare_gradient(p, xp, cosmo, hydro_props);

            /* The particle gradient values are now set.  Do _NOT_
               try to read any particle density variables! */

            /* Prepare the particle for the gradient loop over neighbours
             */
            hydro_reset_gradient(p);
            mhd_reset_gradient(p);

#else
            /* Calculate the time-step for passing to hydro_prepare_force, used
             * for the evolution of alpha factors (i.e. those involved in the
             * artificial viscosity and thermal conduction terms) */
            const double time_base = e->time_base;
            const integertime_t ti_current = e->ti_current;
            double dt_alpha, dt_therm;

            if (with_cosmology) {
              const integertime_t ti_step = get_integer_timestep(p->time_bin);
              const integertime_t ti_begin =
                  get_integer_time_begin(ti_current - 1, p->time_bin);

              dt_alpha =
                  cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
              dt_therm = cosmology_get_therm_kick_factor(cosmo, ti_begin,
                                                         ti_begin + ti_step);
            } else {
              dt_alpha = get_timestep(p->time_bin, time_base);
              dt_therm = get_timestep(p->time_bin, time_base);
            }

            /* As of here, particle force variables will be set. */

            /* Compute variables required for the force loop */
            hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor,
                                dt_alpha, dt_therm);
            chemistry_prepare_force(p, xp, cosmo, dt_alpha, dt_therm,
                                    e->chemistry);
            mhd_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
            timestep_limiter_prepare_force(p, xp);
            rt_prepare_force(p);
            rt_timestep_prepare_force(p);

            /* The particle force values are now set.  Do _NOT_
               try to read any particle density variables! */

            /* Prepare the particle for the force loop over neighbours */
            hydro_reset_acceleration(p);
            mhd_reset_acceleration(p);

#endif /* EXTRA_HYDRO_LOOP */

            if (with_rt) {
#ifdef SWIFT_RT_DEBUG_CHECKS
              rt_debugging_check_nr_subcycles(p, e->rt_props);
#endif
              rt_reset_part(p, cosmo);
            }

            /* Check if h_max has increased */
            h_max = max(h_max, p->h);
            h_max_active = max(h_max_active, p->h);

            /* Ok, we are done with this particle */
            continue;
          }

          /* Normal case: Use Newton-Raphson to get a better value of h */

          /* Avoid floating point exception from f_prime = 0 */
          h_new = h_old - f / (f_prime + FLT_MIN);

          /* Be verbose about the particles that struggle to converge */
          if (num_reruns > max_smoothing_iter - 10) {

            message(
                "Smoothing length convergence problem: iter=%d p->id=%lld "
                "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
                "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
                num_reruns, p->id, h_init, h_old, h_new, f, f_prime, n_sum,
                n_target, left[i], right[i]);
          }

#ifdef SWIFT_DEBUG_CHECKS
          if (((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old)) &&
              (h_old < 0.999f * hydro_props->h_max))
            error(
                "Smoothing length correction not going in the right direction");
#endif

          /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
          h_new = min(h_new, 2.f * h_old);
          h_new = max(h_new, 0.5f * h_old);

          /* Verify that we are actually progrssing towards the answer */
          h_new = max(h_new, left[i]);
          h_new = min(h_new, right[i]);
        }

        /* Check whether the particle has an inappropriate smoothing length
         */
        if (fabsf(h_new - h_old) > eps * h_old) {

          /* Ok, correct then */

          /* Case where we have been oscillating around the solution */
          if ((h_new == left[i] && h_old == right[i]) ||
              (h_old == left[i] && h_new == right[i])) {

            /* Bisect the remaining interval */
            p->h = pow_inv_dimension(
                0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

          } else {

            /* Normal case */
            p->h = h_new;
          }

          /* If within the allowed range, try again */
          if (p->h < hydro_h_max && p->h > hydro_h_min) {

            /* Flag for another round of fun */
            pid[redo] = pid[i];
            h_0[redo] = h_0[i];
            left[redo] = left[i];
            right[redo] = right[i];
            redo += 1;

            /* Re-initialise everything */
            hydro_init_part(p, hs);
            adaptive_softening_init_part(p);
            mhd_init_part(p);
            chemistry_init_part(p, chemistry);
            star_formation_init_part(p, star_formation);
            tracers_after_init(p, xp, e->internal_units, e->physical_constants,
                               with_cosmology, e->cosmology,
                               e->hydro_properties, e->cooling_func, e->time);
            rt_init_part(p);

            /* Off we go ! */
            continue;

          } else if (p->h <= hydro_h_min) {

            /* Ok, this particle is a lost cause... */
            p->h = hydro_h_min;

          } else if (p->h >= hydro_h_max) {

            /* Ok, this particle is a lost cause... */
            p->h = hydro_h_max;

            /* Do some damage control if no neighbours at all were found */
            if (has_no_neighbours) {
              ghost_stats_no_ngb_hydro_converged(&c->ghost_statistics);
              hydro_part_has_no_neighbours(p, xp, cosmo);
              mhd_part_has_no_neighbours(p, xp, cosmo);
              chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
              star_formation_part_has_no_neighbours(p, xp, star_formation,
                                                    cosmo);
              rt_part_has_no_neighbours(p);
            }

          } else {
            error(
                "Fundamental problem with the smoothing length iteration "
                "logic.");
          }
        }

        /* We now have a particle whose smoothing length has converged */

        /* Set the correct depth */
        cell_set_part_h_depth(p, c);

        /* Check if h_max has increased */
        h_max = max(h_max, p->h);
        h_max_active = max(h_max_active, p->h);

        ghost_stats_converged_hydro(&c->ghost_statistics, p);

        /* Update gravitational softening (in adaptive softening case) */
        if (p->gpart)
          gravity_update_softening(p->gpart, p, e->gravity_properties);

#ifdef EXTRA_HYDRO_LOOP

        /* As of here, particle gradient variables will be set. */
        /* The force variables are set in the extra ghost. */

        /* Compute variables required for the gradient loop */
        hydro_prepare_gradient(p, xp, cosmo, hydro_props, pressure_floor);
        mhd_prepare_gradient(p, xp, cosmo, hydro_props);

        /* The particle gradient values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the gradient loop over neighbours */
        hydro_reset_gradient(p);
        mhd_reset_gradient(p);

#else

        /* Calculate the time-step for passing to hydro_prepare_force, used
         * for the evolution of alpha factors (i.e. those involved in the
         * artificial viscosity and thermal conduction terms) */
        const double time_base = e->time_base;
        const integertime_t ti_current = e->ti_current;
        double dt_alpha, dt_therm;

        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(p->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(ti_current - 1, p->time_bin);

          dt_alpha =
              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
          dt_therm = cosmology_get_therm_kick_factor(cosmo, ti_begin,
                                                     ti_begin + ti_step);
        } else {
          dt_alpha = get_timestep(p->time_bin, time_base);
          dt_therm = get_timestep(p->time_bin, time_base);
        }

        /* As of here, particle force variables will be set. */

        /* Compute variables required for the force loop */
        hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, dt_alpha,
                            dt_therm);
        chemistry_prepare_force(p, xp, cosmo, dt_alpha, dt_therm, e->chemistry);
        mhd_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
        timestep_limiter_prepare_force(p, xp);
        rt_prepare_force(p);
        rt_timestep_prepare_force(p);

        /* The particle force values are now set.  Do _NOT_
           try to read any particle density variables! */

        /* Prepare the particle for the force loop over neighbours */
        hydro_reset_acceleration(p);
        mhd_reset_acceleration(p);

#endif /* EXTRA_HYDRO_LOOP */

        if (with_rt) {
#ifdef SWIFT_RT_DEBUG_CHECKS
          rt_debugging_check_nr_subcycles(p, e->rt_props);
#endif
          rt_reset_part(p, cosmo);
        }
      }

      /* We now need to treat the particles whose smoothing length had not
       * converged again */

      /* Re-set the counter for the next loop (potentially). */
      count = redo;
      if (count > 0) {

        /* Climb up the cell hierarchy. */
        for (struct cell *finger = c; finger != NULL; finger = finger->parent) {

          /* Run through this cell's density interactions. */
          for (struct link *l = finger->hydro.density; l != NULL; l = l->next) {

#ifdef SWIFT_DEBUG_CHECKS
            if (l->t->ti_run < r->e->ti_current)
              error("Density task should have been run.");
#endif

            /* Self-interaction? */
            if (l->t->type == task_type_self) {
              runner_dosub_self_subset_density(r, finger, parts, pid, count, 1);
            }

            /* Otherwise, pair interaction? */
            else if (l->t->type == task_type_pair) {

              /* Left or right? */
              if (l->t->ci == finger)
                runner_dosub_pair_subset_density(r, finger, parts, pid, count,
                                                 l->t->cj, 1);
              else
                runner_dosub_pair_subset_density(r, finger, parts, pid, count,
                                                 l->t->ci, 1);
            } else {
#ifdef SWIFT_DEBUG_CHECKS
              error("Invalid sub-type!");
#endif
            }
          }
        }
      }
    }

    if (count) {
      warning(
          "Smoothing length failed to converge for the following gas "
          "particles:");
      for (int i = 0; i < count; i++) {
        struct part *p = &parts[pid[i]];
        warning("ID: %lld, h: %g, wcount: %g", p->id, p->h, p->density.wcount);
      }

      error("Smoothing length failed to converge on %i particles.", count);
    }

    /* Be clean */
    free(left);
    free(right);
    free(pid);
    free(h_0);

#ifdef SWIFT_DEBUG_CHECKS
    /* Update h_max */
    atomic_max_f(&c->hydro.h_max, h_max);
    atomic_max_f(&c->hydro.h_max_active, h_max_active);

    for (int i = offset; i < c->hydro.count; i += ntasks) {
      const struct part *p = &c->hydro.parts[i];
      const float h = c->hydro.parts[i].h;
      if (part_is_inhibited(p, e)) continue;

      if (h > c->hydro.h_max)
        error("Particle has h larger than h_max (id=%lld)", p->id);
      if (part_is_active(p, e) && h > c->hydro.h_max_active)
        error("Active particle has h larger than h_max_active (id=%lld)",
              p->id);
    }
#endif
  }

  /* Update h_max */
  atomic_max_f(&c->hydro.h_max, h_max);
  atomic_max_f(&c->hydro.h_max_active, h_max_active);

  /* The ghost may not always be at the top level.
   * Therefore we need to update h_max between the super- and top-levels
   * Note: We can check the 0th entry as they are either all NULL or none */
  if (c->hydro.ghost[0]) {
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
      atomic_max_f(&tmp->hydro.h_max, h_max);
      atomic_max_f(&tmp->hydro.h_max_active, h_max_active);
    }
  }

  if (timer) TIMER_TOC(timer_do_ghost);
}

/**
 * @brief Intermediate task after the injection to compute the total
 * photon emission rates experienced by the hydro particles.
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_rt_ghost1(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;
  const int with_cosmology = (e->policy & engine_policy_cosmology);
  const struct cosmology *cosmo = e->cosmology;
  int count = c->hydro.count;

  /* Anything to do here? */
  if (count == 0) return;
  if (!cell_is_rt_active(c, e)) return;

  TIMER_TIC;

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

    for (int pid = 0; pid < count; pid++) {
      struct part *restrict p = &(c->hydro.parts[pid]);

      /* Skip inhibited parts */
      if (part_is_inhibited(p, e)) continue;

      /* Skip inactive parts */
      if (!part_is_rt_active(p, e)) continue;

      /* First reset everything that needs to be reset for the following
       * subcycle */
      const integertime_t ti_current_subcycle = e->ti_current_subcycle;
      const integertime_t ti_step =
          get_integer_timestep(p->rt_time_data.time_bin);
      const integertime_t ti_begin = get_integer_time_begin(
          ti_current_subcycle + 1, p->rt_time_data.time_bin);
      const integertime_t ti_end = ti_begin + ti_step;

      const float dt =
          rt_part_dt(ti_begin, ti_end, e->time_base, with_cosmology, cosmo);

      rt_reset_part_each_subcycle(p, cosmo, dt);

      /* Now finish up injection */
      rt_finalise_injection(p, e->rt_props);
    }
  }

  if (timer) TIMER_TOC(timer_do_rt_ghost1);
}

/**
 * @brief Intermediate task after the gradient computation to finish
 *        up the gradients and prepare for the transport step
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_rt_ghost2(struct runner *r, struct cell *c, int timer) {

  const struct engine *e = r->e;
  int count = c->hydro.count;
  const struct cosmology *cosmo = e->cosmology;

  /* Anything to do here? */
  if (count == 0) return;
  if (!cell_is_rt_active(c, e)) return;

  TIMER_TIC;

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

    for (int pid = 0; pid < count; pid++) {
      struct part *restrict p = &(c->hydro.parts[pid]);

      /* Skip inhibited parts */
      if (part_is_inhibited(p, e)) continue;

      /* Skip inactive parts */
      if (!part_is_rt_active(p, e)) continue;

      rt_end_gradient(p, cosmo);
    }
  }

  if (timer) TIMER_TOC(timer_do_rt_ghost2);
}

/**
 * @brief Intermediate task after the density to check that the smoothing
 * lengths are correct, finish density calculations
 * and calculate accretion rates for the particle swallowing step
 *
 * @param r The runner thread.
 * @param c The cell.
 * @param timer Are we timing this ?
 */
void runner_do_sinks_density_ghost(struct runner *r, struct cell *c,
                                   int timer) {

  struct sink *restrict sinks = c->sinks.parts;
  const struct engine *e = r->e;
  const struct cosmology *cosmo = e->cosmology;
  const int with_cosmology = e->policy & engine_policy_cosmology;
  const float sinks_h_max = e->hydro_properties->h_max;
  const float sinks_h_min = e->hydro_properties->h_min;
  const float eps = e->sink_properties->h_tolerance;
  const float sinks_eta_dim = pow_dimension(e->sink_properties->eta_neighbours);
  const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
  int redo = 0, scount = 0;

  /* Running value of the maximal smoothing length */
  float h_max = c->sinks.h_max;
  float h_max_active = 0.f;

  TIMER_TIC;

  /* Anything to do here? */
  if (c->sinks.count == 0) return;
  if (!cell_is_active_sinks(c, e)) return;

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

        /* Update h_max */
        h_max = max(h_max, c->progeny[k]->sinks.h_max);
        h_max_active = max(h_max_active, c->progeny[k]->sinks.h_max_active);
      }
    }
  } else {

    /* Init the list of active particles that have to be updated. */
    int *sid = NULL;
    float *h_0 = NULL;
    float *left = NULL;
    float *right = NULL;
    if ((sid = (int *)malloc(sizeof(int) * c->sinks.count)) == NULL)
      error("Can't allocate memory for sid.");
    if ((h_0 = (float *)malloc(sizeof(float) * c->sinks.count)) == NULL)
      error("Can't allocate memory for h_0.");
    if ((left = (float *)malloc(sizeof(float) * c->sinks.count)) == NULL)
      error("Can't allocate memory for left.");
    if ((right = (float *)malloc(sizeof(float) * c->sinks.count)) == NULL)
      error("Can't allocate memory for right.");
    for (int k = 0; k < c->sinks.count; k++)
      if (sink_is_active(&sinks[k], e)) {
        sid[scount] = k;
        h_0[scount] = sinks[k].h;
        left[scount] = 0.f;
        right[scount] = sinks_h_max;
        ++scount;
      }

    if (e->sink_properties->use_fixed_r_cut) {
      /* If we're using a fixed cutoff rather than a smoothing length, just
       * finish up the density task and leave sp->h untouched. */

      /* Loop over the active sinks in this cell. */
      for (int i = 0; i < scount; i++) {

        /* Get a direct pointer on the part. */
        struct sink *sp = &sinks[sid[i]];

#ifdef SWIFT_DEBUG_CHECKS
        /* Is this part within the timestep? */
        if (!sink_is_active(sp, e)) error("Ghost applied to inactive particle");
#endif

        /* Finish the density calculation */
        sink_end_density(sp, cosmo);

        /* Set these variables to the fixed cutoff radius for the rest of the
         * ghost task */
        h_max = sp->h;
        h_max_active = sp->h;
      }

    } else {
      /* Otherwise we need to iterate to update the smoothing lengths */

      /* While there are particles that need to be updated... */
      for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
           num_reruns++) {

        /* Reset the redo-count. */
        redo = 0;

        /* Loop over the remaining active parts in this cell. */
        for (int i = 0; i < scount; i++) {

          /* Get a direct pointer on the part. */
          struct sink *sp = &sinks[sid[i]];

#ifdef SWIFT_DEBUG_CHECKS
          /* Is this part within the timestep? */
          if (!sink_is_active(sp, e))
            error("Ghost applied to inactive particle");
#endif

          /* Get some useful values */
          const float h_init = h_0[i];
          const float h_old = sp->h;
          const float h_old_dim = pow_dimension(h_old);
          const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);

          float h_new;
          int has_no_neighbours = 0;

          if (sp->density.wcount <
              1.e-5 * kernel_root) { /* No neighbours case */

            /* Flag that there were no neighbours */
            has_no_neighbours = 1;

            /* Double h and try again */
            h_new = 2.f * h_old;

          } else {

            /* Finish the density calculation */
            sink_end_density(sp, cosmo);

            /* Compute one step of the Newton-Raphson scheme */
            const float n_sum = sp->density.wcount * h_old_dim;
            const float n_target = sinks_eta_dim;
            const float f = n_sum - n_target;
            const float f_prime =
                sp->density.wcount_dh * h_old_dim +
                hydro_dimension * sp->density.wcount * h_old_dim_minus_one;

            /* Improve the bisection bounds */
            if (n_sum < n_target)
              left[i] = max(left[i], h_old);
            else if (n_sum > n_target)
              right[i] = min(right[i], h_old);

#ifdef SWIFT_DEBUG_CHECKS
            /* Check the validity of the left and right bounds */
            if (left[i] > right[i])
              error("Invalid left (%e) and right (%e)", left[i], right[i]);
#endif

            /* Skip if h is already h_max and we don't have enough neighbours
             */
            /* Same if we are below h_min */
            if (((sp->h >= sinks_h_max) && (f < 0.f)) ||
                ((sp->h <= sinks_h_min) && (f > 0.f))) {

              /* Check if h_max has increased */
              h_max = max(h_max, sp->h);
              h_max_active = max(h_max_active, sp->h);

              /* Ok, we are done with this particle */
              continue;
            }

            /* Normal case: Use Newton-Raphson to get a better value of h */

            /* Avoid floating point exception from f_prime = 0 */
            h_new = h_old - f / (f_prime + FLT_MIN);

            /* Be verbose about the particles that struggle to converge */
            if (num_reruns > max_smoothing_iter - 10) {

              message(
                  "Smoothing length convergence problem: iter=%d p->id=%lld "
                  "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
                  "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
                  num_reruns, sp->id, h_init, h_old, h_new, f, f_prime, n_sum,
                  n_target, left[i], right[i]);
            }

            /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
            h_new = min(h_new, 2.f * h_old);
            h_new = max(h_new, 0.5f * h_old);

            /* Verify that we are actually progrssing towards the answer */
            h_new = max(h_new, left[i]);
            h_new = min(h_new, right[i]);
          }

          /* Check whether the particle has an inappropriate smoothing length
           */
          if (fabsf(h_new - h_old) > eps * h_old) {

            /* Ok, correct then */

            /* Case where we have been oscillating around the solution */
            if ((h_new == left[i] && h_old == right[i]) ||
                (h_old == left[i] && h_new == right[i])) {

              /* Bisect the remaining interval */
              sp->h = pow_inv_dimension(
                  0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));

            } else {

              /* Normal case */
              sp->h = h_new;
            }

            /* If below the absolute maximum, try again */
            if (sp->h < sinks_h_max && sp->h > sinks_h_min) {

              /* Flag for another round of fun */
              sid[redo] = sid[i];
              h_0[redo] = h_0[i];
              left[redo] = left[i];
              right[redo] = right[i];
              redo += 1;

              /* Re-initialise everything */
              sink_init_sink(sp);

              /* Off we go ! */
              continue;

            } else if (sp->h <= sinks_h_min) {

              /* Ok, this particle is a lost cause... */
              sp->h = sinks_h_min;

            } else if (sp->h >= sinks_h_max) {

              /* Ok, this particle is a lost cause... */
              sp->h = sinks_h_max;

              /* Do some damage control if no neighbours at all were found */
              if (has_no_neighbours) {
                sinks_sink_has_no_neighbours(sp, cosmo);
              }

            } else {
              error(
                  "Fundamental problem with the smoothing length iteration "
                  "logic.");
            }
          }

          /* We now have a particle whose smoothing length has converged */

          /* Set the correct depth */
          cell_set_sink_h_depth(sp, c);

          /* Check if h_max has increased */
          h_max = max(h_max, sp->h);
          h_max_active = max(h_max_active, sp->h);
        }

        /* We now need to treat the particles whose smoothing length had not
         * converged again */

        /* Re-set the counter for the next loop (potentially). */
        scount = redo;
        if (scount > 0) {

          /* Climb up the cell hierarchy. */
          for (struct cell *finger = c; finger != NULL;
               finger = finger->parent) {

            /* Run through this cell's density interactions. */
            for (struct link *l = finger->sinks.density; l != NULL;
                 l = l->next) {

#ifdef SWIFT_DEBUG_CHECKS
              if (l->t->ti_run < r->e->ti_current)
                error("Density task should have been run.");
#endif

              /* Self-interaction? */
              if (l->t->type == task_type_self) {
                runner_dosub_subset_sinks_density(r, finger, sinks, sid, scount,
                                                  NULL, 1);
              }

              /* Otherwise, pair interaction? */
              else if (l->t->type == task_type_pair) {

                /* Left or right? */
                if (l->t->ci == finger)
                  runner_dosub_subset_sinks_density(r, finger, sinks, sid,
                                                    scount, l->t->cj, 1);
                else
                  runner_dosub_subset_sinks_density(r, finger, sinks, sid,
                                                    scount, l->t->ci, 1);
              } else {
#ifdef SWIFT_DEBUG_CHECKS
                error("Invalid sub-type!");
#endif
              }
            }
          }
        }
      }

      if (scount) {
        warning(
            "Smoothing length failed to converge for the following sink "
            "particles:");
        for (int i = 0; i < scount; i++) {
          struct sink *sp = &sinks[sid[i]];
          warning("ID: %lld, h: %g, wcount: %g", sp->id, sp->h,
                  sp->density.wcount);
        }

        error("Smoothing length failed to converge on %i particles.", scount);
      }

      /* Be clean */
      free(left);
      free(right);
      free(sid);
      free(h_0);
    }

    /* We need one more quick loop over the sinks to run prepare_swallow */
    for (int i = 0; i < c->sinks.count; i++) {

      /* Get a direct pointer on the part. */
      struct sink *sp = &sinks[i];

      if (sink_is_active(sp, e)) {

        /* Get particle time-step */
        double dt;
        if (with_cosmology) {
          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
          const integertime_t ti_begin =
              get_integer_time_begin(e->ti_current - 1, sp->time_bin);

          dt = cosmology_get_delta_time(e->cosmology, ti_begin,
                                        ti_begin + ti_step);
        } else {
          dt = get_timestep(sp->time_bin, e->time_base);
        }

        /* Calculate the accretion rate and accreted mass this timestep, for use
         * in swallow loop */
        sink_prepare_swallow(sp, e->sink_properties, e->physical_constants,
                             e->cosmology, e->cooling_func, e->entropy_floor,
                             e->time, with_cosmology, dt, e->ti_current);
      }
    }
  }

  /* Update h_max */
  c->sinks.h_max = h_max;
  c->sinks.h_max_active = h_max_active;

#ifdef SWIFT_DEBUG_CHECKS
  for (int i = 0; i < c->sinks.count; ++i) {
    const struct sink *sp = &c->sinks.parts[i];
    const float h = c->sinks.parts[i].h;
    if (sink_is_inhibited(sp, e)) continue;

    if (h > c->sinks.h_max)
      error("Particle has h larger than h_max (id=%lld)", sp->id);
    if (sink_is_active(sp, e) && h > c->sinks.h_max_active)
      error("Active particle has h larger than h_max_active (id=%lld)", sp->id);
  }
#endif

  /* The ghost may not always be at the top level.
   * Therefore we need to update h_max between the super- and top-levels */
  if (c->sinks.density_ghost) {
    for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
      atomic_max_f(&tmp->sinks.h_max, h_max);
      atomic_max_f(&tmp->sinks.h_max_active, h_max_active);
    }
  }

  if (timer) TIMER_TOC(timer_do_sinks_ghost);
}
