diff options
-rw-r--r-- | param.ccl | 7 | ||||
-rw-r--r-- | src/Hydro_Analysis.c | 51 |
2 files changed, 48 insertions, 10 deletions
@@ -4,6 +4,13 @@ BOOLEAN Hydro_Analysis_comp_rho_max "Look for the value and location of the maxi { } "false" +# the following is useful eg. when using cell-centering where there might not +# be a grid point at the centre of the star and initially several equally valid +# maxima might be clustered around it. +BOOLEAN Hydro_Analysis_average_multiple_maxima_locations "when finding more than one global maximum location, average position and use result" STEERABLE = always +{ +} "false" + BOOLEAN Hydro_Analysis_rho_max_loc_only_positive_x "Restrict location search for density maximum to positive values of x" { } "false" 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]); */ } - |