From 695c3f5b78149e18eeb0cc74f5ae971f5ff9cd32 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 18 Jul 2002 17:40:38 +0000 Subject: changes to compute Jacobian coefficients properly git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@632 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/patch/ghost_zone.cc | 1 + src/patch/ghost_zone.hh | 125 +++++++++++++++---------- src/patch/make.code.defn | 3 +- src/patch/patch.cc | 57 +++++++++--- src/patch/patch.hh | 38 +++++++- src/patch/patch_system.cc | 229 +++++++++++++++++++++++----------------------- src/patch/patch_system.hh | 27 +++++- 7 files changed, 292 insertions(+), 188 deletions(-) (limited to 'src') diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc index 72d9738..939d345 100644 --- a/src/patch/ghost_zone.cc +++ b/src/patch/ghost_zone.cc @@ -585,6 +585,7 @@ void interpatch_ghost_zone::compute_Jacobian bool want_min_par_corner = true, bool want_non_corner = true, bool want_max_par_corner = true) + const { verify_caller_wants_all_of_ghost_zone(want_min_par_corner, want_non_corner, diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh index 408c162..5db2c52 100644 --- a/src/patch/ghost_zone.hh +++ b/src/patch/ghost_zone.hh @@ -213,10 +213,8 @@ // // -// N.b. const qualifiers on member functions of ghost_zone and its derived -// classes refer to the union of the class functions and the underlying -// gridfn data. That is, anything which changes the gridfn data is taken -// to be non-const. +// N.b. const qualifiers in ghost_zone and its derived classes refer to +// the underlying gridfn data. // // forward declarations @@ -268,6 +266,7 @@ public: bool want_min_par_corner = true, bool want_noncorner = true, bool want_max_par_corner = true) + const = 0; // @@ -331,18 +330,6 @@ public: bool is_min() const { return my_edge().is_min(); } bool is_rho() const { return my_edge().is_rho(); } - // adjacent ghost zones to our min/max corners - const ghost_zone& min_par_adjacent_ghost_zone() const - { - return my_patch() - .ghost_zone_on_edge( my_edge().min_par_adjacent_edge() ); - } - const ghost_zone& max_par_adjacent_ghost_zone() const - { - return my_patch() - .ghost_zone_on_edge( my_edge().max_par_adjacent_edge() ); - } - // min/max iperp of the ghost zone int min_iperp() const { @@ -373,14 +360,31 @@ public: virtual int min_ipar(int iperp) const = 0; virtual int max_ipar(int iperp) const = 0; - // assert() that ghost zone is fully setup: - // defined here ==> no-op - // symmetry ghost zone ==> unchanged ==> no-op - // interpatch ghost zone ==> check consistency of this and the - // other patch's ghost zones and - // patch_interp objects - virtual void assert_fully_setup() const { } + // point membership predicate + bool is_in_ghost_zone(int iperp, int ipar) + const + { + // n.b. don't test ipar until we're sure iperp is in range! + return (iperp >= min_iperp()) && (iperp <= max_iperp()) + && (ipar >= min_ipar(iperp)) + && (ipar <= max_ipar(iperp)); + } + // adjacent ghost zones to our min/max corners + const ghost_zone& min_par_adjacent_ghost_zone() const + { + return my_patch() + .ghost_zone_on_edge( my_edge().min_par_adjacent_edge() ); + } + const ghost_zone& max_par_adjacent_ghost_zone() const + { + return my_patch() + .ghost_zone_on_edge( my_edge().max_par_adjacent_edge() ); + } + + // + // ***** constructor, finish setup, destructor ***** + // protected: // ... values for is_interpatch_in constructor argument static const bool ghost_zone_is_symmetry = false; @@ -396,6 +400,14 @@ protected: is_interpatch_(is_interpatch_in) { } public: + // assert() that ghost zone is fully setup: + // defined here ==> no-op + // symmetry ghost zone ==> unchanged ==> no-op + // interpatch ghost zone ==> check consistency of this and the + // other patch's ghost zones and + // patch_interp objects + virtual void assert_fully_setup() const { } + // destructor must be virtual to allow destruction // of derived classes via ptr/ref to this class virtual ~ghost_zone() { } @@ -493,6 +505,7 @@ public: bool want_min_par_corner = true, bool want_noncorner = true, bool want_max_par_corner = true) + const { } // to which patch/edge do the y points in this Jacobian row belong? @@ -544,6 +557,11 @@ public: int min_ipar(int iperp) const { return extreme_min_ipar(); } int max_ipar(int iperp) const { return extreme_max_ipar(); } + // + // ***** constructors, destructor ***** + // +public: + // constructor for mirror-symmetry ghost zone symmetry_ghost_zone(const patch_edge& my_edge_in); @@ -637,7 +655,8 @@ public: void compute_Jacobian(int ghosted_min_gfn, int ghosted_max_gfn, bool want_min_par_corner = true, bool want_noncorner = true, - bool want_max_par_corner = true); + bool want_max_par_corner = true) + const; // to which patch/edge do the y points in this Jacobian row belong? patch& Jacobian_y_patch() const { return other_patch_; } @@ -676,10 +695,22 @@ public: return (*Jacobian_buffer_)(y_iperp, x_ipar, y_ipar_m); } +private: + // helper function for synchronize() and compute_Jacobian(): + // verify (no-op if ok, error_exit() if not) that caller wants + // to operate on the entire ghost zone, since our implementations + // of synchronize() and compute_Jacobian() only support this case + void verify_caller_wants_all_of_ghost_zone(bool want_min_par_corner, + bool want_non_corner, + bool want_max_par_corner) + const; + + // // ***** low-level client interface ***** // +public: // basic connectivity info patch& other_patch() const { return other_patch_; } const patch_edge& other_edge() const { return other_edge_; } @@ -688,29 +719,6 @@ public: // and patch_interp objects void assert_fully_setup() const; - // - // constructor, finish setup, destructor - // - interpatch_ghost_zone(const patch_edge& my_edge_in, - const patch_edge& other_edge_in, - int N_overlap_points); - // finish setup (requires adjacent-side ghost_zone objects - // to exist, though not to have finish_setup() called): - // - setup par coordinate mapping information - // - setup interpolator data pointers & result buffer - // - create patch_interp object to interpolate from *other* patch - void finish_setup(int interp_handle, int interp_par_table_handle); - ~interpatch_ghost_zone(); - -private: - // verify (no-op if ok, error_exit() if not) that caller wants - // to operate on the entire ghost zone, since our implementations - // of synchronize() and compute_Jacobian() only support this case - void verify_caller_wants_all_of_ghost_zone(bool want_min_par_corner, - bool want_non_corner, - bool want_max_par_corner) - const; - // min/max ipar of the ghost zone for specified iperp // with possibly "triangular" corners depending on the type // (symmetry vs interpatch) of the adjacent ghost zones @@ -725,6 +733,23 @@ private: return other_iperp_->map(iperp); } + // + // ***** constructor, finish setup, destructor ***** + // +public: + interpatch_ghost_zone(const patch_edge& my_edge_in, + const patch_edge& other_edge_in, + int N_overlap_points); + + // finish setup (requires adjacent-side ghost_zone objects + // to exist, though not to have finish_setup() called): + // - setup par coordinate mapping information + // - setup interpolator data pointers & result buffer + // - create patch_interp object to interpolate from *other* patch + void finish_setup(int interp_handle, int interp_par_table_handle); + + ~interpatch_ghost_zone(); + private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator @@ -786,13 +811,13 @@ private: // ------------------- // partial gridfn at y // - int Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_; + mutable int Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_; // other patch's y ipar posn for a Jacobian row // ... subscripts are (oiperp, ipar) - jtutil::array2d* Jacobian_y_ipar_posn_; + mutable jtutil::array2d* Jacobian_y_ipar_posn_; // Jacobian values // ... subscripts are (y_iperp, x_ipar, y_ipar_m) - jtutil::array3d* Jacobian_buffer_; + mutable jtutil::array3d* Jacobian_buffer_; }; diff --git a/src/patch/make.code.defn b/src/patch/make.code.defn index ad0b9ac..2919722 100644 --- a/src/patch/make.code.defn +++ b/src/patch/make.code.defn @@ -10,7 +10,8 @@ SRCS = coords.cc \ ghost_zone.cc \ patch_interp.cc \ patch_info.cc \ - patch_system.cc + patch_system.cc \ + Jacobian.cc # test_patch_system.cc # Subdirectories containing source files diff --git a/src/patch/patch.cc b/src/patch/patch.cc index dfa4369..1002a60 100644 --- a/src/patch/patch.cc +++ b/src/patch/patch.cc @@ -8,8 +8,8 @@ // x_patch::x_patch // y_patch::y_patch // -// patch::minmax_ang_ghost_zone // patch::ghost_zone_on_edge +// patch::corner_ghost_zone_containing_point // patch::interpatch_ghost_zone_on_edge // patch::create_mirror_symmetry_ghost_zone // patch::create_periodic_symmetry_ghost_zone @@ -156,29 +156,58 @@ y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in, //****************************************************************************** // -// This function returns a reference to the specified ghost zone -// of this patch. +// This function returns a reference to the ghost zone on a specified +// edge, after first assert()ing that the edge belongs to this patch. // -ghost_zone& patch::minmax_ang_ghost_zone(bool want_min, bool want_rho) +// N.b. This function can't be inline in "patch.hh" because it needs +// member functions of class patch_edge, which comes after class patch +// in our #include order. +// +ghost_zone& patch::ghost_zone_on_edge(const patch_edge& e) const { -return want_min ? (want_rho ? min_rho_ghost_zone() - : min_sigma_ghost_zone()) - : (want_rho ? max_rho_ghost_zone() - : max_sigma_ghost_zone()); +assert(& e.my_patch() == this); +return minmax_ang_ghost_zone(e.is_min(), e.is_rho()); } //****************************************************************************** // -// This function returns a reference to the specified ghost zone -// of this patch. +// This function determines which of the two adjacent ghost zones meeting +// at a specified corner, contains a specified point. If the point isn't +// in either ghost zone, an abort(0) is done. // -ghost_zone& patch::ghost_zone_on_edge(const patch_edge &edge) -const +// Arguments: +// {rho,sigma}_is_min = Specify the corner (and implicitly the ghost zones). +// irho,isigma = Specify the point. +// +// Results: +// This function returns (a reference to) the desired ghost zone. +ghost_zone& patch::corner_ghost_zone_containing_point + (bool rho_is_min, bool sigma_is_min, + int irho, int isigma) + const { -assert(& edge.my_patch() == this); -return minmax_ang_ghost_zone(edge.is_min(), edge.is_rho()); +ghost_zone& rho_gz = minmax_rho_ghost_zone( rho_is_min); +ghost_zone& sigma_gz = minmax_sigma_ghost_zone(sigma_is_min); + +const patch_edge& rho_edge = rho_gz.my_edge(); +const patch_edge& sigma_edge = sigma_gz.my_edge(); + +const int rho_iperp = rho_edge.iperp_of_irho_isigma(irho, isigma); +const int rho_ipar = rho_edge. ipar_of_irho_isigma(irho, isigma); +const int sigma_iperp = sigma_edge.iperp_of_irho_isigma(irho, isigma); +const int sigma_ipar = sigma_edge. ipar_of_irho_isigma(irho, isigma); + +const bool is_in_rho_ghost_zone + = rho_gz.is_in_ghost_zone( rho_iperp, rho_ipar); +const bool is_in_sigma_ghost_zone + = sigma_gz.is_in_ghost_zone(sigma_iperp, sigma_ipar); + +// check that point is in exactly one ghost zone +assert(is_in_rho_ghost_zone ^ is_in_sigma_ghost_zone); + +return is_in_rho_ghost_zone ? rho_gz : sigma_gz; } //****************************************************************************** diff --git a/src/patch/patch.hh b/src/patch/patch.hh index b58882d..6325dfe 100644 --- a/src/patch/patch.hh +++ b/src/patch/patch.hh @@ -312,11 +312,41 @@ public: assert(max_sigma_ghost_zone_ != NULL); return *max_sigma_ghost_zone_; } - ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho) const; - ghost_zone& ghost_zone_on_edge(const patch_edge &edge) const; - // get ghost zone on specified edge, verify that it is indeed - // interpatch, static_cast<> to interpatch_ghost_zone ptr + // ... these fns are defined explicitly (unlike for grid:: stuff) + // because they're actually used by Jacobian code in ../gr/ + ghost_zone& minmax_rho_ghost_zone(bool want_min) + const + { + return want_min ? min_rho_ghost_zone() + : max_rho_ghost_zone(); + } + ghost_zone& minmax_sigma_ghost_zone(bool want_min) + const + { + return want_min ? min_sigma_ghost_zone() + : max_sigma_ghost_zone(); + } + + ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho) + const + { + return want_rho ? minmax_rho_ghost_zone(want_min) + : minmax_sigma_ghost_zone(want_min); + } + + ghost_zone& ghost_zone_on_edge(const patch_edge &e) const; + + // which of the two ghost zones at a specifieid corner, + // contains a specified point? + ghost_zone& corner_ghost_zone_containing_point + (bool rho_is_min, bool sigma_is_min, // specifies corner + int irho, int isigma) // specifies point + const; + + // get ghost zone on specified edge, + // assert() that it is indeed interpatch, + // static_cast<> to interpatch_ghost_zone ptr interpatch_ghost_zone& interpatch_ghost_zone_on_edge (const patch_edge &e) const; diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index ec80133..dcc54c4 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -1151,6 +1151,7 @@ void patch_system::synchronize(int ghosted_min_gfn_to_sync, // void patch_system::compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync, int ghosted_max_gfn_to_sync) + const { for (int pn = 0 ; pn < N_patches() ; ++pn) { @@ -1227,89 +1228,89 @@ return xgz.Jacobian(x_iperp, x_ipar, y_ipar_m); // specified (unknown-grid) radius gridfn. // void patch_system::print_unknown_gridfn -(bool ghosted_flag, int unknown_gfn, - bool print_xyz_flag, bool radius_is_ghosted_flag, - int unknown_radius_gfn, - const char output_file_name[], bool want_ghost_zones) -const + (bool ghosted_flag, int unknown_gfn, + bool print_xyz_flag, bool radius_is_ghosted_flag, + int unknown_radius_gfn, + const char output_file_name[], bool want_ghost_zones) + const { if (want_ghost_zones && !ghosted_flag) -then error_exit(PANIC_EXIT, + then error_exit(PANIC_EXIT, "***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n" " can't have want_ghost_zones && !ghosted_flag !\n" , - unknown_gfn); /*NOTREACHED*/ + unknown_gfn); /*NOTREACHED*/ if (want_ghost_zones && print_xyz_flag && !radius_is_ghosted_flag) -then error_exit(PANIC_EXIT, + then error_exit(PANIC_EXIT, "***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n" " can't have want_ghost_zones && print_xyz_flag\n" " && !radius_is_ghosted_flag!\n" " unknown_radius_gfn=%d\n" , - unknown_gfn, - unknown_radius_gfn); /*NOTREACHED*/ + unknown_gfn, + unknown_radius_gfn); /*NOTREACHED*/ FILE *output_fp = std::fopen(output_file_name, "w"); if (output_fp == NULL) -then error_exit(ERROR_EXIT, + then error_exit(ERROR_EXIT, "***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n" " can't open output file \"%s\"\n!" , - unknown_gfn, - output_file_name); /*NOTREACHED*/ - -for (int pn = 0 ; pn < N_patches() ; ++pn) -{ -const patch& p = ith_patch(pn); - -fprintf(output_fp, "### %s patch\n", p.name()); -fprintf(output_fp, "# %s_gfn=%d\n", - (ghosted_flag ? "ghosted" : "nominal"), unknown_gfn); -fprintf(output_fp, "# dpx = %s\n", p.name_of_dpx()); -fprintf(output_fp, "# dpy = %s\n", p.name_of_dpy()); -fprintf(output_fp, "#\n"); -fprintf(output_fp, - print_xyz_flag - ? "# dpx\tdpy\tgridfn\tglobal_x\tglobal_y\tglobal_z\n" - : "# dpx\tdpy\tgridfn\n"); + unknown_gfn, + output_file_name); /*NOTREACHED*/ - for (int irho = p.effective_min_irho(want_ghost_zones) ; - irho <= p.effective_max_irho(want_ghost_zones) ; - ++irho) - { - for (int isigma = p.effective_min_isigma(want_ghost_zones) ; - isigma <= p.effective_max_isigma(want_ghost_zones) ; - ++isigma) + for (int pn = 0 ; pn < N_patches() ; ++pn) { - const fp rho = p.rho_of_irho(irho); - const fp sigma = p.sigma_of_isigma(isigma); - const fp dpx = p.dpx_of_rho_sigma(rho, sigma); - const fp dpy = p.dpy_of_rho_sigma(rho, sigma); + const patch& p = ith_patch(pn); + + fprintf(output_fp, "### %s patch\n", p.name()); + fprintf(output_fp, "# %s_gfn=%d\n", + (ghosted_flag ? "ghosted" : "nominal"), unknown_gfn); + fprintf(output_fp, "# dpx = %s\n", p.name_of_dpx()); + fprintf(output_fp, "# dpy = %s\n", p.name_of_dpy()); + fprintf(output_fp, "#\n"); fprintf(output_fp, - "%g\t%g\t%#.15g", - dpx, dpy, p.unknown_gridfn(ghosted_flag, - unknown_gfn, irho,isigma)); - if (print_xyz_flag) - then { - const fp r = p.unknown_gridfn(radius_is_ghosted_flag, - unknown_radius_gfn, - irho,isigma); - fp local_x, local_y, local_z; - p.xyz_of_r_rho_sigma(r,rho,sigma, - local_x,local_y,local_z); - const fp global_x = global_x_of_local_x(local_x); - const fp global_y = global_y_of_local_y(local_y); - const fp global_z = global_z_of_local_z(local_z); + print_xyz_flag + ? "# dpx\tdpy\tgridfn\tglobal_x\tglobal_y\tglobal_z\n" + : "# dpx\tdpy\tgridfn\n"); + + for (int irho = p.effective_min_irho(want_ghost_zones) ; + irho <= p.effective_max_irho(want_ghost_zones) ; + ++irho) + { + for (int isigma = p.effective_min_isigma(want_ghost_zones) ; + isigma <= p.effective_max_isigma(want_ghost_zones) ; + ++isigma) + { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + const fp dpx = p.dpx_of_rho_sigma(rho, sigma); + const fp dpy = p.dpy_of_rho_sigma(rho, sigma); fprintf(output_fp, - "\t%#g\t%#g\t%#g", - global_x, global_y, global_z); + "%g\t%g\t%#.15g", + dpx, dpy, p.unknown_gridfn(ghosted_flag, + unknown_gfn, irho,isigma)); + if (print_xyz_flag) + then { + const fp r = p.unknown_gridfn(radius_is_ghosted_flag, + unknown_radius_gfn, + irho,isigma); + fp local_x, local_y, local_z; + p.xyz_of_r_rho_sigma(r,rho,sigma, + local_x,local_y,local_z); + const fp global_x = global_x_of_local_x(local_x); + const fp global_y = global_y_of_local_y(local_y); + const fp global_z = global_z_of_local_z(local_z); + fprintf(output_fp, + "\t%#g\t%#g\t%#g", + global_x, global_y, global_z); + } + fprintf(output_fp, "\n"); + } + fprintf(output_fp, "\n"); } fprintf(output_fp, "\n"); } - fprintf(output_fp, "\n"); - } -fprintf(output_fp, "\n"); -} std::fclose(output_fp); } @@ -1329,89 +1330,89 @@ void patch_system::read_unknown_gridfn(bool ghosted_flag, int unknown_gfn, bool want_ghost_zones) { if (want_ghost_zones && !ghosted_flag) -then error_exit(PANIC_EXIT, + then error_exit(PANIC_EXIT, "***** patch_system::read_unknown_gridfn(unknown_gfn=%d):\n" " can't have want_ghost_zones && !ghosted_flag !\n" , - unknown_gfn); /*NOTREACHED*/ + unknown_gfn); /*NOTREACHED*/ FILE *input_fp = std::fopen(input_file_name, "r"); if (input_fp == NULL) -then error_exit(ERROR_EXIT, + then error_exit(ERROR_EXIT, "***** patch_system::read_unknown_gridfn(unknown_gfn=%d):\n" " can't open input file \"%s\"\n!" , - unknown_gfn, - input_file_name); /*NOTREACHED*/ + unknown_gfn, + input_file_name); /*NOTREACHED*/ int line_number = 1; -for (int pn = 0 ; pn < N_patches() ; ++pn) -{ -patch& p = ith_patch(pn); + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); -for (int irho = p.effective_min_irho(want_ghost_zones) ; - irho <= p.effective_max_irho(want_ghost_zones) ; - ++irho) -{ -for (int isigma = p.effective_min_isigma(want_ghost_zones) ; - isigma <= p.effective_max_isigma(want_ghost_zones) ; - ++isigma) -{ -const fp rho = p.rho_of_irho(irho); -const fp sigma = p.sigma_of_isigma(isigma); -const fp dpx = p.dpx_of_rho_sigma(rho, sigma); -const fp dpy = p.dpy_of_rho_sigma(rho, sigma); - -const int N_buffer = 100; -char buffer[N_buffer]; -// read/discard comments and blank lines - do + for (int irho = p.effective_min_irho(want_ghost_zones) ; + irho <= p.effective_max_irho(want_ghost_zones) ; + ++irho) { - if (std::fgets(buffer, N_buffer, input_fp) == NULL) - then error_exit(ERROR_EXIT, + for (int isigma = p.effective_min_isigma(want_ghost_zones) ; + isigma <= p.effective_max_isigma(want_ghost_zones) ; + ++isigma) + { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + const fp dpx = p.dpx_of_rho_sigma(rho, sigma); + const fp dpy = p.dpy_of_rho_sigma(rho, sigma); + + const int N_buffer = 100; + char buffer[N_buffer]; + // read/discard comments and blank lines + do + { + if (std::fgets(buffer, N_buffer, input_fp) == NULL) + then error_exit(ERROR_EXIT, "***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n" " I/O error or unexpected end-of-file on input!\n" " at irho=%d of [%d,%d], isigma=%d of [%d,%d]\n" " dpx=%g dpy=%g\n" , - p.name(), unknown_gfn, - irho, p.effective_min_irho(want_ghost_zones), - p.effective_max_irho(want_ghost_zones), - isigma, - p.effective_min_isigma(want_ghost_zones), - p.effective_max_isigma(want_ghost_zones), - dpx, dpy); /*NOTREACHED*/ - ++line_number; - } while ((buffer[0] == '#') || (buffer[0] == '\n')); - -double read_dpx, read_dpy, read_gridfn_value; -if (sscanf(buffer, "%lf %lf %lf", - &read_dpx, &read_dpy, &read_gridfn_value) != 3) - then error_exit(ERROR_EXIT, + p.name(), unknown_gfn, + irho, p.effective_min_irho(want_ghost_zones), + p.effective_max_irho(want_ghost_zones), + isigma, + p.effective_min_isigma(want_ghost_zones), + p.effective_max_isigma(want_ghost_zones), + dpx, dpy); /*NOTREACHED*/ + ++line_number; + } while ((buffer[0] == '#') || (buffer[0] == '\n')); + + double read_dpx, read_dpy, read_gridfn_value; + if (sscanf(buffer, "%lf %lf %lf", + &read_dpx, &read_dpy, &read_gridfn_value) != 3) + then error_exit(ERROR_EXIT, "***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n" " bad input data at input line %d!\n" , - p.name(), unknown_gfn, - line_number); /*NOTREACHED*/ -if (! ( jtutil::fuzzy::EQ(read_dpx,dpx) - && jtutil::fuzzy::EQ(read_dpy,dpy) ) ) - then error_exit(ERROR_EXIT, + p.name(), unknown_gfn, + line_number); /*NOTREACHED*/ + if (! ( jtutil::fuzzy::EQ(read_dpx,dpx) + && jtutil::fuzzy::EQ(read_dpy,dpy) ) ) + then error_exit(ERROR_EXIT, "***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n" " wrong (dpx,dpy) at input line %d!\n" " expected (%g,%g)\n" " read (%g,%g)\n" , - p.name(), unknown_gfn, - line_number, - dpx, dpy, - read_dpx, read_dpy); /*NOTREACHED*/ + p.name(), unknown_gfn, + line_number, + dpx, dpy, + read_dpx, read_dpy); /*NOTREACHED*/ -p.unknown_gridfn(ghosted_flag, - unknown_gfn, irho,isigma) = read_gridfn_value; -} -} + p.unknown_gridfn(ghosted_flag, + unknown_gfn, irho,isigma) = read_gridfn_value; + } + } -} + } std::fclose(input_fp); } diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh index a4777f0..a1b4227 100644 --- a/src/patch/patch_system.hh +++ b/src/patch/patch_system.hh @@ -27,7 +27,8 @@ //****************************************************************************** // -// a patch_system object describes a system of interlinked patches. +// A patch_system object describes a system of interlinked patches. +// Its const qualifiers refer to the gridfn data. // class patch_system @@ -245,13 +246,26 @@ public: void synchronize(int ghosted_min_gfn_to_sync, int ghosted_max_gfn_to_sync); + // ... do this for all ghosted gridfns + void synchronize() + { synchronize(ghosted_min_gfn(), ghosted_max_gfn()); } + // // compute Jacobian of synchronize() into internal buffers, // taking into account synchronize()'s full 3-phase algorithm // void compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync, - int ghosted_max_gfn_to_sync); + int ghosted_max_gfn_to_sync) + const; + + // ... do this for all ghosted gridfns + void compute_synchronize_Jacobian() + const + { + compute_synchronize_Jacobian(ghosted_min_gfn(), + ghosted_max_gfn()); + } // // The following functions access the Jacobian computed by @@ -266,7 +280,8 @@ public: const { return xgz.Jacobian_y_patch(); } const patch_edge& synchronize_Jacobian_y_edge (const ghost_zone& xgz) - const { return xgz.Jacobian_y_edge(); } + const + { return xgz.Jacobian_y_edge(); } // what is the [min,max] range of m for a Jacobian row? int synchronize_Jacobian_min_y_ipar_m(const ghost_zone& xgz) @@ -283,7 +298,8 @@ public: // what is the posn value of the y points in this Jacobian row? int synchronize_Jacobian_y_ipar_posn(const ghost_zone& xgz, - int x_iperp, int x_ipar) const; + int x_iperp, int x_ipar) + const; // what is the Jacobian // partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar) @@ -294,7 +310,8 @@ public: // y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar) // taking into account synchronize()'s full 3-phase algorithm fp synchronize_Jacobian(const ghost_zone& xgz, int x_iperp, int x_ipar, - int y_ipar_m) const; + int y_ipar_m) + const; // -- cgit v1.2.3