#include #include #include #include "loopcontrol.h" void Hydro_Analysis_PrepareReduction(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS printf("local\n"); if (Hydro_Analysis_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); } } } // 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"); }