diff options
| -rw-r--r-- | README | 1 | ||||
| -rw-r--r-- | interface.ccl | 4 | ||||
| -rw-r--r-- | param.ccl | 5 | ||||
| -rw-r--r-- | schedule.ccl | 20 | ||||
| -rw-r--r-- | src/make.code.defn | 2 | ||||
| -rw-r--r-- | src/patch/ghost_zone.cc | 14 | ||||
| -rw-r--r-- | src/patch/ghost_zone.hh | 160 | ||||
| -rw-r--r-- | src/patch/test_patch_system.cc | 190 |
8 files changed, 306 insertions, 90 deletions
@@ -66,3 +66,4 @@ Portability This thorn is known to compile and run successfully on the following systems: * Red Hat GNU/Linux, gcc/g++ version 2.96 20000731 (Red Hat Linux 7.1 2.96-98) +* OpenBSD, gcc/g++ version 2.95.3 20010125 (prerelease) diff --git a/interface.ccl b/interface.ccl index 3cb9303..6a8c611 100644 --- a/interface.ccl +++ b/interface.ccl @@ -2,5 +2,5 @@ # $Header$ implements: AHFinder -inherits: grid -inherits: einstein +##inherits: grid +##inherits: einstein @@ -172,3 +172,8 @@ real NP_Jacobian__perturbation_amplitude \ { (0.0:* :: "any real number > 0" } 1.0e-6 + +string Jacobian_file_name "name of the Jacobian output file" +{ +.+ :: "any valid file name" +} "Jacobian.dat" diff --git a/schedule.ccl b/schedule.ccl index afe9ce2..f413ed0 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -8,14 +8,14 @@ # ==> quick-n-dirty hack for now: just schedule at POSTINITIAL # -schedule AHFinderDirect_driver at CCTK_POSTINITIAL -{ -lang: C -options: global -} "find apparent horizon(s)" - -##schedule test_patch_system at CCTK_POSTINITIAL +##schedule AHFinderDirect_driver at CCTK_POSTINITIAL ##{ -##LANG: C -##OPTIONS: global -##} "test driver to verify that src/util/ and src/jtutil/ code works properly" +##lang: C +##options: global +##} "find apparent horizon(s)" + +schedule test_patch_system at CCTK_POSTINITIAL +{ +LANG: C +OPTIONS: global +} "test driver to verify that src/util/ and src/jtutil/ code works properly" diff --git a/src/make.code.defn b/src/make.code.defn index b60ca69..cc06be4 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -5,4 +5,4 @@ SRCS = # Subdirectories containing source files -SUBDIRS = jtutil util gr +SUBDIRS = jtutil util ## gr diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc index 8516241..3197e2e 100644 --- a/src/patch/ghost_zone.cc +++ b/src/patch/ghost_zone.cc @@ -175,7 +175,7 @@ interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in, min_ipar_used_(NULL), max_ipar_used_(NULL), other_par_(NULL), interp_result_buffer_(NULL), - Jacobian_oipar_posn_(NULL), Jacobian_buffer_(NULL) // no comma + Jacobian_y_ipar_posn_(NULL), Jacobian_buffer_(NULL) // no comma { // // verify that we have the expected relationships between @@ -330,7 +330,7 @@ other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(), interpatch_ghost_zone::~interpatch_ghost_zone() { delete Jacobian_buffer_; -delete Jacobian_oipar_posn_; +delete Jacobian_y_ipar_posn_; delete interp_result_buffer_; delete other_par_; delete max_ipar_used_; @@ -589,18 +589,18 @@ assert(other_patch_interp_ != NULL); other_patch_interp_->verify_Jacobian_sparsity_pattern_ok(); -other_patch_interp_->molecule_minmax_ipar_m(min_oipar_m_, max_oipar_m_); +other_patch_interp_->molecule_minmax_ipar_m(min_y_ipar_m_, max_y_ipar_m_); -if (Jacobian_oipar_posn_ == NULL) - then Jacobian_oipar_posn_ = new jtutil::array2d<CCTK_INT> +if (Jacobian_y_ipar_posn_ == NULL) + then Jacobian_y_ipar_posn_ = new jtutil::array2d<CCTK_INT> (min_other_iperp_, max_other_iperp_, extreme_min_ipar_, extreme_max_ipar_); -other_patch_interp_->molecule_posn(*Jacobian_oipar_posn_); +other_patch_interp_->molecule_posn(*Jacobian_y_ipar_posn_); if (Jacobian_buffer_ == NULL) then Jacobian_buffer_ = new jtutil::array3d<fp> (min_other_iperp_, max_other_iperp_, extreme_min_ipar_, extreme_max_ipar_, - min_oipar_m_, max_oipar_m_); + min_y_ipar_m_, max_y_ipar_m_); other_patch_interp_->Jacobian(*Jacobian_buffer_); } diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh index 70b7f3b..1cf1517 100644 --- a/src/patch/ghost_zone.hh +++ b/src/patch/ghost_zone.hh @@ -249,6 +249,11 @@ public: // operation into internal buffers; the following functions // provide access to that Jacobian. // + // n.b. terminology is + // partial gridfn at x + // ------------------- + // partial gridfn at y + // // FIXME: should these be moved out into a separate Jacobian // object/class? // @@ -262,46 +267,44 @@ public: // The API in the remaining functions implicitly assumes that // the Jacobian is independent of ghosted_gfn , and also that // the structure of the Jacobian is "simple", in that the set of - // points (with nonzero Jacobian values) in a single row of the + // y points (with nonzero Jacobian values) in a single row of the // Jacobian matrix (i.e. the set of points on which a single // ghost-zone point depends), - // - lies entirely within a single patch - // - has the same number of points for every point in this ghost zone - // - has a single oiperp value (= iperp value in their patch) - // (depending on our iperp, of course) - // - have a contiguous interval of oipar (= ipar in their patch) - // (depending on our iperp and ipar, of course), whose - // size is + // - lies entirely within a single y patch + // - has the same number of points for every x point in this ghost zone + // - has a single yiperp value (depending on our iperp, of course) + // - have a contiguous interval of yipar (depending on our iperp + //and ipar, of course), whose size is // [or can be taken to be without an unreasonable // amount of zero-padding] // independent of our iperp and ipar; we parameterize this - // interval as oipar = posn+m where posn is determined by + // interval as yipar = posn+m where posn is determined by // our iperp and ipar, and m has a fixed range independent // of our iperp and ipar // - // to which patch/edge do the points in this Jacobian row belong? - virtual const patch& Jacobian_patch() const = 0; - virtual const patch_edge& Jacobian_edge() const = 0; + // to which patch/edge do the y points in this Jacobian row belong? + virtual patch& Jacobian_y_patch() const = 0; + virtual const patch_edge& Jacobian_y_edge() const = 0; // what is the [min,max] range of m for this ghost zone? - virtual int Jacobian_min_oipar_m() const = 0; - virtual int Jacobian_max_oipar_m() const = 0; + virtual int Jacobian_y_min_ipar_m() const = 0; + virtual int Jacobian_y_max_ipar_m() const = 0; - // what is the oiperp of the Jacobian points (= iperp in their patch)? - virtual int Jacobian_oiperp(int iperp) const = 0; + // what is the iperp of the Jacobian y points in their (y) patch? + virtual int Jacobian_y_iperp(int x_iperp) const = 0; - // what is the posn value of the points in this Jacobian row? - virtual int Jacobian_oipar_posn(int iperp, int ipar) const = 0; + // what is the posn value of the y points in this Jacobian row? + virtual int Jacobian_y_ipar_posn(int x_iperp, int x_ipar) const = 0; // what is the Jacobian - // partial synchronize gridfn(ghosted_gfn, iperp, ipar) - // ------------------------------------------------------------ - // partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m) + // partial synchronize() px.gridfn(ghosted_gfn, x_iperp, x_ipar) + // ------------------------------------------------------------- + // partial py.gridfn(ghosted_gfn, y_iperp, y_posn+y_ipar_m) // where - // oiperp = Jacobian_oiperp(iperp) - // posn = Jacobian_oipar_posn(iperp, ipar) - virtual fp Jacobian(int iperp, int ipar, int ipar_m) const = 0; + // y_iperp = Jacobian_y_iperp(x_iperp) + // y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar) + virtual fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const = 0; public: @@ -466,6 +469,11 @@ public: // // ***** Jacobian of synchronize() ***** // + // n.b. terminology is + // partial gridfn at x + // ------------------- + // partial gridfn at y + // // allocate internal buffers, compute Jacobian // ... this function is a no-op in this class @@ -475,30 +483,31 @@ public: bool want_max_par_corner) { } - // to which patch/edge do the points in this Jacobian row belong? - const patch& Jacobian_patch() const { return symmetry_patch_; } - const patch_edge& Jacobian_edge() const { return symmetry_edge_; } + // to which patch/edge do the y points in this Jacobian row belong? + patch& Jacobian_y_patch() const { return symmetry_patch_; } + const patch_edge& Jacobian_y_edge() const { return symmetry_edge_; } // what is the [min,max] range of m for this ghost zone? - int Jacobian_min_oipar_m() const { return 0; } - int Jacobian_max_oipar_m() const { return 0; } + int Jacobian_y_min_ipar_m() const { return 0; } + int Jacobian_y_max_ipar_m() const { return 0; } // what is the oiperp of the Jacobian points (= iperp in their patch)? - virtual int Jacobian_oiperp(int iperp) const - { return iperp_map_of_iperp(iperp); } + virtual int Jacobian_y_iperp(int x_iperp) const + { return iperp_map_of_iperp(x_iperp); } // what is the posn value of the points in this Jacobian row? - int Jacobian_oipar_posn(int iperp, int ipar) const - { return ipar_map_of_ipar(ipar); } + int Jacobian_y_ipar_posn(int x_iperp, int x_ipar) const + { return ipar_map_of_ipar(x_ipar); } // what is the Jacobian - // partial synchronize gridfn(ghosted_gfn, iperp, ipar) - // ------------------------------------------------------------ - // partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m) + // partial synchronize() px.gridfn(ghosted_gfn, x_iperp, x_ipar) + // ------------------------------------------------------------- + // partial py.gridfn(ghosted_gfn, y_iperp, y_posn+y_ipar_m) // where - // oiperp = Jacobian_oiperp(iperp) - // posn = Jacobian_oipar_posn(iperp, ipar) - fp Jacobian(int iperp, int ipar, int ipar_m) const { return 1.0; } + // y_iperp = Jacobian_y_iperp(x_iperp) + // y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar) + fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const + { return (y_ipar_m == 0) ? 1.0 : 0.0; } // @@ -506,7 +515,7 @@ public: // // symmetry-map coordinates - const patch& symmetry_patch() const { return symmetry_patch_; } + patch& symmetry_patch() const { return symmetry_patch_; } const patch_edge& symmetry_edge() const { return symmetry_edge_; } int iperp_map_of_iperp(int iperp) const { return iperp_map_->map(iperp); } @@ -544,7 +553,7 @@ private: private: // symmetry mapping --> interior of which patch? which edge? - const patch& symmetry_patch_; + patch& symmetry_patch_; const patch_edge& symmetry_edge_; // symmetry mappings for (iperp,ipar) @@ -602,6 +611,11 @@ public: // // ***** Jacobian of synchronize() ***** // + // n.b. terminology is + // partial gridfn at x + // ------------------- + // partial gridfn at y + // // allocate internal buffers, compute Jacobian // @@ -613,36 +627,43 @@ public: bool want_noncorner, bool want_max_par_corner); - // to which patch/edge do the points in this Jacobian row belong? - const patch& Jacobian_patch() const { return other_patch_; } - const patch_edge& Jacobian_edge() const { return other_edge_; } + // to which patch/edge do the y points in this Jacobian row belong? + patch& Jacobian_y_patch() const { return other_patch_; } + const patch_edge& Jacobian_y_edge() const { return other_edge_; } // what is the [min,max] range of m for this ghost zone? - int Jacobian_min_oipar_m() const { return min_oipar_m_; } - int Jacobian_max_oipar_m() const { return max_oipar_m_; } + int Jacobian_y_min_ipar_m() const { return min_y_ipar_m_; } + int Jacobian_y_max_ipar_m() const { return max_y_ipar_m_; } - // what is the oiperp of the Jacobian points (= iperp in their patch)? - int Jacobian_oiperp(int iperp) const { return other_iperp(iperp); } + // what is the iperp of the Jacobian y points in their (y) patch? + // ... the ipar row of grid points is actually the same, so + // we just have to translate x_iperp to the y patch's coordinates + int Jacobian_y_iperp(int x_iperp) const { return other_iperp(x_iperp); } - // what is the posn value of the points in this Jacobian row? - int Jacobian_oipar_posn(int iperp, int ipar) const + // what is the posn value of the y points in this Jacobian row? + int Jacobian_y_ipar_posn(int x_iperp, int x_ipar) const { - assert(Jacobian_oipar_posn_ != NULL); - const int oiperp = other_iperp(iperp); - return (*Jacobian_oipar_posn_)(oiperp, ipar); + assert(Jacobian_y_ipar_posn_ != NULL); + const int y_iperp = Jacobian_y_iperp(x_iperp); + return (*Jacobian_y_ipar_posn_)(y_iperp, x_ipar); } // what is the Jacobian - // partial synchronize gridfn(ghosted_gfn, iperp, ipar) - // ------------------------------------------------------------ - // partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m) + // partial synchronize() px.gridfn(ghosted_gfn, x_iperp, x_ipar) + // ------------------------------------------------------------- + // partial py.gridfn(ghosted_gfn, y_iperp, y_posn+y_ipar_m) // where - // posn = Jacobian_oipar_posn(iperp, ipar) - fp Jacobian(int iperp, int ipar, int ipar_m) const + // y_iperp = Jacobian_y_iperp(x_iperp) + // y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar) + fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const { assert(Jacobian_buffer_ != NULL); - const int oiperp = other_iperp(iperp); - return (*Jacobian_buffer_)(oiperp, ipar, ipar_m); + 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); + return (*Jacobian_buffer_)(y_iperp, x_ipar, y_ipar_m); + } + else return 0.0; } // @@ -650,10 +671,8 @@ public: // // basic connectivity info - const patch& other_patch() const - { return other_patch_; } - const patch_edge& other_edge() const - { return other_edge_; } + patch& other_patch() const { return other_patch_; } + const patch_edge& other_edge() const { return other_edge_; } // check consistency of this and the other patch's ghost zones // and patch_interp objects @@ -754,13 +773,18 @@ private: // // stuff computed by compute_Jacobian() // - int min_oipar_m_, max_oipar_m_; + // n.b. terminology is + // partial gridfn at x + // ------------------- + // partial gridfn at y + // + int min_y_ipar_m_, max_y_ipar_m_; - // other patch's ipar posn for a Jacobian row + // other patch's y ipar posn for a Jacobian row // ... subscripts are (oiperp, ipar) - jtutil::array2d<CCTK_INT>* Jacobian_oipar_posn_; + jtutil::array2d<CCTK_INT>* Jacobian_y_ipar_posn_; // Jacobian values - // ... subscripts are (oiperp, ipar, oipar_m) + // ... subscripts are (y_iperp, x_ipar, y_ipar_m) jtutil::array3d<fp>* Jacobian_buffer_; }; diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc index 16d89a4..5143f26 100644 --- a/src/patch/test_patch_system.cc +++ b/src/patch/test_patch_system.cc @@ -98,12 +98,18 @@ extern "C" void test_patch_system(CCTK_ARGUMENTS); namespace { +void test_synchronize_Jacobians(patch_system& ps, + int test_gfn, int NP_test_gfn, + fp perturbation_amplitude, + const char Jacobian_file_name[]); + void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones); void setup_fn_rho_sigma(patch_system& ps, int ghosted_gfn); void finite_diff(patch_system& ps, int ghosted_gfn_src, int gfn_dst, int which_derivs); void analytic_derivs(patch_system& ps, int gfn_dst, int which_derivs); + void gridfn_minus(patch_system& ps, int gfn_x, int gfn_y, int gfn_dst); void ghosted_gridfn_minus(patch_system& ps, @@ -169,28 +175,40 @@ if (STRING_EQUAL(which_test, "gridfn")) setup_sym_fn_xyz(ps, test_fn_gfn, true); ps.print_ghosted_gridfn(test_fn_gfn, "test_fn.dat"); } -else if (STRING_EQUAL(which_test, "read-gridfn")) + +else if (STRING_EQUAL(which_test, "read gridfn")) then { ps.read_ghosted_gridfn(test_fn_gfn, "test_fn.dat"); ps.print_ghosted_gridfn(test_fn_gfn, "test_fn2.dat"); } + else if (STRING_EQUAL(which_test, "synchronize")) then { setup_sym_fn_xyz(ps, test_fn_gfn, false); ps.print_ghosted_gridfn(test_fn_gfn, "test_fn_init.dat"); + setup_sym_fn_xyz(ps, test_fn_copy_gfn, true); ps.print_ghosted_gridfn(test_fn_copy_gfn, "test_fn_copy.dat"); + ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn); ps.print_ghosted_gridfn(test_fn_gfn, "test_fn_sync.dat"); + ghosted_gridfn_minus(ps, test_fn_gfn, test_fn_copy_gfn, ghosted_error_gfn); ps.print_ghosted_gridfn(ghosted_error_gfn, "ghosted_error.dat"); } + +else if (STRING_EQUAL(which_test, "synchronize Jacobian")) + then test_synchronize_Jacobians(ps, + test_fn_copy_gfn, test_fn_gfn, + NP_Jacobian__perturbation_amplitude, + Jacobian_file_name); + else if (STRING_EQUAL(which_test, "derivatives")) then { setup_fn_rho_sigma(ps, test_fn_gfn); - ps.print_gridfn(nominal_error_gfn, "test_fn.dat"); + ps.print_gridfn(test_fn_gfn, "test_fn.dat"); finite_diff(ps, test_fn_gfn, FD_derivs_gfn, which_derivs); ps.print_gridfn(FD_derivs_gfn, "FD_derivs.dat"); analytic_derivs(ps, analytic_derivs_gfn, which_derivs); @@ -200,6 +218,7 @@ else if (STRING_EQUAL(which_test, "derivatives")) nominal_error_gfn); ps.print_gridfn(nominal_error_gfn, "nominal_error.dat"); } + else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown which_test=\"%s\"!", which_test); /*NOTREACHED*/ @@ -212,6 +231,173 @@ CCTK_VInfo(CCTK_THORNSTRING, "destroying patch system"); //****************************************************************************** // +// This function tests the computation of the Jacobian of the +// patch_system::synchronize() operation. In outline, it does the following: +// +// set up a test function in test_gfn on the nominal grid +// 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 +// { +// compute the synchronize() Jacobian for this ghost zone +// (via pgz.compute_Jacobian() et al) +// for each point x in (p,pgz) +// { +// for each point y in gz.Jacobian_y_patch() +// { +// 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 +// } +// } +// } +// +// 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. +// Jacobian_file_name = The name of the output file to which both Jacobians +// should be written. +// +namespace { +void test_synchronize_Jacobians(patch_system& ps, + int test_gfn, int NP_test_gfn, + fp perturbation_amplitude, + const char Jacobian_file_name[]) +{ +setup_sym_fn_xyz(ps, test_gfn, false); +ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn); +ps.print_ghosted_gridfn(test_fn_copy_gfn, "test_fn_copy.dat"); + +setup_sym_fn_xyz(ps, NP_test_gfn, false); + +FILE *fileptr = fopen(Jacobian_file_name, "w"); +if (fileptr == NULL) + then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, +"test_synchronize_Jacobians(): can't open output plot file \"%s\"!", + Jacobian_file_name); /*NOTREACHED*/ +fprintf(fileptr, "# column 1 = x patch number\n"); +fprintf(fileptr, "# column 2 = x ghost_zone is_min()\n"); +fprintf(fileptr, "# column 3 = x ghost_zone is_rho()\n"); +fprintf(fileptr, "# column 4 = x_iperp\n"); +fprintf(fileptr, "# column 5 = x_ipar\n"); +fprintf(fileptr, "# column 6 = x_irho\n"); +fprintf(fileptr, "# column 7 = x_isigma\n"); +fprintf(fileptr, "# column 8 = y patch number\n"); +fprintf(fileptr, "# column 9 = y ghost_zone is_min()\n"); +fprintf(fileptr, "# column 10 = y ghost_zone is_rho()\n"); +fprintf(fileptr, "# column 11 = y_iperp\n"); +fprintf(fileptr, "# column 12 = y_ipar\n"); +fprintf(fileptr, "# column 13 = y_irho\n"); +fprintf(fileptr, "# column 14 = y_isigma\n"); +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) + { + patch& p = ps.ith_patch(pn); + + // n.b. these loops must use _int_ variables for the loop + // to terminate! + for (int want_min = false ; want_min <= true ; ++want_min) + { + 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); + + pgz.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(); + + //*** for each point x in (p,pgz) + for (int x_iperp = pgz.min_iperp() ; + x_iperp <= pgz.max_iperp() ; + ++x_iperp) + { + // FIXME FIXME -- this includes corners when it shouldn't + for (int x_ipar = pgz.ghost_zone_min_ipar() ; + x_ipar <= pgz.ghost_zone_max_ipar() ; + ++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); + + 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_ipar_posn + = pgz.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_isigma = q.min_isigma() ; + y_isigma <= q.max_isigma() ; + ++y_isigma) + { + // 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) + += 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) ) + / perturbation_amplitude; + q.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); + + // print the results + if (m_in_molecule || (Jacobian != NP_Jacobian)) + then 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)); + } + } + } + } + } + } + + } + +fclose(fileptr); +} + } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// // This function sets up the ghosted test function for the gridfn and // synchronize tests, symmetrizing the test function to match the // symmetry of the patch system. |