Rewrite puncture tracker

git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/PunctureTracker/trunk@6 a2659f00-0f4f-0410-9214-a4596bbb8a4f

9 files changed, 962 insertions, 347 deletions

diff --git a/README b/README
index 14af522..dfaeac4 100644
--- a/README
+++ b/README
@@ -1,11 +1,10 @@
-CVS info   : $Header$
-
 Cactus Code Thorn PunctureTracker
 Thorn Author(s)     : Michael Koppitz <koppitz@aei.mpg.de>
-Thorn Maintainer(s) : Michael Koppitz <koppitz@aei.mpg.de>
+                      Erik Schnetter <schnetter@cct.lsu.edu>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@cct.lsu.edu>
 --------------------------------------------------------------------------
 
 Purpose of the thorn:
 
-time integrate shift at location of puncture to follow location of puncture
-
+Time integrate the shift at the location of the punctures to follow the
+location of the punctures.
diff --git a/interface.ccl b/interface.ccl
index 8c79d3f..beac81a 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,20 +1,17 @@
 # Interface definition for thorn PunctureTracker
-# $Header$
 
 IMPLEMENTS: PunctureTracker
 
-INHERITS: ADMBase, TwoPunctures,SphericalSurface
+INHERITS: ADMBase CarpetRegrid2
 
 
-CCTK_REAL pt_loc[pt_num_tracked] TYPE=scalar  timelevels=2
-{
- pt_loc_x, pt_loc_y, pt_loc_z
-
-} "location of punctures"
 
-CCTK_REAL pt_shift[pt_num_tracked] TYPE=scalar
+CCTK_REAL pt_loc[10] TYPE=scalar
 {
- pt_shiftx, pt_shifty, pt_shiftz
-
-} "shift at location of punctures"
+  pt_loc_t pt_loc_x pt_loc_y pt_loc_z
+} "Location of punctures"
 
