From a770662a2f1becd61ab39302aa5636b7a9771427 Mon Sep 17 00:00:00 2001 From: baiotti Date: Wed, 12 Oct 2005 15:39:42 +0000 Subject: - make results independent from the number of processors (not yet finished) - change some parameter names, so that, to me, they are more clear (debatable) - add another way of setting the integration surface position - move constant divisions out of the loops - modify info messages git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/ADMMass/trunk@4 54511f98-0e4f-0410-826e-eb8b393f5a1e --- README | 2 +- interface.ccl | 25 ++++++- param.ccl | 36 ++++++---- schedule.ccl | 1 + src/surface_integral.c | 191 ++++++++++++++++++++++++++++++++++--------------- src/volume_integral.c | 82 ++++++++++++++------- 6 files changed, 238 insertions(+), 99 deletions(-) diff --git a/README b/README index 9c4e62b..883b9c8 100644 --- a/README +++ b/README @@ -4,7 +4,7 @@ This thorn can be used to calculate different approximations of the ADM mass on a finite grid. Warning: This thorn is still _very_ experimental. Do not expect it to work, - especially not the way you might think it is. + especially not the way you might think it should. Copyright --------- diff --git a/interface.ccl b/interface.ccl index bb2660d..05bed81 100644 --- a/interface.ccl +++ b/interface.ccl @@ -6,19 +6,30 @@ inherits: ADMBase ADMMacros StaticConformal SpaceMask USES INCLUDE: SpaceMask.h +CCTK_INT FUNCTION GetDomainSpecification \ + (CCTK_INT IN size, \ + CCTK_REAL OUT ARRAY physical_min, \ + CCTK_REAL OUT ARRAY physical_max, \ + CCTK_REAL OUT ARRAY interior_min, \ + CCTK_REAL OUT ARRAY interior_max, \ + CCTK_REAL OUT ARRAY exterior_min, \ + CCTK_REAL OUT ARRAY exterior_max, \ + CCTK_REAL OUT ARRAY spacing) +REQUIRES FUNCTION GetDomainSpecification + CCTK_INT ADMMass_LoopCounterG type = SCALAR { ADMMass_LoopCounter } "ADMMass LoopCounter" -CCTK_REAL ADMMass_Masses[ADMMass_number] type = SCALAR +CCTK_REAL ADMMass_Masses[ADMMass_number] type = SCALAR tags='checkpoint="no"' { ADMMass_SurfaceMass ADMMass_SurfaceMass_Lapse ADMMass_VolumeMass } "ADMMass Scalars" -CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar"' +CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar" checkpoint="no"' { ADMMass_SurfaceMass_GF ADMMass_VolumeMass_pot_x @@ -26,3 +37,13 @@ CCTK_REAL ADMMass_GFs type = GF Timelevels = 1 tags='Prolongation="none" tensort ADMMass_VolumeMass_pot_z ADMMass_VolumeMass_GF } "ADMMass gridfunctions" + +CCTK_REAL ADMMass_box type = scalar tags='checkpoint="no"' +{ + ADMMass_box_x_min + ADMMass_box_x_max + ADMMass_box_y_min + ADMMass_box_y_max + ADMMass_box_z_min + ADMMass_box_z_max +} "Physical coordinates of the surface on which the integral is computed" diff --git a/param.ccl b/param.ccl index 51a01e2..f2337a7 100644 --- a/param.ccl +++ b/param.ccl @@ -6,59 +6,69 @@ CCTK_INT ADMMass_number "number of measurements" STEERABLE=ALWAYS 0: :: "0 or positive" } 1 -CCTK_REAL ADMMass_x_min[100] "x_min" STEERABLE=ALWAYS +CCTK_REAL ADMMass_x_min[100] "x position of the leftmost yz-plane for the integration box" STEERABLE=ALWAYS { : :: "anything" } -100.0 -CCTK_REAL ADMMass_x_max[100] "x_max" STEERABLE=ALWAYS +CCTK_REAL ADMMass_x_max[100] "x position of the righttmost yz-plane for the integration box" STEERABLE=ALWAYS { : :: "anything" } 100.0 -CCTK_REAL ADMMass_y_min[100] "y_min" STEERABLE=ALWAYS +CCTK_REAL ADMMass_y_min[100] "y position of the leftmost xz-plane for the integration box" STEERABLE=ALWAYS { : :: "anything" } -100.0 -CCTK_REAL ADMMass_y_max[100] "y_max" STEERABLE=ALWAYS +CCTK_REAL ADMMass_y_max[100] "y position of the rightmost xz-plane for the integration box" STEERABLE=ALWAYS { : :: "anything" } 100.0 -CCTK_REAL ADMMass_z_min[100] "z_min" STEERABLE=ALWAYS +CCTK_REAL ADMMass_z_min[100] "z position of the leftmost xy-plane for the integration box" STEERABLE=ALWAYS { : :: "anything" } -100.0 -CCTK_REAL ADMMass_z_max[100] "z_max" STEERABLE=ALWAYS +CCTK_REAL ADMMass_z_max[100] "z position of the rightmost xy-plane for the integration box" STEERABLE=ALWAYS { : :: "anything" } 100.0 -CCTK_REAL ADMMass_x_pos[100] "x_pos" STEERABLE=ALWAYS +CCTK_REAL ADMMass_x_pos[100] "x position of the center of the integration box" STEERABLE=ALWAYS { : :: "anything" } 0.0 -CCTK_REAL ADMMass_y_pos[100] "y_pos" STEERABLE=ALWAYS +CCTK_REAL ADMMass_y_pos[100] "y position of the center of the integration box" STEERABLE=ALWAYS { : :: "anything" } 0.0 -CCTK_REAL ADMMass_z_pos[100] "z_pos" STEERABLE=ALWAYS +CCTK_REAL ADMMass_z_pos[100] "z position of the center of the integration box" STEERABLE=ALWAYS { : :: "anything" } 0.0 -CCTK_REAL ADMMass_boundary_distance[100] "distance to boundary" STEERABLE=ALWAYS +CCTK_REAL ADMMass_distance_from_grid_boundary[100] "distance between the grid boundaries and the surface of the integration box" STEERABLE=ALWAYS { : :: "<=0 for disable, positive otherwise" } -1.0 -CCTK_REAL ADMMass_volume_radius[100] "radius of volume extraction" STEERABLE=ALWAYS + +CCTK_REAL ADMMass_surface_distance[100] "distance between the above-defined center of the integration (cubic) box and its surface" STEERABLE=ALWAYS +{ + : :: "<=0 for disable, positive otherwise" +} -1.0 + +CCTK_REAL ADMMass_volume_radius[100] "radius of the sphere inside which the volume integral is computed" STEERABLE=ALWAYS { : :: "<=0 for disable, positive otherwise" } -1.0 -BOOLEAN ADMMass_use_boundary_distance_as_volume_radius "Use ADMMass_boundary_distance insteat of ADMMass_volume_radius" STEERABLE=ALWAYS +BOOLEAN ADMMass_use_surface_distance_as_volume_radius "Use ADMMass_surface_distance instead of ADMMass_volume_radius" STEERABLE=ALWAYS { } "yes" -BOOLEAN ADMMass_Excise_Horizons "Do not calculate the volume integral inside AHs" STEERABLE=ALWAYS +BOOLEAN ADMMass_use_all_volume_as_volume_radius "Use ADMMass_distance_from_grid_boundary instead of ADMMass_volume_radius" STEERABLE=ALWAYS +{ +} "no" + +BOOLEAN ADMMass_Excise_Horizons "Should we calculate the volume integral inside AHs?" STEERABLE=ALWAYS { } "no" diff --git a/schedule.ccl b/schedule.ccl index 9196694..437ebeb 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -13,6 +13,7 @@ schedule ADMMass_InitLoopCounter AT POSTSTEP AFTER OutsideMask_UpdateMask schedule GROUP ADMMass AT POSTSTEP AFTER ADMMass_InitLoopCounter WHILE ADMMass::ADMMass_LoopCounter { + STORAGE:ADMMass_box } "ADMMass loop" schedule ADMMass_Loop IN ADMMass diff --git a/src/surface_integral.c b/src/surface_integral.c index d84347d..7927a06 100644 --- a/src/surface_integral.c +++ b/src/surface_integral.c @@ -1,20 +1,27 @@ #include +#include +#include +#include #include #include #include -CCTK_INT find_closest(const cGH *cctkGH, CCTK_INT *cctk_lsh, - CCTK_REAL *x, CCTK_REAL x_min, CCTK_INT dir) +#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 + +CCTK_INT find_closest(const cGH *cctkGH, CCTK_INT *cctk_lsh, CCTK_REAL *cctk_delta_space, + CCTK_INT ghost, CCTK_REAL *coord, CCTK_REAL coord_min, CCTK_INT dir) { CCTK_INT i, ijk, min_i = -1; CCTK_REAL min = 1.e100; - for(i=3; i 0.0) - { - x_min = ADMMass_x_pos[*ADMMass_LoopCounter] - - ADMMass_boundary_distance[*ADMMass_LoopCounter]; - x_max = ADMMass_x_pos[*ADMMass_LoopCounter] + - ADMMass_boundary_distance[*ADMMass_LoopCounter]; - y_min = ADMMass_y_pos[*ADMMass_LoopCounter] - - ADMMass_boundary_distance[*ADMMass_LoopCounter]; - y_max = ADMMass_y_pos[*ADMMass_LoopCounter] + - ADMMass_boundary_distance[*ADMMass_LoopCounter]; - z_min = ADMMass_z_pos[*ADMMass_LoopCounter] - - ADMMass_boundary_distance[*ADMMass_LoopCounter]; - z_max = ADMMass_z_pos[*ADMMass_LoopCounter] + - ADMMass_boundary_distance[*ADMMass_LoopCounter]; - } + if (ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter] > 0.0) + { + if (strcmp(CCTK_ParameterValString("type", "cartgrid3d"), "coordbase")) + { + CCTK_WARN(0,"this thorns used with the ADMMass_distance_from_grid_boundar parameter requires to set coordinates through coordbase."); + } + + /* Find the physical coordinates of the boundaries */ + ierr = GetDomainSpecification + (3, physical_min, physical_max, interior_min, interior_max, + exterior_min, exterior_max, spacing); + if (ierr) + CCTK_VWarn ( 0, __LINE__, __FILE__, "CartGrid3D", "error returned from function GetDomainSpecification"); + *ADMMass_box_x_min = physical_min[0] + ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; + *ADMMass_box_x_max = physical_max[0] - ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; + *ADMMass_box_y_min = physical_min[1] + ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; + *ADMMass_box_y_max = physical_max[1] - ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; + *ADMMass_box_z_min = physical_min[2] + ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; + *ADMMass_box_z_max = physical_max[2] - ADMMass_distance_from_grid_boundary[*ADMMass_LoopCounter]; + } + else if (ADMMass_surface_distance[*ADMMass_LoopCounter] > 0.0) + { + *ADMMass_box_x_min = ADMMass_x_pos[*ADMMass_LoopCounter] - + ADMMass_surface_distance[*ADMMass_LoopCounter]; + *ADMMass_box_x_max = ADMMass_x_pos[*ADMMass_LoopCounter] + + ADMMass_surface_distance[*ADMMass_LoopCounter]; + *ADMMass_box_y_min = ADMMass_y_pos[*ADMMass_LoopCounter] - + ADMMass_surface_distance[*ADMMass_LoopCounter]; + *ADMMass_box_y_max = ADMMass_y_pos[*ADMMass_LoopCounter] + + ADMMass_surface_distance[*ADMMass_LoopCounter]; + *ADMMass_box_z_min = ADMMass_z_pos[*ADMMass_LoopCounter] - + ADMMass_surface_distance[*ADMMass_LoopCounter]; + *ADMMass_box_z_max = ADMMass_z_pos[*ADMMass_LoopCounter] + + ADMMass_surface_distance[*ADMMass_LoopCounter]; + } else - { - x_min = ADMMass_x_min[*ADMMass_LoopCounter]; - x_max = ADMMass_x_max[*ADMMass_LoopCounter]; - y_min = ADMMass_y_min[*ADMMass_LoopCounter]; - y_max = ADMMass_y_max[*ADMMass_LoopCounter]; - z_min = ADMMass_z_min[*ADMMass_LoopCounter]; - z_max = ADMMass_z_max[*ADMMass_LoopCounter]; - } - x_min_i = find_closest(cctkGH, cctk_lsh, x, x_min, 0); - x_max_i = find_closest(cctkGH, cctk_lsh, x, x_max, 0); - y_min_j = find_closest(cctkGH, cctk_lsh, y, y_min, 1); - y_max_j = find_closest(cctkGH, cctk_lsh, y, y_max, 1); - z_min_k = find_closest(cctkGH, cctk_lsh, z, z_min, 2); - z_max_k = find_closest(cctkGH, cctk_lsh, z, z_max, 2); + { + *ADMMass_box_x_min = ADMMass_x_min[*ADMMass_LoopCounter]; + *ADMMass_box_x_max = ADMMass_x_max[*ADMMass_LoopCounter]; + *ADMMass_box_y_min = ADMMass_y_min[*ADMMass_LoopCounter]; + *ADMMass_box_y_max = ADMMass_y_max[*ADMMass_LoopCounter]; + *ADMMass_box_z_min = ADMMass_z_min[*ADMMass_LoopCounter]; + *ADMMass_box_z_max = ADMMass_z_max[*ADMMass_LoopCounter]; + } + + if (CCTK_IsThornActive("PUGH")) + { + const void *ghost_ptr = CCTK_ParameterGet("ghost_size","PUGH",NULL); + assert( ghost_ptr != NULL ); + ghost = *(const int *)ghost_ptr; + } + else /* carpet */ + { + const void *ghost_ptr = CCTK_ParameterGet("ghost_size","Carpet",NULL); + assert( ghost_ptr != NULL ); + ghost = *(const int *)ghost_ptr; + } + + ijk = CCTK_GFINDEX3D(cctkGH, cctk_ubnd[0]-cctk_lbnd[0], cctk_ubnd[1]-cctk_lbnd[1],cctk_ubnd[2]-cctk_lbnd[2] ); + + if ( (x[0] < *ADMMass_box_x_min) && (*ADMMass_box_x_min < x[ijk]) ) + x_min_i = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, x, *ADMMass_box_x_min, 0); + + if ( (x[0] < *ADMMass_box_x_max) && (*ADMMass_box_x_max < x[ijk]) ) + x_max_i = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, x, *ADMMass_box_x_max, 0); + + if ( (y[0] < *ADMMass_box_y_min) && (*ADMMass_box_y_min < y[ijk]) ) + y_min_j = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, y, *ADMMass_box_y_min, 1); + + if ( (y[0] < *ADMMass_box_y_max) && (*ADMMass_box_y_max < y[ijk]) ) + y_max_j = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, y, *ADMMass_box_y_max, 1); + + if ( (z[0] < *ADMMass_box_z_min) && (*ADMMass_box_z_min < z[ijk]) ) + z_min_k = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, z, *ADMMass_box_z_min, 2); + + if ( (z[0] < *ADMMass_box_z_max) && (*ADMMass_box_z_max < z[ijk]) ) + z_max_k = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, z, *ADMMass_box_z_max, 2); #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h" #include "CactusEinstein/ADMMacros/src/macro/DG_declare.h" #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h" - for(i=3; i= x_min_i) && + if ((i >= x_min_i) && (i <= x_max_i) && (j >= y_min_j) && (j <= y_max_j) && @@ -149,11 +218,12 @@ void ADMMass_Surface(CCTK_ARGUMENTS) ADMMass_SurfaceMass_GF[ijk] += u[ti][tj] * u[tk][tl] * ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk]; - ADMMass_SurfaceMass_GF[ijk] *= sqrt(DETG_DETG) / 16 / 3.14159265; - local_sum += ADMMass_SurfaceMass_GF[ijk]; + ADMMass_SurfaceMass_GF[ijk] *= sqrt(DETG_DETG); + /* local_sum += ADMMass_SurfaceMass_GF[ijk];*/ } } } + #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/DG_undefine.h" #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_undefine.h" @@ -171,18 +241,23 @@ void ADMMass_Surface_Global(CCTK_ARGUMENTS) reduction_handle = CCTK_ReductionHandle("sum"); if (reduction_handle < 0) - CCTK_WARN(0, "Unable to ge reduction handle."); + CCTK_WARN(0, "Unable to get reduction handle."); if (CCTK_Reduce(cctkGH, -1, reduction_handle, 1, CCTK_VARIABLE_REAL, &ADMMass_SurfaceMass[*ADMMass_LoopCounter], 1, CCTK_VarIndex("ADMMass::ADMMass_SurfaceMass_GF"))) CCTK_WARN(0, "Error while reducing ADMMass_SurfaceMass_GF"); - + + ADMMass_SurfaceMass[*ADMMass_LoopCounter] /= 16.0*PI; + CCTK_VInfo(CCTK_THORNSTRING, - "ADM mass, surface, detector %d: %g %g\n", *ADMMass_LoopCounter, - ADMMass_boundary_distance[*ADMMass_LoopCounter], - ADMMass_SurfaceMass[*ADMMass_LoopCounter]); + "detector %d: ADM mass %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n", + *ADMMass_LoopCounter, + ADMMass_SurfaceMass[*ADMMass_LoopCounter], + *ADMMass_box_x_min,*ADMMass_box_x_max, + *ADMMass_box_y_min,*ADMMass_box_y_max, + *ADMMass_box_z_min,*ADMMass_box_z_max); } void ADMMass_Surface_Lapse(CCTK_ARGUMENTS) @@ -191,15 +266,13 @@ void ADMMass_Surface_Lapse(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT i,j,k, ijk; - CCTK_REAL local_sum = 0.0; - + for(i=3; i #include #include +#include #include #include @@ -16,10 +16,10 @@ void ADMMass_Volume(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT i,j,k, ijk; + CCTK_INT i,j,k, ijk, ghost; CCTK_REAL detg, idetg; CCTK_REAL u[3][3], dg[3][3][3]; - CCTK_REAL local_sum = 0.0; + /* CCTK_REAL local_sum = 0.0;*/ CCTK_REAL radius; CCTK_INT mask_descriptor, state_descriptor_outside; @@ -32,26 +32,40 @@ void ADMMass_Volume(CCTK_ARGUMENTS) const CCTK_INT dj = cctk_lsh[0]; const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1]; - if (ADMMass_use_boundary_distance_as_volume_radius && + if (ADMMass_use_surface_distance_as_volume_radius && (ADMMass_volume_radius[*ADMMass_LoopCounter] < 0.0)) - radius = ADMMass_boundary_distance[*ADMMass_LoopCounter]; + radius = ADMMass_surface_distance[*ADMMass_LoopCounter]; else radius = ADMMass_volume_radius[*ADMMass_LoopCounter]; - if (radius <= 0.0) + if ((radius <= 0.0)&&(!ADMMass_use_all_volume_as_volume_radius)) { CCTK_WARN(1, "radius < 0 / not set, not calculating " "the volume integral to get the ADM mass."); return; } + if (CCTK_IsThornActive("PUGH")) + { + const void *ghost_ptr = CCTK_ParameterGet("ghost_size","PUGH",NULL); + assert( ghost_ptr != NULL ); + ghost = *(const int *)ghost_ptr; + } + else /* carpet */ + { + const void *ghost_ptr = CCTK_ParameterGet("ghost_size","Carpet",NULL); + assert( ghost_ptr != NULL ); + ghost = *(const int *)ghost_ptr; + } + #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h" #include "CactusEinstein/ADMMacros/src/macro/DG_declare.h" #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h" - for(i=3; i