add distance calculation from Whisky

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

5 files changed, 56 insertions, 58 deletions

diff --git a/interface.ccl b/interface.ccl
index 060ec29..5d6ff08 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -13,3 +13,5 @@ CCTK_REAL Hydro_Analysis_rho_max_loc[3] type = SCALAR tags='checkpoint="no"' "co
 CCTK_REAL Hydro_Analysis_rho_center_volume_weighted_gf[3] type = GF Timelevels = 3 tags='checkpoint="no" tensortypealias="vector"' "temporary GF to obtain the coordinate location of the volume weightes center of mass via a reduction"
 
 CCTK_REAL Hydro_Analysis_rho_center_volume_weighted[3] type = SCALAR tags='checkpoint="no"' "coordinate location of the volume weightes center of mass"
+
+CCTK_REAL Hydro_Analysis_rho_max_origin_distance  type = SCALAR tags='checkpoint="no"' "proper distance between the maximum of the density and the origin (along a straight coordinate line)"
diff --git a/param.ccl b/param.ccl
index 3ab7385..b30d7e8 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,6 +1,6 @@
 # Parameter definitions for thorn Hydro_Analysis
 
-BOOLEAN Hydro_Analysis_rho_max_search "Look for the value and location of the maximum of rho"
+BOOLEAN Hydro_Analysis_comp_rho_max "Look for the value and location of the maximum of rho"
 {
 } "false"
 
@@ -9,7 +9,16 @@ BOOLEAN Hydro_Analysis_rho_max_loc_only_positive_x "Restrict location search for
 } "false"
 
 
-BOOLEAN Hydro_Analysis_weighted_center_of_mass "Look for the location of the center of mass"
+BOOLEAN Hydro_Analysis_comp_vol_weighted_center_of_mass "Look for the location of the volume-weighted center of mass"
 {
 } "false"
 
+BOOLEAN Hydro_Analysis_comp_rho_max_origin_distance "Look for the proper distance between the maximum of the density and the origin (along a straight coordinate line)"
+{
+} "false"
+
+CCTK_INT Hydro_Analysis_rho_max_origin_distance_npoints "Number of points along the straight line for measuring proper distance"
+{
+  1:* :: "Any positive number"
+} 100
+
diff --git a/schedule.ccl b/schedule.ccl
index b9fba6e..f3f11a9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,13 +1,24 @@
 # Schedule definitions for thorn HydroAnalysis
 
-STORAGE: Hydro_Analysis_rho_max
-STORAGE: Hydro_Analysis_rho_max_loc
-STORAGE: Hydro_Analysis_rho_center_volume_weighted
+if (Hydro_Analysis_comp_rho_max || Hydro_Analysis_comp_rho_max_origin_distance)
+{
+  STORAGE: Hydro_Analysis_rho_max
+  STORAGE: Hydro_Analysis_rho_max_loc
+}
+if (Hydro_Analysis_comp_vol_weighted_center_of_mass)
+{
+  STORAGE: Hydro_Analysis_rho_sum
+  STORAGE: Hydro_Analysis_rho_center_volume_weighted
   STORAGE: Hydro_Analysis_rho_center_volume_weighted_gf[3]
+}
+if (Hydro_Analysis_comp_rho_max_origin_distance)
+{
+  STORAGE: Hydro_Analysis_rho_max_origin_distance
+}
+
 
 schedule group Hydro_Analysis AT POSTSTEP
 {
-  STORAGE: Hydro_Analysis_rho_sum
 } "Group for Hydro_Analysis routines"
 
 schedule Hydro_Analysis_PrepareReduction IN Hydro_Analysis
@@ -22,7 +33,7 @@ schedule Hydro_Analysis_Reduction IN Hydro_Analysis AFTER Hydro_Analysis_Prepare
   OPTIONS: global
 } "Compute the global reduction results"
 
-if (Hydro_Analysis_rho_max_search)
+if (Hydro_Analysis_comp_rho_max || Hydro_Analysis_comp_rho_max_origin_distance)
 {
   schedule Hydro_Analysis_LocationSearch IN Hydro_Analysis AFTER Hydro_Analysis_Reduction
   {
@@ -31,3 +42,11 @@ if (Hydro_Analysis_rho_max_search)
   } "Look for the location of the maximum density"
 }
 
+if (Hydro_Analysis_comp_rho_max_origin_distance)
+{
+  schedule Hydro_Analysis_FindSeparation IN Hydro_Analysis AFTER Hydro_Analysis_LocationSearch
+  {
+    LANG: Fortran
+    OPTIONS: global
+  } "Compute the proper distance between the maximum of the density and the origin (along a straight coordinate line)"
+}
diff --git a/src/Hydro_Analysis.c b/src/Hydro_Analysis.c
index 5fef8a1..7d182ea 100644
--- a/src/Hydro_Analysis.c
+++ b/src/Hydro_Analysis.c
@@ -4,13 +4,13 @@
 
 #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
-  printf("local\n");
 
-  if (Hydro_Analysis_weighted_center_of_mass)
+  if (Hydro_Analysis_comp_vol_weighted_center_of_mass)
   {
     #pragma omp parallel
     {
@@ -30,12 +30,11 @@ void Hydro_Analysis_PrepareReduction(CCTK_ARGUMENTS)
   }
 }
 
-// Hydro_Analysis_rho_center_volume_weighted
+/* Reduce grid functions */
 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)
@@ -44,68 +43,37 @@ void Hydro_Analysis_Reduction(CCTK_ARGUMENTS)
   if (handle_sum < 0)
     CCTK_WARN(0, "Unable to get reduction handle 'sum'.");
 
-  if (Hydro_Analysis_rho_max_search)
+  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_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
-  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");
-}
-
diff --git a/src/make.code.defn b/src/make.code.defn
index cb708ae..8372f8a 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn Hydro_Analysis
 
 # Source files in this directory
-SRCS = Hydro_Analysis.c
+SRCS = Hydro_Analysis.c Whisky_Separation.F90
 
 # Subdirectories containing source files
 SUBDIRS =