some bugfixes

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

2 files changed, 140 insertions, 43 deletions

diff --git a/param.ccl b/param.ccl
index 61914c7..cfa152d 100644
--- a/param.ccl
+++ b/param.ccl
@@ -40,14 +40,22 @@ REAL pt_initial_z[10] "Initial z coordinate positions of BHs"
   :  :: "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
+0:42 :: "" 
+} 0 
 
 
 SHARES: TwoPunctures
 
 USES REAL par_b 
-
-
+ 
diff --git a/src/punc_track.c b/src/punc_track.c
index 7fb26b9..69cb754 100644
--- a/src/punc_track.c
+++ b/src/punc_track.c
@@ -7,9 +7,15 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 #include "util_Table.h"
-
+#include <math.h>
 #define REFLEVEL round(log10((CCTK_REAL)(cctkGH->cctk_levfac[0]))/log10(2.0))
 
+#define ABS(x) ((x)>=0?x:-x)
+
+double pt_loc_old_x[10];
+double pt_loc_old_y[10];
+double pt_loc_old_z[10];
+
 void PunctureTracker (CCTK_ARGUMENTS);
 void PunctureTracker_init (CCTK_ARGUMENTS);
 
@@ -37,10 +43,23 @@ void PunctureTracker_init (CCTK_ARGUMENTS)
     pt_loc_z_p[0] = 0.0;
     pt_loc_z_p[1] = 0.0;
   }
-  if(pt_verbose>0)
-    CCTK_INFO("done initializing PunctureTracker");
+  for (i=0; i< pt_num_tracked; i++) {
+    pt_loc_old_x[i] = pt_loc_x_p[i];
+    pt_loc_old_y[i] = pt_loc_y_p[i];
+    pt_loc_old_z[i] = pt_loc_z_p[i];
+  
+    sf_centroid_x[pt_which_surface_to_take[i]]=pt_loc_old_x[i];
+    sf_centroid_y[pt_which_surface_to_take[i]]=pt_loc_old_y[i];      
+    sf_centroid_z[pt_which_surface_to_take[i]]=pt_loc_old_z[i];      
+
+    if(pt_verbose>0){
+      printf("set puncture %d to location=(%g,%g,%g)\n",i,pt_loc_old_x[i],pt_loc_old_y[i],pt_loc_old_z[i]);
+    }
+  }
+      CCTK_INFO("done initializing PunctureTracker");
 }
 
+
 void PunctureTracker (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_PARAMETERS;
@@ -51,9 +70,9 @@ void PunctureTracker (CCTK_ARGUMENTS)
   //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));
+  CCTK_REAL pt_betax[10];
+  CCTK_REAL pt_betay[10];
+  CCTK_REAL pt_betaz[10];
 
   int         operator_handle, param_table_handle, coordsys_handle;
   int         input_array_indices[3];
@@ -83,12 +102,7 @@ void PunctureTracker (CCTK_ARGUMENTS)
   
   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;
@@ -116,16 +130,15 @@ void PunctureTracker (CCTK_ARGUMENTS)
   
 
   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];
+    if(pt_verbose>4)
+      printf("at location of puncture %i is \
+              (beta,betay,betaz) = (%f,%f,%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){
@@ -133,32 +146,108 @@ void update_punc_loc(CCTK_ARGUMENTS){
   DECLARE_CCTK_PARAMETERS;
   DECLARE_CCTK_ARGUMENTS;
 
-  CCTK_REAL pt_dt;
-  int n;
-
+  CCTK_REAL pt_dt,pt_x_diff,pt_y_diff,pt_z_diff;
+  int n,ierr;
+  int regrid_now = 0;
+  char parname[200];
+  char value[200];
+  
   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);
-    }
+    //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;    
+  
+    //update
     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];
+    //info
+    if(pt_verbose>0) CCTK_VInfo(CCTK_THORNSTRING,
+      "Puncture %d is at (x,y,z)=(%f,%f,%f)",n,pt_loc_x[n],pt_loc_y[n],pt_loc_z[n]);
+    if(pt_verbose > 2){
+      printf("     pt_dt      is %f\n",pt_dt);
+      printf("     pt_loc_x_p is %f\n",pt_loc_x_p[n]);
+      printf("     pt_shiftx  is %f\n",pt_shiftx[n]);
+      fflush(stdout);
+    }
+  
+
+    //have we moved enough to regrid now?
+    printf("regrid_now is %i",regrid_now);
+    if((pt_regrid)&&(cctkGH->cctk_iteration>1)){
+      if(!regrid_now){
+        pt_x_diff = ABS(pt_loc_old_x[n] - pt_loc_x[n]);
+        pt_y_diff = ABS(pt_loc_old_y[n] - pt_loc_y[n]);
+        pt_z_diff = ABS(pt_loc_old_z[n] - pt_loc_z[n]);
+
+        if(pt_verbose > 2){
+          printf("     pt_x_diff      is %f\n",pt_x_diff);
+          printf("     pt_y_diff      is %f\n",pt_y_diff);
+          printf("     pt_z_diff      is %f\n",pt_z_diff);
+          printf("     pt_reg_buf : %f\n",pt_regrid_buffer);
+        }
+        
+        if((pt_x_diff > pt_regrid_buffer)||
+           (pt_y_diff > pt_regrid_buffer)||
+           (pt_z_diff > pt_regrid_buffer)){
+          regrid_now = 1;
+        }
+      }
+    }
   }
+  
+  if(regrid_now){
+    regrid_now = 0;
+
+    if(pt_verbose>0) CCTK_INFO("regridding now");
+    for (n=0; n< pt_num_tracked; n++) {
+      ierr = CCTK_ParameterSet ("regrid_every","CarpetReGrid2","1");
+      assert (ierr>=0);
+
+      ierr = sprintf (parname, "position_x_%d", n+1);
+      assert (ierr>=0);
+      ierr = sprintf (value, "%f", pt_loc_x[n]);
+      assert (ierr>=0);
+      ierr = CCTK_ParameterSet(parname,"CarpetRegrid2",value);
+
+      //info
+      if(pt_verbose>3) printf("%c was set to %c",parname,value);
+
+      ierr = sprintf (parname, "position_y_%d", n+1);
+      assert (ierr>=0 && ierr<sizeof parname);
+      ierr = sprintf (value, "%f", pt_loc_y[n]);
+      assert (ierr>=0 && ierr<sizeof value);
+      ierr = CCTK_ParameterSet(parname,"CarpetRegrid2",value);
+
+      //info
+      if(pt_verbose>3) printf("%c was set to %c",parname,value);
+
+      ierr = sprintf (parname, "position_z_%d", n+1);
+      assert (ierr>=0 && ierr<sizeof parname);
+      ierr = sprintf (value, "%f", pt_loc_z[n]);
+      assert (ierr>=0 && ierr<sizeof value);
+      ierr = CCTK_ParameterSet(parname,"CarpetRegrid2",value);
+
+      //info
+      if(pt_verbose>3) printf("%c was set to %c",parname,value);        
+
+      //remember this location
+      pt_loc_old_x[n] = pt_loc_x[n];
+      pt_loc_old_y[n] = pt_loc_y[n];
+      pt_loc_old_z[n] = pt_loc_z[n];
+    }
+  }
+  else{
+    ierr = CCTK_ParameterSet ("regrid_every","CarpetReGrid2","-1");
+    if(pt_verbose>0) CCTK_INFO("not regridding now");
+    regrid_now = 0;           //just to make sure
+  }
+
+  if(pt_verbose>0)
+    CCTK_INFO ("done");
 }