diff options
Diffstat (limited to 'src/punc_track.c')
-rw-r--r-- | src/punc_track.c | 170 |
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; - } -} |