diff options
-rw-r--r-- | schedule.ccl | 2 | ||||
-rw-r--r-- | src/setup.cc | 92 |
2 files changed, 65 insertions, 29 deletions
diff --git a/schedule.ccl b/schedule.ccl index 53827c4..89aec32 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -16,11 +16,13 @@ STORAGE: sf_maxreflevel SCHEDULE SphericalSurface_SetupRes AT basegrid BEFORE SphericalSurface_Setup AFTER SpatialCoordinates AFTER CorrectCoordinates { LANG: C + OPTIONS: GLOBAL LOOP-LOCAL } "Set surface resolution automatically" SCHEDULE SphericalSurface_Setup AT basegrid { LANG: C + OPTIONS: GLOBAL } "Calculate surface coordinate descriptors" diff --git a/src/setup.cc b/src/setup.cc index 897efb3..fa22c5d 100644 --- a/src/setup.cc +++ b/src/setup.cc @@ -164,15 +164,15 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS) if (verbose) { - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_ntheta = %d\n", n, sf_ntheta[n]); - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_nphi = %d\n", n, sf_nphi[n]); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_ntheta = %d", n, sf_ntheta[n]); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_nphi = %d", n, sf_nphi[n]); } } if (verbose) { - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that fully contains this surface: sf_maxreflevel = %d\n", n, sf_maxreflevel[n]); - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that overlaps with this surface : sf_minreflevel = %d\n", n, sf_minreflevel[n]); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that fully contains this surface: sf_maxreflevel = %d", n, sf_maxreflevel[n]); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that overlaps with this surface : sf_minreflevel = %d", n, sf_minreflevel[n]); } // coordinates in the theta direction @@ -249,6 +249,14 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) static bool first_call = true; + int min_reduction_handle = CCTK_ReductionArrayHandle("minimum"); + int max_reduction_handle = CCTK_ReductionArrayHandle("maximum"); + + if (min_reduction_handle < 0) + CCTK_WARN(0, "Cannot get reduction handle for minimum operation."); + if (max_reduction_handle < 0) + CCTK_WARN(0, "Cannot get reduction handle for maximum operation."); + for (int n=0; n<nsurfaces; ++n) { @@ -278,7 +286,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) } if (verbose) { - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: myradius = %f\n", double(my_radius)); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: myradius = %f", double(my_radius)); } // get theta/phi delta-spacing by looking at Cartesian grid @@ -293,27 +301,29 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) bool is_overlapping = false; CCTK_REAL cart_radius = 0; // radius of Cart. grid CCTK_REAL cart_origin[3]; // origin of Cart. grid - + + //CCTK_ORIGIN_SPACE(0) + (cctk_lbnd[0]) * CCTK_DELTA_SPACE(0) - (CCTK_ORIGIN_SPACE(0) + (cctk_lbnd[0] + cctk_lsh[0]) * CCTK_DELTA_SPACE(0)) + // calculate radius of Cartesian grid - cart_origin[0] = (x[l[1]] - x[l[0]]) / 2.0; - cart_origin[1] = (y[l[1]] - y[l[0]]) / 2.0; - cart_origin[2] = (z[l[1]] - z[l[0]]) / 2.0; + cart_origin[0] = 0.5 * cctk_lsh[0] * CCTK_DELTA_SPACE(0); //(x[l[1]] - x[l[0]]) / 2.0; + cart_origin[1] = 0.5 * cctk_lsh[1] * CCTK_DELTA_SPACE(1); //(y[l[1]] - y[l[0]]) / 2.0; + cart_origin[2] = 0.5 * cctk_lsh[2] * CCTK_DELTA_SPACE(2); //(z[l[1]] - z[l[0]]) / 2.0; cart_radius = sqrt(cart_origin[0]*cart_origin[0] + cart_origin[1]*cart_origin[1] + cart_origin[2]*cart_origin[2]); if (verbose) { - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: cart_radius = %f, cart_origin = %f, %f, %f\n", double(cart_radius), double(cart_origin[0]), double(cart_origin[1]), double(cart_origin[2])); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: cart_radius = %f, cart_origin = %f, %f, %f", double(cart_radius), double(cart_origin[0]), double(cart_origin[1]), double(cart_origin[2])); } // set origin of Cartesian grid - cart_origin[0] += x[l[0]]; - cart_origin[1] += y[l[0]]; - cart_origin[2] += z[l[0]]; + cart_origin[0] += CCTK_ORIGIN_SPACE(0) + cctk_lbnd[0] * CCTK_DELTA_SPACE(0); //x[l[0]]; + cart_origin[1] += CCTK_ORIGIN_SPACE(1) + cctk_lbnd[1] * CCTK_DELTA_SPACE(1); //y[l[0]]; + cart_origin[2] += CCTK_ORIGIN_SPACE(2) + cctk_lbnd[2] * CCTK_DELTA_SPACE(2); //z[l[0]]; // norm of displacement vector between origin of Cart. grid and origin of SphericalSurface CCTK_REAL dist = sqrt( (cart_origin[0]-origin_x[n]) * (cart_origin[0]-origin_x[n]) - + (cart_origin[1]-origin_y[n]) * (cart_origin[1]-origin_y[n]) - + (cart_origin[2]-origin_z[n]) * (cart_origin[2]-origin_z[n])); + + (cart_origin[1]-origin_y[n]) * (cart_origin[1]-origin_y[n]) + + (cart_origin[2]-origin_z[n]) * (cart_origin[2]-origin_z[n])); // does our spherical shell overlap with the processor's local grid patch? if (my_radius+cart_radius >= dist && // Cart. grid and SphericalSurface overlap @@ -350,11 +360,13 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) { // find minimum value of "sf_delta_phi/theta" over all processors CCTK_REAL dummy = 0; - MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + //MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + CCTK_ReduceLocScalar(cctkGH, -1, min_reduction_handle, &sf_delta_theta_estimate[n], &dummy, CCTK_VARIABLE_REAL); sf_delta_theta_estimate[n] = dummy; dummy = 0; - MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + //MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + CCTK_ReduceLocScalar(cctkGH, -1, min_reduction_handle, &sf_delta_phi_estimate[n], &dummy, CCTK_VARIABLE_REAL); sf_delta_phi_estimate[n] = dummy; /* @@ -379,7 +391,8 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) // find maximum value of "sf_minreflevel" over all processors int idummy; - MPI_Allreduce(&sf_minreflevel[n], &idummy, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + //MPI_Allreduce(&sf_minreflevel[n], &idummy, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + CCTK_ReduceLocScalar(cctkGH, -1, max_reduction_handle, &sf_minreflevel[n], &idummy, CCTK_VARIABLE_INT); sf_minreflevel[n] = idummy; @@ -398,12 +411,12 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) if (verbose) { //CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: calculated max_radius = %.12f\n", max_radius); - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: buffer_radius = %.12f\n", double (sqrt( np*np*CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) - + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) - + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ))); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: buffer_radius = %.12f", double (sqrt( np*np*CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) + + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) + + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ))); } - + /* unsigned long k[8]; k[0] = CCTK_GFINDEX3D(cctkGH, 0, 0, 0); @@ -414,12 +427,30 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) k[5] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, 0, cctk_lsh[2]-1); k[6] = CCTK_GFINDEX3D(cctkGH, 0, cctk_lsh[1]-1, cctk_lsh[2]-1); k[7] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1); - + */ + + unsigned long kx[8], ky[8], kz[8]; + + kx[0] = 0; ky[0] = 0; kz[0] = 0; + kx[1] = cctk_lsh[0]-1, ky[1] = 0; kz[1] = 0; + kx[2] = 0; ky[2] = cctk_lsh[1]-1; kz[2] = 0; + kx[3] = cctk_lsh[0]-1; ky[3] = cctk_lsh[1]-1; kz[3] = 0; + kx[4] = 0; ky[4] = 0; kz[4] = cctk_lsh[2]-1; + kx[5] = cctk_lsh[0]-1; ky[5] = 0; kz[5] = cctk_lsh[2]-1; + kx[6] = 0; ky[6] = cctk_lsh[1]-1; kz[6] = cctk_lsh[2]-1; + kx[7] = cctk_lsh[0]-1; ky[7] = cctk_lsh[1]-1; kz[7] = cctk_lsh[2]-1; + + int i; bool is_contained = false; for (i=0; i < 8; i++) { - if (max_radius < sqrt(x[k[i]]*x[k[i]] + y[k[i]]*y[k[i]] + z[k[i]]*z[k[i]])) + //if (max_radius < sqrt(x[k[i]]*x[k[i]] + y[k[i]]*y[k[i]] + z[k[i]]*z[k[i]])) + // is_contained = true; + CCTK_REAL xx = CCTK_ORIGIN_SPACE(0) + (cctk_lbnd[0] + kx[i]) * CCTK_DELTA_SPACE(0); + CCTK_REAL yy = CCTK_ORIGIN_SPACE(1) + (cctk_lbnd[1] + ky[i]) * CCTK_DELTA_SPACE(1); + CCTK_REAL zz = CCTK_ORIGIN_SPACE(2) + (cctk_lbnd[2] + kz[i]) * CCTK_DELTA_SPACE(2); + if (max_radius < sqrt(xx*xx + yy*yy + zz*zz)) is_contained = true; } @@ -448,7 +479,8 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) // find minimum value of "sf_maxreflevel" over all processors - MPI_Allreduce(&max_rl, &idummy, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + //MPI_Allreduce(&max_rl, &idummy, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + CCTK_ReduceLocScalar(cctkGH, -1, min_reduction_handle, &max_rl, &idummy, CCTK_VARIABLE_INT); if (idummy != 0 || first_call) sf_maxreflevel[n] = idummy; @@ -462,17 +494,19 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS) { // find maximum value of "sf_delta_phi/theta" over all processors CCTK_REAL dummy = 0; - MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + //MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + CCTK_ReduceLocScalar(cctkGH, -1, max_reduction_handle, &sf_delta_theta_estimate[n], &dummy, CCTK_VARIABLE_REAL); sf_delta_theta_estimate[n] = dummy; dummy = 0; - MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + //MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + CCTK_ReduceLocScalar(cctkGH, -1, max_reduction_handle, &sf_delta_phi_estimate[n], &dummy, CCTK_VARIABLE_REAL); sf_delta_phi_estimate[n] = dummy; } if (verbose) { - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_maxreflevel[%d] = %d, sf_minreflevel[%d] = %d \n", n, sf_maxreflevel[n], n, sf_minreflevel[n]); - CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_delta_theta[%d] = %.6f, sf_delta_phi[%d] = %.6f, \n", n, sf_delta_theta_estimate[n], n, sf_delta_phi_estimate[n]); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_maxreflevel[%d] = %d, sf_minreflevel[%d] = %d", n, sf_maxreflevel[n], n, sf_minreflevel[n]); + CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_delta_theta[%d] = %.6f, sf_delta_phi[%d] = %.6f\n", n, sf_delta_theta_estimate[n], n, sf_delta_phi_estimate[n]); } } |