diff options
| -rw-r--r-- | src/gr/horizon_Jacobian.cc | 35 | ||||
| -rw-r--r-- | src/patch/patch_system.cc | 123 | ||||
| -rw-r--r-- | src/patch/patch_system.hh | 81 |
3 files changed, 185 insertions, 54 deletions
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index d4d11f9..0dd7a31 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -269,17 +269,11 @@ CCTK_VInfo(CCTK_THORNSTRING, ++x_isigma) { // - // compute the Jacobian row for this grid point, i.e. + // compute the main Jacobian terms for this grid point, i.e. // partial H(this point x, Jacobian row II) // --------------------------------------------- // partial h(other points y, Jacobian column JJ) // - // FIXME FIXME: we still have to take into account the - // position-dependence of the coefficients, - // cf the difference between J[3H(h)] and J[2H(h)]; - // here we're sort of computing the former, but - // Newton's method really wants the latter - // // Jacobian row index const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma); @@ -397,15 +391,24 @@ const patch_edge& xme = xmgz.my_edge(); const int xm_iperp = xme.iperp_of_irho_isigma(xm_irho, xm_isigma); const int xm_ipar = xme. ipar_of_irho_isigma(xm_irho, xm_isigma); +// FIXME: this won't change from one call to another +// ==> it would be more efficient to reuse the same buffer +// across multiple calls on this function +int global_min_m, global_max_m; +ps.synchronize_Jacobian_global_minmax_y_ipar_m(global_min_m, global_max_m); +jtutil::array1d<fp> Jacobian_buffer(global_min_m, global_max_m); + // on what other points ym does this molecule point xm depend // via the patch_system::synchronize() operation? -const patch& ymp = ps.synchronize_Jacobian_y_patch(xmgz); -const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz); -const int min_ym_ipar_m = ps.synchronize_Jacobian_min_y_ipar_m(xmgz); -const int max_ym_ipar_m = ps.synchronize_Jacobian_max_y_ipar_m(xmgz); -const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp); -const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn - (xmgz, xm_iperp, xm_ipar); +const patch& ymp = ps.synchronize_Jacobian_y_patch(xmgz, xm_iperp,xm_ipar); +const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz, xm_iperp,xm_ipar); +const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp,xm_ipar); +const int ym_ipar_posn + = ps.synchronize_Jacobian_y_ipar_posn(xmgz, xm_iperp,xm_ipar); +int min_ym_ipar_m, max_ym_ipar_m; +ps.synchronize_Jacobian_minmax_y_ipar_m(xmgz, xm_iperp, xm_ipar, + min_ym_ipar_m, max_ym_ipar_m); +ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar, Jacobian_buffer); // add the Jacobian contributions from the ym points for (int ym_ipar_m = min_ym_ipar_m ; @@ -416,9 +419,7 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn const int ym_irho = yme. irho_of_iperp_ipar(ym_iperp,ym_ipar); const int ym_isigma = yme.isigma_of_iperp_ipar(ym_iperp,ym_ipar); const int JJ = ps.gpn_of_patch_irho_isigma(ymp, ym_irho, ym_isigma); - const fp sync_Jac = ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar, - ym_ipar_m); - Jac(x_II,JJ) += mol*sync_Jac; + Jac(x_II,JJ) += mol*Jacobian_buffer(ym_ipar_m); } } } diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index e03e8bd..4f409f0 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -30,7 +30,12 @@ // // patch_system::synchronize // patch_system::compute_synchronize_Jacobian +// patch_system::synchronize_Jacobian_global_minmax_y_ipar_m +// patch_system::synchronize_Jacobian_y_patch +// patch_system::synchronize_Jacobian_y_edge +// patch_system::synchronize_Jacobian_y_iperp // patch_system::synchronize_Jacobian_y_ipar_posn +// patch_system::synchronize_Jacobian_minmax_y_ipar_m // patch_system::synchronize_Jacobian // @@ -39,6 +44,7 @@ #include <cstring> #include <vector> #include <assert.h> +#include <limits.h> using std::printf; using std::strcmp; @@ -1405,6 +1411,8 @@ void patch_system::synchronize(int ghosted_min_gfn_to_sync, // This function computes the Jacobian of synchronize() into internal // buffers, taking into account synchronize()'s full 3-phase algorithm. // +// It also properly sets global_{min,max}_m_ . +// // FIXME FIXME: at the moment we ignore the 3-phase algorithm and just // pass all the Jacobian operations down to the ghost zones // @@ -1412,6 +1420,8 @@ void patch_system::compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync, int ghosted_max_gfn_to_sync) const { +global_min_m_ = +INT_MAX; +global_max_m_ = -INT_MAX; for (int pn = 0 ; pn < N_patches() ; ++pn) { patch& p = ith_patch(pn); @@ -1424,6 +1434,11 @@ void patch_system::compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync, ghost_zone& gz = p.minmax_ang_ghost_zone(want_min, want_rho); gz.compute_Jacobian(ghosted_min_gfn_to_sync, ghosted_max_gfn_to_sync); + + global_min_m_ = jtutil::min(global_min_m_, + gz.Jacobian_min_y_ipar_m()); + global_max_m_ = jtutil::max(global_max_m_, + gz.Jacobian_max_y_ipar_m()); } } } @@ -1433,6 +1448,71 @@ void patch_system::compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync, // // Given that compute_synchronize_Jacobian() has been called, this +// function computes the global min/max m over all ghost zone points. +// This is useful for sizing the buffer for synchronize_Jacobian(). +// +void patch_system::synchronize_Jacobian_global_minmax_y_ipar_m + (int& min_m, int& max_m) + const +{ +min_m = global_min_m_; +max_m = global_max_m_; +} + +//****************************************************************************** + +// +// Given that compute_synchronize_Jacobian() has been called, this +// function computes the patch containing the y points in a Jacobian row. +// +// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just +// pass the operation down to the ghost zone +// +patch& patch_system::synchronize_Jacobian_y_patch(const ghost_zone& xgz, + int x_iperp, int x_ipar) + const +{ +return xgz.Jacobian_y_patch(); +} + +//****************************************************************************** + +// +// Given that compute_synchronize_Jacobian() has been called, this +// function computes the patch_edge containing the y points in a +// Jacobian row. +// +// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just +// pass the operation down to the ghost zone +// +const patch_edge& patch_system::synchronize_Jacobian_y_edge + (const ghost_zone& xgz, + int x_iperp, int x_ipar) + const +{ +return xgz.Jacobian_y_edge(); +} + +//****************************************************************************** + +// +// Given that compute_synchronize_Jacobian() has been called, this +// function computes the iperp value of the y points in a Jacobian row. +// +// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just +// pass the operation down to the ghost zones +// +int patch_system::synchronize_Jacobian_y_iperp(const ghost_zone& xgz, + int x_iperp, int x_ipar) + const +{ +return xgz.Jacobian_y_iperp(x_iperp); +} + +//****************************************************************************** + +// +// Given that compute_synchronize_Jacobian() has been called, this // function computes the posn value of the y points in a Jacobian row. // // FIXME FIXME: at the moment we ignore the 3-phase algorithm and just @@ -1449,22 +1529,47 @@ return xgz.Jacobian_y_ipar_posn(x_iperp, x_ipar); // // Given that compute_synchronize_Jacobian() has been called, this -// function computes the Jacobian +// function computes the min/max m values of the y points in a Jacobian row. +// +// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just +// pass all the operation down to the ghost zones +// +void patch_system::synchronize_Jacobian_minmax_y_ipar_m(const ghost_zone& xgz, + int x_iperp, int x_ipar, + int& min_m, int& max_m) + const +{ +min_m = xgz.Jacobian_min_y_ipar_m(); +max_m = xgz.Jacobian_max_y_ipar_m(); +} + +//****************************************************************************** + +// +// Given that compute_synchronize_Jacobian() has been called, this +// function stores the Jacobian // partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar) // ------------------------------------------------------------- // partial gridfn(ghosted_gfn, py, y_iperp, y_posn+y_ipar_m) -// where -// y_iperp = synchronize_Jacobian_y_iperp(xgz, x_iperp) -// y_posn = synchronize_Jacobian_y_ipar_posn(xgz, x_iperp, x_ipar) -// taking into account synchronize()'s full 3-phase algorithm +// (taking into account synchronize()'s full 3-phase algorithm) +// in the caller-supplied buffer +// Jacobian_buffer(m) +// for each m , where +// y_iperp = Jacobian_y_iperp(x_iperp) +// y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar) // // FIXME FIXME: at the moment we ignore the 3-phase algorithm and just // pass all the operation down to the ghost zones // -fp patch_system::synchronize_Jacobian - (const ghost_zone& xgz, int x_iperp, int x_ipar, - int y_ipar_m) +void patch_system::synchronize_Jacobian(const ghost_zone& xgz, + int x_iperp, int x_ipar, + jtutil::array1d<fp>& Jacobian_buffer) const { -return xgz.Jacobian(x_iperp, x_ipar, y_ipar_m); +const int min_m = xgz.Jacobian_min_y_ipar_m(); +const int max_m = xgz.Jacobian_max_y_ipar_m(); + for (int m = min_m ; m <= max_m ; ++m) + { + Jacobian_buffer(m) = xgz.Jacobian(x_iperp, x_ipar, m); + } } diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh index d7b39bd..44f8304 100644 --- a/src/patch/patch_system.hh +++ b/src/patch/patch_system.hh @@ -28,7 +28,12 @@ // // A patch_system object describes a system of interlinked patches. -// Its const qualifiers refer to the gridfn data. +// +// Its const qualifiers refer (only) to the gridfn data. Notably, this +// means that synchronize() is a non-const function (it modifies gridfn +// data), while synchronize_Jacobian() et al are const functions (they +// don't modify gridfn data) even though they may update other internal +// state in the patch_system object and its subobjects. // class patch_system @@ -158,48 +163,65 @@ public: // // The following functions access the Jacobian computed by - // compute_synchronize_Jacobian() . Note this API has the - // same implicit assumptions on the Jacobian structure documented - // in the comments in "ghost_zone.hh" immediately following - // ghost_zone::compute_Jacobian() . + // compute_synchronize_Jacobian() . Note this API is rather + // different than that of ghost_zone::comute_Jacobian() et al: + // here we must take into account synchronize()'s full 3-phase + // algorithm, and this may lead to a more general Jacobian + // structure. + // + // This API still implicitly assumes that the Jacobian is + // independent of ghosted_gfn , and that the set of 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 y patch + // - has a single yiperp value + // - have a contiguous interval of yipar; we parameterize this + // interval as yipar = posn+m // - // to which patch/edge do the y points in a Jacobian row belong? - patch& synchronize_Jacobian_y_patch(const ghost_zone& xgz) - const - { return xgz.Jacobian_y_patch(); } - const patch_edge& synchronize_Jacobian_y_edge (const ghost_zone& xgz) - const - { return xgz.Jacobian_y_edge(); } + // what are the global min/max m over all ghost zone points? + // (this is useful for sizing the buffer for synchronize_Jacobian()) + void synchronize_Jacobian_global_minmax_y_ipar_m + (int& min_m, int& max_m) + const; - // what is the [min,max] range of m for a Jacobian row? - int synchronize_Jacobian_min_y_ipar_m(const ghost_zone& xgz) - const - { return xgz.Jacobian_min_y_ipar_m(); } - int synchronize_Jacobian_max_y_ipar_m(const ghost_zone& xgz) - const - { return xgz.Jacobian_max_y_ipar_m(); } + // to which patch/edge do the y points in a Jacobian row belong? + patch& synchronize_Jacobian_y_patch(const ghost_zone& xgz, + int x_iperp, int x_ipar) + const; + const patch_edge& synchronize_Jacobian_y_edge(const ghost_zone& xgz, + int x_iperp, int x_ipar) + const; // what is the iperp of the Jacobian y points in their (y) patch? - int synchronize_Jacobian_y_iperp(const ghost_zone& xgz, int x_iperp) - const - { return xgz.Jacobian_y_iperp(x_iperp); } + int synchronize_Jacobian_y_iperp(const ghost_zone& xgz, + int x_iperp, int x_ipar) + const; - // what is the posn value of the y points in this Jacobian row? + // what are the posn and min/max m for a given x point? int synchronize_Jacobian_y_ipar_posn(const ghost_zone& xgz, int x_iperp, int x_ipar) const; + void synchronize_Jacobian_minmax_y_ipar_m(const ghost_zone& xgz, + int x_iperp, int x_ipar, + int& min_m, int& max_m) + const; - // what is the Jacobian +public: + // store the Jacobian // partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar) // ------------------------------------------------------------- // partial gridfn(ghosted_gfn, py, y_iperp, y_posn+y_ipar_m) - // where + // (taking into account synchronize()'s full 3-phase algorithm) + // in the caller-supplied buffer + // Jacobian_buffer(m) + // for each m , where // y_iperp = Jacobian_y_iperp(x_iperp) // 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) + void synchronize_Jacobian(const ghost_zone& xgz, + int x_iperp, int x_ipar, + jtutil::array1d<fp>& Jacobian_buffer) const; @@ -440,4 +462,7 @@ private: // ... patches point into these, but we own the storage blocks fp* gridfn_storage_; fp* ghosted_gridfn_storage_; + + // min/max m over all ghost zone points + mutable int global_min_m_, global_max_m_; }; |