+CCTK_REAL pt_loc_p[10] TYPE=scalar
+{
+  pt_loc_t_p pt_loc_x_p pt_loc_y_p pt_loc_z_p
+} "Previous location of punctures"
diff --git a/par/qc0-punctrac.par b/par/qc0-punctrac.par
new file mode 100644
index 0000000..19ba635
--- /dev/null
+++ b/par/qc0-punctrac.par
@@ -0,0 +1,638 @@
+Cactus::cctk_run_title = "QC-0"
+
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+Cactus::cctk_timer_output          = "full"
+
+Cactus::terminate       = "time"
+#Cactus::cctk_final_time = 1.0
+#Cactus::cctk_final_time = 10.0
+Cactus::cctk_final_time = 100.0
+
+
+
+ActiveThorns = "IOUtil"
+
+IO::out_dir = $parfile
+
+
+
+ActiveThorns = "ManualTermination"
+
+ManualTermination::max_walltime          = @WALLTIME_HOURS@   # hours
+ManualTermination::on_remaining_walltime = 60                 # minutes
+
+ManualTermination::termination_from_file   = yes
+ManualTermination::create_termination_file = yes
+ManualTermination::termination_file        = "../TERMINATE"
+
+
+
+ActiveThorns = "AEILocalInterp"
+
+ActiveThorns = "Fortran"
+
+ActiveThorns = "LocalInterp"
+
+ActiveThorns = "MPIClock"
+
+#ActiveThorns = "NaNCatcher"
+
+ActiveThorns = "Slab"
+
+ActiveThorns = "TGRtensor"
+
+
+
+ActiveThorns = "SphericalSurface"
+
+SphericalSurface::nsurfaces = 6
+SphericalSurface::maxntheta = 39
+SphericalSurface::maxnphi   = 76
+
+SphericalSurface::ntheta      [0] = 39
+SphericalSurface::nphi        [0] = 76
+SphericalSurface::nghoststheta[0] = 2
+SphericalSurface::nghostsphi  [0] = 2
+
+SphericalSurface::ntheta      [1] = 39
+SphericalSurface::nphi        [1] = 76
+SphericalSurface::nghoststheta[1] = 2
+SphericalSurface::nghostsphi  [1] = 2
+
+SphericalSurface::ntheta      [2] = 39
+SphericalSurface::nphi        [2] = 76
+SphericalSurface::nghoststheta[2] = 2
+SphericalSurface::nghostsphi  [2] = 2
+
+SphericalSurface::ntheta      [3] = 39
+SphericalSurface::nphi        [3] = 76
+SphericalSurface::nghoststheta[3] = 2
+SphericalSurface::nghostsphi  [3] = 2
+
+SphericalSurface::ntheta      [4] = 39
+SphericalSurface::nphi        [4] = 76
+SphericalSurface::nghoststheta[4] = 2
+SphericalSurface::nghostsphi  [4] = 2
+
+SphericalSurface::ntheta      [5] = 39
+SphericalSurface::nphi        [5] = 76
+SphericalSurface::nghoststheta[5] = 2
+SphericalSurface::nghostsphi  [5] = 2
+
+
+
+ActiveThorns = "SummationByParts"
+
+SummationByParts::order = 4
+
+
+
+ActiveThorns = "InitBase"
+
+
+
+ActiveThorns = "Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab"
+
+Carpet::verbose           = no
+Carpet::schedule_barriers = no
+Carpet::veryverbose       = no
+CarpetLib::output_bboxes  = no
+
+Carpet::domain_from_coordbase = yes
+Carpet::max_refinement_levels = 10
+
+driver::ghost_size       = 3
+Carpet::use_buffer_zones = yes
+
+Carpet::prolongation_order_space = 5
+Carpet::prolongation_order_time  = 2
+
+Carpet::convergence_level = 0
+
+Carpet::init_3_timelevels = yes
+
+Carpet::poison_new_timelevels = yes
+CarpetLib::poison_new_memory  = yes
+
+Carpet::output_timers_every      = 512
+CarpetLib::print_timestats_every = 512
+CarpetLib::timestat_timer        = "rdtsc"
+CarpetLib::print_memstats_every  = 512
+
+
+
+ActiveThorns = "NaNChecker"
+
+NaNChecker::check_every     = 512
+NaNChecker::action_if_found = "just warn"
+NaNChecker::check_vars      = "
+        ADM_BSSN::ADM_BSSN_phi
+        ADM_BSSN::ADM_BSSN_metric
+        ADM_BSSN::ADM_BSSN_curv
+        ADM_BSSN::ADM_BSSN_K
+        ADM_BSSN::ADM_BSSN_gamma
+        ADMBase::lapse
+        ADMBase::shift
+        ADM_BSSN::ADM_BSSN_B
+"
+
+
+
+ActiveThorns = "Boundary CartGrid3D CoordBase ReflectionSymmetry RotatingSymmetry180 SymBase"
+
+CoordBase::domainsize = "minmax"
+
+CoordBase::xmin =    0.00
+CoordBase::ymin = -128.00
+CoordBase::zmin =    0.00
+CoordBase::xmax = +128.00
+CoordBase::ymax = +128.00
+CoordBase::zmax = +128.00
+CoordBase::dx   =    2.56
+CoordBase::dy   =    2.56
+CoordBase::dz   =    2.56
+
+CoordBase::boundary_size_x_lower     = 3
+CoordBase::boundary_size_z_lower     = 3
+CoordBase::boundary_shiftout_x_lower = 1
+CoordBase::boundary_shiftout_z_lower = 1
+
+CartGrid3D::type = "coordbase"
+
+ReflectionSymmetry::reflection_z   = yes
+ReflectionSymmetry::avoid_origin_z = no
+
+
+
+ActiveThorns = "CarpetMask"
+
+CarpetMask::verbose = yes
+
+CarpetMask::excluded_surface       [0] = 0
+CarpetMask::excluded_surface_factor[0] = 1.0
+
+CarpetMask::excluded_surface       [1] = 1
+CarpetMask::excluded_surface_factor[1] = 1.0
+
+CarpetMask::excluded_surface       [2] = 2
+CarpetMask::excluded_surface_factor[2] = 1.0
+
+
+
+ActiveThorns = "CarpetRegrid2"
+
+CarpetRegrid2::regrid_every = 128
+
+CarpetRegrid2::symmetry_rotating180 = yes
+
+CarpetRegrid2::num_centres = 2
+
+CarpetRegrid2::num_levels_1         =  7
+CarpetRegrid2::position_x_1         = +1.168642873
+CarpetRegrid2::radius_1[ 1]         =  64.0 # 1.28
+CarpetRegrid2::radius_1[ 2]         =  16.0 # 0.64
+CarpetRegrid2::radius_1[ 3]         =   8.0 # 0.32
+CarpetRegrid2::radius_1[ 4]         =   4.0 # 0.16
+CarpetRegrid2::radius_1[ 5]         =   2.0 # 0.08
+CarpetRegrid2::radius_1[ 6]         =   1.0 # 0.04
+CarpetRegrid2::movement_threshold_1 =   0.16
+
+CarpetRegrid2::num_levels_2         =  7
+CarpetRegrid2::position_x_2         = -1.168642873
+CarpetRegrid2::radius_2[ 1]         =  64.0 # 1.28
+CarpetRegrid2::radius_2[ 2]         =  16.0 # 0.64
+CarpetRegrid2::radius_2[ 3]         =   8.0 # 0.32
+CarpetRegrid2::radius_2[ 4]         =   4.0 # 0.16
+CarpetRegrid2::radius_2[ 5]         =   2.0 # 0.08
+CarpetRegrid2::radius_2[ 6]         =   1.0 # 0.04
+CarpetRegrid2::movement_threshold_2 =   0.16
+
+
+
+ActiveThorns = "MoL Time"
+
+MoL::ODE_Method             = "RK4"
+MoL::MoL_Intermediate_Steps = 4
+MoL::MoL_Num_Scratch_Levels = 1
+
+Time::dtfac = 0.40
+
+
+
+ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal TmunuBase"
+
+ADMMacros::spatial_order = 4
+
+
+
+ActiveThorns = "TwoPunctures"
+
+ADMBase::metric_type = "physical"
+
+ADMBase::initial_data  = "twopunctures"
+ADMBase::initial_lapse = "twopunctures-averaged"
+ADMBase::initial_shift = "zero"
+
+TwoPunctures::par_b          =  1.168642873
+TwoPunctures::par_m_plus     =  0.453
+TwoPunctures::par_m_minus    =  0.453
+TwoPunctures::par_P_plus [1] = +0.3331917498
+TwoPunctures::par_P_minus[1] = -0.3331917498
+
+#TODO# TwoPunctures::grid_setup_method = "evaluation"
+
+#TwoPunctures::TP_epsilon = 1.0e-4
+TwoPunctures::TP_epsilon = 1.0e-2
+TwoPunctures::TP_Tiny    = 1.0e-2
+
+TwoPunctures::verbose = yes
+
+
+
+ActiveThorns = "BSSN_MoL"
+
+ADMBase::evolution_method = "ADM_BSSN"
+
+ADM_BSSN::timelevels   = 3
+ADM_BSSN::stencil_size = 3
+ADM_BSSN::advection    = "upwind4"
+ADM_BSSN::bound        = "newrad"
+
+ADMBase::lapse_evolution_method = "1+log"
+ADM_BSSN::lapsesource           = "straight"
+ADM_BSSN::harmonic_f            = 2.0
+ADM_BSSN::lapse_advection_coeff = 1.0
+ADM_BSSN::force_lapse_positive  = yes
+
+ADMBase::shift_evolution_method        = "gamma0"
+ADM_BSSN::ShiftGammaCoeff              = 0.75
+ADM_BSSN::BetaDriver                   = 0.5
+ADM_BSSN::gamma_driver_advection_coeff = 1.0
+ADM_BSSN::ApplyShiftBoundary           = yes
+
+
+
+ActiveThorns = "Dissipation"
+
+Dissipation::order = 5
+Dissipation::vars  = "
+        ADM_BSSN::ADM_BSSN_phi
+        ADM_BSSN::ADM_BSSN_metric
+        ADM_BSSN::ADM_BSSN_curv
+        ADM_BSSN::ADM_BSSN_K
+        ADM_BSSN::ADM_BSSN_gamma
+        ADMBase::lapse
+        ADMBase::shift
+        ADM_BSSN::ADM_BSSN_B
+"
+
+
+
+ActiveThorns = "ADMConstraints"
+
+ADMConstraints::constraints_persist    = yes
+ADMConstraints::constraints_timelevels = 3
+
+
+
+ActiveThorns = "AHFinderDirect"
+
+AHFinderDirect::find_every = 128
+
+AHFinderDirect::run_at_CCTK_POST_RECOVER_VARIABLES = no
+
+AHFinderDirect::move_origins            = yes
+AHFinderDirect::reshape_while_moving    = yes
+AHFinderDirect::predict_origin_movement = yes
+
+AHFinderDirect::geometry_interpolator_name = "Lagrange polynomial interpolation"
+AHFinderDirect::geometry_interpolator_pars = "order=4"
+AHFinderDirect::surface_interpolator_name  = "Lagrange polynomial interpolation"
+AHFinderDirect::surface_interpolator_pars  = "order=4"
+
+AHFinderDirect::output_h_every = 0
+
+AHFinderDirect::N_horizons = 6
+
+AHFinderDirect::origin_x                             [1] = +1.168642873
+AHFinderDirect::initial_guess__coord_sphere__x_center[1] = +1.168642873
+AHFinderDirect::initial_guess__coord_sphere__radius  [1] =  0.25
+AHFinderDirect::which_surface_to_store_info          [1] = 0
+AHFinderDirect::reset_horizon_after_not_finding      [1] = no
+AHFinderDirect::dont_find_after_individual           [1] = 24000 # t=30.0
+
+AHFinderDirect::origin_x                             [2] = -1.168642873
+AHFinderDirect::initial_guess__coord_sphere__x_center[2] = -1.168642873
+AHFinderDirect::initial_guess__coord_sphere__radius  [2] =  0.25
+AHFinderDirect::which_surface_to_store_info          [2] = 1
+AHFinderDirect::reset_horizon_after_not_finding      [2] = no
+AHFinderDirect::dont_find_after_individual           [2] = 24000
+
+AHFinderDirect::initial_guess__coord_sphere__radius  [3] = 1.0
+AHFinderDirect::which_surface_to_store_info          [3] = 2
+AHFinderDirect::reset_horizon_after_not_finding      [3] = no
+AHFinderDirect::find_after_individual                [3] = 12000 # t=15.0
+
+AHFinderDirect::surface_definition                   [4] = "expansion product"
+AHFinderDirect::surface_selection                    [4] = "areal radius"
+AHFinderDirect::desired_value                        [4] = 50.0
+AHFinderDirect::initial_guess__coord_sphere__radius  [4] = 50.0
+AHFinderDirect::which_surface_to_store_info          [4] = 3
+AHFinderDirect::reset_horizon_after_not_finding      [4] = no
+
+AHFinderDirect::depends_on                           [5] = 1
+AHFinderDirect::desired_value_offset                 [5] = 0.001
+AHFinderDirect::which_surface_to_store_info          [5] = 4
+AHFinderDirect::reset_horizon_after_not_finding      [5] = no
+AHFinderDirect::dont_find_after_individual           [5] = 24000 # t=30.0
+
+AHFinderDirect::depends_on                           [6] = 3
+AHFinderDirect::desired_value_offset                 [6] = 0.001
+AHFinderDirect::which_surface_to_store_info          [6] = 5
+AHFinderDirect::reset_horizon_after_not_finding      [6] = no
+AHFinderDirect::find_after_individual                [6] = 12000 # t=15.0
+
+
+
+ActiveThorns = "IsolatedHorizon"
+
+IsolatedHorizon::verbose              = yes
+IsolatedHorizon::interpolator         = "Lagrange polynomial interpolation"
+IsolatedHorizon::interpolator_options = "order=4"
+IsolatedHorizon::spatial_order        = 6
+
+IsolatedHorizon::num_horizons = 6
+
+IsolatedHorizon::surface_index  [0] = 0
+IsolatedHorizon::companion_index[0] = 4
+
+IsolatedHorizon::surface_index  [1] = 1
+
+IsolatedHorizon::surface_index  [2] = 2
+
+IsolatedHorizon::surface_index  [3] = 3
+IsolatedHorizon::companion_index[3] = 5
+
+IsolatedHorizon::surface_index  [4] = 4
+
+IsolatedHorizon::surface_index  [5] = 5
+
+
+
+ActiveThorns = "ProperDistance ProperTime"
+
+ProperDistance::number_geodesics      = 512
+ProperDistance::direction             = "x"
+ProperDistance::plane                 = "xy"
+ProperDistance::opening_angle         = 180.0
+ProperDistance::step_size             =   0.04
+ProperDistance::horizon_number        = 1
+ProperDistance::interpolation_order   = 3
+ProperDistance::integration_method    = "rk4"
+ProperDistance::eps                   = 0.00001
+ProperDistance::max_proper_distance   = 15
+ProperDistance::calc_every            = 128
+ProperDistance::use_second_horizon    = yes
+ProperDistance::second_horizon_number = 2
+
+
+
+ActiveThorns = "PsiKadelia"
+
+PsiKadelia::psikadelia_persists     = yes
+PsiKadelia::calc_every              = 256   # only on levels 0 and 1
+PsiKadelia::ricci_prolongation_type = "none"
+PsiKadelia::weyl_timelevels         = 3
+PsiKadelia::psif_vec                = "standard-radial"
+
+
+
+ActiveThorns = "PunctureTracker"
+
+PunctureTracker::verbose = yes
+
+PunctureTracker::track    [0] = yes
+PunctureTracker::initial_x[0] = +1.168642873
+
+PunctureTracker::track    [1] = yes
+PunctureTracker::initial_x[1] = -1.168642873
+
+
+
+ActiveThorns = "SphericalHarmonics"
+
+SphericalHarmonics::grid_type                = "cart3d"
+SphericalHarmonics::interp_integration_order = 4
+SphericalHarmonics::InterpPointsTheta        = 41
+SphericalHarmonics::InterpPointsPhi          = 76
+
+SphericalHarmonics::number_of_radii = 5
+SphericalHarmonics::ex_radii[0] = 30
+SphericalHarmonics::ex_radii[1] = 40
+SphericalHarmonics::ex_radii[2] = 50
+SphericalHarmonics::ex_radii[3] = 60
+SphericalHarmonics::ex_radii[4] = 70
+
+SphericalHarmonics::number_of_vars = 10
+SphericalHarmonics::vars[0]           = "PsiKadelia::psi0re"
+SphericalHarmonics::vars[1]           = "PsiKadelia::psi0im"
+SphericalHarmonics::vars[2]           = "PsiKadelia::psi1re"
+SphericalHarmonics::vars[3]           = "PsiKadelia::psi1im"
+SphericalHarmonics::vars[4]           = "PsiKadelia::psi2re"
+SphericalHarmonics::vars[5]           = "PsiKadelia::psi2im"
+SphericalHarmonics::vars[6]           = "PsiKadelia::psi3re"
+SphericalHarmonics::vars[7]           = "PsiKadelia::psi3im"
+SphericalHarmonics::vars[8]           = "PsiKadelia::psi4re"
+SphericalHarmonics::vars[9]           = "PsiKadelia::psi4im"
+SphericalHarmonics::SH_spin_weight[0] = +2
+SphericalHarmonics::SH_spin_weight[1] = +2
+SphericalHarmonics::SH_spin_weight[2] = +1
+SphericalHarmonics::SH_spin_weight[3] = +1
+SphericalHarmonics::SH_spin_weight[4] =  0
+SphericalHarmonics::SH_spin_weight[5] =  0
+SphericalHarmonics::SH_spin_weight[6] = -1
+SphericalHarmonics::SH_spin_weight[7] = -1
+SphericalHarmonics::SH_spin_weight[8] = -2
+SphericalHarmonics::SH_spin_weight[9] = -2
+
+
+
+ActiveThorns = "WaveExtract"
+
+WaveExtract::out_every               = 128
+WaveExtract::maximum_detector_number = 5
+WaveExtract::switch_output_format    = 100
+WaveExtract::rsch2_computation       = "average Schwarzschild metric"
+WaveExtract::l_mode                  = 8
+WaveExtract::m_mode                  = 8
+WaveExtract::detector_radius[0]      = 30
+WaveExtract::detector_radius[1]      = 40
+WaveExtract::detector_radius[2]      = 50
+WaveExtract::detector_radius[3]      = 60
+WaveExtract::detector_radius[4]      = 70
+
+
+
+ActiveThorns = "CarpetIOBasic"
+
+IOBasic::outInfo_every      = 128
+IOBasic::outInfo_reductions = "norm2"
+IOBasic::outInfo_vars       = "
+        Carpet::timing
+        ADMConstraints::ham
+        IsolatedHorizon::ih_spin[0]
+        IsolatedHorizon::ih_radius[0]
+"
+
+
+
+ActiveThorns = "CarpetIOScalar"
+
+IOScalar::one_file_per_group = yes
+
+IOScalar::outScalar_every = 128
+IOScalar::outScalar_vars  = "
+        CarpetReduce::weight
+        ADMBase::metric
+        ADMBase::curv
+        ADMBase::lapse
+        ADMBase::shift
+        ADMConstraints::hamiltonian
+        ADMConstraints::momentum
+        SphericalSurface::sf_radius
+        IsolatedHorizon::ih_shapes
+        IsolatedHorizon::ih_coordinates
+        IsolatedHorizon::ih_tetrad_l
+        IsolatedHorizon::ih_tetrad_n
+        IsolatedHorizon::ih_tetrad_m
+        IsolatedHorizon::ih_newman_penrose
+        IsolatedHorizon::ih_weyl_scalars
+        IsolatedHorizon::ih_twometric
+        IsolatedHorizon::ih_killing_vector
+        IsolatedHorizon::ih_killed_twometric
+        IsolatedHorizon::ih_invariant_coordinates
+        IsolatedHorizon::ih_multipole_moments
+        IsolatedHorizon::ih_3determinant
+        PsiKadelia::WeylComponents
+        PunctureTracker::pt_loc
+        SphericalHarmonics::decomposed_vars
+"
+
+
+
+ActiveThorns = "CarpetIOASCII"
+
+IOASCII::one_file_per_group = yes
+
+IOASCII::output_symmetry_points = no
+IOASCII::out3D_ghosts           = no
+
+IOASCII::out0D_every = 128
+IOASCII::out0D_vars  = "
+        Carpet::timing
+        CarpetReduce::weight
+        ADMBase::metric
+        ADMBase::curv
+        ADMBase::lapse
+        ADMBase::shift
+        ADMConstraints::hamiltonian
+        ADMConstraints::momentum
+        SphericalSurface::sf_active
+        SphericalSurface::sf_valid
+        SphericalSurface::sf_info
+        SphericalSurface::sf_radius
+        SphericalSurface::sf_origin
+        SphericalSurface::sf_coordinate_descriptors
+        IsolatedHorizon::ih_shapes
+        IsolatedHorizon::ih_state
+        IsolatedHorizon::ih_grid_int
+        IsolatedHorizon::ih_grid_real
+        IsolatedHorizon::ih_shapes
+        IsolatedHorizon::ih_coordinates
+        IsolatedHorizon::ih_tetrad_l
+        IsolatedHorizon::ih_tetrad_n
+        IsolatedHorizon::ih_tetrad_m
+        IsolatedHorizon::ih_newman_penrose
+        IsolatedHorizon::ih_weyl_scalars
+        IsolatedHorizon::ih_twometric
+        IsolatedHorizon::ih_killing_vector
+        IsolatedHorizon::ih_killed_twometric
+        IsolatedHorizon::ih_scalars
+        IsolatedHorizon::ih_invariant_coordinates
+        IsolatedHorizon::ih_multipole_moments
+        IsolatedHorizon::ih_3determinant
+        PsiKadelia::WeylComponents
+        SphericalHarmonics::decomposed_vars
+        ProperDistance::pdistance
+        ProperTime::ptime
+"
+
+IOASCII::out1D_every = 128
+IOASCII::out1D_vars  = "
+        CarpetReduce::weight
+        ADMBase::metric
+        ADMBase::curv
+        ADMBase::lapse
+        ADMBase::shift
+        ADMConstraints::hamiltonian
+        ADMConstraints::momentum
+        SphericalSurface::sf_radius
+        IsolatedHorizon::ih_shapes
+        IsolatedHorizon::ih_weyl_scalars
+        IsolatedHorizon::ih_killing_vector
+        IsolatedHorizon::ih_killed_twometric
+        IsolatedHorizon::ih_3determinant
+        PsiKadelia::WeylComponents
+        SphericalHarmonics::decomposed_vars
+"
+
+IOASCII::out2D_every = 128
+IOASCII::out2D_vars  = "
+        SphericalSurface::sf_radius
+        SphericalHarmonics::decomposed_vars
+"
+
+
+
+Activethorns = "CarpetIOHDF5"
+
+IOHDF5::out_every          = 512
+IOHDF5::one_file_per_group = yes
+IOHDF5::compression_level  = 1
+IOHDF5::out_vars  = "
+        CarpetReduce::weight
+        ADMBase::metric
+        ADMBase::curv
+        ADMBase::lapse
+        ADMBase::shift
+        ADMConstraints::hamiltonian
+        ADMConstraints::momentum
+        PsiKadelia::WeylComponents
+"
+
+IOHDF5::checkpoint          = yes
+IO::checkpoint_dir          = $parfile
+IO::checkpoint_ID           = yes
+IO::checkpoint_every        = 512
+IO::checkpoint_on_terminate = yes
+
+IO::recover     = "autoprobe"
+IO::recover_dir = $parfile
+
+
+
+ActiveThorns = "Formaline"
+
+Formaline::collect_metadata = no
+
+#Formaline::send_as_rdf     = yes
+#Formaline::rdf_hostname[0] = "numrel07.cct.lsu.edu"
+#Formaline::rdf_port    [0] = 24000
+
+
+
+ActiveThorns = "TimerReport"
+
+TimerReport::out_every    = 512
+TimerReport::out_filename = "TimerReport"
diff --git a/param.ccl b/param.ccl
index cfa152d..1475ebb 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,61 +1,24 @@
 # Parameter definitions for thorn PunctureTracker
