From f088625096574f7347e3792b5dc54d51edb38391 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Thu, 27 Jun 2024 13:24:30 +0000 Subject: [PATCH 01/14] PunctureTracker: remove unsued files --- PunctureTracker/src/make.code.defn | 2 +- PunctureTracker/src/paramcheck.cc | 21 --------------------- 2 files changed, 1 insertion(+), 22 deletions(-) delete mode 100644 PunctureTracker/src/paramcheck.cc diff --git a/PunctureTracker/src/make.code.defn b/PunctureTracker/src/make.code.defn index 942e3e9c..a980b052 100644 --- a/PunctureTracker/src/make.code.defn +++ b/PunctureTracker/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn PunctureTracker -*-Makefile-*- # Source files in this directory -SRCS = puncture_tracker.cc paramcheck.cc +SRCS = puncture_tracker.cc # Subdirectories containing source files SUBDIRS = diff --git a/PunctureTracker/src/paramcheck.cc b/PunctureTracker/src/paramcheck.cc deleted file mode 100644 index 1a82b43b..00000000 --- a/PunctureTracker/src/paramcheck.cc +++ /dev/null @@ -1,21 +0,0 @@ -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Functions.h" -#include "cctk_Parameters.h" - -extern "C" void PunctureTracker_ParamCheck(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - // do this only on processor 0 - if (CCTK_MyProc(cctkGH) != 0) - return; - - // for (int i = 0; i < 10; ++i) { - // if (which_surface_to_store_info[i] != -1) { - // if (which_surface_to_store_info[i] >= nsurfaces) - // CCTK_PARAMWARN("You assigned a greater surface index than there are " - // "spherical surfaces!"); - // } - // } -} From 2d76c71289dc7bf133f94289bc6233616ed399ff Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Thu, 27 Jun 2024 16:27:58 +0000 Subject: [PATCH 02/14] PunctureTracker: beautify PunctureTracker_Track --- PunctureTracker/src/puncture_tracker.cc | 108 +++++++++--------------- 1 file changed, 38 insertions(+), 70 deletions(-) diff --git a/PunctureTracker/src/puncture_tracker.cc b/PunctureTracker/src/puncture_tracker.cc index 2ad9e9a5..35386114 100644 --- a/PunctureTracker/src/puncture_tracker.cc +++ b/PunctureTracker/src/puncture_tracker.cc @@ -100,84 +100,60 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { } // Interpolate - - // Dimensions - const int dim = 3; - - // Interpolation operator - const int operator_handle = - CCTK_InterpHandle("CarpetX"); - if (operator_handle < 0) { - CCTK_WARN(CCTK_WARN_ALERT, "Can't get interpolation handle"); - return; - } // TODO: Not used by CarpetX Interpolator - - // Interpolation parameter table - int ierr; - int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); - if (param_table_handle < 0) - CCTK_VERROR("Can't create parameter table: %d", param_table_handle); - if ((ierr = Util_TableSetInt(param_table_handle, interp_order, "order")) < 0) - CCTK_VERROR("Can't set order in parameter table: %d", ierr); - { - // Interpolation coordinate system: Not used in CarpetX_DriverInterpolate - const int coordsys_handle = 0; - // const int coordsys_handle = CCTK_CoordSystemHandle("cart3d"); - // if (coordsys_handle < 0) { - // CCTK_WARN(CCTK_WARN_ALERT, "Can't get coordinate system handle"); - // goto label_free_param_table; - // } - // Only processor 0 interpolates - const int num_points = CCTK_MyProc(cctkGH) == 0 ? max_num_tracked : 0; + const CCTK_INT nPoints = CCTK_MyProc(cctkGH) == 0 ? max_num_tracked : 0; // Interpolation coordinates - assert(dim == 3); - CCTK_POINTER_TO_CONST interp_coords[dim]; - interp_coords[0] = pt_loc_x; - interp_coords[1] = pt_loc_y; - interp_coords[2] = pt_loc_z; - - // Number of interpolation variables - int const num_vars = 3; + const void *interpCoords[Loop::dim] = {pt_loc_x, pt_loc_y, pt_loc_z}; // Interpolated variables - assert(num_vars == 3); - int input_array_indices[3]; - input_array_indices[0] = CCTK_VarIndex("ADMBaseX::betax"); - input_array_indices[1] = CCTK_VarIndex("ADMBaseX::betay"); - input_array_indices[2] = CCTK_VarIndex("ADMBaseX::betaz"); - - // Interpolation result types: Not used by CarpetX DriverInterp - assert(num_vars == 3); - CCTK_INT output_array_type_codes[3]; - output_array_type_codes[0] = CCTK_VARIABLE_REAL; - output_array_type_codes[1] = CCTK_VARIABLE_REAL; - output_array_type_codes[2] = CCTK_VARIABLE_REAL; + const CCTK_INT nInputArrays = 3; + const CCTK_INT inputArrayIndices[3] = {CCTK_VarIndex("ADMBaseX::betax"), + CCTK_VarIndex("ADMBaseX::betay"), + CCTK_VarIndex("ADMBaseX::betaz")}; // Interpolation result CCTK_REAL pt_betax[max_num_tracked]; CCTK_REAL pt_betay[max_num_tracked]; CCTK_REAL pt_betaz[max_num_tracked]; + CCTK_POINTER outputArrays[3] = {pt_betax, pt_betay, pt_betaz}; - assert(num_vars == 3); - CCTK_POINTER output_arrays[3]; - output_arrays[0] = pt_betax; - output_arrays[1] = pt_betay; - output_arrays[2] = pt_betaz; + /* DriverInterpolate arguments that aren't currently used */ + const int coordSystemHandle = 0; + const CCTK_INT interpCoordsTypeCode = 0; + const CCTK_INT outputArrayTypes[1] = {0}; + + const int interpHandle = CCTK_InterpHandle("CarpetX"); + if (interpHandle < 0) { + CCTK_WARN(CCTK_WARN_ALERT, "Can't get interpolation handle"); + return; + } - // Interpolate int ierr; - ierr = DriverInterpolate( - cctkGH, dim, operator_handle, param_table_handle, coordsys_handle, - num_points, CCTK_VARIABLE_REAL, interp_coords, num_vars, - input_array_indices, num_vars, output_array_type_codes, output_arrays); + + int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); + if (paramTableHandle < 0) { + CCTK_VERROR("Can't create parameter table: %d", paramTableHandle); + } + + if ((ierr = Util_TableSetInt(paramTableHandle, interp_order, "order")) < + 0) { + CCTK_VERROR("Can't set order in parameter table: %d", ierr); + } + + // Interpolate + ierr = DriverInterpolate(cctkGH, Loop::dim, interpHandle, paramTableHandle, + coordSystemHandle, nPoints, interpCoordsTypeCode, + interpCoords, nInputArrays, inputArrayIndices, + nInputArrays, outputArrayTypes, outputArrays); + if (ierr < 0) { CCTK_WARN(CCTK_WARN_ALERT, "Interpolation error"); - goto label_free_param_table; } + Util_TableDestroy(paramTableHandle); + if (CCTK_MyProc(cctkGH) == 0) { // Some more output @@ -208,7 +184,6 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { } // Time evolution - for (int n = 0; n < max_num_tracked; ++n) { if (track[n]) { const CCTK_REAL dt = pt_loc_t[n] - pt_t_prev[n]; @@ -225,10 +200,8 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { } } - // Broadcast result - - CCTK_REAL loc_global[6 * max_num_tracked]; /* 3 components for location, 3 - components for velocity */ + // Broadcast result: 3 components for location, 3 components for velocity + CCTK_REAL loc_global[6 * max_num_tracked]; if (CCTK_MyProc(cctkGH) == 0) { for (int n = 0; n < max_num_tracked; ++n) { loc_global[n] = pt_loc_x[n]; @@ -261,11 +234,6 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { position_z[i] = pt_loc_z[i]; } } -// Done - -// Poor man's exception handling -label_free_param_table: - Util_TableDestroy(param_table_handle); } using namespace Arith; From 9a341185b23a5a2b2a79cfc9d7bb57219ad75542 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Thu, 27 Jun 2024 16:29:35 +0000 Subject: [PATCH 03/14] Multipole: remove redundant nOutputArrays in interpolator --- Multipole/src/surface.cxx | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/Multipole/src/surface.cxx b/Multipole/src/surface.cxx index 93ddadb0..d35a0220 100644 --- a/Multipole/src/surface.cxx +++ b/Multipole/src/surface.cxx @@ -20,7 +20,6 @@ void Surface::interpolate(CCTK_ARGUMENTS, int realFieldIndex, const void *interpCoords[Loop::dim] = {x_.data(), y_.data(), z_.data()}; CCTK_INT nInputArrays = imagFieldIndex == -1 ? 1 : 2; - CCTK_INT nOutputArrays = imagFieldIndex == -1 ? 1 : 2; const CCTK_INT inputArrayIndices[2] = {realFieldIndex, imagFieldIndex}; @@ -32,22 +31,27 @@ void Surface::interpolate(CCTK_ARGUMENTS, int realFieldIndex, const CCTK_INT interpCoordsTypeCode = 0; const CCTK_INT outputArrayTypes[1] = {0}; - int interpHandle = CCTK_InterpHandle("CarpetX"); + const int interpHandle = CCTK_InterpHandle("CarpetX"); if (interpHandle < 0) { CCTK_VERROR("Could not obtain interpolator handle for built-in 'CarpetX' " "interpolator: %d", interpHandle); } + int ierr; + // Interpolation parameter table int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); - int ierr = Util_TableSetFromString(paramTableHandle, interpolator_pars); + if ((ierr = Util_TableSetFromString(paramTableHandle, interpolator_pars)) < + 0) { + CCTK_VERROR("Can't set pars in parameter table: %d", ierr); + } ierr = DriverInterpolate(cctkGH, Loop::dim, interpHandle, paramTableHandle, coordSystemHandle, nPoints, interpCoordsTypeCode, interpCoords, nInputArrays, inputArrayIndices, - nOutputArrays, outputArrayTypes, outputArrays); + nInputArrays, outputArrayTypes, outputArrays); if (ierr < 0) { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, From d0fff9ad6b8da443467dfd01ff1938592f328c3e Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Thu, 27 Jun 2024 17:22:35 +0000 Subject: [PATCH 04/14] PunctureTracker: rename puncture_tracker.cxx --- PunctureTracker/src/make.code.defn | 4 ++-- .../src/{puncture_tracker.cc => puncture_tracker.cxx} | 0 2 files changed, 2 insertions(+), 2 deletions(-) rename PunctureTracker/src/{puncture_tracker.cc => puncture_tracker.cxx} (100%) diff --git a/PunctureTracker/src/make.code.defn b/PunctureTracker/src/make.code.defn index a980b052..20573869 100644 --- a/PunctureTracker/src/make.code.defn +++ b/PunctureTracker/src/make.code.defn @@ -1,7 +1,7 @@ -# Main make.code.defn file for thorn PunctureTracker -*-Makefile-*- +# Main make.code.defn file for thorn PunctureTracker # Source files in this directory -SRCS = puncture_tracker.cc +SRCS = puncture_tracker.cxx # Subdirectories containing source files SUBDIRS = diff --git a/PunctureTracker/src/puncture_tracker.cc b/PunctureTracker/src/puncture_tracker.cxx similarity index 100% rename from PunctureTracker/src/puncture_tracker.cc rename to PunctureTracker/src/puncture_tracker.cxx From 6a0ce41c8d892181dd4272015475b7c05a9b1972 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Thu, 27 Jun 2024 17:26:57 +0000 Subject: [PATCH 05/14] PunctureTracker: rename using namespace Loop --- PunctureTracker/src/puncture_tracker.cxx | 26 +++++++----------------- 1 file changed, 7 insertions(+), 19 deletions(-) diff --git a/PunctureTracker/src/puncture_tracker.cxx b/PunctureTracker/src/puncture_tracker.cxx index 35386114..c2eb1904 100644 --- a/PunctureTracker/src/puncture_tracker.cxx +++ b/PunctureTracker/src/puncture_tracker.cxx @@ -12,8 +12,6 @@ #include namespace PunctureTracker { -using namespace std; -using namespace Loop; const int max_num_tracked = 10; @@ -239,20 +237,9 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { using namespace Arith; extern "C" void PunctureTracker_CheckShift(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_PunctureTracker_CheckShift; + DECLARE_CCTK_ARGUMENTSX_PunctureTracker_CheckShift; DECLARE_CCTK_PARAMETERS; - const int dim = 3; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout(cctkGH, indextype); - - const GF3D2 betax_(layout, betax); - const GF3D2 betay_(layout, betay); - const GF3D2 betaz_(layout, betaz); - - const GridDescBaseDevice grid(cctkGH); - const int level = ilogb(CCTK_REAL(cctk_levfac[0])); const int finest_lvl = num_levels[0] - 1; @@ -260,17 +247,18 @@ extern "C" void PunctureTracker_CheckShift(CCTK_ARGUMENTS) { for (int n = 0; n < max_num_tracked; ++n) { if (track[n]) { - const vect loc_vec = {pt_loc_x[n], pt_loc_y[n], - pt_loc_z[n]}; + const vect loc_vec = {pt_loc_x[n], pt_loc_y[n], + pt_loc_z[n]}; grid.loop_all_device<0, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_DEVICE( + const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (maximum(abs(p.X - loc_vec)) <= (1 / pow(2, level))) { printf("Shift at level %d near puncture #%d is {%g, %g, %g} at " "coords {%g, %g, %g}.", - level, n, betax_(p.I), betay_(p.I), betaz_(p.I), p.x, - p.y, p.z); + level, n, betax(p.I), betay(p.I), betaz(p.I), p.x, p.y, + p.z); } }); } From b89ee4152aafdc36816cd5d296d3797b49f12380 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Thu, 27 Jun 2024 17:28:34 +0000 Subject: [PATCH 06/14] PunctureTracker: rename using namespace Arith --- PunctureTracker/src/puncture_tracker.cxx | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/PunctureTracker/src/puncture_tracker.cxx b/PunctureTracker/src/puncture_tracker.cxx index c2eb1904..aca41750 100644 --- a/PunctureTracker/src/puncture_tracker.cxx +++ b/PunctureTracker/src/puncture_tracker.cxx @@ -234,8 +234,6 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { } } -using namespace Arith; - extern "C" void PunctureTracker_CheckShift(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_PunctureTracker_CheckShift; DECLARE_CCTK_PARAMETERS; @@ -246,9 +244,8 @@ extern "C" void PunctureTracker_CheckShift(CCTK_ARGUMENTS) { if (level == finest_lvl) { for (int n = 0; n < max_num_tracked; ++n) { if (track[n]) { - - const vect loc_vec = {pt_loc_x[n], pt_loc_y[n], - pt_loc_z[n]}; + const Arith::vect loc_vec = { + pt_loc_x[n], pt_loc_y[n], pt_loc_z[n]}; grid.loop_all_device<0, 0, 0>( grid.nghostzones, From 02df55ee012d78de6e32542434d41f435a8e44ca Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 01:07:34 +0000 Subject: [PATCH 07/14] PunctureTracker: add class PunctureContainer --- PunctureTracker/schedule.ccl | 12 +- PunctureTracker/src/make.code.defn | 2 +- PunctureTracker/src/puncture.cxx | 113 +++++++++++ PunctureTracker/src/puncture.hxx | 57 ++++++ PunctureTracker/src/puncture_tracker.cxx | 238 +++++++++-------------- 5 files changed, 270 insertions(+), 152 deletions(-) create mode 100644 PunctureTracker/src/puncture.cxx create mode 100644 PunctureTracker/src/puncture.hxx diff --git a/PunctureTracker/schedule.ccl b/PunctureTracker/schedule.ccl index efd8783b..a9a075a4 100644 --- a/PunctureTracker/schedule.ccl +++ b/PunctureTracker/schedule.ccl @@ -1,24 +1,30 @@ # Schedule definitions for thorn PunctureTracker -SCHEDULE PunctureTracker_Init AT initial +SCHEDULE PunctureTracker_Setup AT initial { LANG: C + OPTIONS: GLOBAL WRITES: pt_loc WRITES: pt_vel WRITES: BoxInBox::positions - OPTIONS: GLOBAL } "Calculate initial location of punctures" SCHEDULE PunctureTracker_Track AT evol AFTER ODESolvers_Solve { LANG: C + OPTIONS: GLOBAL READS: ADMBaseX::shift(everywhere) WRITES: pt_loc WRITES: pt_vel WRITES: BoxInBox::positions - OPTIONS: GLOBAL } "Calculate new location of punctures" +SCHEDULE PunctureTracker_Finalize AT terminate BEFORE driver_terminate +{ + LANG: C + OPTIONS: GLOBAL +} "Free PunctureContainer instance and stuff" + if (read_shifts) { SCHEDULE PunctureTracker_CheckShift AT evol BEFORE PunctureTracker_Track diff --git a/PunctureTracker/src/make.code.defn b/PunctureTracker/src/make.code.defn index 20573869..ddb76be6 100644 --- a/PunctureTracker/src/make.code.defn +++ b/PunctureTracker/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn PunctureTracker # Source files in this directory -SRCS = puncture_tracker.cxx +SRCS = puncture.cxx puncture_tracker.cxx # Subdirectories containing source files SUBDIRS = diff --git a/PunctureTracker/src/puncture.cxx b/PunctureTracker/src/puncture.cxx new file mode 100644 index 00000000..ece9d94c --- /dev/null +++ b/PunctureTracker/src/puncture.cxx @@ -0,0 +1,113 @@ +#include "puncture.hxx" + +#include + +namespace PunctureTracker { + +void PunctureContainer::updatePreviousTime(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTS; + for (size_t n = 0; n < time_.size(); ++n) { + previousTime_[n] = time_[n]; + time_[n] = cctk_time; + } +} + +void PunctureContainer::interpolate(CCTK_ARGUMENTS) { + DECLARE_CCTK_PARAMETERS; + + // Only processor 0 interpolates + const CCTK_INT nPoints = CCTK_MyProc(cctkGH) == 0 ? numPunctures_ : 0; + + // Interpolation coordinates + const void *interpCoords[Loop::dim] = { + location_[0].data(), location_[1].data(), location_[2].data()}; + + // Interpolated variables + const CCTK_INT nInputArrays = 3; + const CCTK_INT inputArrayIndices[3] = {CCTK_VarIndex("ADMBaseX::betax"), + CCTK_VarIndex("ADMBaseX::betay"), + CCTK_VarIndex("ADMBaseX::betaz")}; + + CCTK_POINTER outputArrays[3] = {beta_[0].data(), beta_[1].data(), + beta_[2].data()}; + + /* DriverInterpolate arguments that aren't currently used */ + const int coordSystemHandle = 0; + const CCTK_INT interpCoordsTypeCode = 0; + const CCTK_INT outputArrayTypes[1] = {0}; + + const int interpHandle = CCTK_InterpHandle("CarpetX"); + if (interpHandle < 0) { + CCTK_WARN(CCTK_WARN_ALERT, "Can't get interpolation handle"); + return; + } + + int ierr; + + int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); + if (paramTableHandle < 0) { + CCTK_VERROR("Can't create parameter table: %d", paramTableHandle); + } + + if ((ierr = Util_TableSetInt(paramTableHandle, interp_order, "order")) < 0) { + CCTK_VERROR("Can't set order in parameter table: %d", ierr); + } + + // Interpolate + ierr = DriverInterpolate(cctkGH, Loop::dim, interpHandle, paramTableHandle, + coordSystemHandle, nPoints, interpCoordsTypeCode, + interpCoords, nInputArrays, inputArrayIndices, + nInputArrays, outputArrayTypes, outputArrays); + + if (ierr < 0) { + CCTK_WARN(CCTK_WARN_ALERT, "Interpolation error"); + } + + Util_TableDestroy(paramTableHandle); +} + +void PunctureContainer::evolve(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTS; + + if (CCTK_MyProc(cctkGH) == 0) { + // First order time integrator + // Michael Koppitz says this works... + // if it doesn't, we can make it second order accurate + for (size_t n = 0; n < time_.size(); ++n) { + const CCTK_REAL dt = time_[n] - previousTime_[n]; + for (int i = 0; i < Loop::dim; ++i) { + location_[i][n] += dt * (-beta_[i][n]); + velocity_[i][n] = -beta_[i][n]; + } + } + } +} + +void PunctureContainer::broadcast(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTS; + + // 3 components for location, 3 components for velocity + const CCTK_INT bufferSize = 6 * numPunctures_; + + CCTK_REAL buffer[bufferSize]; + + if (CCTK_MyProc(cctkGH) == 0) { + for (int i = 0; i < Loop::dim; ++i) { + for (int n = 0; n < numPunctures_; ++n) { + buffer[i * numPunctures_ + n] = location_[i][n]; + buffer[(i + Loop::dim) * numPunctures_ + n] = velocity_[i][n]; + } + } + } + + MPI_Bcast(buffer, bufferSize, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + for (int i = 0; i < Loop::dim; ++i) { + for (int n = 0; n < numPunctures_; ++n) { + location_[i][n] = buffer[i * numPunctures_ + n]; + velocity_[i][n] = buffer[(i + Loop::dim) * numPunctures_ + n]; + } + } +} + +} // namespace PunctureTracker diff --git a/PunctureTracker/src/puncture.hxx b/PunctureTracker/src/puncture.hxx new file mode 100644 index 00000000..586c86e5 --- /dev/null +++ b/PunctureTracker/src/puncture.hxx @@ -0,0 +1,57 @@ +/* puncture.hxx */ +/* (c) Liwei Ji 06/2024 */ + +#ifndef PUNCTURETRACKER_PUNCTURE_HXX +#define PUNCTURETRACKER_PUNCTURE_HXX + +#include + +#include + +#include + +namespace PunctureTracker { + +class PunctureContainer { +public: + PunctureContainer() {} + + virtual ~PunctureContainer() = default; // Use default for trivial destructors + + void setNumPunctures() { numPunctures_ = location_[0].size(); } + + CCTK_INT getNumPunctures() { return numPunctures_; } + + std::vector &getTime() { return time_; } + + std::vector &getPreviousTime() { return previousTime_; } + + std::array, Loop::dim> &getLocation() { + return location_; + } + + std::array, Loop::dim> &getVelocity() { + return velocity_; + } + + std::array, Loop::dim> &getBeta() { return beta_; } + + void updatePreviousTime(CCTK_ARGUMENTS); + + void interpolate(CCTK_ARGUMENTS); + + void evolve(CCTK_ARGUMENTS); + + void broadcast(CCTK_ARGUMENTS); + +private: + CCTK_INT numPunctures_; + std::vector time_, previousTime_; + std::array, Loop::dim> location_; + std::array, Loop::dim> velocity_; + std::array, Loop::dim> beta_; +}; + +} // namespace PunctureTracker + +#endif // #ifndef PUNCTURETRACKER_PUNCTURE_HXX diff --git a/PunctureTracker/src/puncture_tracker.cxx b/PunctureTracker/src/puncture_tracker.cxx index aca41750..74b6d90b 100644 --- a/PunctureTracker/src/puncture_tracker.cxx +++ b/PunctureTracker/src/puncture_tracker.cxx @@ -1,8 +1,10 @@ +#include "puncture.hxx" + #include #include #include #include -#include + #include #include @@ -13,10 +15,12 @@ namespace PunctureTracker { +static PunctureContainer *g_punctures = nullptr; + const int max_num_tracked = 10; -extern "C" void PunctureTracker_Init(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_PunctureTracker_Init; +extern "C" void PunctureTracker_Setup(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTS_PunctureTracker_Setup; DECLARE_CCTK_PARAMETERS; if (verbose) { @@ -34,7 +38,6 @@ extern "C" void PunctureTracker_Init(CCTK_ARGUMENTS) { pt_vel_y[n] = 0.0; pt_vel_z[n] = 0.0; } else { - // Initialise to some sensible but unimportant values pt_loc_t[n] = 0.0; pt_loc_x[n] = 0.0; pt_loc_y[n] = 0.0; @@ -46,191 +49,130 @@ extern "C" void PunctureTracker_Init(CCTK_ARGUMENTS) { } } + // Initialize PunctureContainer + if (g_punctures == nullptr) { + g_punctures = new PunctureContainer(); + } + + for (int n = 0; n < max_num_tracked; ++n) { + if (track[n]) { + g_punctures->getTime().push_back(cctk_time); + g_punctures->getLocation()[0].push_back(initial_x[n]); + g_punctures->getLocation()[1].push_back(initial_y[n]); + g_punctures->getLocation()[2].push_back(initial_z[n]); + g_punctures->getVelocity()[0].push_back(0.0); + g_punctures->getVelocity()[1].push_back(0.0); + g_punctures->getVelocity()[2].push_back(0.0); + } + } + const int nPunctures = g_punctures->getLocation()[0].size(); + g_punctures->getPreviousTime().resize(nPunctures); + g_punctures->getBeta()[0].resize(nPunctures); + g_punctures->getBeta()[1].resize(nPunctures); + g_punctures->getBeta()[2].resize(nPunctures); + g_punctures->setNumPunctures(); + assert(g_punctures->getNumPunctures() == nPunctures); + // enabled if refinement regions should follow the punctures if (track_boxes) { - const int max_num_regions = 2; - for (int i = 0; i < max_num_regions; i++) { - CCTK_VINFO("Writing punc coords to box %d.", i); - position_x[i] = pt_loc_x[i]; - position_y[i] = pt_loc_y[i]; - position_z[i] = pt_loc_z[i]; + const std::array, Loop::dim> &location = + g_punctures->getLocation(); + for (size_t i = 0; i < location[0].size(); ++i) { + CCTK_VINFO("Writing punc coords to box %zu.", i); + position_x[i] = location[0][i]; + position_y[i] = location[1][i]; + position_z[i] = location[2][i]; } } } +extern "C" void PunctureTracker_Finalize(CCTK_ARGUMENTS) { + delete g_punctures; + g_punctures = nullptr; +} + extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS_PunctureTracker_Track; DECLARE_CCTK_PARAMETERS; - // Do not track while setting up initial data; - // time interpolation may fail - + // Do not track while setting up initial data; time interpolation may fail if (cctk_iteration == 0) { return; } // Some output - if (verbose) { CCTK_INFO("Tracking punctures..."); } + const std::array, Loop::dim> &location = + g_punctures->getLocation(); + if (verbose) { - for (int n = 0; n < max_num_tracked; ++n) { + for (size_t n = 0; n < location[0].size(); ++n) { if (track[n]) { - CCTK_VINFO("Puncture #%d is at (%g,%g,%g)", n, double(pt_loc_x[n]), - double(pt_loc_y[n]), double(pt_loc_z[n])); + CCTK_VINFO("Puncture #%zu is at (%g,%g,%g)", n, double(location[0][n]), + double(location[1][n]), double(location[2][n])); } } } // Manual time level cycling - CCTK_REAL pt_t_prev[max_num_tracked]; - - for (int n = 0; n < max_num_tracked; ++n) { - if (track[n]) { - pt_t_prev[n] = pt_loc_t[n]; - pt_loc_t[n] = cctk_time; - pt_vel_t[n] = cctk_time; - } else { - pt_t_prev[n] = 0.0; - } - } + g_punctures->updatePreviousTime(CCTK_PASS_CTOC); // Interpolate - { - // Only processor 0 interpolates - const CCTK_INT nPoints = CCTK_MyProc(cctkGH) == 0 ? max_num_tracked : 0; - - // Interpolation coordinates - const void *interpCoords[Loop::dim] = {pt_loc_x, pt_loc_y, pt_loc_z}; - - // Interpolated variables - const CCTK_INT nInputArrays = 3; - const CCTK_INT inputArrayIndices[3] = {CCTK_VarIndex("ADMBaseX::betax"), - CCTK_VarIndex("ADMBaseX::betay"), - CCTK_VarIndex("ADMBaseX::betaz")}; - - // Interpolation result - CCTK_REAL pt_betax[max_num_tracked]; - CCTK_REAL pt_betay[max_num_tracked]; - CCTK_REAL pt_betaz[max_num_tracked]; - CCTK_POINTER outputArrays[3] = {pt_betax, pt_betay, pt_betaz}; - - /* DriverInterpolate arguments that aren't currently used */ - const int coordSystemHandle = 0; - const CCTK_INT interpCoordsTypeCode = 0; - const CCTK_INT outputArrayTypes[1] = {0}; - - const int interpHandle = CCTK_InterpHandle("CarpetX"); - if (interpHandle < 0) { - CCTK_WARN(CCTK_WARN_ALERT, "Can't get interpolation handle"); - return; - } - - int ierr; - - int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); - if (paramTableHandle < 0) { - CCTK_VERROR("Can't create parameter table: %d", paramTableHandle); - } - - if ((ierr = Util_TableSetInt(paramTableHandle, interp_order, "order")) < - 0) { - CCTK_VERROR("Can't set order in parameter table: %d", ierr); - } - - // Interpolate - ierr = DriverInterpolate(cctkGH, Loop::dim, interpHandle, paramTableHandle, - coordSystemHandle, nPoints, interpCoordsTypeCode, - interpCoords, nInputArrays, inputArrayIndices, - nInputArrays, outputArrayTypes, outputArrays); - - if (ierr < 0) { - CCTK_WARN(CCTK_WARN_ALERT, "Interpolation error"); - } + g_punctures->interpolate(CCTK_PASS_CTOC); - Util_TableDestroy(paramTableHandle); + if (CCTK_MyProc(cctkGH) == 0) { + const std::array, Loop::dim> &beta = + g_punctures->getBeta(); - if (CCTK_MyProc(cctkGH) == 0) { - - // Some more output - - if (verbose) { - for (int n = 0; n < max_num_tracked; ++n) { - if (track[n]) { - CCTK_VINFO("Shift at puncture #%d is at (%g,%g,%g)", n, - double(pt_betax[n]), double(pt_betay[n]), - double(pt_betaz[n])); - } - } + // More output + if (verbose) { + for (size_t n = 0; n < beta[0].size(); ++n) { + CCTK_VINFO("Shift at puncture #%zu is at (%g,%g,%g)", n, + double(beta[0][n]), double(beta[1][n]), double(beta[2][n])); } + } - // Check for NaNs and large shift components - for (int n = 0; n < max_num_tracked; ++n) { - if (track[n]) { - CCTK_REAL norm = sqrt(pow(pt_betax[n], 2) + pow(pt_betay[n], 2) + - pow(pt_betaz[n], 2)); - - if (!CCTK_isfinite(norm) || norm > shift_limit) { - CCTK_VERROR("Shift at puncture #%d is (%g,%g,%g). This likely " - "indicates an error in the simulation.", - n, double(pt_betax[n]), double(pt_betay[n]), - double(pt_betaz[n])); - } - } - } + // Check for NaNs and large shift components + for (size_t n = 0; n < beta[0].size(); ++n) { + CCTK_REAL norm = + sqrt(pow(beta[0][n], 2) + pow(beta[1][n], 2) + pow(beta[2][n], 2)); - // Time evolution - for (int n = 0; n < max_num_tracked; ++n) { - if (track[n]) { - const CCTK_REAL dt = pt_loc_t[n] - pt_t_prev[n]; - // First order time integrator - // Michael Koppitz says this works... - // if it doesn't, we can make it second order accurate - pt_loc_x[n] += dt * (-pt_betax[n]); - pt_loc_y[n] += dt * (-pt_betay[n]); - pt_loc_z[n] += dt * (-pt_betaz[n]); - pt_vel_x[n] = -pt_betax[n]; - pt_vel_y[n] = -pt_betay[n]; - pt_vel_z[n] = -pt_betaz[n]; - } + if (!CCTK_isfinite(norm) || norm > shift_limit) { + CCTK_VERROR("Shift at puncture #%zu is (%g,%g,%g). This likely " + "indicates an error in the simulation.", + n, double(beta[0][n]), double(beta[1][n]), + double(beta[2][n])); } } + } - // Broadcast result: 3 components for location, 3 components for velocity - CCTK_REAL loc_global[6 * max_num_tracked]; - if (CCTK_MyProc(cctkGH) == 0) { - for (int n = 0; n < max_num_tracked; ++n) { - loc_global[n] = pt_loc_x[n]; - loc_global[max_num_tracked + n] = pt_loc_y[n]; - loc_global[2 * max_num_tracked + n] = pt_loc_z[n]; - loc_global[3 * max_num_tracked + n] = pt_vel_x[n]; - loc_global[4 * max_num_tracked + n] = pt_vel_y[n]; - loc_global[5 * max_num_tracked + n] = pt_vel_z[n]; - } - } + // Time evolution + g_punctures->evolve(CCTK_PASS_CTOC); - // CarpetX doesn't register reduction handle, here's a quick fix. - MPI_Bcast(loc_global, 6 * max_num_tracked, MPI_DOUBLE, 0, MPI_COMM_WORLD); + // Broadcast result: 3 components for location, 3 components for velocity + g_punctures->broadcast(CCTK_PASS_CTOC); - for (int n = 0; n < max_num_tracked; ++n) { - pt_loc_x[n] = loc_global[n]; - pt_loc_y[n] = loc_global[max_num_tracked + n]; - pt_loc_z[n] = loc_global[2 * max_num_tracked + n]; - pt_vel_x[n] = loc_global[3 * max_num_tracked + n]; - pt_vel_y[n] = loc_global[4 * max_num_tracked + n]; - pt_vel_z[n] = loc_global[5 * max_num_tracked + n]; + if (track_boxes) { + for (size_t i = 0; i < location[0].size(); ++i) { + position_x[i] = location[0][i]; + position_y[i] = location[1][i]; + position_z[i] = location[2][i]; } } - if (track_boxes) { - const int max_num_regions = 2; - for (int i = 0; i < max_num_regions; i++) { - position_x[i] = pt_loc_x[i]; - position_y[i] = pt_loc_y[i]; - position_z[i] = pt_loc_z[i]; - } + // Write to pt_loc_foo and pt_vel_foo + const std::array, Loop::dim> &velocity = + g_punctures->getVelocity(); + for (size_t i = 0; i < location[0].size(); ++i) { + pt_loc_x[i] = location[0][i]; + pt_loc_y[i] = location[1][i]; + pt_loc_z[i] = location[2][i]; + pt_vel_x[i] = velocity[0][i]; + pt_vel_y[i] = velocity[1][i]; + pt_vel_z[i] = velocity[2][i]; } } From e0db5a80729297ca5e3d6321c3d9c3dc6cf2813c Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 01:30:31 +0000 Subject: [PATCH 08/14] PunctureTracker: beautify member function broadcast --- PunctureTracker/src/puncture.cxx | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/PunctureTracker/src/puncture.cxx b/PunctureTracker/src/puncture.cxx index ece9d94c..cc054948 100644 --- a/PunctureTracker/src/puncture.cxx +++ b/PunctureTracker/src/puncture.cxx @@ -67,8 +67,6 @@ void PunctureContainer::interpolate(CCTK_ARGUMENTS) { } void PunctureContainer::evolve(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS; - if (CCTK_MyProc(cctkGH) == 0) { // First order time integrator // Michael Koppitz says this works... @@ -84,12 +82,11 @@ void PunctureContainer::evolve(CCTK_ARGUMENTS) { } void PunctureContainer::broadcast(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS; + const CCTK_INT numComponents = 6; + const CCTK_INT bufferSize = numComponents * numPunctures_; // 3 components for location, 3 components for velocity - const CCTK_INT bufferSize = 6 * numPunctures_; - - CCTK_REAL buffer[bufferSize]; + std::vector buffer(bufferSize); if (CCTK_MyProc(cctkGH) == 0) { for (int i = 0; i < Loop::dim; ++i) { @@ -100,7 +97,12 @@ void PunctureContainer::broadcast(CCTK_ARGUMENTS) { } } - MPI_Bcast(buffer, bufferSize, MPI_DOUBLE, 0, MPI_COMM_WORLD); + int mpiError = + MPI_Bcast(buffer.data(), bufferSize, MPI_DOUBLE, 0, MPI_COMM_WORLD); + if (mpiError != MPI_SUCCESS) { + CCTK_VINFO("MPI_Bcast failed with error code %d", mpiError); + MPI_Abort(MPI_COMM_WORLD, mpiError); + } for (int i = 0; i < Loop::dim; ++i) { for (int n = 0; n < numPunctures_; ++n) { From cd0534106f4cc2cb7f2a530d30699905cda27cd1 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 01:37:56 +0000 Subject: [PATCH 09/14] PunctureTracker: beautify member function evolve --- PunctureTracker/src/puncture.cxx | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/PunctureTracker/src/puncture.cxx b/PunctureTracker/src/puncture.cxx index cc054948..f5eca306 100644 --- a/PunctureTracker/src/puncture.cxx +++ b/PunctureTracker/src/puncture.cxx @@ -71,7 +71,7 @@ void PunctureContainer::evolve(CCTK_ARGUMENTS) { // First order time integrator // Michael Koppitz says this works... // if it doesn't, we can make it second order accurate - for (size_t n = 0; n < time_.size(); ++n) { + for (int n = 0; n < numPunctures_; ++n) { const CCTK_REAL dt = time_[n] - previousTime_[n]; for (int i = 0; i < Loop::dim; ++i) { location_[i][n] += dt * (-beta_[i][n]); @@ -83,10 +83,8 @@ void PunctureContainer::evolve(CCTK_ARGUMENTS) { void PunctureContainer::broadcast(CCTK_ARGUMENTS) { const CCTK_INT numComponents = 6; - const CCTK_INT bufferSize = numComponents * numPunctures_; - // 3 components for location, 3 components for velocity - std::vector buffer(bufferSize); + std::vector buffer(numComponents * numPunctures_); if (CCTK_MyProc(cctkGH) == 0) { for (int i = 0; i < Loop::dim; ++i) { @@ -98,7 +96,7 @@ void PunctureContainer::broadcast(CCTK_ARGUMENTS) { } int mpiError = - MPI_Bcast(buffer.data(), bufferSize, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(buffer.data(), buffer.size(), MPI_DOUBLE, 0, MPI_COMM_WORLD); if (mpiError != MPI_SUCCESS) { CCTK_VINFO("MPI_Bcast failed with error code %d", mpiError); MPI_Abort(MPI_COMM_WORLD, mpiError); From ea90c3de44f8e8e0594e31243d0443cf6596b970 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 01:46:41 +0000 Subject: [PATCH 10/14] PunctureTracker: beautify member function interpolate --- PunctureTracker/src/puncture.cxx | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/PunctureTracker/src/puncture.cxx b/PunctureTracker/src/puncture.cxx index f5eca306..6df4cb22 100644 --- a/PunctureTracker/src/puncture.cxx +++ b/PunctureTracker/src/puncture.cxx @@ -6,7 +6,7 @@ namespace PunctureTracker { void PunctureContainer::updatePreviousTime(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; - for (size_t n = 0; n < time_.size(); ++n) { + for (int n = 0; n < numPunctures_; ++n) { previousTime_[n] = time_[n]; time_[n] = cctk_time; } @@ -16,7 +16,7 @@ void PunctureContainer::interpolate(CCTK_ARGUMENTS) { DECLARE_CCTK_PARAMETERS; // Only processor 0 interpolates - const CCTK_INT nPoints = CCTK_MyProc(cctkGH) == 0 ? numPunctures_ : 0; + const CCTK_INT nPoints = (CCTK_MyProc(cctkGH) == 0) ? numPunctures_ : 0; // Interpolation coordinates const void *interpCoords[Loop::dim] = { @@ -24,14 +24,14 @@ void PunctureContainer::interpolate(CCTK_ARGUMENTS) { // Interpolated variables const CCTK_INT nInputArrays = 3; - const CCTK_INT inputArrayIndices[3] = {CCTK_VarIndex("ADMBaseX::betax"), - CCTK_VarIndex("ADMBaseX::betay"), - CCTK_VarIndex("ADMBaseX::betaz")}; + const CCTK_INT inputArrayIndices[nInputArrays] = { + CCTK_VarIndex("ADMBaseX::betax"), CCTK_VarIndex("ADMBaseX::betay"), + CCTK_VarIndex("ADMBaseX::betaz")}; - CCTK_POINTER outputArrays[3] = {beta_[0].data(), beta_[1].data(), - beta_[2].data()}; + CCTK_POINTER outputArrays[nInputArrays] = {beta_[0].data(), beta_[1].data(), + beta_[2].data()}; - /* DriverInterpolate arguments that aren't currently used */ + // DriverInterpolate arguments that aren't currently used const int coordSystemHandle = 0; const CCTK_INT interpCoordsTypeCode = 0; const CCTK_INT outputArrayTypes[1] = {0}; @@ -42,18 +42,19 @@ void PunctureContainer::interpolate(CCTK_ARGUMENTS) { return; } - int ierr; - - int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); + // Create parameter table for interpolation + const int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); if (paramTableHandle < 0) { CCTK_VERROR("Can't create parameter table: %d", paramTableHandle); } - if ((ierr = Util_TableSetInt(paramTableHandle, interp_order, "order")) < 0) { + // Set interpolation order in the parameter table + int ierr = Util_TableSetInt(paramTableHandle, interp_order, "order"); + if (ierr < 0) { CCTK_VERROR("Can't set order in parameter table: %d", ierr); } - // Interpolate + // Perform the interpolation ierr = DriverInterpolate(cctkGH, Loop::dim, interpHandle, paramTableHandle, coordSystemHandle, nPoints, interpCoordsTypeCode, interpCoords, nInputArrays, inputArrayIndices, @@ -63,6 +64,7 @@ void PunctureContainer::interpolate(CCTK_ARGUMENTS) { CCTK_WARN(CCTK_WARN_ALERT, "Interpolation error"); } + // Destroy the parameter table Util_TableDestroy(paramTableHandle); } From dddf97c4032b3b937e6dee641e274d5e276259c7 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 02:46:34 +0000 Subject: [PATCH 11/14] PunctureTracker: remove unused function --- PunctureTracker/param.ccl | 4 ---- PunctureTracker/schedule.ccl | 13 ----------- PunctureTracker/src/puncture_tracker.cxx | 29 ------------------------ 3 files changed, 46 deletions(-) diff --git a/PunctureTracker/param.ccl b/PunctureTracker/param.ccl index 62754105..e82efe50 100644 --- a/PunctureTracker/param.ccl +++ b/PunctureTracker/param.ccl @@ -4,10 +4,6 @@ BOOLEAN verbose "speak up?" STEERABLE=always { } "no" -BOOLEAN read_shifts "Enable to write nearby shifts to stdout" STEERABLE=always -{ -} "no" - BOOLEAN track[10] "Track this puncture" { } "no" diff --git a/PunctureTracker/schedule.ccl b/PunctureTracker/schedule.ccl index a9a075a4..3433339a 100644 --- a/PunctureTracker/schedule.ccl +++ b/PunctureTracker/schedule.ccl @@ -24,16 +24,3 @@ SCHEDULE PunctureTracker_Finalize AT terminate BEFORE driver_terminate LANG: C OPTIONS: GLOBAL } "Free PunctureContainer instance and stuff" - -if (read_shifts) -{ - SCHEDULE PunctureTracker_CheckShift AT evol BEFORE PunctureTracker_Track - { - LANG: C - READS: ADMBaseX::shift(interior) - READS: BoxInBox::num_levels - READS: pt_loc - READS: CoordinatesX::vertex_coords(interior) - OPTIONS: LOCAL - } "Outputs nearby shift quantities to stdout" -} diff --git a/PunctureTracker/src/puncture_tracker.cxx b/PunctureTracker/src/puncture_tracker.cxx index 74b6d90b..36f2c748 100644 --- a/PunctureTracker/src/puncture_tracker.cxx +++ b/PunctureTracker/src/puncture_tracker.cxx @@ -176,33 +176,4 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { } } -extern "C" void PunctureTracker_CheckShift(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_PunctureTracker_CheckShift; - DECLARE_CCTK_PARAMETERS; - - const int level = ilogb(CCTK_REAL(cctk_levfac[0])); - const int finest_lvl = num_levels[0] - 1; - - if (level == finest_lvl) { - for (int n = 0; n < max_num_tracked; ++n) { - if (track[n]) { - const Arith::vect loc_vec = { - pt_loc_x[n], pt_loc_y[n], pt_loc_z[n]}; - - grid.loop_all_device<0, 0, 0>( - grid.nghostzones, - [=] CCTK_DEVICE( - const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - if (maximum(abs(p.X - loc_vec)) <= (1 / pow(2, level))) { - printf("Shift at level %d near puncture #%d is {%g, %g, %g} at " - "coords {%g, %g, %g}.", - level, n, betax(p.I), betay(p.I), betaz(p.I), p.x, p.y, - p.z); - } - }); - } - } - } -} - } // namespace PunctureTracker From 963c05bc906640595c98e5a3e8fea1279cdda25e Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 02:48:44 +0000 Subject: [PATCH 12/14] PunctureTracker: remove unused parameters --- PunctureTracker/param.ccl | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/PunctureTracker/param.ccl b/PunctureTracker/param.ccl index e82efe50..97c63939 100644 --- a/PunctureTracker/param.ccl +++ b/PunctureTracker/param.ccl @@ -27,32 +27,11 @@ REAL initial_z[10] "Initial z coordinate positions of punctures" *:* :: "" } 0.0 -INT modify_puncture[2] "Punctures to use for modification criteria" -{ - -1:9 :: "One of the tracking punctures or negative for no modification" -} -1 - -REAL modify_distance "Modify levels when the distance is less than this" -{ - 0.0:* :: "zero or positive" -} 0.0 - -INT new_reflevel_number[2] "The new number of refinement levels" -{ - -1:* :: "Negative for no change" -} -1 - INT interp_order "Shift interpolation order" { 0:* :: "" } 1 -INT which_surface_to_store_info[10] "A spherical surface index where we can store the puncture location" -{ - -1 :: "don't store puncture location" - 0:* :: "any spherical surface index" -} -1 - REAL shift_limit "Shift components must be less than this in magnitude" { 0.0:* :: "zero or positive" From 0f20e78b35270dd8c06890d91ee72853cd16308f Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 03:23:20 +0000 Subject: [PATCH 13/14] PunctureTracker: beautify PunctureTracker_Track --- PunctureTracker/src/puncture_tracker.cxx | 45 ++++++++++++------------ 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/PunctureTracker/src/puncture_tracker.cxx b/PunctureTracker/src/puncture_tracker.cxx index 36f2c748..accc9ff2 100644 --- a/PunctureTracker/src/puncture_tracker.cxx +++ b/PunctureTracker/src/puncture_tracker.cxx @@ -65,7 +65,7 @@ extern "C" void PunctureTracker_Setup(CCTK_ARGUMENTS) { g_punctures->getVelocity()[2].push_back(0.0); } } - const int nPunctures = g_punctures->getLocation()[0].size(); + const int nPunctures = g_punctures->getTime().size(); g_punctures->getPreviousTime().resize(nPunctures); g_punctures->getBeta()[0].resize(nPunctures); g_punctures->getBeta()[1].resize(nPunctures); @@ -77,8 +77,8 @@ extern "C" void PunctureTracker_Setup(CCTK_ARGUMENTS) { if (track_boxes) { const std::array, Loop::dim> &location = g_punctures->getLocation(); - for (size_t i = 0; i < location[0].size(); ++i) { - CCTK_VINFO("Writing punc coords to box %zu.", i); + for (int i = 0; i < nPunctures; ++i) { + CCTK_VINFO("Writing punc coords to box %d.", i); position_x[i] = location[0][i]; position_y[i] = location[1][i]; position_z[i] = location[2][i]; @@ -105,15 +105,16 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { CCTK_INFO("Tracking punctures..."); } + const int nPunctures = g_punctures->getNumPunctures(); const std::array, Loop::dim> &location = g_punctures->getLocation(); + const std::array, Loop::dim> &velocity = + g_punctures->getVelocity(); if (verbose) { - for (size_t n = 0; n < location[0].size(); ++n) { - if (track[n]) { - CCTK_VINFO("Puncture #%zu is at (%g,%g,%g)", n, double(location[0][n]), - double(location[1][n]), double(location[2][n])); - } + for (int n = 0; n < nPunctures; ++n) { + CCTK_VINFO("Puncture #%d is at (%g,%g,%g)", n, double(location[0][n]), + double(location[1][n]), double(location[2][n])); } } @@ -129,19 +130,19 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { // More output if (verbose) { - for (size_t n = 0; n < beta[0].size(); ++n) { - CCTK_VINFO("Shift at puncture #%zu is at (%g,%g,%g)", n, + for (int n = 0; n < nPunctures; ++n) { + CCTK_VINFO("Shift at puncture #%d is at (%g,%g,%g)", n, double(beta[0][n]), double(beta[1][n]), double(beta[2][n])); } } // Check for NaNs and large shift components - for (size_t n = 0; n < beta[0].size(); ++n) { + for (int n = 0; n < nPunctures; ++n) { CCTK_REAL norm = sqrt(pow(beta[0][n], 2) + pow(beta[1][n], 2) + pow(beta[2][n], 2)); if (!CCTK_isfinite(norm) || norm > shift_limit) { - CCTK_VERROR("Shift at puncture #%zu is (%g,%g,%g). This likely " + CCTK_VERROR("Shift at puncture #%d is (%g,%g,%g). This likely " "indicates an error in the simulation.", n, double(beta[0][n]), double(beta[1][n]), double(beta[2][n])); @@ -155,18 +156,8 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { // Broadcast result: 3 components for location, 3 components for velocity g_punctures->broadcast(CCTK_PASS_CTOC); - if (track_boxes) { - for (size_t i = 0; i < location[0].size(); ++i) { - position_x[i] = location[0][i]; - position_y[i] = location[1][i]; - position_z[i] = location[2][i]; - } - } - // Write to pt_loc_foo and pt_vel_foo - const std::array, Loop::dim> &velocity = - g_punctures->getVelocity(); - for (size_t i = 0; i < location[0].size(); ++i) { + for (int i = 0; i < nPunctures; ++i) { pt_loc_x[i] = location[0][i]; pt_loc_y[i] = location[1][i]; pt_loc_z[i] = location[2][i]; @@ -174,6 +165,14 @@ extern "C" void PunctureTracker_Track(CCTK_ARGUMENTS) { pt_vel_y[i] = velocity[1][i]; pt_vel_z[i] = velocity[2][i]; } + + 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]; + } + } } } // namespace PunctureTracker From 1d6ee7b268ebe151a0060c13f4ed7cc58457ae54 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Fri, 28 Jun 2024 03:29:06 +0000 Subject: [PATCH 14/14] PunctureTracker: add copy right to puncture.cxx --- PunctureTracker/src/puncture.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/PunctureTracker/src/puncture.cxx b/PunctureTracker/src/puncture.cxx index 6df4cb22..403b4b2b 100644 --- a/PunctureTracker/src/puncture.cxx +++ b/PunctureTracker/src/puncture.cxx @@ -1,3 +1,6 @@ +/* puncture.cxx */ +/* (c) Liwei Ji 06/2024 */ + #include "puncture.hxx" #include