From 0f9bd2adc8e1a09fc5f549fc83c9b96724e56b0e Mon Sep 17 00:00:00 2001 From: koppitz Date: Tue, 6 Mar 2007 13:04:52 +0000 Subject: some bugfixes git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/PunctureTracker/trunk@4 a2659f00-0f4f-0410-9214-a4596bbb8a4f --- src/punc_track.c | 167 ++++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 128 insertions(+), 39 deletions(-) (limited to 'src') 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 #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; i1) - 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=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