-# $Header$
 
-INT pt_verbose "speak up?"
+BOOLEAN verbose "speak up?" STEERABLE=always
 {
-0:10 :: "eloquence"
-} 1
-
-INT pt_num_tracked "Number of punctures that should be tracked"
-{
-  0:10 :: ""
-} 0
-
-# not to be used 
-# because the puncture has to be tracked every timestep
-#INT pt_every "how often to calculate new position"
-#{
-#  0:*::""
-#} 0
+} "no"
 
-BOOLEAN pt_track_punctures "do it?"
+BOOLEAN track[10] "Track this puncture"
 {
 } "no"
 
-REAL pt_initial_x[10] "Initial x coordinate positions of BHs"
-#This overwrites parameter par_b from Twopunctures if nonzero
+REAL initial_x[10] "Initial x coordinate positions of punctures"
 {
-  :  :: "Anything goes"
+  *:* :: ""
 } 0.0
 
-REAL pt_initial_y[10] "Initial y coordinate positions of BHs"
-#This overwrites parameter par_b from Twopunctures if nonzero
+REAL initial_y[10] "Initial y coordinate positions of punctures"
 {
-  :  :: "Anything goes"
+  *:* :: ""
 } 0.0
 
-REAL pt_initial_z[10] "Initial z coordinate positions of BHs"
-#This overwrites parameter par_b from Twopunctures if nonzero
+REAL initial_z[10] "Initial z coordinate positions of punctures"
 {
-  :  :: "Anything goes"
+  *:* :: ""
 } 0.0
