#include <math.h>

#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>

#include "loopcontrol.h"

void Hydro_Analysis_Init(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  if (Hydro_Analysis_comp_rho_max || Hydro_Analysis_comp_rho_max_origin_distance)
    *Hydro_Analysis_rho_max = -1.;
}

/* Populate grid functions to be reduced in Hydro_Analysis_Reduction */
void Hydro_Analysis_PrepareReduction(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  if (Hydro_Analysis_comp_vol_weighted_center_of_mass)
  {
    #pragma omp parallel
    {
      LC_LOOP3(local_reduction, i,j,k,
               cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[2],
               cctk_lsh[0] - cctk_nghostzones[0], cctk_lsh[1] - cctk_nghostzones[1],
               cctk_lsh[2] - cctk_nghostzones[2], cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
      {
        CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
        CCTK_INT size= cctk_lsh[0] *  cctk_lsh[1] * cctk_lsh[2];
        Hydro_Analysis_rho_center_volume_weighted_gf[i3D+0*size] = rho[i3D] * x[i3D];
        Hydro_Analysis_rho_center_volume_weighted_gf[i3D+1*size] = rho[i3D] * y[i3D];
        Hydro_Analysis_rho_center_volume_weighted_gf[i3D+2*size] = rho[i3D] * z[i3D];
      }
      LC_ENDLOOP3(local_reduction);
    }
  }
}

/* Reduce grid functions */
void Hydro_Analysis_Reduction(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT handle_max = CCTK_ReductionHandle("maximum");
  if (handle_max < 0)
    CCTK_WARN(0, "Unable to get reduction handle 'max'.");
  CCTK_INT handle_sum = CCTK_ReductionHandle("sum");
  if (handle_sum < 0)
    CCTK_WARN(0, "Unable to get reduction handle 'sum'.");

  if (Hydro_Analysis_comp_rho_max)
  {
    if (CCTK_Reduce(cctkGH, -1, handle_max, 1, CCTK_VARIABLE_REAL,
                    Hydro_Analysis_rho_max, 1,
                    CCTK_VarIndex("HydroBase::rho")))
      CCTK_WARN(0, "Error while reducing HydroBase::rho");
  }
  if (Hydro_Analysis_comp_vol_weighted_center_of_mass)
  {
    if (CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
                    Hydro_Analysis_rho_sum, 1,
                    CCTK_VarIndex("HydroBase::rho")))
      CCTK_WARN(0, "Error while reducing HydroBase::rho");
    if (CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
                    &(Hydro_Analysis_rho_center_volume_weighted[0]), 1,
                    CCTK_VarIndex(
                      "Hydro_Analysis::Hydro_Analysis_rho_center_volume_weighted_gf[0]")) ||
        CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
                    &(Hydro_Analysis_rho_center_volume_weighted[1]), 1,
                    CCTK_VarIndex(
                      "Hydro_Analysis::Hydro_Analysis_rho_center_volume_weighted_gf[1]")) ||
        CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
                    &(Hydro_Analysis_rho_center_volume_weighted[2]), 1,
                    CCTK_VarIndex(
                      "Hydro_Analysis::Hydro_Analysis_rho_center_volume_weighted_gf[2]")))
      CCTK_WARN(0, "Error while reducing Hydro_Analysis_rho_center_volume_weighted_gf");
  }
}

/* Look for the location of the density maximum */
/* RH: this is scheduled global loop-local but does MPI calls. It will hang if
 * some processors have more components than others. */
void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  int set_maximum_location;

  /* Initialized MPI-local quantities */
  CCTK_REAL local_rho_max_loc[4], level_rho_max_loc[4];
  local_rho_max_loc[0] = 0.0;
  local_rho_max_loc[1] = 0.0;
  local_rho_max_loc[2] = 0.0;
  local_rho_max_loc[3] = 0.0; // Information if it was found at all on this level

  /* Look for the location of the global maximum.
   * This algorithm will have problems when that occurs at more than
   * one location. */
  #pragma omp parallel
  {
    LC_LOOP3(location_reduction, i,j,k,
             cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[2],
             cctk_lsh[0] - cctk_nghostzones[0], cctk_lsh[1] - cctk_nghostzones[1],
             cctk_lsh[2] - cctk_nghostzones[2], cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
    {
      CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
      if (rho[i3D] == *Hydro_Analysis_rho_max &&
          (!Hydro_Analysis_rho_max_loc_only_positive_x || x[i3D] >= 0.0))
      {
        #pragma omp critical
        {
          local_rho_max_loc[0] += x[i3D];
          local_rho_max_loc[1] += y[i3D];
          local_rho_max_loc[2] += z[i3D];
          local_rho_max_loc[3] += 1.;
        }
      }
    }
    LC_ENDLOOP3(location_reduction);
  }
  /* After this at least one MPI process should have the location set.
   * The others should have a value of 0.0, so that the sum reduction
   * gives the value on each MPI process. This however is not true sometimes
   * and this is the point where the algorithm can fail. */
  /* RH: this assumption is false if a processor owns more than a single patch
   * */
  CCTK_INT handle_sum = CCTK_LocalArrayReductionHandle("sum");
  if (handle_sum < 0)
     CCTK_WARN(0, "Unable to get reduction handle 'sum'.");
  if (CCTK_ReduceLocArrayToArray1D(cctkGH, -1, handle_sum, &local_rho_max_loc,
                                   level_rho_max_loc, 4, CCTK_VARIABLE_REAL))
    CCTK_WARN(0, "Error while reducing local_rho_max_loc");
  if (round(level_rho_max_loc[3]) > 0.) {
    if (round(level_rho_max_loc[3]) == 1.)
    {
      set_maximum_location = 1;
    } else {
      if (Hydro_Analysis_average_multiple_maxima_locations)
      {
        static int have_warned_about_multiple_maxima = 0;

        level_rho_max_loc[0] /= level_rho_max_loc[3];
        level_rho_max_loc[1] /= level_rho_max_loc[3];
        level_rho_max_loc[2] /= level_rho_max_loc[3];

        set_maximum_location = 1;
      } else {
        CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "Found more than one (%d) identical maximum, not setting anything.",
                   (int)round(level_rho_max_loc[3]));
        }
        set_maximum_location = 0;
      }
    }

  } else {
    set_maximum_location = 0;
  }


  if (set_maximum_location)
  {
      Hydro_Analysis_rho_max_loc[0] = level_rho_max_loc[0];
      Hydro_Analysis_rho_max_loc[1] = level_rho_max_loc[1];
      Hydro_Analysis_rho_max_loc[2] = level_rho_max_loc[2];
  }
/*  CCTK_VInfo(CCTK_THORNSTRING, "New location: %g,%g,%g",
    Hydro_Analysis_rho_max_loc[0],
    Hydro_Analysis_rho_max_loc[1],
    Hydro_Analysis_rho_max_loc[2]); */
}