From f87829d2acde4a9913037ef93954905a0d85bf3f Mon Sep 17 00:00:00 2001 From: koppitz Date: Sun, 11 Mar 2007 19:17:23 +0000 Subject: 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 --- src/punc_track.c | 173 +++++++++++++++---------------------------------------- 1 file 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 -#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; i4) - 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=0 && ierr3) printf("%c was set to %c",parname,value); - - ierr = sprintf (parname, "position_z_%d", n+1); - assert (ierr>=0 && ierr=0 && ierr3) 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"); } -- cgit v1.2.3