add forgotten file

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Hydro_Analysis/trunk@27 b6729ddc-ac74-4bd1-908c-9dc7244c52a1

1 files changed, 108 insertions, 0 deletions

diff --git a/src/Hydro_Analysis_Maxima.c b/src/Hydro_Analysis_Maxima.c
new file mode 100644
index 0000000..b3aae2a
--- /dev/null
+++ b/src/Hydro_Analysis_Maxima.c
@@ -0,0 +1,108 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include "loopcontrol.h"
+
+void Hydro_Analysis_PrepareReduction(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  printf("local\n");
+
+    #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);
+  }
+}
+
+// Hydro_Analysis_rho_center_volume_weighted
+void Hydro_Analysis_Reduction(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  printf("global\n");
+
+  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_rho_max_search)
+  {
+    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 (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_WARN(0, "Error while reducing Hydro_Analysis_rho_center_volume_weighted_gf");
+}
+
+void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  printf("search\n");
+
+  /* Initialized MPI-local quantities */
+  CCTK_REAL local_rho_max_loc[3];
+  local_rho_max_loc[0] = 0.0;
+  local_rho_max_loc[1] = 0.0;
+  local_rho_max_loc[2] = 0.0;
+
+  /* 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];
+        }
+      }
+    }
+    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. */
+  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,
+                                   Hydro_Analysis_rho_max_loc, 3, CCTK_VARIABLE_REAL))
+    CCTK_WARN(0, "Error while reducing local_rho_max_loc");
+}
+