aboutsummaryrefslogtreecommitdiff
path: root/src/punc_track.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/punc_track.c')
-rw-r--r--src/punc_track.c170
1 files changed, 0 insertions, 170 deletions
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;
- }
-}