diff options
author | schnetter <schnetter@a2659f00-0f4f-0410-9214-a4596bbb8a4f> | 2007-09-21 15:47:13 +0000 |
---|---|---|
committer | schnetter <schnetter@a2659f00-0f4f-0410-9214-a4596bbb8a4f> | 2007-09-21 15:47:13 +0000 |
commit | 83b6898036608da56fe70d8654019a38f9dfca81 (patch) | |
tree | b3ff961b41bd093527446478637d8397194afac7 /src | |
parent | f87829d2acde4a9913037ef93954905a0d85bf3f (diff) |
Rewrite puncture tracker
git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/PunctureTracker/trunk@6 a2659f00-0f4f-0410-9214-a4596bbb8a4f
Diffstat (limited to 'src')
-rw-r--r-- | src/adjust.c | 90 | ||||
-rw-r--r-- | src/make.code.defn | 5 | ||||
-rw-r--r-- | src/punc_track.c | 170 | ||||
-rw-r--r-- | src/puncture_tracker.c | 285 |
4 files changed, 287 insertions, 263 deletions
diff --git a/src/adjust.c b/src/adjust.c deleted file mode 100644 index 43a2df0..0000000 --- a/src/adjust.c +++ /dev/null @@ -1,90 +0,0 @@ -/* $Header$ */ - -#include <assert.h> -#include <stdio.h> - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - - - -void -BHTracker_Adjust (CCTK_ARGUMENTS); - -static void -setpar (char const * name, int index, char const * thorn, CCTK_REAL value); - - - -void -BHTracker_Adjust (CCTK_ARGUMENTS) -{ - DECLARE_CCTK_PARAMETERS; - DECLARE_CCTK_ARGUMENTS; - - int n; - int si, ri; - char name[100], val[100]; - int icnt, ierr; - - CCTK_REAL xpos,ypos,zpos; - - if ((cctk_iteration-1) < track_first) return; - if ((cctk_iteration-1) % track_every != 0) return; - - - - CCTK_INFO ("Tracking..."); - - for (n=0; n<num_tracked; ++n) { - - si = surface_index[n]; - assert (si>=0 && si<nsurfaces); - - ri = region_index[n]; - assert (ri>=0 && ri<num_offsets); - - if (sf_valid[si] > 0) { - - CCTK_VInfo (CCTK_THORNSTRING, - " Object #%d has moved to [%g,%g,%g]", n, - (double)sf_centroid_x[si], (double)sf_centroid_y[si], - (double)sf_centroid_z[si]); - - /* calculate new positions */ - xpos = (sf_centroid_x[si] - initial_x[si]); - ypos = (sf_centroid_y[si] - initial_y[si]); - zpos = (sf_centroid_z[si] - initial_z[si]); - - setpar ("offsetx", ri, "CarpetRegrid", xpos); - setpar ("offsety", ri, "CarpetRegrid", ypos); - setpar ("offsetz", ri, "CarpetRegrid", zpos); - - } else { - - CCTK_VInfo (CCTK_THORNSTRING, - " Object #%d cannot be detected at this time", n); - - } - } -} - - - -static void -setpar (char const * const name, int const index, char const * const thorn, - CCTK_REAL const value) -{ - char nambuf[100], valbuf[100]; - int icnt, ierr; - - icnt = snprintf (nambuf, sizeof nambuf, "%s[%d]", name, index); - assert (icnt>=0 && icnt<sizeof nambuf); - - icnt = snprintf (valbuf, sizeof valbuf, "%g", value); - assert (icnt>=0 && icnt<sizeof valbuf); - - ierr = CCTK_ParameterSet (nambuf, thorn, valbuf); - assert (! ierr); -} diff --git a/src/make.code.defn b/src/make.code.defn index e0e5f26..31baa1c 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -1,8 +1,7 @@ -# Main make.code.defn file for thorn PunctureTracker -# $Header$ +# Main make.code.defn file for thorn PunctureTracker -*-Makefile-*- # Source files in this directory -SRCS = punc_track.c +SRCS = puncture_tracker.c # Subdirectories containing source files SUBDIRS = 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; - } -} diff --git a/src/puncture_tracker.c b/src/puncture_tracker.c new file mode 100644 index 0000000..d3da795 --- /dev/null +++ b/src/puncture_tracker.c @@ -0,0 +1,285 @@ +#include <assert.h> +#include <stdio.h> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_Table.h" + + + +static int const max_num_tracked = 10; + + + +void +PunctureTracker_Init (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + CCTK_INFO ("Initializing PunctureTracker"); + } + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + pt_loc_t[n] = cctk_time; + pt_loc_x[n] = initial_x[n]; + pt_loc_y[n] = initial_y[n]; + pt_loc_z[n] = initial_z[n]; + } else { + // Initialise to some sensible but unimportant values + pt_loc_t[n] = 0.0; + pt_loc_x[n] = 0.0; + pt_loc_y[n] = 0.0; + pt_loc_z[n] = 0.0; + } + pt_loc_t_p[n] = 0.0; + pt_loc_x_p[n] = 0.0; + pt_loc_y_p[n] = 0.0; + pt_loc_z_p[n] = 0.0; + } +} + + + +void +PunctureTracker_Track (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + // Do not track while setting up initial data; + // time interpolation may fail + + if (cctk_iteration == 0) { + return; + } + + // Some output + + if (verbose) { + CCTK_INFO ("Tracking punctures..."); + } + + if (verbose) { + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + CCTK_VInfo (CCTK_THORNSTRING, + "Puncture #%d is at (%g,%g,%g)", + n, + (double)pt_loc_x[n], + (double)pt_loc_y[n], + (double)pt_loc_z[n]); + } + } + } + + // Manual time level cycling + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + pt_loc_t_p[n] = pt_loc_t[n]; + pt_loc_x_p[n] = pt_loc_x[n]; + pt_loc_y_p[n] = pt_loc_y[n]; + pt_loc_z_p[n] = pt_loc_z[n]; + + pt_loc_t[n] = cctk_time; + } + } + + // Interpolate + + // Dimensions + int const dim = 3; + + // Interpolation operator + int const operator_handle = + CCTK_InterpHandle ("Lagrange polynomial interpolation"); + if (operator_handle < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can’t get interpolation handle"); + goto label_return; + } + + // Interpolation parameter table + int const param_table_handle = Util_TableCreateFromString ("order=4"); + if (param_table_handle < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can't create parameter table"); + goto label_return; + } + + { + + // Interpolation coordinate system + int const coordsys_handle = CCTK_CoordSystemHandle ("cart3d"); + if (coordsys_handle < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can't get coordinate system handle"); + goto label_free_param_table; + } + + // Only processor 0 interpolates + int const num_points = CCTK_MyProc(cctkGH) == 0 ? max_num_tracked : 0; + + // Interpolation coordinates + assert (dim == 3); + CCTK_POINTER_TO_CONST interp_coords[3]; + interp_coords[0] = pt_loc_x_p; + interp_coords[1] = pt_loc_y_p; + interp_coords[2] = pt_loc_z_p; + + // Number of interpolation variables + int const num_vars = 3; + + // Interpolated variables + assert (num_vars == 3); + int input_array_indices[3]; + input_array_indices[0] = CCTK_VarIndex ("ADMBase::betax"); + input_array_indices[1] = CCTK_VarIndex ("ADMBase::betay"); + input_array_indices[2] = CCTK_VarIndex ("ADMBase::betaz"); + + // Interpolation result types + assert (num_vars == 3); + CCTK_INT output_array_type_codes[3]; + output_array_type_codes[0] = CCTK_VARIABLE_REAL; + output_array_type_codes[1] = CCTK_VARIABLE_REAL; + output_array_type_codes[2] = CCTK_VARIABLE_REAL; + + // Interpolation result + CCTK_REAL pt_betax[max_num_tracked]; + CCTK_REAL pt_betay[max_num_tracked]; + CCTK_REAL pt_betaz[max_num_tracked]; + + assert (num_vars == 3); + CCTK_POINTER output_arrays[3]; + output_arrays[0] = pt_betax; + output_arrays[1] = pt_betay; + output_arrays[2] = pt_betaz; + + // Interpolate + int const ierr = CCTK_InterpGridArrays + (cctkGH, dim, + operator_handle, param_table_handle, coordsys_handle, + num_points, + CCTK_VARIABLE_REAL, + interp_coords, + num_vars, input_array_indices, + num_vars, output_array_type_codes, output_arrays); + if (ierr < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Interpolation error"); + goto label_free_param_table; + } + + if (CCTK_MyProc(cctkGH) == 0) { + + // Some more output + + if (verbose && CCTK_MyProc(cctkGH) == 0) { + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + CCTK_VInfo (CCTK_THORNSTRING, + "Shift at puncture #%d is at (%g,%g,%g)", + n, + (double)pt_betax[n], + (double)pt_betay[n], + (double)pt_betaz[n]); + } + } + } + + // Time evolution + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + CCTK_REAL const dt = pt_loc_t[n] - pt_loc_t_p[n]; + // First order time integrator + // Michael Koppitz says this works... + // if it doesn't, we can make it second order accurate + pt_loc_x[n] = pt_loc_x_p[n] + dt * (- pt_betax[n]); + pt_loc_y[n] = pt_loc_y_p[n] + dt * (- pt_betay[n]); + pt_loc_z[n] = pt_loc_z_p[n] + dt * (- pt_betaz[n]); + } + } + + } + + // Broadcast result + + CCTK_REAL loc_local[3*max_num_tracked]; + if (CCTK_MyProc(cctkGH) == 0) { + for (int n = 0; n < max_num_tracked; ++ n) { + loc_local[ n] = pt_loc_x[n]; + loc_local[ max_num_tracked+n] = pt_loc_y[n]; + loc_local[2*max_num_tracked+n] = pt_loc_z[n]; + } + } else { + for (int n = 0; n < max_num_tracked; ++ n) { + loc_local[ n] = 0.0; + loc_local[ max_num_tracked+n] = 0.0; + loc_local[2*max_num_tracked+n] = 0.0; + } + } + + CCTK_REAL loc_global[3*max_num_tracked]; + + int const handle_sum = CCTK_ReductionHandle ("sum"); + if (handle_sum < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can't get redunction handle"); + goto label_free_param_table; + } + + int const ierr2 = CCTK_ReduceLocArrayToArray1D + (cctkGH, -1, handle_sum, + loc_local, loc_global, 3*max_num_tracked, CCTK_VARIABLE_REAL); + if (ierr2 < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Reduction error"); + goto label_free_param_table; + } + + for (int n = 0; n < max_num_tracked; ++ n) { + pt_loc_x[n] = loc_global[ n]; + pt_loc_y[n] = loc_global[ max_num_tracked+n]; + pt_loc_z[n] = loc_global[2*max_num_tracked+n]; + } + + } + + // Done + + // Poor man's exception handling + label_free_param_table: + Util_TableDestroy (param_table_handle); + + label_return: + return; +} + + + +void +PunctureTracker_SetPositions (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Setting position of refined region from puncture #%d to (%g,%g,%g)", + n, + (double)pt_loc_x[n], + (double)pt_loc_y[n], + (double)pt_loc_z[n]); + } + + // Set position in CarpetRegrid2 + position_x[n] = pt_loc_x[n]; + position_y[n] = pt_loc_y[n]; + position_z[n] = pt_loc_z[n]; + + } + } +} |