diff --git a/AHFinderDirect/schedule.ccl b/AHFinderDirect/schedule.ccl index e3d0de09..2dfd5b2a 100644 --- a/AHFinderDirect/schedule.ccl +++ b/AHFinderDirect/schedule.ccl @@ -6,12 +6,23 @@ storage: ahmask[1] # # setup # +schedule AHFinderDirect_init at CCTK_WRAGH +{ + lang: C + options: global + WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius +} "initialise the horizon grid arrays" + schedule AHFinderDirect_setup at CCTK_BASEGRID \ after SpatialCoordinates { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags + READS: AHFinderDirect::ah_origin, ah_radius WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius } "setup data structures" schedule AHFinderDirect_import_mask at CCTK_POSTREGRIDINITIAL \ @@ -38,6 +49,7 @@ schedule AHFinderDirect_recover at CCTK_POST_RECOVER_VARIABLES { lang: C options: global + READS: AHFinderDirect::ah_flags, ah_origin, ah_radius } "import horizon data from Cactus variables" # @@ -49,6 +61,7 @@ if (run_at_CCTK_ANALYSIS != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags } "find apparent horizon(s) after this time step" } if (run_at_CCTK_POSTSTEP != 0) @@ -57,6 +70,7 @@ if (run_at_CCTK_POSTSTEP != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags } "find apparent horizon(s) after this time step" } if (run_at_CCTK_POSTINITIAL != 0) @@ -65,6 +79,7 @@ if (run_at_CCTK_POSTINITIAL != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags } "find apparent horizon(s) after this time step" } if (run_at_CCTK_POSTPOSTINITIAL != 0) @@ -73,6 +88,7 @@ if (run_at_CCTK_POSTPOSTINITIAL != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags } "find apparent horizon(s) after this time step" } if (run_at_CCTK_POST_RECOVER_VARIABLES != 0) @@ -82,6 +98,7 @@ if (run_at_CCTK_POST_RECOVER_VARIABLES != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags } "find apparent horizon(s) after this time step" } @@ -116,6 +133,18 @@ if (run_at_CCTK_ANALYSIS != 0) { lang: C options: global + READS: SphericalSurface::sf_shape_descriptors + READS: SphericalSurface::sf_coordinate_descriptors + READS: SphericalSurface::sf_active + READS: SphericalSurface::sf_valid + READS: SphericalSurface::sf_info + READS: SphericalSurface::sf_origin + READS: SphericalSurface::sf_radius + WRITES: SphericalSurface::sf_active + WRITES: SphericalSurface::sf_valid + WRITES: SphericalSurface::sf_info + WRITES: SphericalSurface::sf_origin + WRITES: SphericalSurface::sf_radius } "store apparent horizon(s) into spherical surface(s)" schedule AHFinderDirect_save at CCTK_ANALYSIS \ @@ -123,6 +152,9 @@ if (run_at_CCTK_ANALYSIS != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius } "save apparent horizon(s) into Cactus variables" if (which_horizon_to_announce_centroid != 0) @@ -158,6 +190,18 @@ if (run_at_CCTK_POSTSTEP != 0) { lang: C options: global + READS: SphericalSurface::sf_shape_descriptors + READS: SphericalSurface::sf_coordinate_descriptors + READS: SphericalSurface::sf_active + READS: SphericalSurface::sf_valid + READS: SphericalSurface::sf_info + READS: SphericalSurface::sf_origin + READS: SphericalSurface::sf_radius + WRITES: SphericalSurface::sf_active + WRITES: SphericalSurface::sf_valid + WRITES: SphericalSurface::sf_info + WRITES: SphericalSurface::sf_origin + WRITES: SphericalSurface::sf_radius } "store apparent horizon(s) into spherical surface(s)" schedule AHFinderDirect_save at CCTK_POSTSTEP \ @@ -165,6 +209,9 @@ if (run_at_CCTK_POSTSTEP != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius } "save apparent horizon(s) into Cactus variables" if (which_horizon_to_announce_centroid != 0) @@ -200,6 +247,18 @@ if (run_at_CCTK_POSTINITIAL != 0) { lang: C options: global + READS: SphericalSurface::sf_shape_descriptors + READS: SphericalSurface::sf_coordinate_descriptors + READS: SphericalSurface::sf_active + READS: SphericalSurface::sf_valid + READS: SphericalSurface::sf_info + READS: SphericalSurface::sf_origin + READS: SphericalSurface::sf_radius + WRITES: SphericalSurface::sf_active + WRITES: SphericalSurface::sf_valid + WRITES: SphericalSurface::sf_info + WRITES: SphericalSurface::sf_origin + WRITES: SphericalSurface::sf_radius } "store apparent horizon(s) into spherical surface(s)" schedule AHFinderDirect_save at CCTK_POSTINITIAL \ @@ -207,6 +266,9 @@ if (run_at_CCTK_POSTINITIAL != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius } "save apparent horizon(s) into Cactus variables" if (which_horizon_to_announce_centroid != 0) @@ -243,6 +305,18 @@ if (run_at_CCTK_POSTPOSTINITIAL != 0) { lang: C options: global + READS: SphericalSurface::sf_shape_descriptors + READS: SphericalSurface::sf_coordinate_descriptors + READS: SphericalSurface::sf_active + READS: SphericalSurface::sf_valid + READS: SphericalSurface::sf_info + READS: SphericalSurface::sf_origin + READS: SphericalSurface::sf_radius + WRITES: SphericalSurface::sf_active + WRITES: SphericalSurface::sf_valid + WRITES: SphericalSurface::sf_info + WRITES: SphericalSurface::sf_origin + WRITES: SphericalSurface::sf_radius } "store apparent horizon(s) into spherical surface(s)" schedule AHFinderDirect_save at CCTK_POSTPOSTINITIAL \ @@ -250,6 +324,9 @@ if (run_at_CCTK_POSTPOSTINITIAL != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius } "save apparent horizon(s) into Cactus variables" if (which_horizon_to_announce_centroid != 0) @@ -285,6 +362,18 @@ if (run_at_CCTK_POST_RECOVER_VARIABLES != 0) { lang: C options: global + READS: SphericalSurface::sf_shape_descriptors + READS: SphericalSurface::sf_coordinate_descriptors + READS: SphericalSurface::sf_active + READS: SphericalSurface::sf_valid + READS: SphericalSurface::sf_info + READS: SphericalSurface::sf_origin + READS: SphericalSurface::sf_radius + WRITES: SphericalSurface::sf_active + WRITES: SphericalSurface::sf_valid + WRITES: SphericalSurface::sf_info + WRITES: SphericalSurface::sf_origin + WRITES: SphericalSurface::sf_radius } "store apparent horizon(s) into spherical surface(s)" schedule AHFinderDirect_save at CCTK_POST_RECOVER_VARIABLES \ @@ -292,6 +381,9 @@ if (run_at_CCTK_POST_RECOVER_VARIABLES != 0) { lang: C options: global + READS: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_centroid, ah_flags + WRITES: AHFinderDirect::ah_origin, ah_radius } "save apparent horizon(s) into Cactus variables" if (which_horizon_to_announce_centroid != 0) diff --git a/AHFinderDirect/src/driver/setup.cc b/AHFinderDirect/src/driver/setup.cc index 0d12033f..db532735 100644 --- a/AHFinderDirect/src/driver/setup.cc +++ b/AHFinderDirect/src/driver/setup.cc @@ -223,6 +223,62 @@ extern struct state state; // This function is called by the Cactus scheduler to set up all our // persistent data structures. (These are stored in struct state .) // +// +// AHFinderDirect_setup() below declares ah_radius, ah_origin, ah_centroid and +// ah_flags as output, but it only fills them the first time it runs. The +// static already_ran guard makes every later call a no-op. CarpetX +// re-traverses CCTK_BASEGRID after each regrid, and it poisons a routine's +// write-only outputs immediately before calling it, so those later calls used +// to leave the whole of ah_radius (and the centroids) as nans and abort in +// valid.cxx. The BASEGRID entry is therefore declared as READS *and* WRITES, +// which is what it really does, a partial, preserving update, and CarpetX +// then leaves the arrays alone. That only works if they are already valid the +// first time AHFinderDirect_setup runs, which is what this routine is for. +// +// The values match what AHFinderDirect_setup would write for a horizon that +// has not been found yet. +// +extern "C" + void AHFinderDirect_init(CCTK_ARGUMENTS) +{ +DECLARE_CCTK_ARGUMENTS_AHFinderDirect_init +DECLARE_CCTK_PARAMETERS + +for (int n = 0; n < N_horizons; ++n) { + ah_origin_x[n] = 0.0; + ah_origin_y[n] = 0.0; + ah_origin_z[n] = 0.0; + + ah_centroid_x[n] = 0.0; + ah_centroid_y[n] = 0.0; + ah_centroid_z[n] = 0.0; + ah_centroid_t[n] = 0.0; + ah_centroid_x_p[n] = 0.0; + ah_centroid_y_p[n] = 0.0; + ah_centroid_z_p[n] = 0.0; + ah_centroid_t_p[n] = 0.0; + + ah_initial_find_flag[n] = 0; + ah_really_initial_find_flag[n] = 0; + ah_search_flag[n] = 0; + ah_found_flag[n] = 0; + ah_centroid_valid[n] = 0; + ah_centroid_valid_p[n] = 0; + ah_centroid_iteration[n] = -1; + ah_centroid_iteration_p[n] = -1; + + // the whole array is poisoned, so the whole array has to be set, not only + // the (N_zones_per_right_angle+1)^2 * N_patches points that are in use + const int nzones = max_N_zones_per_right_angle + 1; + for (int pn = 0; pn < 6; ++pn) + for (int j = 0; j < nzones; ++j) + for (int i = 0; i < nzones; ++i) + ah_radius[i + nzones * (j + nzones * (pn + 6 * n))] = 0.0; + } +} + +//****************************************************************************** + extern "C" void AHFinderDirect_setup(CCTK_ARGUMENTS) { diff --git a/PunctureTracker/src/puncture_tracker.cxx b/PunctureTracker/src/puncture_tracker.cxx index b76f9020..5931d081 100644 --- a/PunctureTracker/src/puncture_tracker.cxx +++ b/PunctureTracker/src/puncture_tracker.cxx @@ -7,6 +7,7 @@ #include +#include #include #include #include @@ -19,6 +20,72 @@ static PunctureContainer *g_punctures = nullptr; const int max_num_tracked = 10; +// `BoxInBox::positions` is a vector grid scalar with one element per +// refinement region, and both `PunctureTracker_Setup` and +// `PunctureTracker_Track` declare `WRITES: BoxInBox::positions`. A WRITES +// clause names the whole vector group, so CarpetX poisons every region +// immediately before the routine and checks every region immediately after it. +// Writing only the tracked punctures leaves the remaining regions as nans +// and aborts in `valid.cxx`. Write all of them: the tracked punctures where +// there is one, and BoxInBox's own parameters (i.e. what `BoxInBox_Init` put +// there) everywhere else. This also caps the write at the number of regions +// that actually exist, which the bare `n < nPunctures` loops did not do. + +namespace { + +const int max_num_boxes = 3; // BoxInBox::max_num_regions + +int get_num_boxes() { + const int gi = CCTK_GroupIndex("BoxInBox::positions"); + if (gi < 0) + CCTK_VERROR("Could not find group BoxInBox::positions"); + cGroup gdata; + const int ierr = CCTK_GroupData(gi, &gdata); + if (ierr != 0) + CCTK_VERROR("Could not query group BoxInBox::positions"); + if (gdata.vectorlength != max_num_boxes) + CCTK_VERROR("BoxInBox::positions has %d regions, but PunctureTracker knows " + "how to restore %d of them", + gdata.vectorlength, max_num_boxes); + return gdata.vectorlength; +} + +// BoxInBox's `position_[xyz]_N` are private parameters, so they cannot be +// pulled in with `SHARES: BoxInBox` / `USES`. We read them at runtime instead. +CCTK_REAL get_box_position_param(const char *const component, const int box) { + char name[32]; + snprintf(name, sizeof name, "position_%s_%d", component, box + 1); + int type = -1; + const void *const value = CCTK_ParameterGet(name, "BoxInBox", &type); + if (!value) + CCTK_VERROR("Could not read parameter BoxInBox::%s", name); + if (type != PARAMETER_REAL) + CCTK_VERROR("Parameter BoxInBox::%s is not a real", name); + return *static_cast(value); +} + +void set_box_positions(CCTK_REAL *restrict const position_x, + CCTK_REAL *restrict const position_y, + CCTK_REAL *restrict const position_z, + const int num_tracked_boxes, + const std::array, Loop::dim> + &location) { + const int num_boxes = get_num_boxes(); + for (int n = 0; n < num_boxes; ++n) { + if (n < num_tracked_boxes) { + position_x[n] = location[0][n]; + position_y[n] = location[1][n]; + position_z[n] = location[2][n]; + } else { + position_x[n] = get_box_position_param("x", n); + position_y[n] = get_box_position_param("y", n); + position_z[n] = get_box_position_param("z", n); + } + } +} + +} // namespace + extern "C" void PunctureTracker_Init(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS_PunctureTracker_Init; DECLARE_CCTK_PARAMETERS; @@ -79,17 +146,17 @@ extern "C" void PunctureTracker_Setup(CCTK_ARGUMENTS) { g_punctures->setNumPunctures(); assert(g_punctures->getNumPunctures() == nPunctures); - // enabled if refinement regions should follow the punctures - if (track_boxes) { - const std::array, Loop::dim> &location = - g_punctures->getLocation(); - for (int n = 0; n < nPunctures; ++n) { + // enabled if refinement regions should follow the punctures. The regions + // that do not follow a puncture are still written, see `set_box_positions` + const std::array, Loop::dim> &location = + g_punctures->getLocation(); + const int num_tracked_boxes = + track_boxes ? std::min(nPunctures, max_num_boxes) : 0; + if (verbose) + for (int n = 0; n < num_tracked_boxes; ++n) CCTK_VINFO("Writing punc coords to box %d.", n); - position_x[n] = location[0][n]; - position_y[n] = location[1][n]; - position_z[n] = location[2][n]; - } - } + set_box_positions(position_x, position_y, position_z, num_tracked_boxes, + location); } extern "C" void PunctureTracker_Finalize(CCTK_ARGUMENTS) { @@ -163,25 +230,26 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { // Broadcast result: 3 components for location, 3 components for velocity g_punctures->broadcast(CCTK_PASS_CTOC); - // Write to pt_loc_foo and pt_vel_foo - for (int i = 0; i < nPunctures; ++i) { - pt_loc_t[i] = time[i]; - pt_loc_x[i] = location[0][i]; - pt_loc_y[i] = location[1][i]; - pt_loc_z[i] = location[2][i]; - pt_vel_t[i] = time[i]; - pt_vel_x[i] = velocity[0][i]; - pt_vel_y[i] = velocity[1][i]; - pt_vel_z[i] = velocity[2][i]; + // Write to pt_loc_foo and pt_vel_foo. `pt_loc` and `pt_vel` are + // max_num_tracked-element vector groups and this routine declares both as + // WRITES, so the driver poisons all of them before the call and checks all + // of them after it. The slots beyond the tracked punctures have to be given + // the same zeros PunctureTracker_Init writes for an untracked puncture. + for (int i = 0; i < max_num_tracked; ++i) { + const bool tracked = i < nPunctures; + pt_loc_t[i] = tracked ? time[i] : 0.0; + pt_loc_x[i] = tracked ? location[0][i] : 0.0; + pt_loc_y[i] = tracked ? location[1][i] : 0.0; + pt_loc_z[i] = tracked ? location[2][i] : 0.0; + pt_vel_t[i] = tracked ? time[i] : 0.0; + pt_vel_x[i] = tracked ? velocity[0][i] : 0.0; + pt_vel_y[i] = tracked ? velocity[1][i] : 0.0; + pt_vel_z[i] = tracked ? velocity[2][i] : 0.0; } - if (track_boxes) { - for (int i = 0; i < nPunctures; ++i) { - position_x[i] = location[0][i]; - position_y[i] = location[1][i]; - position_z[i] = location[2][i]; - } - } + set_box_positions(position_x, position_y, position_z, + track_boxes ? std::min(nPunctures, max_num_boxes) : 0, + location); } } // namespace PunctureTracker diff --git a/SphericalSurface/schedule.ccl b/SphericalSurface/schedule.ccl index 8a1bc770..9447e70b 100644 --- a/SphericalSurface/schedule.ccl +++ b/SphericalSurface/schedule.ccl @@ -53,10 +53,20 @@ SCHEDULE GROUP SphericalSurface_HasBeenSet AT basegrid +# Same clauses as the basegrid entry above. Without them CarpetX nulls +# cctkGH->data[] for every undeclared variable, and this routine segfaults on +# the first sf_ write. SCHEDULE SphericalSurface_Set AT poststep BEFORE SphericalSurface_HasBeenSet { LANG: C OPTIONS: global + READS: SphericalSurface::sf_shape_descriptors + READS: SphericalSurface::sf_coordinate_descriptors + WRITES: SphericalSurface::sf_active(everywhere) + WRITES: SphericalSurface::sf_valid(everywhere) + WRITES: SphericalSurface::sf_info(everywhere) + WRITES: SphericalSurface::sf_origin(everywhere) + WRITES: SphericalSurface::sf_radius(everywhere) } "Set surface radii" SCHEDULE GROUP SphericalSurface_HasBeenSet AT poststep diff --git a/SphericalSurface/src/radius.c b/SphericalSurface/src/radius.c index aeafa023..c79bb069 100644 --- a/SphericalSurface/src/radius.c +++ b/SphericalSurface/src/radius.c @@ -146,6 +146,55 @@ void SphericalSurface_Set (CCTK_ARGUMENTS) sf_origin_y[n] = origin_y[n]; sf_origin_z[n] = origin_z[n]; + } else { + + //This routine declares all of the surface variables as unconditional + // output, so the driver is free to poison them before every call. + // Surfaces that are neither spherical nor elliptic are not set here, so + // give them the same "no surface" state that SphericalSurface_Setup + // uses. + + sf_active[n] = 0; + sf_valid[n] = 0; + + sf_area[n] = 0.0; + + sf_mean_radius[n] = 0.0; + + sf_centroid_x[n] = 0.0; + sf_centroid_y[n] = 0.0; + sf_centroid_z[n] = 0.0; + + sf_quadrupole_xx[n] = 0.0; + sf_quadrupole_xy[n] = 0.0; + sf_quadrupole_xz[n] = 0.0; + sf_quadrupole_yy[n] = 0.0; + sf_quadrupole_yz[n] = 0.0; + sf_quadrupole_zz[n] = 0.0; + + sf_min_radius[n] = 0.0; + sf_max_radius[n] = 0.0; + + sf_min_x[n] = 0.0; + sf_min_y[n] = 0.0; + sf_min_z[n] = 0.0; + sf_max_x[n] = 0.0; + sf_max_y[n] = 0.0; + sf_max_z[n] = 0.0; + + // the whole array is poisoned, so the whole array has to be set, not + // only the sf_ntheta * sf_nphi points that are in use + for (j=0; j