diff options
Diffstat (limited to 'src/patch')
| -rw-r--r-- | src/patch/coords.cc | 11 | ||||
| -rw-r--r-- | src/patch/fd_grid.cc | 10 | ||||
| -rw-r--r-- | src/patch/ghost_zone.cc | 15 | ||||
| -rw-r--r-- | src/patch/ghost_zone.hh | 8 | ||||
| -rw-r--r-- | src/patch/grid.cc | 10 | ||||
| -rw-r--r-- | src/patch/make.code.defn | 2 | ||||
| -rw-r--r-- | src/patch/patch.cc | 12 | ||||
| -rw-r--r-- | src/patch/patch.hh | 15 | ||||
| -rw-r--r-- | src/patch/patch_interp.cc | 548 | ||||
| -rw-r--r-- | src/patch/patch_interp.hh | 75 | ||||
| -rw-r--r-- | src/patch/patch_system.cc | 12 | ||||
| -rw-r--r-- | src/patch/test_coords.cc | 6 | ||||
| -rw-r--r-- | src/patch/test_coords2.cc | 6 | ||||
| -rw-r--r-- | src/patch/test_patch_system.cc | 12 |
14 files changed, 386 insertions, 356 deletions
diff --git a/src/patch/coords.cc b/src/patch/coords.cc index af32552..b529a1c 100644 --- a/src/patch/coords.cc +++ b/src/patch/coords.cc @@ -23,18 +23,17 @@ #include <assert.h> #include <limits.h> -#include "jt/stdc.h" -#include "jt/util.hh" - -#include "../config.hh" -#include "coords.hh" - +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" using jtutil::error_exit; using jtutil::arctan_xy; using jtutil::signum; using jtutil::pow2; using jtutil::hypot3; +#include "coords.hh" + //****************************************************************************** //****************************************************************************** //****************************************************************************** diff --git a/src/patch/fd_grid.cc b/src/patch/fd_grid.cc index 10ac5b5..d0ca181 100644 --- a/src/patch/fd_grid.cc +++ b/src/patch/fd_grid.cc @@ -9,13 +9,13 @@ #include <assert.h> #include <math.h> -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../config.hh" #include "coords.hh" #include "grid.hh" #include "fd_grid.hh" diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc index d686460..ecb397e 100644 --- a/src/patch/ghost_zone.cc +++ b/src/patch/ghost_zone.cc @@ -24,14 +24,15 @@ #include "cctk.h" -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/cpm_map.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" + +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../config.hh" #include "coords.hh" #include "grid.hh" #include "fd_grid.hh" @@ -390,6 +391,8 @@ const int other_min_iperp = std::min(other_iperp(min_iperp()), other_iperp(max_iperp())); const int other_max_iperp = std::max(other_iperp(min_iperp()), other_iperp(max_iperp())); + + // // compute extreme range of ipar that we'll use, // and set up arrays giving actual [min,max] ipar that we'll use diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh index 53c4bd4..1728095 100644 --- a/src/patch/ghost_zone.hh +++ b/src/patch/ghost_zone.hh @@ -37,10 +37,10 @@ // ghost zone or zones) based on either the patch system's symmetry // or interpolation from a neighboring patch. ghost_zone is an abstract // base class, from which we derive two classes: -// * A symmetry_ghost_zone object describes a ghost zone which -// is a (discrete) symmetry of spacetime, either mirror-image or -// periodic. Such an object knows how to fill in ghost-zone gridfn -// data from the "other side" of the symmetry. +// * A symmetry_ghost_zone object describes a ghost zone which is a +// (discrete) symmetry of spacetime, either mirror-image or periodic. +// Such an object knows how to fill in ghost-zone gridfn data from +// the "other side" of the symmetry. // * An interpatch_ghost_zone object describes a ghost zone which // overlaps another patch. Such an object knows how to get ghost // zone gridfn data from the other patch. More accurately, it gets diff --git a/src/patch/grid.cc b/src/patch/grid.cc index 0e96b78..2189d5f 100644 --- a/src/patch/grid.cc +++ b/src/patch/grid.cc @@ -12,12 +12,12 @@ #include <assert.h> #include <math.h> -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/linear_map.hh" -#include "../config.hh" #include "coords.hh" #include "grid.hh" diff --git a/src/patch/make.code.defn b/src/patch/make.code.defn index a3fb1a7..1929236 100644 --- a/src/patch/make.code.defn +++ b/src/patch/make.code.defn @@ -10,7 +10,7 @@ SRCS = coords.cc \ ghost_zone.cc \ patch_interp.cc \ patch_system.cc \ - ## test_patch_system.cc + test_patch_system.cc # Subdirectories containing source files SUBDIRS = diff --git a/src/patch/patch.cc b/src/patch/patch.cc index e396c52..8035a8e 100644 --- a/src/patch/patch.cc +++ b/src/patch/patch.cc @@ -30,14 +30,14 @@ using std::printf; #include "cctk.h" -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/cpm_map.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../config.hh" #include "coords.hh" #include "grid.hh" #include "fd_grid.hh" diff --git a/src/patch/patch.hh b/src/patch/patch.hh index a7d64d7..6f812e5 100644 --- a/src/patch/patch.hh +++ b/src/patch/patch.hh @@ -18,7 +18,7 @@ // "jt/util.hh" // "jt/array.hh" // "jt/linear_map.hh" -// "../config.hh" +// "config.hh" // "coords.hh" // "grid.hh" // "fd_grid.hh" @@ -76,13 +76,12 @@ // | | <--- p.patch_interp <--- q.interpatch_ghost_zone <--> | | // +-----+ +-----+ // -// Because of the mutual pointers, we can't easily construct (say) -// p's interpatch_ghost_zone until after q itself has been constructed, -// and vice versa. Moreover, the patch_interp:: constructor -// needs the adjacent-side ghost_zone objects to already exist, -// and it needs to know the iperp range of the interpolation region, -// which can only be computed from the adjacent-patch interpatch_ghost_zone -// object. +// Because of the mutual pointers, we can't easily construct (say) p's +// interpatch_ghost_zone until after q itself has been constructed, and +// vice versa. Moreover, the patch_interp:: constructor needs the +// adjacent-side ghost_zone objects to already exist, and it needs to +// know the iperp range of the interpolation region, which can only be +// computed from the adjacent-patch interpatch_ghost_zone object. // // The solution adopted here is to use a 3-phase algorithm, ultimately // driven by the patch_system constructor: diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index 5a3a03a..3e08add 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -4,11 +4,13 @@ // // patch_interp::patch_interp // patch_interp::~patch_interp -// patch_interp::[min,max]_ipar_for_gridfn_data +// patch_interp::compute_[min,max]_ipar // patch_interp::interpolate +// patch_interp::query_interpolator // patch_interp::verify_Jacobian_sparsity_pattern_ok -// patch_interp::interp_molecule_minmax_mipar -// patch_interp::interp_molecule_posn +// patch_interp::molecule_minmax_m_ipar +// patch_interp::molecule_posn +// patch_interp::Jacobian // #include <stdio.h> @@ -18,14 +20,14 @@ #include "util_Table.h" #include "cctk.h" -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/cpm_map.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../config.hh" #include "coords.hh" #include "grid.hh" #include "fd_grid.hh" @@ -70,8 +72,7 @@ patch_interp::patch_interp(const patch_edge& my_edge_in, gridfn_coord_delta_ (my_edge().par_map().delta_fp()), gridfn_type_codes_ (min_gfn_, max_gfn_), gridfn_data_ptrs_ (min_gfn_, max_gfn_), - interp_data_buffer_ptrs_(min_gfn_, max_gfn_), - Jacobian_ptrs_ (min_gfn_, max_gfn_) // no comma + interp_data_buffer_ptrs_(min_gfn_, max_gfn_) // no comma { int status; @@ -102,9 +103,9 @@ if (status < 0) then error_exit(ERROR_EXIT, "***** patch_interp::patch_interp():\n" " can't set gridfn stride in interpolator parmameter table!\n" -" table handle=%d error status=%d\n" +" error status=%d\n" , - interp_par_table_handle_, status); /*NOTREACHED*/ + status); /*NOTREACHED*/ // // set up the gridfn type codes for the interpolator @@ -158,31 +159,17 @@ 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 stores +// 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) . // -// If Jacobian_buffer != NULL, then ... -// this function also stores 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+mipar) -// in -// (*Jacobian_buffer)(iperp, parindex, mipar) -// where we implicitly assume the Jacobian to be independent of -// ghosted_gfn, and where -// posn = molecule_posn_buffer(iperp, parindex) -// is the molecule position as returned by interp_molecule_posn() . -// -// This function calls error_exit() if any interpolation point is out -// of range or any other error is detected. -// void patch_interp::interpolate(int ghosted_min_gfn_to_interp, int ghosted_max_gfn_to_interp, - jtutil::array3d<fp>& interp_data_buffer, - jtutil::array3d<fp>* Jacobian_buffer = NULL) + jtutil::array3d<fp>& data_buffer) const { +int status; + const int N_dims = 1; const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp, ghosted_max_gfn_to_interp); @@ -193,46 +180,6 @@ const CCTK_INT *gridfn_type_codes_ptr = & gridfn_type_codes_(ghosted_min_gfn_to_interp); // -// we condition all the Jacobian-query stuff on this flag to ensure that: -// - there is always a valid "first gridfn" for the Jacobian query -// - we don't have to worry about deleting table entries that we never -// set due to there being no iperp values in this interpolation region -// -const bool Jacobian_flag - = ( (Jacobian_buffer != NULL) - && (N_gridfns != 0) - && (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0) ); - -int status, status1, status2; - -// -// if we're doing a Jacobian query, -// set Jacobian stride info in parameter table -// -if (Jacobian_flag) - then { - const int Jacobian_interp_point_stride - = 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() }; - status2 = Util_TableSetIntArray(interp_par_table_handle_, - N_dims, Jacobian_m_strides, - "Jacobian_m_strides"); - - if ((status1 < 0) || (status2 < 0)) - then error_exit(ERROR_EXIT, -"***** patch_interp::interpolate():\n" -" can't set Jacobian stride info in interpolator parmameter table!\n" -" table handle=%d error status1=%d status2=%d\n" -, - interp_par_table_handle_, status1, status2); - /*NOTREACHED*/ - } - -// // do the interpolations at each iperp // for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) @@ -261,16 +208,15 @@ if (Jacobian_flag) = my_edge(). irho_of_iperp_ipar(iperp, min_ipar()); const int start_isigma = my_edge().isigma_of_iperp_ipar(iperp, min_ipar()); - gridfn_data_ptrs_(gfn) + gridfn_data_ptrs_(ghosted_gfn) = static_cast<const void *>( & my_patch() .ghosted_gridfn(ghosted_gfn, start_irho, start_isigma) ); - interp_data_buffer_ptrs_(gfn) + interp_data_buffer_ptrs_(ghosted_gfn) = static_cast<void *>( - & interp_data_buffer(ghosted_gfn, - iperp,min_parindex) + & data_buffer(ghosted_gfn, iperp,min_parindex) ); } const void* const* const gridfn_data @@ -279,46 +225,7 @@ if (Jacobian_flag) = & interp_data_buffer_ptrs_(ghosted_min_gfn_to_interp); // - // if we're doing a Jacobian query, - // set up the query itself in the parameter table - // - if (Jacobian_flag) - then { - // - // since (we assume) the Jacobian to be independent of - // ghosted_gfn, we only query ghosted_min_gfn_to_interp; - // we pass null pointers for all other ghosted_gfns - // - for (int ghosted_gfn = ghosted_min_gfn_to_interp ; - ghosted_gfn <= ghosted_max_gfn_to_interp ; - ++ghosted_gfn) - { - Jacobian_ptrs_(ghosted_gfn) = NULL; - } - Jacobian_ptrs_(ghosted_min_gfn_to_interp) - = static_cast<CCTK_POINTER>( - & (*Jacobian_buffer)(iperp, min_parindex, 0) - ); - status = Util_TableSetPointerArray - (interp_par_table_handle_, - N_gridfns, - & Jacobian_ptrs_(ghosted_min_gfn_to_interp), - "Jacobian_pointer"); - - if (status < 0) - then error_exit(ERROR_EXIT, -"***** patch_interp::interpolate():\n" -" can't set Jacobian pointer in interpolator parmameter table!\n" -" table handle=%d error status=%d\n" -" at iperp=%d of [%d,%d]!\n" -, - interp_par_table_handle_, status, - iperp, min_iperp(), max_iperp()); - } - - // - // interpolate (and optionally query the Jacobian) - // for all the parindex values at this iperp + // interpolate for all the parindex values at this iperp // status = CCTK_InterpLocalUniform(N_dims, interp_handle_, @@ -343,30 +250,70 @@ if (Jacobian_flag) status, iperp, min_iperp(), max_iperp()); /*NOTREACHED*/ } +} + +//****************************************************************************** // -// if we're doing a Jacobian query, -// delete query entries from the parameter table -// (so future interpolator calls don't mistakenly redo this query) +// 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. // -if (Jacobian_flag) - then { - status = Util_TableDeleteKey(interp_par_table_handle_, - "Jacobian_pointer"); - status1 = Util_TableDeleteKey(interp_par_table_handle_, - "Jacobian_interp_point_stride"); - status2 = Util_TableDeleteKey(interp_par_table_handle_, - "Jacobian_m_strides"); - if ((status < 0) || (status1 < 0) || (status2 < 0)) - then error_exit(ERROR_EXIT, -"***** patch_interp::interpolat():\n" -" can't delete Jacobian query entries\n" -" from interpolator parmameter table!\n" -" table handle=%d error status=%d status1=%d status2=%d\n" -, - interp_par_table_handle_, status, status1, status2); - /*NOTREACHED*/ - } +// Results: +// This function returns the interpolator status code. +// +int patch_interp::query_interpolator(int N_interp_points = 0, + const void *const *interp_coords = NULL) + const +{ +#ifdef NOT_YET +const bool interp_points_flag = (iperp != INT_SENTINAL_VALUE); + +xxxxx +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 interp_coords_type_code = CCTK_VARIABLE_REAL; +const int N_input_arrays = 1; +const CCTK_INT input_array_dims[N_dims] = { 0 }; +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 }; + + + // 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<const void *>(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 } //****************************************************************************** @@ -391,34 +338,12 @@ void patch_interp::verify_Jacobian_sparsity_pattern_ok() int status, status1, status2, status3; // -// do the actual interpolator query -// ... explicit arguments all set to no-op values -// so we only do the query, not any actual interpolation +// query the interpolator +// ... we don't need to set up anything in the parameter table: the +// Jacobian metainfo is returned by the interpolator on every call, +// even without any explicit request // -const int N_dims = 1; -const int N_interp_points = 0; -const int interp_coords_type_code = CCTK_VARIABLE_REAL; -const int N_input_arrays = 0; -const CCTK_INT* input_array_dims = NULL; -const CCTK_INT* input_array_type_codes = NULL; -const void* const input_arrays = NULL; -const int N_output_arrays = 0; -const CCTK_INT output_array_type_codes = NULL; -void* const 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); +status = query_interpolator(); if (status < 0) then error_exit(ERROR_EXIT, "***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n" @@ -429,7 +354,8 @@ if (status < 0) // // get Jacobian-sparsity query results from parameter table // -CCTK_INT MSS_is_fn_of_interp_coords, MSS_is_fn_of_which_output; +CCTK_INT MSS_is_fn_of_interp_coords, MSS_is_fn_of_input_array_values; +CCTK_INT Jacobian_is_fn_of_input_array_values; status1 = Util_TableGetInt(interp_par_table_handle_, &MSS_is_fn_of_interp_coords, "MSS_is_fn_of_interp_coords"); @@ -444,10 +370,8 @@ if ((status1 < 0) || (status2 < 0) || (status3 < 0)) "***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n" " can't get Jacobian sparsity info\n" " from interpolator parmameter table!\n" -" table handle=%d\n" " error status1=%d status2=%d status3=%d\n" , - interp_par_table_handle_, status1, status2, status3); /*NOTREACHED*/ // @@ -459,10 +383,10 @@ if ( MSS_is_fn_of_interp_coords || MSS_is_fn_of_input_array_values "***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n" " implementation restriction: we only grok Jacobians with\n" " fixed-sized hypercube-shaped molecules, independent of\n" -" the interpolation coordinates and the floating-point values -" MSS_is_fn_of_interp_coords=(int)%d\n" -" MSS_is_fn_of_input_array_values=(int)%d\n" -" Jacobian_is_fn_of_input_array_values=(int)%d\n" +" the interpolation coordinates and the floating-point values!\n" +" MSS_is_fn_of_interp_coords=(int)%d (we only grok 0)\n" +" MSS_is_fn_of_input_array_values=(int)%d (we only grok 0)\n" +" Jacobian_is_fn_of_input_array_values=(int)%d (we only grok 0)\n" , MSS_is_fn_of_interp_coords, MSS_is_fn_of_input_array_values, @@ -475,7 +399,10 @@ if ( MSS_is_fn_of_interp_coords || MSS_is_fn_of_input_array_values // This function queries the interpolator to get the [min,max] m ipar // coordinates of the interpolation molecules. // -void patch_interp::interp_molecule_minmax_mipar(int& min_mipar, int& max_mipar) +// (This API implicitly assumes that the Jacobian sparsity is one which +// is "ok" as verified by verify_Jacobian_sparsity_pattern_ok() .) +// +void patch_interp::molecule_minmax_m_ipar(int& min_m_ipar, int& max_m_ipar) const { int status, status1, status2; @@ -484,50 +411,25 @@ int status, status1, status2; // set molecule min/max m query entries in parameter table // status1 = Util_TableSetInt(interp_par_table_handle_, - 0, "interp_molecule_min_m"); + 0, "molecule_min_m"); status2 = Util_TableSetInt(interp_par_table_handle_, - 0, "interp_molecule_max_m"); -if ((status1 < 0) || (status2 < 0) || (status1 < 0) || (status2 < 0)) + 0, "molecule_max_m"); +if ((status1 < 0) || (status2 < 0)) then error_exit(ERROR_EXIT, -"***** patch_interp::interp_molecule_minmax_mipar():\n" +"***** patch_interp::molecule_minmax_m_ipar():\n" " can't set molecule min/max m queries\n" " in interpolator parmameter table!\n" -" table handle=%d error status1=%d status2=%d\n" +" error status1=%d status2=%d\n" , - interp_par_table_handle_, status1, status2); /*NOTREACHED*/ + status1, status2); /*NOTREACHED*/ // -// do the actual interpolator query -// ... explicit arguments all set to no-op values -// so we only do the query, not any actual interpolation +// query the interpolator // -const int N_dims = 1; -const int N_interp_points = 0; -const int interp_coords_type_code = CCTK_VARIABLE_REAL; -const int N_input_arrays = 0; -const CCTK_INT* input_array_dims = NULL; -const CCTK_INT* input_array_type_codes = NULL; -const void* const input_arrays = NULL; -const int N_output_arrays = 0; -const CCTK_INT output_array_type_codes = NULL; -void* const 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); +status = query_interpolator(); if (status < 0) then error_exit(ERROR_EXIT, -"***** patch_interp::interp_molecule_minmax_mipar():\n" +"***** patch_interp::molecule_minmax_m_ipar():\n" " error return %d from interpolator query!\n" , status); /*NOTREACHED*/ @@ -536,51 +438,53 @@ if (status < 0) // get molecule min/max m query results from parameter table // ... must go through temp variables for CCTK_INT to int type conversion // -CCTK_INT interp_molecule_min_m, interp_molecule_max_m; +CCTK_INT molecule_min_m, molecule_max_m; status1 = Util_TableGetIntArray(interp_par_table_handle_, - N_dims, &interp_molecule_min_m, - "interp_molecule_min_m"); + N_dims, &molecule_min_m, + "molecule_min_m"); status2 = Util_TableGetIntArray(interp_par_table_handle_, - N_dims, &interp_molecule_max_m, - "interp_molecule_max_m"); + N_dims, &molecule_max_m, + "molecule_max_m"); if ((status1 < 0) || (status2 < 0)) then error_exit(ERROR_EXIT, -"***** patch_interp::interp_molecule_maxmax_mipar():\n" +"***** patch_interp::molecule_maxmax_m_ipar():\n" " can't get molecule min/max m query results\n" " from interpolator parmameter table!\n" -" table handle=%d error status1=%d status2=%d\n" +" error status1=%d status2=%d\n" , - interp_par_table_handle_, status1, status2); /*NOTREACHED*/ -min_mipar = interp_molecule_min_m; -max_mipar = interp_molecule_min_m; + status1, status2); /*NOTREACHED*/ +min_m_ipar = molecule_min_m; +max_m_ipar = molecule_min_m; // // delete Jacobian-sparsity-pattern query entries from the parameter table // (so future interpolator calls don't mistakenly redo this query) // status1 = Util_TableDeleteKey(interp_par_table_handle_, - "interp_molecule_min_m"); + "molecule_min_m"); status2 = Util_TableDeleteKey(interp_par_table_handle_, - "interp_molecule_max_m"); + "molecule_max_m"); if ((status1 < 0) || (status2 < 0)) then error_exit(ERROR_EXIT, "***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n" " can't delete molecule min/max m queries\n" " from interpolator parmameter table!\n" -" table handle=%d error status1=%d status2=%d\n" +" error status1=%d status2=%d\n" , - interp_par_table_handle_, status1, status2); /*NOTREACHED*/ + status1, status2); /*NOTREACHED*/ } //****************************************************************************** // -// This function queries our interpolator to find out the molecule ipar -// positions (which we implicitly assume to be independent of ghosted_gfn), -// and stores these in molecule_posn_buffer(iperp, parindex) . +// This function queries the interpolator at each iperp to find out the +// molecule ipar positions (which we implicitly assume to be independent +// of ghosted_gfn), and stores these in posn_buffer(iperp, parindex) . +// +// (This API implicitly assumes that the Jacobian sparsity is one which +// is "ok" as verified by verify_Jacobian_sparsity_pattern_ok() .) // -void patch_interp::interp_molecule_posn - (jtutil::array2d<CCTK_INT>& molecule_posn_buffer) +void patch_interp::molecule_posn(jtutil::array2d<CCTK_INT>& posn_buffer) const { const int N_dims = 1; @@ -593,17 +497,13 @@ int status; // set up the molecule-position query in the parameter table CCCTK_POINTER molecule_posn_ptrs[N_dims] - = { - static_cast<CCTK_POINTER>( - & molecule_posn_buffer(iperp, min_parindex) - ) - }; + = { static_cast<CCTK_POINTER>(& posn_buffer(iperp, min_parindex)) }; status = Util_TableSetPointerArray(interp_par_table_handle_, N_dims, molecule_posn_ptrs, - "interp_molecule_positions"); + "molecule_positions"); if (status < 0) then error_exit(ERROR_EXIT, -"***** patch_interp::interp_molecule_posn():\n" +"***** patch_interp::molecule_posn():\n" " can't set molecule position query\n" " in interpolator parmameter table at iperp=%d of [%d,%d]!\n" " error status=%d\n" @@ -614,40 +514,17 @@ int status; // interpolation-point coordinates const CCTK_INT N_interp_points = jtutil::how_many_in_range(min_parindex, max_parindex); - const int interp_coords_type_code = CCTK_VARIABLE_REAL; const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex); const void* const interp_coords[N_dims] = { static_cast<const void *>(interp_coords_ptr) }; // query the interpolator to get the molecule positions - // for all parindex at this iperp - // ... N_{input,output}_arrays set to zero - // so we only do the query, not any actual interpolation - const int N_input_arrays = 0; - const CCTK_INT* input_array_dims = NULL; - const CCTK_INT* input_array_type_codes = NULL; - const void* const input_arrays = NULL; - const int N_output_arrays = 0; - const CCTK_INT output_array_type_codes = NULL; - void* const 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); + // 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::interp_molecule_posn():\n" +"***** patch_interp::molecule_posn():\n" " error return %d from interpolator query at iperp=%d of [%d,%d]!\n" , status, iperp, min_iperp(), max_iperp()); @@ -655,20 +532,163 @@ int status; } // +// if we actually did any queries, // delete molecule-position query entry from the parameter table // (so future interpolator calls don't mistakenly redo this query) // if (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0) then { status = Util_TableDeleteKey(interp_par_table_handle_, - "interp_molecule_positions"); + "molecule_positions"); if (status < 0) then error_exit(ERROR_EXIT, -"***** patch_interp::interp_molecule_posn():\n" -" can't delete molecule position query +"***** patch_interp::molecule_posn():\n" +" can't delete molecule position query" " from interpolator parmameter table!\n" -" table handle=%d error status=%d\n" +" error status=%d\n" , - interp_par_table_handle_, status); /*NOTREACHED*/ + status); /*NOTREACHED*/ } } + +//***************************************************************************** + +// +// 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) +// +// This function stores the Jacobian in +// Jacobian_buffer(iperp, parindex, m_ipar) +// where we implicitly assume the Jacobian to be independent of +// ghosted_gfn[*], and where +// posn = posn_buffer(iperp, parindex) +// is the molecule position as returned by molecule_posn() . +// +// [*] We actually do the queries for min_gfn_ (as computed in our +// constructor). +// +// This function calls error_exit() if any interpolation point is out +// of range or any other error is detected. +// +// (This API implicitly assumes that the Jacobian sparsity is one which +// is "ok" as verified by verify_Jacobian_sparsity_pattern_ok() .) +// +void patch_interp::Jacobian(jtutil::array3d<fp>& Jacobian_buffer) + const +{ +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; + + +// +// set Jacobian stride info in parameter table +// +const int Jacobian_interp_point_stride + = 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 +status2 = Util_TableSetIntArray(interp_par_table_handle_, + N_dims, Jacobian_m_strides, + "Jacobian_m_strides"); +if ((status1 < 0) || (status2 < 0)) + then error_exit(ERROR_EXIT, +"***** patch_interp::Jacobian():\n" +" can't set Jacobian stride info in interpolator parmameter table!\n" +" error status1=%d status2=%d\n" +, + status1, status2); /*NOTREACHED*/ + + +// +// query the Jacobians at each iperp +// + 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<const void *>(interp_coords_ptr) }; + + // + // set up the Jacobian query in the parameter table + // + CCTK_POINTER const Jacobian_ptrs[N_gridfns] + = { static_cast<CCTK_POINTER>( + & Jacobian_buffer(iperp, min_parindex, 0) + ) }; + status = Util_TableSetPointerArray(interp_par_table_handle_, + N_gridfns, Jacobian_ptrs, + "Jacobian_pointer"); + if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_interp::Jacobian():\n" +" can't set Jacobian pointer in interpolator parmameter table!\n" +" error status=%d at iperp=%d of [%d,%d]!\n" +, + status, iperp, min_iperp(), max_iperp()); + /*NOTREACHED*/ + + // + // 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*/ + } + +// +// delete the query entries from the parameter table +// (so future interpolator calls don't mistakenly redo this query) +// +status = (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0) + ? Util_TableDeleteKey(interp_par_table_handle_, + "Jacobian_pointer") + : 0; +status1 = Util_TableDeleteKey(interp_par_table_handle_, + "Jacobian_interp_point_stride"); +status2 = Util_TableDeleteKey(interp_par_table_handle_, + "Jacobian_m_strides"); +if ((status < 0) || (status1 < 0) || (status2 < 0)) + then error_exit(ERROR_EXIT, +"***** patch_interp::interpolat():\n" +" can't delete Jacobian query entries\n" +" from interpolator parmameter table!\n" +" error status=%d status1=%d status2=%d\n" +, + status, status1, status2); /*NOTREACHED*/ +} diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index 9e944be..4bd402b 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -42,10 +42,10 @@ // // -// The way our coordinates are constructed, any two adjacent patches +// The way the patch coordnates are constructed, any two adjacent patches // share a common (perpendicular) coordinate. Thus we only have to do -// 1-dimensional interpolation here (in the parallel direction). -// In other words, for each iperp we do interpolation in par. +// 1-dimensional interpolation here (in the parallel direction). In +// other words, for each iperp we interpolate in par. // // In general we interpolate each gridfn at a number of distinct par // for each iperp; the integer "parindex" indexes these points. We @@ -121,23 +121,9 @@ public: // at all the coordinates specified when we were constructed, // store interpolated data in // interp_data_buffer(ghosted_gfn, iperp, parindex) - // - // if Jacobian_buffer != NULL , then ... - // also store 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+mipar) - // in - // (*Jacobian_buffer)(iperp, parindex, mipar) - // where we implicitly assume the Jacobian to be independent - // of ghosted_gfn, and where - // posn = molecule_posn_buffer(iperp, parindex) - // as returned by interp_molecule_posn() void interpolate(int ghosted_min_gfn_to_interp, int ghosted_max_gfn_to_interp, - jtutil::array3d<fp>& interp_data_buffer, - jtutil::array3d<fp>* Jacobian_buffer = NULL); + jtutil::array3d<fp>& data_buffer) const; // verify (no-op if ok, error_exit() if not) that interpolator @@ -145,20 +131,35 @@ public: // means molecules are fixed-sized hypercubes, with size/shape // independent of interpolation coordinates and the floating-point // values in the input arrays - void verify_Jacobian_sparsity_pattern_ok() const + void verify_Jacobian_sparsity_pattern_ok() const; + + // + // the remaining functions in the main client interface all + // implicitly assume that the Jacobian sparsity pattern is "ok" + // as verified by verify_Jacobian_sparsity_pattern_ok() + // // get [min,max] m ipar coordinates of interpolation molecules - void interp_molecule_minmax_mipar(int& min_mipar, int& max_mipar) - const; + void molecule_minmax_m_ipar(int& min_m_ipar, int& max_m_ipar) const; // get interpolation molecule ipar positions in // molecule_posn_buffer(iperp, parindex) // ... array type is CCTK_INT so we can pass by reference // to interpolator - void interp_molecule_posn - (jtutil::array2d<CCTK_INT>& molecule_posn_buffer) - const; + void molecule_posn(jtutil::array2d<CCTK_INT>& posn_buffer) const; + // 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) + // store Jacobian in + // Jacobian_buffer(iperp, parindex, m_ipar) + // where we implicitly assume the Jacobian to be independent of + // ghosted_gfn, and where + // posn = posn_buffer(iperp, parindex) + // as returned by molecule_posn() + void Jacobian(jtutil::array3d<fp>& Jacobian_buffer) const; // // ***** internal functions ***** @@ -180,6 +181,14 @@ private: 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) + const; + // // ***** constructor, destructor, et al ***** @@ -205,8 +214,9 @@ public: // interp_par_table_handle_in // = Cactus handle to a Cactus key/value table giving // parameters (eg order) for the interpolation operator. - // This table is not modified by any actions of this - // class. + // This class internally clones this table and modifies + // the clone, so the original table is not modified by + // any actions of this class. // // Note that this constructor looks at what type of ghost zones // are on the adjacent edges to decide how to handle the corners: @@ -216,8 +226,11 @@ public: // * if an adjacent ghost zone is an interpatch ghost zone, // we can't assume it's already been interpatch-interpolated // by the time we're called ==> we don't use data from it - // Thus the constructor requires that the adjacent-side ghost_zone - // objects already exist! + // Thus the adjacent-side ghost_zone objects must already exist! + // + // This constructor also requires that this patch's gridfns already + // exist, since we size various arrays based on the patch's min/max + // ghosted gfn. // patch_interp(const patch_edge& my_edge_in, int min_iperp_in, int max_iperp_in, @@ -282,7 +295,7 @@ private: // (par) origin and delta values of the gridfn data const fp gridfn_coord_origin_, gridfn_coord_delta_; - // gridfn type codes for interpolator + // array of gridfn type codes for interpolator // ... must be CCTK_INT so we can pass by reference to interpolator // ... values set in ctor body, const thereafter // ... index is (gfn) @@ -299,8 +312,4 @@ private: // ... we do *not* own the pointed-to data! // ... index is (gfn) mutable jtutil::array1d<void*> interp_data_buffer_ptrs_; - - // array of pointers to where Jacobian info should be stored - // if we do Jacobian queries - mutable jtutil::array1d<void*> Jacobian_ptrs_; }; diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index d4fb1c1..2cbf953 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -34,14 +34,14 @@ using std::strcmp; #include "cctk.h" -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/cpm_map.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../config.hh" #include "coords.hh" #include "grid.hh" #include "fd_grid.hh" diff --git a/src/patch/test_coords.cc b/src/patch/test_coords.cc index 4f0f117..bf50599 100644 --- a/src/patch/test_coords.cc +++ b/src/patch/test_coords.cc @@ -6,10 +6,10 @@ #include <math.h> #include <string.h> -#include "jt/stdc.h" -#include "jt/util.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" -#include "../config.hh" #include "coords.hh" using jtutil::error_exit; diff --git a/src/patch/test_coords2.cc b/src/patch/test_coords2.cc index 3dbde2e..526ec06 100644 --- a/src/patch/test_coords2.cc +++ b/src/patch/test_coords2.cc @@ -6,10 +6,10 @@ #include <math.h> #include <string.h> -#include "jt/stdc.h" -#include "jt/util.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" -#include "../config.hh" #include "coords.hh" using jtutil::fuzzy; diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc index 02138c6..16d89a4 100644 --- a/src/patch/test_patch_system.cc +++ b/src/patch/test_patch_system.cc @@ -33,14 +33,14 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" -#include "jt/stdc.h" -#include "jt/util.hh" -#include "jt/array.hh" -#include "jt/cpm_map.hh" -#include "jt/linear_map.hh" +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" using jtutil::error_exit; -#include "../config.hh" #include "coords.hh" #include "grid.hh" #include "fd_grid.hh" |