Merge pull request #28 from EinsteinToolkit/lwji/puncture-tracker

Beautify PunctureTracker by introducing class PunctureContainer

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,

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"
@@ -31,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"

PunctureTracker/schedule.ccl

@@ -1,33 +1,26 @@
 # 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"
 
-if (read_shifts)
+SCHEDULE PunctureTracker_Finalize AT terminate BEFORE driver_terminate
 {
-	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"
-}
+  LANG: C
+  OPTIONS: GLOBAL
+} "Free PunctureContainer instance and stuff"

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 paramcheck.cc
+SRCS = puncture.cxx puncture_tracker.cxx
 
 # Subdirectories containing source files
 SUBDIRS = 

PunctureTracker/src/puncture.cxx

@@ -0,0 +1,118 @@
+/* puncture.cxx */
+/* (c) Liwei Ji 06/2024 */
+
+#include "puncture.hxx"
+
+#include <util_Table.h>
+
+namespace PunctureTracker {
+
+void PunctureContainer::updatePreviousTime(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTS;
+  for (int n = 0; n < numPunctures_; ++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[nInputArrays] = {
+      CCTK_VarIndex("ADMBaseX::betax"), CCTK_VarIndex("ADMBaseX::betay"),
+      CCTK_VarIndex("ADMBaseX::betaz")};
+
+  CCTK_POINTER outputArrays[nInputArrays] = {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;
+  }
+
+  // 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);
+  }
+
+  // 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);
+  }
+
+  // Perform the interpolation
+  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");
+  }
+
+  // Destroy the parameter table
+  Util_TableDestroy(paramTableHandle);
+}
+
+void PunctureContainer::evolve(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 (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]);
+        velocity_[i][n] = -beta_[i][n];
+      }
+    }
+  }
+}
+
+void PunctureContainer::broadcast(CCTK_ARGUMENTS) {
+  const CCTK_INT numComponents = 6;
+  // 3 components for location, 3 components for velocity
+  std::vector<CCTK_REAL> buffer(numComponents * numPunctures_);
+
+  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];
+      }
+    }
+  }
+
+  int mpiError =
+      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);
+  }
+
+  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

PunctureTracker/src/puncture.hxx

@@ -0,0 +1,57 @@
+/* puncture.hxx */
+/* (c) Liwei Ji 06/2024 */
+
+#ifndef PUNCTURETRACKER_PUNCTURE_HXX
+#define PUNCTURETRACKER_PUNCTURE_HXX
+
+#include <loop_device.hxx>
+
+#include <cctk.h>
+
+#include <vector>
+
+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<CCTK_REAL> &getTime() { return time_; }
+
+  std::vector<CCTK_REAL> &getPreviousTime() { return previousTime_; }
+
+  std::array<std::vector<CCTK_REAL>, Loop::dim> &getLocation() {
+    return location_;
+  }
+
+  std::array<std::vector<CCTK_REAL>, Loop::dim> &getVelocity() {
+    return velocity_;
+  }
+
+  std::array<std::vector<CCTK_REAL>, 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<CCTK_REAL> time_, previousTime_;
+  std::array<std::vector<CCTK_REAL>, Loop::dim> location_;
+  std::array<std::vector<CCTK_REAL>, Loop::dim> velocity_;
+  std::array<std::vector<CCTK_REAL>, Loop::dim> beta_;
+};
+
+} // namespace PunctureTracker
+
+#endif // #ifndef PUNCTURETRACKER_PUNCTURE_HXX