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

#include "loopcontrol.h"

/* 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 */
void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS