Hydro_Analysis: add option to average the location of

 multiple identical maxima this is useful for cell
 centering where intially there might be several
 points with precisely the same density around the
 center of a star



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

1 files changed, 41 insertions, 10 deletions

diff --git a/src/Hydro_Analysis.c b/src/Hydro_Analysis.c
index adb2e99..bcd3430 100644
--- a/src/Hydro_Analysis.c
+++ b/src/Hydro_Analysis.c
@@ -1,3 +1,5 @@
+#include <math.h>
+
 #include <cctk.h>
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
@@ -81,11 +83,15 @@ void Hydro_Analysis_Reduction(CCTK_ARGUMENTS)
 }
 
 /* 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;
@@ -109,10 +115,10 @@ void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
       {
         #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.;
+          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.;
         }
       }
     }
@@ -122,25 +128,50 @@ void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
    * 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 (fabs(level_rho_max_loc[3]) > 1.e-14) {
-    if (fabs(level_rho_max_loc[3]-1.) < 1.e-14)
+  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];
-    } else
-      CCTK_WARN(1, "Found more than one identical maximum, not setting anything.");
-
   }
 /*  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]); */
 }
-