add sf_active once surface has new location

return to use of spherical surface


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

1 files changed, 45 insertions, 128 deletions

diff --git a/src/punc_track.c b/src/punc_track.c
index 69cb754..2b9bd26 100644
--- a/src/punc_track.c
+++ b/src/punc_track.c
@@ -7,14 +7,8 @@
 #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];
+#define REFLEVEL round(log10((CCTK_REAL)(cctkGH->cctk_levfac[0]))/log10(2.0))
 
 void PunctureTracker (CCTK_ARGUMENTS);
 void PunctureTracker_init (CCTK_ARGUMENTS);
@@ -43,23 +37,10 @@ void PunctureTracker_init (CCTK_ARGUMENTS)
     pt_loc_z_p[0] = 0.0;
     pt_loc_z_p[1] = 0.0;
   }
-  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");
+  if(pt_verbose>0)
+    CCTK_INFO("done initializing PunctureTracker");
 }
 
-
 void PunctureTracker (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_PARAMETERS;
@@ -70,9 +51,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 pt_betax[10];
-  CCTK_REAL pt_betay[10];
-  CCTK_REAL pt_betaz[10];
+  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];
@@ -102,7 +83,12 @@ 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;
@@ -130,15 +116,16 @@ void PunctureTracker (CCTK_ARGUMENTS)
   
 
   for (n=0; n< pt_num_tracked; 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];
+    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){
@@ -146,108 +133,38 @@ void update_punc_loc(CCTK_ARGUMENTS){
   DECLARE_CCTK_PARAMETERS;
   DECLARE_CCTK_ARGUMENTS;
 
-  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];
-  
+  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;    
-  
-    //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;
-
-    //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_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("     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]);
+      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;
 
-    //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;
-        }
-      }
-    }
+    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;
   }
-  
-  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");
 }