From caa53b52471955e21d1fbfd0ae1aca6872316491 Mon Sep 17 00:00:00 2001 From: jthorn Date: Wed, 26 Jun 2002 10:43:13 +0000 Subject: finish (?) Jacobian computation git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@596 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/patch/patch_interp.cc | 214 +++++++++++++++++++--------------------------- src/patch/patch_interp.hh | 37 ++++---- 2 files changed, 110 insertions(+), 141 deletions(-) (limited to 'src') diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index 3e08add..0784f74 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -6,6 +6,7 @@ // patch_interp::~patch_interp // patch_interp::compute_[min,max]_ipar // patch_interp::interpolate +// // patch_interp::query_interpolator // patch_interp::verify_Jacobian_sparsity_pattern_ok // patch_interp::molecule_minmax_m_ipar @@ -36,6 +37,8 @@ using jtutil::error_exit; #include "patch_interp.hh" #include "ghost_zone.hh" +//***************************************************************************** +//***************************************************************************** //***************************************************************************** // @@ -50,8 +53,8 @@ using jtutil::error_exit; // patch_interp::patch_interp(const patch_edge& my_edge_in, int min_iperp_in, int max_iperp_in, - const jtutil::array1d& min_parindex_used_in, - const jtutil::array1d& max_parindex_used_in, + const jtutil::array1d& min_parindex_array_in, + const jtutil::array1d& max_parindex_array_in, const jtutil::array2d& interp_par_in, int interp_handle_in, int interp_par_table_handle_in) : my_patch_(my_edge_in.my_patch()), @@ -61,8 +64,8 @@ patch_interp::patch_interp(const patch_edge& my_edge_in, min_iperp_(min_iperp_in), max_iperp_(max_iperp_in), min_ipar_(compute_min_ipar()), max_ipar_(compute_max_ipar()), - min_parindex_used_(min_parindex_used_in), - max_parindex_used_(max_parindex_used_in), + min_parindex_array_(min_parindex_array_in), + max_parindex_array_(max_parindex_array_in), interp_par_(interp_par_in), interp_handle_(interp_handle_in), interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)), @@ -161,7 +164,7 @@ return max_par_adjacent_ghost_zone.is_symmetry() // This function interpolates the specified range of ghosted gridfns // at all the coordinates specified when we were constructed. (It does // a separate interpolator call for each iperp.) This function stores -// the interpolated data in interp_data_buffer(gfn, iperp,parindex) . +// the interpolated data in data_buffer(gfn, iperp,parindex) . // void patch_interp::interpolate(int ghosted_min_gfn_to_interp, int ghosted_max_gfn_to_interp, @@ -175,7 +178,7 @@ const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp, ghosted_max_gfn_to_interp); const CCTK_INT N_gridfn_data_points = jtutil::how_many_in_range(min_ipar(), max_ipar()); -const int interp_coords_type_code = CCTK_VARIABLE_REAL; // fp == CCTK_REAL +const int interp_coords_type_code = CCTK_VARIABLE_REAL; const CCTK_INT *gridfn_type_codes_ptr = & gridfn_type_codes_(ghosted_min_gfn_to_interp); @@ -184,12 +187,11 @@ const CCTK_INT *gridfn_type_codes_ptr // for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) { - const int min_parindex = min_parindex_used_(iperp); - const int max_parindex = max_parindex_used_(iperp); - // // interpolation-point coordinates // + const int min_parindex = min_parindex_array_(iperp); + const int max_parindex = max_parindex_array_(iperp); const CCTK_INT N_interp_points = jtutil::how_many_in_range(min_parindex, max_parindex); const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex); @@ -252,68 +254,89 @@ const CCTK_INT *gridfn_type_codes_ptr } } +//****************************************************************************** +//****************************************************************************** //****************************************************************************** // // This function queries the interpolator with whatever arguments are -// in the parameter table. It passes the specified interpolation-point -// arguments to the interpolator. It passes no-op values (a single CCTK_REAL -// input and output array, but NULL input and output array pointers) for -// the input/output array arguments to the interpolator. The result is -// that the interpolator should do only the desired query, not any actual -// interpolation operation. +// in the parameter table. It specifies arguments such that no actual +// interpolation is done. +// +// In particular, the following interpolator arguments are set up properly +// based on iperp : +// N_dims +// operator_handle, param_table_handle, +// coord_origin, coord_delta, +// N_interp_points, interp_coords_type_code, interp_coords, +// input_array_dims # specifies the correct par size of the +// # patch interpolation region at this iperp +// The following arguments are set to specify a single input array, but +// with a NULL data pointer so no actual data is used: +// N_input_arrays +// input_array_type_codes, input_arrays +// The following arguments are set to specify a single output array, but +// with a NULL data pointer so no actual data is stored: +// N_output_arrays +// output_array_type_codes, output_arrays +// +// Arguments: +// iperp = (in) Specifies where in the patch interpolation region the +// interpolator query should be done. // // Results: -// This function returns the interpolator status code. +// If the interpolator returns an error status, this function does an +// error_exit() (and does not return to its caller). Otherwise, this +// function returns the interpolator status code. // -int patch_interp::query_interpolator(int N_interp_points = 0, - const void *const *interp_coords = NULL) +int patch_interp::query_interpolator(const char function_name[], int iperp) const { -#ifdef NOT_YET -const bool interp_points_flag = (iperp != INT_SENTINAL_VALUE); - -xxxxx +int status; const int N_dims = 1; -const int N_interp_points = interp_points_flag - ? jtutil::how_many_in_range(min_parindex, - max_parindex) - : 0; -xxxx + +const int min_parindex = min_parindex_array_(iperp); +const int max_parindex = max_parindex_array_(iperp); +const int N_interp_points = jtutil::how_many_in_range(min_parindex, + max_parindex); const int interp_coords_type_code = CCTK_VARIABLE_REAL; +const void* const interp_coords[N_dims] + = { static_cast(& interp_par_(iperp, min_parindex)) }; + const int N_input_arrays = 1; -const CCTK_INT input_array_dims[N_dims] = { 0 }; +const CCTK_INT input_array_dims[N_dims] + = { jtutil::how_many_in_range(min_ipar(), max_ipar()) }; const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL }; const void* const input_arrays[N_input_arrays] = { NULL }; + const int N_output_arrays = 1; const CCTK_INT output_array_type_codes[N_output_arrays] = { CCTK_VARIABLE_REAL }; void* const output_arrays[N_output_arrays] = { NULL }; +status = CCTK_InterpLocalUniform(N_dims, + interp_handle_, interp_par_table_handle_, + &gridfn_coord_origin_, &gridfn_coord_delta_, + N_interp_points, + interp_coords_type_code, + interp_coords, + N_input_arrays, + input_array_dims, + input_array_type_codes, + input_arrays, + N_output_arrays, + output_array_type_codes, + output_arrays); +if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_interp::query_interpolator():\n" +" on behalf of patch_interp::%s()\n" +" error return %d from interpolator query at iperp=%d of [%d,%d]!\n" +, + function_name, + status, iperp, min_iperp(), max_iperp()); /*NOTREACHED*/ - // interpolation-point coordinates - const CCTK_INT N_interp_points - = jtutil::how_many_in_range(min_parindex, max_parindex); - const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex); - const void* const interp_coords[N_dims] - = { static_cast(interp_coords_ptr) }; - -return CCTK_InterpLocalUniform(N_dims, - interp_handle_, interp_par_table_handle_, - &gridfn_coord_origin_, &gridfn_coord_delta_, - N_interp_points, - interp_coords_type_code, - interp_coords, - N_input_arrays, - input_array_dims, - input_array_type_codes, - input_arrays, - N_output_arrays, - output_array_type_codes, - output_arrays); -#else -return 0; -#endif +return status; } //****************************************************************************** @@ -335,7 +358,7 @@ return 0; void patch_interp::verify_Jacobian_sparsity_pattern_ok() const { -int status, status1, status2, status3; +int status1, status2, status3; // // query the interpolator @@ -343,13 +366,7 @@ int status, status1, status2, status3; // Jacobian metainfo is returned by the interpolator on every call, // even without any explicit request // -status = query_interpolator(); -if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n" -" error return %d from interpolator query!\n" -, - status); /*NOTREACHED*/ +query_interpolator("verify_Jacobian_sparsity_pattern_ok"); // // get Jacobian-sparsity query results from parameter table @@ -405,7 +422,8 @@ if ( MSS_is_fn_of_interp_coords || MSS_is_fn_of_input_array_values void patch_interp::molecule_minmax_m_ipar(int& min_m_ipar, int& max_m_ipar) const { -int status, status1, status2; +const int N_dims = 1; +int status1, status2; // // set molecule min/max m query entries in parameter table @@ -426,13 +444,7 @@ if ((status1 < 0) || (status2 < 0)) // // query the interpolator // -status = query_interpolator(); -if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::molecule_minmax_m_ipar():\n" -" error return %d from interpolator query!\n" -, - status); /*NOTREACHED*/ +query_interpolator("molecule_minmax_m_ipar"); // // get molecule min/max m query results from parameter table @@ -492,11 +504,10 @@ int status; for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) { - const int min_parindex = min_parindex_used_(iperp); - const int max_parindex = max_parindex_used_(iperp); + const int min_parindex = min_parindex_array_(iperp); // set up the molecule-position query in the parameter table - CCCTK_POINTER molecule_posn_ptrs[N_dims] + CCTK_POINTER molecule_posn_ptrs[N_dims] = { static_cast(& posn_buffer(iperp, min_parindex)) }; status = Util_TableSetPointerArray(interp_par_table_handle_, N_dims, molecule_posn_ptrs, @@ -511,24 +522,10 @@ int status; iperp, min_iperp(), max_iperp(), status); /*NOTREACHED*/ - // interpolation-point coordinates - const CCTK_INT N_interp_points - = jtutil::how_many_in_range(min_parindex, max_parindex); - const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex); - const void* const interp_coords[N_dims] - = { static_cast(interp_coords_ptr) }; - // query the interpolator to get the molecule positions // for all parindex at this iperp (the interpolator will - // store them directly in the posn_buffer(iperp,*) array - status = query_interpolator(N_interp_points, interp_coords); - if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::molecule_posn():\n" -" error return %d from interpolator query at iperp=%d of [%d,%d]!\n" -, - status, iperp, min_iperp(), max_iperp()); - /*NOTREACHED*/ + // store them directly in the posn_buffer(iperp,*) array) + query_interpolator("molecule_positions", iperp); } // @@ -557,9 +554,9 @@ if (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0) // This function queries the interpolator at each iperp to get the // Jacobian of the interpolated data with respect to this patch's // ghosted gridfns, -// partial interp_data_buffer(ghosted_gfn, iperp, parindex) -// -------------------------------------------------------- -// partial ghosted_gridfn(ghosted_gfn, iperp, posn+m_ipar) +// partial interpolate() data_buffer(ghosted_gfn, iperp, parindex) +// --------------------------------------------------------------- +// partial ghosted_gridfn(ghosted_gfn, iperp, posn+m_ipar) // // This function stores the Jacobian in // Jacobian_buffer(iperp, parindex, m_ipar) @@ -582,10 +579,6 @@ void patch_interp::Jacobian(jtutil::array3d& Jacobian_buffer) { const int N_dims = 1; const int N_gridfns = 1; -const CCTK_INT N_gridfn_data_points - = jtutil::how_many_in_range(min_ipar(), max_ipar()); -const CCTK_INT *gridfn_type_codes_ptr - = & gridfn_type_codes_(ghosted_min_gfn_to_interp); int status, status1, status2; @@ -594,12 +587,12 @@ int status, status1, status2; // set Jacobian stride info in parameter table // const int Jacobian_interp_point_stride - = Jacobian_buffer->subscript_stride_j(); + = Jacobian_buffer.subscript_stride_j(); status1 = Util_TableSetInt(interp_par_table_handle_, Jacobian_interp_point_stride, "Jacobian_interp_point_stride"); const CCTK_INT Jacobian_m_strides[N_dims] - = { Jacobian_buffer->subscript_stride_k() }; // = { 1 } in practice + = { Jacobian_buffer.subscript_stride_k() }; // = { 1 } in practice status2 = Util_TableSetIntArray(interp_par_table_handle_, N_dims, Jacobian_m_strides, "Jacobian_m_strides"); @@ -617,17 +610,7 @@ if ((status1 < 0) || (status2 < 0)) // for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) { - const int min_parindex = min_parindex_used_(iperp); - const int max_parindex = max_parindex_used_(iperp); - - // - // interpolation-point coordinates - // - const CCTK_INT N_interp_points - = jtutil::how_many_in_range(min_parindex, max_parindex); - const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex); - const void* const interp_coords[N_dims] - = { static_cast(interp_coords_ptr) }; + const int min_parindex = min_parindex_array_(iperp); // // set up the Jacobian query in the parameter table @@ -651,24 +634,7 @@ if ((status1 < 0) || (status2 < 0)) // // query the interpolator to get the Jacobian for this iperp // - const CCTK_INT output_array_type_codes[N_gridfns] - = { CCTK_VARIABLE_REAL }; - CCTK_POINTER const output_arrays[N_gridfns] = { NULL }; -#ifdef NOT_YET - status = query_interpolator(N_interp_points, interp_coords, - N_gridfns, - output_array_type_codes, - output_arrays); -#else - status = query_interpolator(N_interp_points, interp_coords); -#endif - if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::Jacobian():\n" -" error return %d from interpolator query at iperp=%d of [%d,%d]!\n" -, - status, iperp, min_iperp(), max_iperp()); - /*NOTREACHED*/ + query_interpolator("Jacobian", iperp); } // diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index 4bd402b..92bc127 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -67,8 +67,8 @@ // region would be described by // [min,max]_iperp_=[10,11] // [min,max]_ipar_=[1,9] -// [min,max]_parindex_used_(10)=[2,5] -// [min,max]_parindex_used_(11)=[2,6] +// [min,max]_parindex_array_(10)=[2,5] +// [min,max]_parindex_array_(11)=[2,6] // interp_par_(10,2) = x[a] // interp_par_(10,3) = x[b] // interp_par_(10,4) = x[c] @@ -120,7 +120,7 @@ public: // interpolate specified range of ghosted gridfns // at all the coordinates specified when we were constructed, // store interpolated data in - // interp_data_buffer(ghosted_gfn, iperp, parindex) + // data_buffer(ghosted_gfn, iperp, parindex) void interpolate(int ghosted_min_gfn_to_interp, int ghosted_max_gfn_to_interp, jtutil::array3d& data_buffer) @@ -150,9 +150,9 @@ public: // get Jacobian of interpolated data with respect to this patch's // ghosted gridfns, - // partial interp_data_buffer(ghosted_gfn, iperp, parindex) - // -------------------------------------------------------- - // partial ghosted_gridfn(ghosted_gfn, iperp, posn+m_ipar) + // partial interpolate() data_buffer(ghosted_gfn, iperp, parindex) + // --------------------------------------------------------------- + // partial ghosted_gridfn(ghosted_gfn, iperp, posn+m_ipar) // store Jacobian in // Jacobian_buffer(iperp, parindex, m_ipar) // where we implicitly assume the Jacobian to be independent of @@ -177,17 +177,20 @@ private: int compute_max_ipar() const; // min/max (iperp,ipar) of the gridfn data we'll use - // (uses precomputed values) + // (uses cached values set up by our constructor) int min_ipar() const { return min_ipar_; } int max_ipar() const { return max_ipar_; } // query the interpolator via the parameter table - // passing the specified arguments to the interpolator, - // and no-op values for all the other interpolator arguments - // ... return interpolator status code - int query_interpolator(int N_interp_points = 0, - const void* const *interp_coords = NULL) + // on behalf of the specified function (name used for error msgs only) + // ... error_exit() if error return from interpolator, + // otherwise return interpolator status code + int query_interpolator(const char function_name[], int iperp) const; + // ... use default iperp for queries that don't care + int query_interpolator(const char function_name[]) + const + { return query_interpolator(function_name, min_iperp_); } // @@ -200,7 +203,7 @@ public: // interpolation region is to belong. // [min,max]_iperp_in = The range of iperp for this interpolation // region - // [min,max]_parindex_used_in(iperp) + // [min,max]_parindex_array_in(iperp) // = [min,max] range of parindex actually used at each iperp. // We keep references to these arrays, so they should have // lifetimes at last as long as that of this object. @@ -234,8 +237,8 @@ public: // patch_interp(const patch_edge& my_edge_in, int min_iperp_in, int max_iperp_in, - const jtutil::array1d& min_parindex_used_in, - const jtutil::array1d& max_parindex_used_in, + const jtutil::array1d& min_parindex_array_in, + const jtutil::array1d& max_parindex_array_in, const jtutil::array2d& interp_par_in, int interp_handle_in, int interp_par_table_handle_in); ~patch_interp(); @@ -269,8 +272,8 @@ private: // ... these are references to arrays passed in to our constructor // ==> we do *not* own them! // ... indices are (iperp) - const jtutil::array1d& min_parindex_used_; - const jtutil::array1d& max_parindex_used_; + const jtutil::array1d& min_parindex_array_; + const jtutil::array1d& max_parindex_array_; // interp_par_(iperp,parindex) // = Gives the par coordinates at which we will interpolate; -- cgit v1.2.3