diff options
Diffstat (limited to 'src/patch')
| -rw-r--r-- | src/patch/patch_interp.cc | 113 | ||||
| -rw-r--r-- | src/patch/patch_interp.hh | 20 |
2 files changed, 117 insertions, 16 deletions
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index c9e170b..e21b32f 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -12,8 +12,14 @@ #include <assert.h> #include <math.h> +// Cactus stuff +#include "util_Table.h" +#include "cctk.h" + +// jtutil:: stuff #include "jt/stdc.h" #include "jt/util.hh" +using jtutil::error_exit; #include "jt/array.hh" #include "jt/cpm_map.hh" #include "jt/linear_map.hh" @@ -72,12 +78,23 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in, interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)), gridfn_coord_origin_(my_edge().par_of_ipar(min_ipar())), gridfn_coord_delta_(my_edge().par_map().delta_fp()), - gridfn_data_ptrs_(ghosted_min_gfn_in, ghosted_max_gfn_in, - min_iperp_in, max_iperp_in), - gridfn_offsets_(ghosted_min_gfn_in, ghosted_max_gfn_in, - min_iperp_in, max_iperp_in), - gridfn_par_stride_(0) // computed later + gridfn_data_ptrs_(min_iperp_in, max_iperp_in, + ghosted_min_gfn_in, ghosted_max_gfn_in), + gridfn_type_codes_(ghosted_min_gfn_in, ghosted_max_gfn_in) { +int status; + +// did interpolator-parameter-table cloning work ok? +status = interp_par_table_handle_; +if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_frontier::patch_frontier():\n" +" can't clone interpolator parmameter table (handle=%d)!\n" +" error status=%d\n" +, + interp_par_table_handle_in, + status); /*NOTREACHED*/ + // number of points to interpolate at each iperp for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) { @@ -85,6 +102,42 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in, = jtutil::how_many_in_range(min_parindex_used_in(iperp), max_parindex_used_in(iperp)); } + +// set up pointers to start of gridfn data we can use for each iperp + for (int gfn = ghosted_min_gfn() ; gfn <= ghosted_max_gfn() ; ++gfn) + { + for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) + { + const int irho = my_edge(). irho_of_iperp_ipar(iperp, min_ipar()); + const int isigma = my_edge().isigma_of_iperp_ipar(iperp, min_ipar()); + gridfn_data_ptrs_(iperp,gfn) + = & my_patch().ghosted_gridfn(gfn, irho,isigma); + } + } + +// set up CCTK_VARIABLE_* type codes for gridfns + for (int gfn = ghosted_min_gfn() ; gfn <= ghosted_max_gfn() ; ++gfn) + { + gridfn_type_codes_(gfn) = CCTK_VARIABLE_REAL; // fp == CCTK_REAL + } + +// +// set up parameter table for interpolator +// +const int N_dims = 1; +const CCTK_INT stride = my_patch().ghosted_iang_stride(my_edge().is_rho()); +status = Util_TableSetIntArray(interp_par_table_handle_, + N_dims, &stride, + "input_array_strides"); +if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_frontier::patch_frontier():\n" +" can't set stride in interpolator parmameter table (handle=%d)!\n" +" error status=%d\n" +, + interp_par_table_handle_, + status); /*NOTREACHED*/ + } //***************************************************************************** @@ -106,12 +159,50 @@ Util_TableDestroy(interp_par_table_handle_); // It calls error_exit() if the point is out of range or any other // interpolator error is detected. // -fp patch_frontier::interpolate(int ghosted_min_gfn, ghosted_max_gfn, - jtutil::array3d<fp>& result_buffer) const +void patch_frontier::interpolate(int ghosted_min_gfn, ghosted_max_gfn, + jtutil::array3d<fp>& result_buffer) const { const int N_dims = 1; +const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp, + ghosted_max_gfn_to_interp); -const int status - = CCTK_InterpLocalUniform(N_dims, - interp_handle_, - interp_par_table_handle_, + for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) + { + const int min_parindex = min_parindex_used_(iperp); + const CCTK_INT N_interp_points = N_interp_points_(iperp); + const void* const interp_coords[N_dims] + = { static_cast<const void *>( + & interp_par_(iperp,min_parindex) + ) }; + const void* const* input_arrays + = static_cast<const void *const *>( + & gridfn_data_ptrs(iperp,ghosted_min_gfn_to_interp) + ); + // FIXME FIXME input_arrays and output_arrays are dubious + const CCTK_INT *gridfn_type_codes_ptr + = & gridfn_type_codes_(ghosted_min_gfn_to_interp); + const int status + = CCTK_InterpLocalUniform(N_dims, + interp_handle_, + interp_par_table_handle_, + gridfn_coord_origin_, + gridfn_coord_delta_, + N_interp_points, + CCTK_VARIABLE_REAL, + interp_coords, + N_gridfns, + &N_interp_points, + gridfn_type_codes_ptr, + input_arrays, + N_gridfns, + gridfn_type_codes_ptr, + output_arrays); + if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_frontier::interpolate():\n" +" error return %d from interpolator at iperp=%d of [%d,%d]!\n" +, + status, iperp, min_iperp(), max_iperp()); + /*NOTREACHED*/ + } +} diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index ba62792..2c92da2 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -152,6 +152,7 @@ private: // min parindex at each iperp // ... this is a reference to an array passed in to our constructor // ==> we do *not* own this! + // ... index is (iperp) const jtutil::array1d<int>& min_parindex_used_; // number of points to interpolate at each iperp @@ -187,6 +188,8 @@ private: fp gridfn_coord_origin_, gridfn_coord_delta_; // + // We use (the 0-origin) parsub as an interpolator point index + // within each interpolation (i.e. at each value of iperp). // The interpolator accesses the gridfn data with the generic // subscripting expression // data[offset + parsub*stride] @@ -196,13 +199,20 @@ private: // Our strategy will be to leave offset unspecified so it // defaults to 0, specify a common stride for all (gfn,iperp), // and specify the appropriate data pointer for each (gfn,iperp). + // Thus data should point to the start of the gridfn data + // we should use for each iperp // - // --> start of gridfn data we can use for each iperp - // subscripts are (gfn,iperp) + // --> start of gridfn data we should use for each iperp + // subscripts are (iperp,gfn) so [gfn] is contiguous at each iperp + // ... n.b. we do *not* own this data! jtutil::array2d<fp *> gridfn_data_ptrs_; - // stride of par coordinate in gridfn data - // (same for all gfn and iperp) - int gridfn_par_stride_; + // CCTK_VARIABLE_ type codes for the gridfns + // ... all are actually CCTK_REAL, but we need the array for + // the interpolator + // ... datatype is CCTK_INT so we can pass this by reference + // to the interpolator + // ... index is (ghosted_gfn) + jtutil::array1d<CCTK_INT> gridfn_type_codes_; }; |