From 04cab2b135932cb205829b5c16e4d11c2d7288f0 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 4 Jul 2002 12:15:44 +0000 Subject: * finish Jacobian support in src/util * test driver for this in src/util/test_patch_system.cc --> IT WORKS!!!!! * also fix comments on #include prerequisites to say ../jtutil/ not jt/ git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@605 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- param.ccl | 11 ++ src/patch/coords.hh | 6 +- src/patch/fd_grid.hh | 17 ++- src/patch/ghost_zone.cc | 15 ++- src/patch/ghost_zone.hh | 110 +++++++++++--------- src/patch/grid.hh | 15 ++- src/patch/make.code.defn | 1 + src/patch/patch.cc | 111 -------------------- src/patch/patch.hh | 63 +----------- src/patch/patch_edge.hh | 17 +-- src/patch/patch_interp.cc | 229 ++++++++++++++++++++--------------------- src/patch/patch_interp.hh | 72 ++++++------- src/patch/patch_system.cc | 1 + src/patch/patch_system.hh | 23 +++-- src/patch/patch_system_info.hh | 13 +++ src/patch/test_patch_system.cc | 189 +++++++++++++++++++++++++--------- 16 files changed, 429 insertions(+), 464 deletions(-) diff --git a/param.ccl b/param.ccl index 3b43c96..628e5d7 100644 --- a/param.ccl +++ b/param.ccl @@ -173,6 +173,17 @@ real NP_Jacobian__perturbation_amplitude \ (0.0:* :: "any real number > 0" } 1.0e-6 +# +# true ==> gives a more thorough test of the Jacobian, +# but makes the test run much slower +# false ==> gives a slightly less thorough test, but runs faster +# +boolean NP_Jacobian__perturb_all_y_patch_points \ + "should we perturb at *all* points in the y patch, or just those with the \ + iperp which is (supposedly) involved in the interpatch interpolation?" +{ +} true; + string Jacobian_file_name "name of the Jacobian output file" { .+ :: "any valid file name" diff --git a/src/patch/coords.hh b/src/patch/coords.hh index c0ae3e8..8543f9c 100644 --- a/src/patch/coords.hh +++ b/src/patch/coords.hh @@ -9,8 +9,10 @@ // // prerequisites: -// "jt/util.hh" -// fp.hh +// +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" // //****************************************************************************** diff --git a/src/patch/fd_grid.hh b/src/patch/fd_grid.hh index e9a88c7..697def4 100644 --- a/src/patch/fd_grid.hh +++ b/src/patch/fd_grid.hh @@ -13,15 +13,14 @@ // prerequisites: // // -// // for M_PI (used by degree/radian conversions) -// "jt/util.hh" // jtutil:: stuff: -// // how_many_in_range(), -// // degrees_of_radians(), radians_of_degrees(), -// "jt/array.hh" -// "jt/linear_map.hh" -// fp.hh -// coords.hh -// grid.hh +// +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" +// "grid.hh" // //****************************************************************************** diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc index ab40101..0b7825e 100644 --- a/src/patch/ghost_zone.cc +++ b/src/patch/ghost_zone.cc @@ -451,11 +451,24 @@ interp_result_buffer_ // // construct the patch_interp object to interpolate from the *other* patch -// +// ... tell it to use corners if and only if the adjacent ghost zones +// are symmetry; cf min_ipar() and max_ipar() functions above, +// and header comments in this file +// +const bool ok_to_use_min_par_corner + = min_par_adjacent_ghost_zone().is_symmetry() + ? true + : false; +const bool ok_to_use_max_par_corner + = max_par_adjacent_ghost_zone().is_symmetry() + ? true + : false; other_patch_interp_ = new patch_interp(other_edge(), min_other_iperp_, max_other_iperp_, *min_ipar_used_, *max_ipar_used_, *other_par_, + ok_to_use_min_par_corner, + ok_to_use_max_par_corner, interp_handle, interp_par_table_handle); } diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh index 6b81ee2..99da6a9 100644 --- a/src/patch/ghost_zone.hh +++ b/src/patch/ghost_zone.hh @@ -12,17 +12,18 @@ // // // -// "jt/util.hh" -// "jt/array.hh" -// "jt/cpm_map.hh" -// "jt/linear_map.hh" -// fp.hh -// coords.hh -// grid.hh -// fd_grid.hh -// patch.hh -// patch_edge.hh -// patch_interp.hh +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/cpm_map.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" +// "grid.hh" +// "fd_grid.hh" +// "patch.hh" +// "patch_edge.hh" +// "patch_interp.hh" // //***************************************************************************** @@ -65,27 +66,37 @@ // corner in the example below. In this case we could get the gridfn // data by either of two distinct interpolation operations (presumably // from two distinct patches), which would in general give slightly -// different results. In some ideal world we might take the average -// of these, but at present we split the corner down its diagonal. -// For the points on the diagonal we make an arbitrary choice; at -// present this is that they belong to (and get their data via) the -// rho ghost zone. -// [For elliptic equations it's important to ensure that -// each grid point get a value from exactly one computation, -// and in particular we can't allow two different ghost -// zones to "own" (assign values) to the same point, since -// then the first value would be "lost" and wouldn't show -// up in the Jacobian matrix.] +// different results. In some ideal world we might do a centered +// interpolation using data from both patches, but this would be +// complicated: +// - it would require a 2-D interpolation +// - it would require bookkeeping for interpolating from multiple +// patches within the same ghost zone, indeed for the same ghost +// zone point +// At present, we follow a simpler approach: we split the corner down +// its diagonal, +// [for the points on the diagonal we make an arbitrary choice; +// at present this is that they belong to (and get their data via) +// the rho ghost zone.] +// and off-center the interpolation as necessary so each ghost-zone +// point gets data solely from the neighboring patch on its own side. // * A corner between a symmetry and an interpatch ghost zone, for // example the +x/-y or -x/+y corners in the example below. In this -// case we need to first do a symmetry operation in the neighboring -// patch so we don't have to off-center the interpolation in the -// non-corner part of the interpatch ghost zone. Then after the -// interpatch interpolation, we need to do a final symmetry operation -// in this patch to set up gridfn data in the corner. -// -// To handle all these cases, we use a 3-phase algorithm to synchronize -// ghost zones: +// case we first do a symmetry operation in the neighboring patch, +// then a fully centered interpolation (using the data just obtained +// from a symmetry operation) to get data in the non-corner part of +// the interpatch ghost zone. After the interpatch interpolation, +// we do a final symmetry operation to get gridfn data in the corner. +// +// In general, then, a ghost zone is rhomboid-shaped: iperp lies in a +// fixed interval, while ipar lies in an interval which may depend on +// iperp. In general, this shape depends on the type (symmetry vs interpatch) +// of the adjacent ghost zones. +// + +// +// To properly handle all the symmetry/interpatch cases described above, +// we use a 3-phase algorithm to synchronize ghost zones: // Phase 1: Fill in gridfn data at all the non-corner points of symmetry // ghost zones, by using the symmetries to get this data from // its "home patch" nominal grids. @@ -171,6 +182,8 @@ // adjacent edges! There are no interpolation regions inside the -x or // -y boundaries, since no interpolation is needed across those boundaries // of this patch. +// The diagonal *** line shows the boundary between the +x and +y ghost +// zones. // // Our 3-phase algorithm described above thus becomes: // Phase 1: Fill in gridfn values at points marked with "s-x" below or @@ -187,21 +200,6 @@ // (yz plane) or -y boundary (xz plane), as described by the // +z patch's -x or -y symmetry_ghost_zone object respectively. // -// The diagonal *** line shows the boundary between the +x and +y ghost -// zones; points there may be interpolated via either of the two possible -// interpatch_ghost_zone objects. -// [At present we require all points in a given ghost zone -// to be interpolated from the same neighboring patch and -// patch_interp object, so must arbitrarily choose one -// of the two neighbors for the diagonal points. In theory -// it would be better to take the average of the two neighbors, -// but in practice this doesn't matter for horizon finding -// or other elliptic stuff (it would matter for stability in -// time evolutions using multipatch finite differencing).] -// -// In general, then, a ghost zone is trapezoid-shaped: iperp lies in a -// fixed interval, while ipar lies in an interval which may depend on iperp. -// //***************************************************************************** @@ -230,7 +228,7 @@ class ghost_zone { public: // - // ***** main client interface ***** + // ***** main high-level client interface ***** // // "synchronize" a ghost zone, i.e. update the ghost-zone values // of the specified gridfns via the appropriate sequence of @@ -255,6 +253,12 @@ public: // FIXME: should these be moved out into a separate Jacobian // object/class? // + // Note that this function just computes the Jacobian of this + // ghost zone's synchronize() operation -- it does *NOT* take + // into account the 3-phase synchronization algorithm described + // in the header comments for this file. (That's done by + // patch_system::synchronize_ghost_zones_Jacobian() .) + // // n.b. terminology is // partial gridfn at x // ------------------- @@ -365,7 +369,7 @@ public: { return my_edge().max_ipar_with_corners(); } // actual min/max ipar in the ghost zone at a particular iperp - // (only defined in derived classes) + // (may depend on type of the adjacent ghost zones) virtual int min_ipar(int iperp) const = 0; virtual int max_ipar(int iperp) const = 0; @@ -463,7 +467,7 @@ class symmetry_ghost_zone { public: // - // ***** main client interface ***** + // ***** main high-level client interface ***** // // "synchronize" a ghost zone, i.e. update the ghost-zone values // of the specified gridfns via the appropriate symmetry operations @@ -601,7 +605,7 @@ class interpatch_ghost_zone { public: // - // ***** main client interface ***** + // ***** main high-level client interface ***** // // "synchronize" a ghost zone, i.e. update the ghost-zone // values of the specified gridfns via the appropriate @@ -666,6 +670,9 @@ public: fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const { assert(Jacobian_buffer_ != NULL); + // we're actually only supposed to be called with y_ipar_m + // in range, but as a backup we validate this here and return + // 0.0 for out of range if ((y_ipar_m >= min_y_ipar_m_) && (y_ipar_m <= max_y_ipar_m_)) then { const int y_iperp = Jacobian_y_iperp(x_iperp); @@ -710,8 +717,9 @@ private: const; // min/max ipar of the ghost zone for specified iperp - // taking into account our "triangular" corners - // (cf. the example at the start of this file) + // with possibly "triangular" corners depending on the type + // (symmetry vs interpatch) of the adjacent ghost zones + // (cf. comments & example at the start of this file) int min_ipar(int iperp) const; int max_ipar(int iperp) const; diff --git a/src/patch/grid.hh b/src/patch/grid.hh index 5d63698..aadec38 100644 --- a/src/patch/grid.hh +++ b/src/patch/grid.hh @@ -9,14 +9,13 @@ // prerequisites: // // -// // for M_PI (used by degree/radian conversions) -// "jt/util.hh" // jtutil:: stuff: -// // how_many_in_range(), -// // degrees_of_radians(), radians_of_degrees(), -// "jt/array.hh" -// "jt/linear_map.hh" -// fp.hh -// coords.hh +// +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" // //***************************************************************************** diff --git a/src/patch/make.code.defn b/src/patch/make.code.defn index 1929236..148b738 100644 --- a/src/patch/make.code.defn +++ b/src/patch/make.code.defn @@ -9,6 +9,7 @@ SRCS = coords.cc \ patch.cc \ ghost_zone.cc \ patch_interp.cc \ + patch_info.cc \ patch_system.cc \ test_patch_system.cc diff --git a/src/patch/patch.cc b/src/patch/patch.cc index 8035a8e..dfa4369 100644 --- a/src/patch/patch.cc +++ b/src/patch/patch.cc @@ -18,9 +18,6 @@ // patch::edge_adjacent_to_patch // patch::assert_all_ghost_zones_fully_setup // -// patch_info::grid_array_pars -// patch_info::grid_pars -// #include #include @@ -413,111 +410,3 @@ max_rho_ghost_zone().assert_fully_setup(); min_sigma_ghost_zone().assert_fully_setup(); max_sigma_ghost_zone().assert_fully_setup(); } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function computes, and returns a reference to, a -// struct grid_arrays::grid_array_pars from the info in a -// struct patch_info and the additional information in the arguments. -// -// The result refers to an internal static buffer in this function; the -// usual caveats about lifetimes/overwriting apply. -// -// Arguments: -// N_ghost_points = Width in grid points of all ghost zones. -// N_extend_points = Number of grid points to extend each patch past -// "just touching" so as to overlap neighboring patches. -// Thus patches overlap by -// N_overlap_points = 2*N_extend_points + 1 -// grid points. For example, with N_extend_points == 2, -// here are the grid points of two neighboring patches: -// x x x x x X X -// | -// O O o o o o o -// Here | marks the "just touching" boundary, -// x and o the grid points before this extension, -// and X and O the extra grid points added by this -// extension. -// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees. -// -const grid_arrays::grid_array_pars& - patch_info::grid_array_pars(int N_ghost_points, int N_extend_points, - fp delta_drho_dsigma) - const -{ -static - struct grid_arrays::grid_array_pars grid_array_pars_buffer; - -grid_array_pars_buffer.min_irho - = jtutil::round::to_integer(min_drho /delta_drho_dsigma); -grid_array_pars_buffer.min_isigma - = jtutil::round::to_integer(min_dsigma/delta_drho_dsigma); -grid_array_pars_buffer.max_irho - = grid_array_pars_buffer.min_irho - + jtutil::round::to_integer( - (max_drho -min_drho ) / delta_drho_dsigma - ); -grid_array_pars_buffer.max_isigma - = grid_array_pars_buffer.min_isigma - + jtutil::round::to_integer( - (max_dsigma-min_dsigma) / delta_drho_dsigma - ); -grid_array_pars_buffer.min_irho -= N_extend_points; -grid_array_pars_buffer.min_isigma -= N_extend_points; -grid_array_pars_buffer.max_irho += N_extend_points; -grid_array_pars_buffer.max_isigma += N_extend_points; - -grid_array_pars_buffer.min_rho_N_ghost_points = N_ghost_points; -grid_array_pars_buffer.max_rho_N_ghost_points = N_ghost_points; -grid_array_pars_buffer.min_sigma_N_ghost_points = N_ghost_points; -grid_array_pars_buffer.max_sigma_N_ghost_points = N_ghost_points; - -return grid_array_pars_buffer; -} - -//****************************************************************************** -// -// -// This function computes, and returns a reference to, a -// struct grid_arrays::grid_pars from the info in a struct patch_info -// and the additional information in the arguments. -// -// The result refers to an internal static buffer in this function; the -// usual caveats about lifetimes/overwriting apply. -// -// Arguments: -// N_extend_points = Number of grid points to extend each patch past -// "just touching" so as to overlap neighboring patches. -// Thus patches overlap by 2*N_extend_points + 1 grid -// points. For example, with N_extend_points == 2, here -// are the grid points of two neighboring patches: -// x x x x x X X -// | -// O O o o o o o -// Here | marks the "just touching" boundary, -// x and o the grid points before this extension, -// and X and O the extra grid points added by this -// extension. -// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees. -// -const grid::grid_pars& patch_info::grid_pars(int N_extend_points, - fp delta_drho_dsigma) - const -{ -static - struct grid::grid_pars grid_pars_buffer; - -const fp extend_drho_dsigma = fp(N_extend_points) * delta_drho_dsigma; - -grid_pars_buffer. min_drho = min_drho - extend_drho_dsigma; -grid_pars_buffer.delta_drho = delta_drho_dsigma; -grid_pars_buffer. max_drho = max_drho + extend_drho_dsigma; -grid_pars_buffer. min_dsigma = min_dsigma - extend_drho_dsigma; -grid_pars_buffer.delta_dsigma = delta_drho_dsigma; -grid_pars_buffer. max_dsigma = max_dsigma + extend_drho_dsigma; - -return grid_pars_buffer; -} diff --git a/src/patch/patch.hh b/src/patch/patch.hh index 6f812e5..b58882d 100644 --- a/src/patch/patch.hh +++ b/src/patch/patch.hh @@ -3,7 +3,6 @@ // // ***** how patch boundaries are handled ***** // patch - abstract base class to describe a coordinate/grid patch -// patch_info - POD struct with summary info for patch construction // // z_patch - derived class for a +/- z patch // x_patch - derived class for a +/- x patch @@ -15,10 +14,11 @@ // // // -// "jt/util.hh" -// "jt/array.hh" -// "jt/linear_map.hh" +// "stdc.h" // "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/linear_map.hh" // "coords.hh" // "grid.hh" // "fd_grid.hh" @@ -424,61 +424,6 @@ private: //***************************************************************************** //***************************************************************************** -// -// This (POD, and hence static-initializable) struct gives a minimal -// set of information which varies from one patch to another. -// -// The member functions allow a patch to be constructed from this struct -// together with the additional information provided by their arguments. -// This doesn't allow the most general possible type of patch (eg it -// constrains all ghost zones to have the same width), but does cover -// the common casees. -// -// Arguments for member functions: -// N_ghost_points = Width in grid points of all ghost zones. -// N_extend_points = Number of grid points to extend each patch past -// "just touching" so as to overlap neighboring patches. -// Thus patches overlap by -// N_overlap_points = 2*N_extend_points + 1 -// grid points. For example, with N_extend_points == 2, -// here are the grid points of two neighboring patches: -// x x x x x X X -// | -// O O o o o o o -// Here | marks the "just touching" boundary, -// x and o the grid points before this extension, -// and X and O the extra grid points added by this -// extension. -// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees. -// -struct patch_info - { - const char* name; - bool is_plus; - char ctype; - fp min_drho, max_drho; - fp min_dsigma, max_dsigma; - - // compute and return reference to struct grid_arrays::grid_array_pars - // ... result refers to internal static buffer; - // the usual caveats about lifetimes/overwriting apply - const grid_arrays::grid_array_pars& - grid_array_pars(int N_ghost_points, int N_extend_points, - fp delta_drho_dsigma) - const; - - // compute and return reference to struct grid::grid_pars - // ... result refers to internal static buffer; - // the usual caveats about lifetimes/overwriting apply - const grid::grid_pars& grid_pars(int N_extend_points, - fp delta_drho_dsigma) - const; - }; - -//***************************************************************************** -//***************************************************************************** -//***************************************************************************** - // // This class describes a +/- z patch. It doesn't define any new // functions not already present in class patch ; it "just" defines diff --git a/src/patch/patch_edge.hh b/src/patch/patch_edge.hh index 58191b7..6140172 100644 --- a/src/patch/patch_edge.hh +++ b/src/patch/patch_edge.hh @@ -7,14 +7,15 @@ // // // -// "jt/util.hh" -// "jt/array.hh" -// "jt/linear_map.hh" -// fp.hh -// coords.hh -// grid.hh -// fd_grid.hh -// patch.hh +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" +// "grid.hh" +// "fd_grid.hh" +// "patch.hh" // //***************************************************************************** diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc index 0c7181f..c19dd42 100644 --- a/src/patch/patch_interp.cc +++ b/src/patch/patch_interp.cc @@ -4,14 +4,13 @@ // // patch_interp::patch_interp // 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_ipar_m // patch_interp::molecule_posn // patch_interp::Jacobian +// patch_interp::query_interpolator // #include @@ -56,22 +55,28 @@ patch_interp::patch_interp(const patch_edge& my_edge_in, const jtutil::array1d& min_parindex_array_in, const jtutil::array1d& max_parindex_array_in, const jtutil::array2d& interp_par_in, + bool ok_to_use_min_par_corner, + bool ok_to_use_max_par_corner, int interp_handle_in, int interp_par_table_handle_in) : my_patch_(my_edge_in.my_patch()), my_edge_(my_edge_in), min_gfn_(my_patch().ghosted_min_gfn()), max_gfn_(my_patch().ghosted_max_gfn()), min_iperp_(min_iperp_in), max_iperp_(max_iperp_in), - min_ipar_(compute_min_ipar()), - max_ipar_(compute_max_ipar()), + min_ipar_(ok_to_use_min_par_corner + ? my_edge_in.min_ipar_with_corners() + : my_edge_in.min_ipar_without_corners()), + max_ipar_(ok_to_use_max_par_corner + ? my_edge_in.max_ipar_with_corners() + : my_edge_in.max_ipar_without_corners()), 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)), - // note interpolator coordinate origin = par of min_ipar() point + // note interpolator coordinate origin = min_ipar_ // cf. comments in "patch_interp.hh" on gridfn array subscripting - gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar())), + gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar_)), gridfn_coord_delta_ (my_edge().par_map().delta_fp()), gridfn_type_codes_ (min_gfn_, max_gfn_), gridfn_data_ptrs_ (min_gfn_, max_gfn_), @@ -131,35 +136,6 @@ Util_TableDestroy(interp_par_table_handle_); //***************************************************************************** -// -// These functions compute the [min,max] range of ipar in this patch -// from which we may use gridfn data in interpolation. This computation -// depends on our adjacent ghost zones' types as described in the comments -// to our constructor, but needs only my_patch_ and my_edge_ members of -// this object to be valid. -// -int patch_interp::compute_min_ipar() const -{ -const ghost_zone& min_par_adjacent_ghost_zone - = my_patch() - .ghost_zone_on_edge(my_edge().min_par_adjacent_edge()); -return min_par_adjacent_ghost_zone.is_symmetry() - ? my_edge().min_ipar_with_corners() - : my_edge().min_ipar_without_corners(); -} - -int patch_interp::compute_max_ipar() const -{ -const ghost_zone& max_par_adjacent_ghost_zone - = my_patch() - .ghost_zone_on_edge(my_edge().max_par_adjacent_edge()); -return max_par_adjacent_ghost_zone.is_symmetry() - ? my_edge().max_ipar_with_corners() - : my_edge().max_ipar_without_corners(); -} - -//***************************************************************************** - // // This function interpolates the specified range of ghosted gridfns // at all the coordinates specified when we were constructed. (It does @@ -258,89 +234,6 @@ const CCTK_INT *gridfn_type_codes_ptr //****************************************************************************** //****************************************************************************** -// -// This function queries the interpolator with whatever arguments are -// 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: -// 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(const char function_name[], int iperp) - const -{ -int status; -const int N_dims = 1; - -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] - = { 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*/ - -return status; -} - -//****************************************************************************** - // // This function verifies that our interpolator has a Jacobian sparsity // pattern which we grok: at present this means that the molecules are @@ -505,6 +398,7 @@ int status; for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp) { const int min_parindex = min_parindex_array_(iperp); + const int max_parindex = max_parindex_array_(iperp); // set up the molecule-position query in the parameter table CCTK_POINTER molecule_posn_ptrs[N_dims] @@ -523,9 +417,21 @@ int status; status); /*NOTREACHED*/ // 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) + // for all parindex at this iperp; the interpolator will + // store the parindex-min_parindex values (cf comments + // on array subscripting at the start of "patch_interp.hh") + // directly in the posn_buffer(iperp,*) array) query_interpolator("molecule_positions", iperp); + + // convert the molecule positions from parindex-min_ipar + // to parindex values (again, cf comments on array subscripting + // at the start of "patch_interp.hh") + for (int parindex = min_parindex ; + parindex <= max_parindex ; + ++parindex) + { + posn_buffer(iperp, parindex) += min_ipar(); + } } // @@ -658,3 +564,86 @@ if ((status < 0) || (status1 < 0) || (status2 < 0)) , status, status1, status2); /*NOTREACHED*/ } + +//****************************************************************************** + +// +// This function queries the interpolator with whatever arguments are +// 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: +// 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(const char function_name[], int iperp) + const +{ +int status; +const int N_dims = 1; + +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] + = { 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*/ + +return status; +} diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh index b98d39f..5693a4a 100644 --- a/src/patch/patch_interp.hh +++ b/src/patch/patch_interp.hh @@ -1,4 +1,4 @@ -// patch_interp.hh -- patch interpolation region +// patch_interp.hh -- interpolation from a patch // $Id$ // @@ -9,16 +9,15 @@ // // "cctk.h" // -// "jt/util.hh" -// "jt/array.hh" -// "jt/linear_map.hh" -// "jt/interpolate.hh" -// fp.hh -// coords.hh -// grid.hh -// fd_grid.hh -// patch.hh -// patch_edge.hh +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" +// "grid.hh" +// "fd_grid.hh" +// "patch.hh" // //***************************************************************************** @@ -36,16 +35,11 @@ // patch's coordinates into our coordinates. // +// // A patch_interp defines a "patch interpolation region", the region of // its owning patch from which this interpolation will use gridfn data. // -// -// Since all patch_interp member functions are const, a patch_interp -// object is effectively always const whether declared that way or not, -// so there's no harm in always declaring patch_interp objects as const. -// - // // The way the patch coordnates are constructed, any two adjacent patches // share a common (perpendicular) coordinate. Thus we only have to do @@ -94,20 +88,23 @@ // the gridfn data is a bit tricky: The interpolator accesses the // gridfn data using the generic (1-D) subscripting expression // data[offset + i*stride] -// where i is the data array index. However, we don't want to use -// offset, because it has to be supplied in the parameter table as -// an array subscripted by gfn, and so would require changing the -// parameter table for each call on interpolate() (with potentially +// where i is the data array index. However, we'd rather not use +// offset , because it has to be supplied in the parameter table as +// an array subscripted by gfn , and so would require changing the +// parameter table for each call on interpolate() (with potentially // different numbers of gridfns being interpolated). Instead, at each -// iperp we use i = ipar-min_ipar, so the default offset=0 makes the -// subscripting expression zero for ipar = min_ipar. This also makes -// the interpolator's min_i = 0 and max_i be dims-1 (both the defaults), -// so those also don't have to be set in the parameter table. We set -// the interpolator's data coordinate origin to the par coordinate for -// min_ipar, so it correctly maps i --> par. With this strategy we can -// share the interpolator parameter table across all the iperp values, -// and we don't need to modify the parameter table at all after the -// initial setup in our constructor. +// iperp we use i = ipar-min_ipar , so the default offset=0 makes +// the subscripting expression zero for ipar = min_ipar . This also +// makes the interpolator's min_i = 0 and max_i be dims-1 (both +// the defaults), so those also don't have to be set in the parameter +// table either. We set the interpolator's data coordinate origin to +// the par coordinate for min_ipar , so it correctly maps i --> par . +// With this strategy we can share the interpolator parameter table +// across all the iperp values, and we don't need to modify the +// parameter table at all after the initial setup in our constructor. +// However, we do have to adjust the molecule positions in +// patch_interp::molecule_posn() , since the interpolator will return +// i values, while molecule_posn() needs ipar values. // class patch_interp @@ -179,15 +176,7 @@ private: int min_iperp() const { return min_iperp_; } int max_iperp() const { return max_iperp_; } - // compute [min,max] ipar in this patch from which we may use - // gridfn data in interpolation -- depends on adjacent ghost zones' - // types as described in comments to constructor (below), but needs - // only my_patch_ and my_edge_ members of this object to be valid - int compute_min_ipar() const; - int compute_max_ipar() const; - - // min/max (iperp,ipar) of the gridfn data we'll use - // (uses cached values set up by our constructor) + // min/max (iperp,ipar) of the gridfn data to use for interpolation int min_ipar() const { return min_ipar_; } int max_ipar() const { return max_ipar_; } @@ -223,6 +212,9 @@ public: // are unused. We keep a reference to this array, so it // should have a lifetime at last as long as that of this // object. + // ok_to_use_[min,max]_par_corner + // = Boolean flags saying whether or not we should use gridfn + // data from the [min,max]_par corners in the interpolation. // interp_handle_in = Cactus handle to the interpolation operator. // interp_par_table_handle_in // = Cactus handle to a Cactus key/value table giving @@ -250,6 +242,8 @@ public: const jtutil::array1d& min_parindex_array_in, const jtutil::array1d& max_parindex_array_in, const jtutil::array2d& interp_par_in, + bool ok_to_use_min_par_corner, + bool ok_to_use_max_par_corner, int interp_handle_in, int interp_par_table_handle_in); ~patch_interp(); diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index fb30ded..4a3e3d1 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -50,6 +50,7 @@ using jtutil::error_exit; #include "patch_edge.hh" #include "patch_interp.hh" #include "ghost_zone.hh" +#include "patch_info.hh" #include "patch_system.hh" #include "patch_system_info.hh" diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh index e00758f..8fd9541 100644 --- a/src/patch/patch_system.hh +++ b/src/patch/patch_system.hh @@ -9,16 +9,19 @@ // // // -// "jt/util.hh" -// "jt/array.hh" -// "jt/linear_map.hh" -// fp.hh -// coords.hh -// grid.hh -// patch.hh -// patch_edge.hh -// ghost_zone.hh -// patch_interp.hh +// "stdc.h" +// "config.hh" +// "../jtutil/util.hh" +// "../jtutil/array.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" +// "grid.hh" +// "fd_grid.hh" +// "patch.hh" +// "patch_edge.hh" +// "ghost_zone.hh" +// "patch_interp.hh" +// "patch_info.hh" // //****************************************************************************** diff --git a/src/patch/patch_system_info.hh b/src/patch/patch_system_info.hh index 72b0191..dd6cf9f 100644 --- a/src/patch/patch_system_info.hh +++ b/src/patch/patch_system_info.hh @@ -1,6 +1,19 @@ // patch_system_info.hh -- static info describing various types of patch systems // $Id$ +// +// prerequisites: +// +// +// +// "../jtutil/util.hh" +// "../jtutil//array.hh" +// "../jtutil/linear_map.hh" +// "coords.hh" +// "grid.hh" +// "patch_info.hh" +// + //****************************************************************************** namespace patch_system_info diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc index 6371305..dbde22b 100644 --- a/src/patch/test_patch_system.cc +++ b/src/patch/test_patch_system.cc @@ -101,6 +101,7 @@ namespace { void test_synchronize_Jacobians(patch_system& ps, int test_gfn, int NP_test_gfn, fp perturbation_amplitude, + bool perturb_all_y_patch_points, const char Jacobian_file_name[]); void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones); @@ -203,6 +204,7 @@ else if (STRING_EQUAL(which_test, "synchronize Jacobian")) then test_synchronize_Jacobians(ps, test_fn_gfn, test_fn_copy_gfn, NP_Jacobian__perturbation_amplitude, + (NP_Jacobian__perturb_all_y_patch_points != 0), Jacobian_file_name); else if (STRING_EQUAL(which_test, "derivatives")) @@ -238,32 +240,51 @@ CCTK_VInfo(CCTK_THORNSTRING, "destroying patch system"); // synchronize test_gfn // print ths synchronized function // set up the same test function in NP_test_gfn on the nominal grid -// for each patch p and ghost zone pgz +// for each patch xp and ghost zone xgz // { // compute the synchronize() Jacobian for this ghost zone -// (via pgz.compute_Jacobian() et al) -// for each point x in (p,pgz) +// (via xgz.compute_Jacobian() et al) +// for each point x in (p,xgz) // { -// for each point y in gz.Jacobian_y_patch() +// for each point y in xgz.Jacobian_y_patch() // { +// if ( !perturb_all_y_patch_points +// && (y iperp != xgz Jacobian y iperp) ) +// then continue; // *** LOOP CONTROL *** // const fp save_y_gridfn = NP_test_gfn at y // NP_test_gfn at y += perturbation_amplitude // synchronize NP_test_gfn // NP_Jacobian = (NP_test_gfn at x - test_gfn at x) // / perturbation_amplitude // NP_test_gfn at y = save_y_gridfn -// Jacobian = pgz.Jacobian(...) -// if (m in molecule || the Jacobians differ) -// then print Jacobian, NP_Jacobian, -// Jacobian error +// Jacobian = xgz.Jacobian(...) +// print Jacobians to output file +// if (the Jacobians differ) +// then { +// print lots of debugging info to stdout +// goto done; +// } // } // } // } +// done: +// +// Note that this test can be *very* slow to run! Eg. on a 731 MHz PIII +// laptop, a full-sphere patch system with 5 degree resolution, a ghost zone +// iperp width of 2 points, and perturb_all_y_patch_points set to true, +// took around 25 minutes of cpu time, and produced an 20 MB Jacobian +// output file. Even with perturb_all_y_patch_points set to false, it +// took 1 minute 15 seconds cpu time, and produced a 1 MB Jacobian output +// file. // // Arguments: // ps = (in out) THe patch system in/on which to do the computations. // test_gfn, NP_test_gfn = The gfns of two ghosted test gridfns. // perturbation_amplitude = The perturbation amplitude for the NP Jacobian. +// perturb_all_y_patch_points +// = true ==> When computing the NP Jacobian, try perturbing at +// every point in the y patch. This gives +// false ==> Only try perturbing at y_iper == Jacobian_y_iperp. // Jacobian_file_name = The name of the output file to which both Jacobians // should be written. // @@ -271,8 +292,11 @@ namespace { void test_synchronize_Jacobians(patch_system& ps, int test_gfn, int NP_test_gfn, fp perturbation_amplitude, + bool perturb_all_y_patch_points, const char Jacobian_file_name[]) { +CCTK_VInfo(CCTK_THORNSTRING, "testing synchronize() Jacobian..."); + setup_sym_fn_xyz(ps, test_gfn, false); ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn); ps.print_ghosted_gridfn(test_fn_gfn, "test_fn.dat"); @@ -302,10 +326,10 @@ fprintf(fileptr, "# column 15 = Jacobian\n"); fprintf(fileptr, "# column 16 = NP_Jacobian\n"); fprintf(fileptr, "# column 17 = Jacobian error\n"); - //*** for each patch p and ghost zone pgz - for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + //*** for each patch p and ghost zone xgz + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) { - patch& p = ps.ith_patch(pn); + patch& xp = ps.ith_patch(xpn); // n.b. these loops must use _int_ variables for the loop // to terminate! @@ -313,72 +337,144 @@ fprintf(fileptr, "# column 17 = Jacobian error\n"); { for (int want_rho = false ; want_rho <= true ; ++want_rho) { - const patch_edge& pe = p.minmax_ang_patch_edge(want_min, want_rho); - ghost_zone& pgz = p.minmax_ang_ghost_zone(want_min, want_rho); + const patch_edge& xe = xp.minmax_ang_patch_edge(want_min, want_rho); + ghost_zone& xgz = xp.minmax_ang_ghost_zone(want_min, want_rho); + + patch& yp = xgz.Jacobian_y_patch(); + const patch_edge& ye = xgz.Jacobian_y_edge(); - pgz.compute_Jacobian(test_gfn, test_gfn, + CCTK_VInfo(CCTK_THORNSTRING, + " testing x patch %s, edge %s", + xp.name(), xe.name()); + CCTK_VInfo(CCTK_THORNSTRING, + " y patch %s, edge %s", + yp.name(), ye.name()); + + xgz.compute_Jacobian(test_gfn, test_gfn, true, true, true); // want *all* of ghost zone - const int Jacobian_y_min_ipar_m = pgz.Jacobian_y_min_ipar_m(); - const int Jacobian_y_max_ipar_m = pgz.Jacobian_y_max_ipar_m(); + const int Jacobian_y_min_ipar_m = xgz.Jacobian_y_min_ipar_m(); + const int Jacobian_y_max_ipar_m = xgz.Jacobian_y_max_ipar_m(); - //*** for each point x in (p,pgz) - for (int x_iperp = pgz.min_iperp() ; - x_iperp <= pgz.max_iperp() ; + //*** for each point x in (p,xgz) + for (int x_iperp = xgz.min_iperp() ; + x_iperp <= xgz.max_iperp() ; ++x_iperp) { - for (int x_ipar = pgz.min_ipar(x_iperp) ; - x_ipar <= pgz.max_ipar(x_iperp) ; + for (int x_ipar = xgz.min_ipar(x_iperp) ; + x_ipar <= xgz.max_ipar(x_iperp) ; ++x_ipar) { - const int x_irho = pe. irho_of_iperp_ipar(x_iperp, x_ipar); - const int x_isigma = pe.isigma_of_iperp_ipar(x_iperp, x_ipar); + const int x_irho = xe. irho_of_iperp_ipar(x_iperp, x_ipar); + const int x_isigma = xe.isigma_of_iperp_ipar(x_iperp, x_ipar); - patch& q = pgz.Jacobian_y_patch(); - const patch_edge& qe = pgz.Jacobian_y_edge(); - const int Jacobian_y_iperp = pgz.Jacobian_y_iperp(x_iperp); + const int Jacobian_y_iperp = xgz.Jacobian_y_iperp(x_iperp); const int Jacobian_y_ipar_posn - = pgz.Jacobian_y_ipar_posn(x_iperp, x_ipar); + = xgz.Jacobian_y_ipar_posn(x_iperp, x_ipar); //*** for each point y in gz.Jacobian_patch() - for (int y_irho = q.min_irho() ; y_irho <= q.max_irho() ; ++y_irho) + for (int y_irho = yp.min_irho() ; + y_irho <= yp.max_irho() ; + ++y_irho) { - for (int y_isigma = q.min_isigma() ; - y_isigma <= q.max_isigma() ; + for (int y_isigma = yp.min_isigma() ; + y_isigma <= yp.max_isigma() ; ++y_isigma) { + const int y_iperp = ye.iperp_of_irho_isigma(y_irho,y_isigma); + const int y_ipar = ye. ipar_of_irho_isigma(y_irho,y_isigma); + + if ( !perturb_all_y_patch_points && (y_iperp != Jacobian_y_iperp) ) + then continue; // *** LOOP CONTROL *** + // compute the NP Jacobian const fp save_y_gridfn - = q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma); - q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) + = yp.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma); + yp.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) += perturbation_amplitude; ps.synchronize_ghost_zones(NP_test_gfn, NP_test_gfn); const fp NP_Jacobian - = ( p.ghosted_gridfn(NP_test_gfn, x_irho,x_isigma) - - p.ghosted_gridfn( test_gfn, x_irho,x_isigma) ) + = ( xp.ghosted_gridfn(NP_test_gfn, x_irho,x_isigma) + - xp.ghosted_gridfn( test_gfn, x_irho,x_isigma) ) / perturbation_amplitude; - q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) = save_y_gridfn; + yp.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) = save_y_gridfn; // compute the query Jacobian - const int y_iperp = qe.iperp_of_irho_isigma(y_irho,y_isigma); - const int y_ipar = qe. ipar_of_irho_isigma(y_irho,y_isigma); const int y_ipar_m = y_ipar - Jacobian_y_ipar_posn; const bool m_in_molecule = (y_iperp == Jacobian_y_iperp) && (y_ipar_m >= Jacobian_y_min_ipar_m) && (y_ipar_m <= Jacobian_y_max_ipar_m); - const fp Jacobian = pgz.Jacobian(x_iperp, x_ipar, y_ipar_m); + const fp Jacobian = m_in_molecule + ? xgz.Jacobian(x_iperp, x_ipar, y_ipar_m) + : 0.0; // print the results - if (m_in_molecule || (Jacobian != NP_Jacobian)) - then fprintf(fileptr, + const fp error = Jacobian - NP_Jacobian; + fprintf(fileptr, "%d %d %d\t%d %d\t%d %d\t%d %d %d\t%d %d\t%d %d\t%.10g\t%.10g\t%e\n", - p.patch_number(), pe.is_min(), pe.is_rho(), - x_iperp, x_ipar, x_irho, x_isigma, - q.patch_number(), qe.is_min(), qe.is_rho(), - y_iperp, y_ipar, y_irho, y_isigma, - double(Jacobian), double(NP_Jacobian), - double(Jacobian-NP_Jacobian)); + xp.patch_number(), xe.is_min(), xe.is_rho(), + x_iperp, x_ipar, x_irho, x_isigma, + yp.patch_number(), ye.is_min(), ye.is_rho(), + y_iperp, y_ipar, y_irho, y_isigma, + double(Jacobian), double(NP_Jacobian), + double(error)); + + // debugging code in case the Jacobian is wrong :( + if (jtutil::abs(error) > 1.0e-6) + then { + printf("### large Jacobian error!\n"); + + printf("x: p patch %s, edge %s\n", xp.name(), xe.name()); + printf("y: q patch %s, edge %s\n", yp.name(), ye.name()); + + const fp x_rho = xp.rho_of_irho (x_irho); + const fp x_sigma = xp.sigma_of_isigma(x_isigma); + const fp y_rho = yp.rho_of_irho (y_irho); + const fp y_sigma = yp.sigma_of_isigma(y_isigma); + + const fp x_drho = jtutil::degrees_of_radians(x_rho); + const fp x_dsigma = jtutil::degrees_of_radians(x_sigma); + const fp y_drho = jtutil::degrees_of_radians(y_rho); + const fp y_dsigma = jtutil::degrees_of_radians(y_sigma); + + printf( +"x iperp=%d ipar=%d irho=%d isigma=%d drho=%g dsigma=%g\n", + x_iperp, x_ipar, x_irho, x_isigma, x_drho, x_dsigma); + printf( +"y iperp=%d ipar=%d irho=%d isigma=%d drho=%g dsigma=%g\n", + y_iperp, y_ipar, y_irho, y_isigma, y_drho, y_dsigma); + + const fp x_mu = xp.mu_of_rho_sigma (x_rho, x_sigma); + const fp x_nu = xp.nu_of_rho_sigma (x_rho, x_sigma); + const fp x_phi = xp.phi_of_rho_sigma(x_rho, x_sigma); + const fp y_mu = yp.mu_of_rho_sigma (y_rho, y_sigma); + const fp y_nu = yp.nu_of_rho_sigma (y_rho, y_sigma); + const fp y_phi = yp.phi_of_rho_sigma(y_rho, y_sigma); + + const fp x_dmu = jtutil::degrees_of_radians(x_mu); + const fp x_dnu = jtutil::degrees_of_radians(x_nu); + const fp x_dphi = jtutil::degrees_of_radians(x_phi); + const fp y_dmu = jtutil::degrees_of_radians(y_mu); + const fp y_dnu = jtutil::degrees_of_radians(y_nu); + const fp y_dphi = jtutil::degrees_of_radians(y_phi); + + printf("x dmu=%g dnu=%g dphi=%g\n", x_dmu, x_dnu, x_dphi); + printf("y dmu=%g dnu=%g dphi=%g\n", y_dmu, y_dnu, y_dphi); + + printf("Jacobian=%.10g\tNP_Jacobian=%.10g\terror=%e\n", + double(Jacobian), double(NP_Jacobian), + double(error)); + + printf("Jacobian_y_[min,max]_ipar_m=[%d,%d]\n", + Jacobian_y_min_ipar_m, Jacobian_y_max_ipar_m); + printf("Jacobian_y_iperp=%d Jacobian_y_ipar_posn=%d\n", + Jacobian_y_iperp, Jacobian_y_ipar_posn); + printf("y_ipar_m=%d\n", y_ipar_m); + + printf("###\n"); + goto done; + } } } } @@ -388,6 +484,7 @@ fprintf(fileptr, "# column 17 = Jacobian error\n"); } +done: fclose(fileptr); } } -- cgit v1.2.3