add PunctureTracker

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

8 files changed, 495 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..14af522
--- /dev/null
+++ b/README
@@ -0,0 +1,11 @@
+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>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+time integrate shift at location of puncture to follow location of puncture
+
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..f9b451f
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,122 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header$
+%
+%  Thorn documentation in the latex file doc/documentation.tex
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include
+%  pages about variables, parameters and scheduling parsed from the
+%  relevant thorn CCL files.
+%
+%  This template contains guidelines which help to assure that your
+%  documentation will be correctly added to ThornGuides. More
+%  information is available in the Cactus UsersGuide.
+%
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations',
+%     'references', and 'image' names should conform to the following
+%     convention:
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package.
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command.
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as
+%     standard sections.
+%   - For the benefit of our Perl scripts, and for future extensions,
+%     please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header$
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{BHTracker}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+\date{$ $Date$ $}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+   \texttt{BHTracker} tracks black holes or other features as they
+   move over the grid and makes Carpet's refinement regions follow
+   them.
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Overview}
+
+   \texttt{BHTracker} expects that some other mechanism has performed
+   the actual tracking and stored the information about the tracked
+   black hole location in a \texttt{SphericalSurface}.  It can thus
+   track any object, not just black holes.  \texttt{BHTracker} then
+   steers certain parameters of \texttt{CarpetRegrid} to make it move
+   the refined regions.  \texttt{BHTracker} is not intelligent, it is
+   just a glob of glue to make ends meet.
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..8c79d3f
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,20 @@
+# Interface definition for thorn PunctureTracker
+# $Header$
+
+IMPLEMENTS: PunctureTracker
+
+INHERITS: ADMBase, TwoPunctures,SphericalSurface
+
+
+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
+{
+ pt_shiftx, pt_shifty, pt_shiftz
+
+} "shift at location of punctures"
+
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..61914c7
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,53 @@
+# Parameter definitions for thorn PunctureTracker
+# $Header$
+
+INT pt_verbose "speak up?"
+{
+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
+
+BOOLEAN pt_track_punctures "do it?"
+{
+} "no"
+
+REAL pt_initial_x[10] "Initial x coordinate positions of BHs"
+#This overwrites parameter par_b from Twopunctures if nonzero
+{
+  :  :: "Anything goes"
+} 0.0
+
+REAL pt_initial_y[10] "Initial y coordinate positions of BHs"
+#This overwrites parameter par_b from Twopunctures if nonzero
+{
+  :  :: "Anything goes"
+} 0.0
+
+REAL pt_initial_z[10] "Initial z coordinate positions of BHs"
+#This overwrites parameter par_b from Twopunctures if nonzero
+{
+  :  :: "Anything goes"
+} 0.0
+
+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
new file mode 100644
index 0000000..f95f829
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,27 @@
+# Schedule definitions for thorn PunctureTracker
+# $Header$
+
+
+storage: pt_loc[2]
+storage: pt_shift
+
+if(pt_track_punctures){
+
+  SCHEDULE PunctureTracker_init at CCTK_INITIAL
+  {
+    LANG: C
+    OPTIONS: GLOBAL
+  } "calculate new location of punctures"
+
+  SCHEDULE PunctureTracker at CCTK_POSTSTEP
+  {
+    LANG: C
+    OPTIONS: GLOBAL
+  } "calculate new location of punctures"
+  
+  SCHEDULE update_punc_loc at CCTK_POSTSTEP after PunctureTracker
+  {
+    LANG: C
+    OPTIONS: GLOBAL
+  } "update new location of punctures"
+}
diff --git a/src/adjust.c b/src/adjust.c
new file mode 100644
index 0000000..43a2df0
--- /dev/null
+++ b/src/adjust.c
@@ -0,0 +1,90 @@
+/* $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
new file mode 100644
index 0000000..e0e5f26
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn PunctureTracker
+# $Header$
+
+# Source files in this directory
+SRCS = punc_track.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/punc_track.c b/src/punc_track.c
new file mode 100644
index 0000000..d4b83fa
--- /dev/null
+++ b/src/punc_track.c
@@ -0,0 +1,164 @@
+/* $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);
+    }
+    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[n];    
+    sf_centroid_y[pt_which_surface_to_take[n]]=pt_loc_y[n];
+    sf_centroid_z[pt_which_surface_to_take[n]]=pt_loc_z[n];
+  }
+}