-
-BOOLEAN pt_regrid "should we use PT to Regrid?"
-{
-} "no"
-
-REAL pt_regrid_buffer "regrid if new location-old_location > buffer"
-{
-  : ::"anything greater zero"
-} 0.5
-
-INT pt_which_surface_to_take[10] "which surface"
-{  
-0:42 :: "" 
-} 0 
-
-
-SHARES: TwoPunctures
-
-USES REAL par_b 
- 
diff --git a/schedule.ccl b/schedule.ccl
index ab3b151..cce26af 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,27 +1,21 @@
 # Schedule definitions for thorn PunctureTracker
-# $Header$
 
+STORAGE: pt_loc pt_loc_p
 
-storage: pt_loc[2]
-storage: pt_shift
+SCHEDULE PunctureTracker_Init AT initial
+{
+  LANG: C
+  OPTIONS: GLOBAL
+} "Calculate initial location of punctures"
 
-if(pt_track_punctures){
+SCHEDULE PunctureTracker_Track AT evol AFTER MoL_Evolution
+{
+  LANG: C
+  OPTIONS: GLOBAL
+} "Calculate new location of punctures"
 
-  SCHEDULE PunctureTracker_init at CCTK_INITIAL
-  {
-    LANG: C
-    OPTIONS: GLOBAL
-  } "calculate new location of punctures"
-
-  SCHEDULE PunctureTracker at CCTK_POSTSTEP before SphericalSurface_HasBeenSet
-  {
-    LANG: C
-    OPTIONS: GLOBAL
-  } "calculate new location of punctures"
-  
-  SCHEDULE update_punc_loc at CCTK_POSTSTEP after PunctureTracker before SphericalSurface_HasBeenSet
-  {
-    LANG: C
-    OPTIONS: GLOBAL
-  } "update new location of punctures"
-}
+SCHEDULe PunctureTracker_SetPositions AT preregrid
+{
+  LANG: C
+  OPTIONS: global
+} "Set positions of refined regions"
diff --git a/src/adjust.c b/src/adjust.c
deleted file mode 100644
index 43a2df0..0000000
--- a/src/adjust.c
+++ /dev/null
@@ -1,90 +0,0 @@
-/* $Header$ */
-
-#include <assert.h>
-#include <stdio.h>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-
-
-void
-BHTracker_Adjust (CCTK_ARGUMENTS);
-
-static void
-setpar (char const * name, int index, char const * thorn, CCTK_REAL value);
-
-
-
-void
-BHTracker_Adjust (CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_PARAMETERS;
-  DECLARE_CCTK_ARGUMENTS;
-  
-  int n;
-  int si, ri;
-  char name[100], val[100];
-  int icnt, ierr;
-
-  CCTK_REAL xpos,ypos,zpos;
-  
-  if ((cctk_iteration-1) < track_first) return;
-  if ((cctk_iteration-1) % track_every != 0) return;
-
-
-
-  CCTK_INFO ("Tracking...");
-  
-  for (n=0; n<num_tracked; ++n) {
-    
-    si = surface_index[n];
-    assert (si>=0 && si<nsurfaces);
-    
-    ri = region_index[n];
-    assert (ri>=0 && ri<num_offsets);
-    
-    if (sf_valid[si] > 0) {
-      
-      CCTK_VInfo (CCTK_THORNSTRING,
-                  "   Object #%d has moved to [%g,%g,%g]", n,
-                  (double)sf_centroid_x[si], (double)sf_centroid_y[si],
-                  (double)sf_centroid_z[si]);
-
-      /* calculate new positions */
-      xpos = (sf_centroid_x[si] - initial_x[si]);
-      ypos = (sf_centroid_y[si] - initial_y[si]);
-      zpos = (sf_centroid_z[si] - initial_z[si]);
-
-      setpar ("offsetx", ri, "CarpetRegrid", xpos);
-      setpar ("offsety", ri, "CarpetRegrid", ypos);
-      setpar ("offsetz", ri, "CarpetRegrid", zpos);
-      
-    } else {
-      
-      CCTK_VInfo (CCTK_THORNSTRING,
-                  "   Object #%d cannot be detected at this time", n);
-      
-    }
-  }
-}
-
-
-
-static void
-setpar (char const * const name, int const index, char const * const thorn,
-        CCTK_REAL const value)
-{
-  char nambuf[100], valbuf[100];
-  int icnt, ierr;
-  
-  icnt = snprintf (nambuf, sizeof nambuf, "%s[%d]", name, index);
-  assert (icnt>=0 && icnt<sizeof nambuf);
-  
-  icnt = snprintf (valbuf, sizeof valbuf, "%g", value);
-  assert (icnt>=0 && icnt<sizeof valbuf);
-  
-  ierr = CCTK_ParameterSet (nambuf, thorn, valbuf);
-  assert (! ierr);
-}
diff --git a/src/make.code.defn b/src/make.code.defn
index e0e5f26..31baa1c 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,8 +1,7 @@
-# Main make.code.defn file for thorn PunctureTracker
-# $Header$
+# Main make.code.defn file for thorn PunctureTracker -*-Makefile-*-
 
 # Source files in this directory
-SRCS = punc_track.c
+SRCS = puncture_tracker.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/punc_track.c b/src/punc_track.c
deleted file mode 100644
index 2b9bd26..0000000
--- a/src/punc_track.c
+++ /dev/null
@@ -1,170 +0,0 @@
-/* $Header$ */
- 
-#include <assert.h>
-#include <stdio.h>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-#include "util_Table.h"
-
-#define REFLEVEL round(log10((CCTK_REAL)(cctkGH->cctk_levfac[0]))/log10(2.0))
-
-void PunctureTracker (CCTK_ARGUMENTS);
-void PunctureTracker_init (CCTK_ARGUMENTS);
-
-
-void PunctureTracker_init (CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_PARAMETERS;
-  DECLARE_CCTK_ARGUMENTS;
-  int i;
-
-  if(pt_verbose>0)
-    CCTK_INFO("Initializing PunctureTracker");
-
-  if((pt_initial_x[0]!=0.0)||(pt_initial_y[0]!=0.0)||(pt_initial_z[0]!=0.0))
-    for (i=0;i<pt_num_tracked;i++){
-      pt_loc_x_p[i] = pt_initial_x[i];
-      pt_loc_y_p[i] = pt_initial_y[i];
-      pt_loc_z_p[i] = pt_initial_z[i];
-    }
-  else{
-    pt_loc_x_p[0] = par_b;
-    pt_loc_x_p[1] =-par_b;
-    pt_loc_y_p[0] = 0.0;
-    pt_loc_y_p[1] = 0.0;
-    pt_loc_z_p[0] = 0.0;
-    pt_loc_z_p[1] = 0.0;
-  }
-  if(pt_verbose>0)
-    CCTK_INFO("done initializing PunctureTracker");
-}
-
-void PunctureTracker (CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_PARAMETERS;
-  DECLARE_CCTK_ARGUMENTS;
-  
-  int n,ierror,i,numpoints;
-  // portland doesnt like that
-  //CCTK_REAL  pt_betax[pt_num_tracked], pt_betay[pt_num_tracked], pt_betaz[pt_num_tracked];
-  //so I do that
-
-  CCTK_REAL * restrict const pt_betax = malloc(pt_num_tracked*sizeof(CCTK_REAL));
-  CCTK_REAL * restrict const pt_betay = malloc(pt_num_tracked*sizeof(CCTK_REAL));
-  CCTK_REAL * restrict const pt_betaz = malloc(pt_num_tracked*sizeof(CCTK_REAL));
-
-  int         operator_handle, param_table_handle, coordsys_handle;
-  int         input_array_indices[3];
-  const void *interp_coords[3];
-  void       *output_arrays[3];
-
-  const CCTK_INT input_array_type_codes[3] = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL};   
-  const CCTK_INT output_array_type_codes[3]= { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL};
-   
-  if(pt_verbose >0)
-    CCTK_INFO("tracking Punctures");
-  
-  interp_coords[0] = (const void *) pt_loc_x_p;
-  interp_coords[1] = (const void *) pt_loc_y_p;
-  interp_coords[2] = (const void *) pt_loc_z_p;
-  
-  operator_handle     = -1;
-  coordsys_handle     = -1;
-  param_table_handle  = -1;
-  
-  // set interpolation handles
-  operator_handle = CCTK_InterpHandle("Lagrange polynomial interpolation");
-  if (operator_handle < 0) CCTK_WARN(0, "can’t get interpolation handle!");
-
-  param_table_handle = Util_TableCreateFromString("order=4");
-  if (param_table_handle < 0) CCTK_WARN(0, "can't ge parameter table!");
-  
-  coordsys_handle = CCTK_CoordSystemHandle("cart3d");
-  if (coordsys_handle < 0) CCTK_WARN(0,"can't get coordsys handle!");
-  if(pt_verbose>0)
-    for (i=0; i<pt_num_tracked; i++) {
-      CCTK_VInfo(CCTK_THORNSTRING,
-      "Puncture %d is at (x,y,z)=(%f,%f,%f)\n",i,pt_loc_x_p[i],pt_loc_y_p[i],pt_loc_z_p[i]);
-  }
-  
-  output_arrays[0] = (void *) pt_betax;
-  output_arrays[1] = (void *) pt_betay;
-  output_arrays[2] = (void *) pt_betaz;
-
-  input_array_indices[0] = CCTK_VarIndex ("ADMBase::betax");
-  input_array_indices[1] = CCTK_VarIndex ("ADMBase::betay");
-  input_array_indices[2] = CCTK_VarIndex ("ADMBase::betaz");
-
-  //only proc 0 does interpolation
-  numpoints = 0;
-  if (CCTK_MyProc(cctkGH) == 0)
-    numpoints = pt_num_tracked;
-
-  ierror = CCTK_InterpGridArrays(cctkGH,3,
-                                 operator_handle, param_table_handle,coordsys_handle,
-                                 numpoints,
-                                 CCTK_VARIABLE_REAL,
-                                 interp_coords,
-                                 3,input_array_indices,
-                                 3,output_array_type_codes,
-                                  output_arrays);
-  if (ierror<0)
-    CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
-               "error return from interpolator ierror=%d",(int)ierror);
-  
-
-  for (n=0; n< pt_num_tracked; n++) {
-    if(pt_verbose>1)
-      printf("at the location of puncture %i:\n              is pt_betax=%f, pt_betay=%f,pt_betaz=%f \n",n,pt_betax[n],pt_betay[n],pt_betaz[n]);
-        pt_shiftx[n] = pt_betax[n];
-        pt_shifty[n] = pt_betay[n];
-        pt_shiftz[n] = pt_betaz[n];
-  }
-
-    free(pt_betax,pt_betay,pt_betaz);
-  if(pt_verbose>0)
-    CCTK_INFO ("done");
-}
-
-void update_punc_loc(CCTK_ARGUMENTS){
-
-  DECLARE_CCTK_PARAMETERS;
-  DECLARE_CCTK_ARGUMENTS;
-
-  CCTK_REAL pt_dt;
-  int n;
-
-  pt_dt = CCTK_DELTA_TIME;
-
-  for (n=0; n< pt_num_tracked; n++) {
-    //correct small errors
-    //    if (pt_shiftx[n] < 10e-14) pt_shiftx[n] = 0.0;
-    //    if (pt_shifty[n] < 10e-14) pt_shifty[n] = 0.0;
-    //    if (pt_shiftz[n] < 10e-14) pt_shiftz[n] = 0.0;
-    
-    if(pt_verbose > 2){
-      printf("updating puncture %i\n",n);
-      printf("   pt_dt      is %f\n",pt_dt);
-      printf("   pt_loc_x   is %f\n",pt_loc_x[n]);
-      printf("   pt_loc_x_p is %f\n",pt_loc_x_p[n]);
-      printf("   pt_shiftx  is %f\n",pt_shiftx[n]);
-      printf("   pt_dt      is %f\n",pt_dt);
-      fflush(stdout);
-    }
-     
-     printf("WRITING TO SURFACE %d\n",pt_which_surface_to_take[n]);
-    pt_loc_x[n] = pt_loc_x_p[n] - pt_shiftx[n] * pt_dt;
-    pt_loc_y[n] = pt_loc_y_p[n] - pt_shifty[n] * pt_dt;
-    pt_loc_z[n] = pt_loc_z_p[n] - pt_shiftz[n] * pt_dt;
-
-    sf_centroid_x[pt_which_surface_to_take[n]]=pt_loc_x_p[n];    
-    sf_centroid_y[pt_which_surface_to_take[n]]=pt_loc_y_p[n];
-    sf_centroid_z[pt_which_surface_to_take[n]]=pt_loc_z_p[n];
-     
-    sf_valid[pt_which_surface_to_take[n]]=1.0;
-    sf_valid[pt_which_surface_to_take[n]]=1.0;
-    sf_valid[pt_which_surface_to_take[n]]=1.0;
-  }
-}
diff --git a/src/puncture_tracker.c b/src/puncture_tracker.c
new file mode 100644
index 0000000..d3da795
--- /dev/null
+++ b/src/puncture_tracker.c
@@ -0,0 +1,285 @@
+#include <assert.h>
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_Table.h"
+
+
+
+static int const max_num_tracked = 10;
+
+
+
+void
+PunctureTracker_Init (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  if (verbose) {
+    CCTK_INFO ("Initializing PunctureTracker");
+  }
+  
+  for (int n = 0; n < max_num_tracked; ++ n) {
+    if (track[n]) {
+      pt_loc_t[n] = cctk_time;
+      pt_loc_x[n] = initial_x[n];
+      pt_loc_y[n] = initial_y[n];
+      pt_loc_z[n] = initial_z[n];
+    } 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;
+      pt_loc_z[n] = 0.0;
+    }
+    pt_loc_t_p[n] = 0.0;
+    pt_loc_x_p[n] = 0.0;
+    pt_loc_y_p[n] = 0.0;
+    pt_loc_z_p[n] = 0.0;
+  }
+}
+
+
+
+void
+PunctureTracker_Track (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  // 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...");
+  }
+  
+  if (verbose) {
+    for (int n = 0; n < max_num_tracked; ++ n) {
+      if (track[n]) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Puncture #%d is at (%g,%g,%g)",
+                    n,
+                    (double)pt_loc_x[n],
+                    (double)pt_loc_y[n],
+                    (double)pt_loc_z[n]);
+      }
+    }
+  }
+  
+  // Manual time level cycling
+  
+  for (int n = 0; n < max_num_tracked; ++ n) {
+    if (track[n]) {
+      pt_loc_t_p[n] = pt_loc_t[n];
+      pt_loc_x_p[n] = pt_loc_x[n];
+      pt_loc_y_p[n] = pt_loc_y[n];
+      pt_loc_z_p[n] = pt_loc_z[n];
+      
+      pt_loc_t[n] = cctk_time;
+    }
+  }
+  
+  // Interpolate
+  
+  // Dimensions
+  int const dim = 3;
+  
+  // Interpolation operator
+  int const operator_handle =
+    CCTK_InterpHandle ("Lagrange polynomial interpolation");
+  if (operator_handle < 0) {
+    CCTK_WARN (CCTK_WARN_ALERT, "Can’t get interpolation handle");
+    goto label_return;
+  }
+  
+  // Interpolation parameter table
+  int const param_table_handle = Util_TableCreateFromString ("order=4");
+  if (param_table_handle < 0) {
+    CCTK_WARN (CCTK_WARN_ALERT, "Can't create parameter table");
+    goto label_return;
+  }
+  
+  {
+    
+    // Interpolation coordinate system
+    int const 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
+    int const num_points = CCTK_MyProc(cctkGH) == 0 ? max_num_tracked : 0;
+    
+    // Interpolation coordinates
+    assert (dim == 3);
+    CCTK_POINTER_TO_CONST interp_coords[3];
+    interp_coords[0] = pt_loc_x_p;
+    interp_coords[1] = pt_loc_y_p;
+    interp_coords[2] = pt_loc_z_p;
+    
+    // Number of interpolation variables
+    int const num_vars = 3;
+    
+    // Interpolated variables
+    assert (num_vars == 3);
+    int input_array_indices[3];
+    input_array_indices[0] = CCTK_VarIndex ("ADMBase::betax");
+    input_array_indices[1] = CCTK_VarIndex ("ADMBase::betay");
+    input_array_indices[2] = CCTK_VarIndex ("ADMBase::betaz");
+    
+    // Interpolation result types
+    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;
+    
+    // Interpolation result
+    CCTK_REAL pt_betax[max_num_tracked];
+    CCTK_REAL pt_betay[max_num_tracked];
+    CCTK_REAL pt_betaz[max_num_tracked];
+    
+    assert (num_vars == 3);
+    CCTK_POINTER output_arrays[3];
+    output_arrays[0] = pt_betax;
+    output_arrays[1] = pt_betay;
+    output_arrays[2] = pt_betaz;
+    
+    // Interpolate
+    int const ierr = CCTK_InterpGridArrays
+      (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);
+    if (ierr < 0) {
+      CCTK_WARN (CCTK_WARN_ALERT, "Interpolation error");
+      goto label_free_param_table;
+    }
+    
+    if (CCTK_MyProc(cctkGH) == 0) {
+      
+      // Some more output
+      
+      if (verbose && CCTK_MyProc(cctkGH) == 0) {
+        for (int n = 0; n < max_num_tracked; ++ n) {
+          if (track[n]) {
+          CCTK_VInfo (CCTK_THORNSTRING,
+                      "Shift at puncture #%d is at (%g,%g,%g)",
+                      n,
+                      (double)pt_betax[n],
+                      (double)pt_betay[n],
+                      (double)pt_betaz[n]);
+          }
+        }
+      }
+      
+      // Time evolution
+      
+      for (int n = 0; n < max_num_tracked; ++ n) {
+        if (track[n]) {
+          CCTK_REAL const dt = pt_loc_t[n] - pt_loc_t_p[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] = pt_loc_x_p[n] + dt * (- pt_betax[n]);
+          pt_loc_y[n] = pt_loc_y_p[n] + dt * (- pt_betay[n]);
+          pt_loc_z[n] = pt_loc_z_p[n] + dt * (- pt_betaz[n]);
+        }
+      }
+      
+    }
+    
+    // Broadcast result
+    
+    CCTK_REAL loc_local[3*max_num_tracked];
+    if (CCTK_MyProc(cctkGH) == 0) {
+      for (int n = 0; n < max_num_tracked; ++ n) {
+        loc_local[                  n] = pt_loc_x[n];
+        loc_local[  max_num_tracked+n] = pt_loc_y[n];
+        loc_local[2*max_num_tracked+n] = pt_loc_z[n];
+      }
+    } else {
+      for (int n = 0; n < max_num_tracked; ++ n) {
+        loc_local[                  n] = 0.0;
+        loc_local[  max_num_tracked+n] = 0.0;
+        loc_local[2*max_num_tracked+n] = 0.0;
+      }
+    }
+    
+    CCTK_REAL loc_global[3*max_num_tracked];
+    
+    int const handle_sum = CCTK_ReductionHandle ("sum");
+    if (handle_sum < 0) {
+      CCTK_WARN (CCTK_WARN_ALERT, "Can't get redunction handle");
+      goto label_free_param_table;
+    }
+    
+    int const ierr2 = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, handle_sum,
+       loc_local, loc_global, 3*max_num_tracked, CCTK_VARIABLE_REAL);
+    if (ierr2 < 0) {
+      CCTK_WARN (CCTK_WARN_ALERT, "Reduction error");
+      goto label_free_param_table;
+    }
+    
+    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];
+    }
+    
+  }
+  
+  // Done
+  
+  // Poor man's exception handling
+ label_free_param_table:
+  Util_TableDestroy (param_table_handle);
+  
+ label_return:
+  return;
+}
+
+
+
+void
+PunctureTracker_SetPositions (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  for (int n = 0; n < max_num_tracked; ++ n) {
+    if (track[n]) {
+      
+      if (verbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Setting position of refined region from puncture #%d to (%g,%g,%g)",
+                    n,
+                    (double)pt_loc_x[n],
+                    (double)pt_loc_y[n],
+                    (double)pt_loc_z[n]);
+      }
+      
+      // Set position in CarpetRegrid2
+      position_x[n] = pt_loc_x[n];
+      position_y[n] = pt_loc_y[n];
+      position_z[n] = pt_loc_z[n];
+      
+    }
+  }
+}