diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-18 17:42:03 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-18 17:42:03 +0000 |
commit | 07784c803efe3a065851d5e0e273cac82c12f0cb (patch) | |
tree | 60bff0605f54ecf0f9906df412947acbb6bf0d5b /src | |
parent | 695c3f5b78149e18eeb0cc74f5ae971f5ff9cd32 (diff) |
move Jacobian.* to ../elliptic/
update other files to compute more of Jacobian coeffs
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@633 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src')
-rw-r--r-- | src/gr/AHFinderDirect.hh | 5 | ||||
-rw-r--r-- | src/gr/Jacobian.hh | 34 | ||||
-rw-r--r-- | src/gr/cg.hh | 4 | ||||
-rw-r--r-- | src/gr/driver.cc | 170 | ||||
-rw-r--r-- | src/gr/horizon_Jacobian.cc | 186 | ||||
-rw-r--r-- | src/gr/horizon_function.cc | 5 |
6 files changed, 177 insertions, 227 deletions
diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index 8bd0225..8e5fc8c 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -72,6 +72,7 @@ void horizon_function(patch_system& ps, bool Jacobian_flag); // horizon_Jacobian.cc -class Jacobian; -void horizon_Jacobian(const patch_system& ps, +Jacobian& create_Jacobian(const patch_system& ps, + const char Jacobian_type[]); +void horizon_Jacobian(patch_system& ps, Jacobian& Jac); diff --git a/src/gr/Jacobian.hh b/src/gr/Jacobian.hh deleted file mode 100644 index e49187a..0000000 --- a/src/gr/Jacobian.hh +++ /dev/null @@ -1,34 +0,0 @@ -// Jacobian.hh -- data structures for the Jacobian matrix -// $Id$ - -// -// prerequisites: -// - -// -// A Jacobian object stores the (a) Jacobian matrix of H with respect to h. -// Jacobian is an abstract base class, from which we derive a separate -// concrete class for each type of (storage scheme for the) Jacobian matrix. -// -// A row/column index of the Jacobian (denoted II/JJ) is a linearization -// of the patch system's (patch,irho,isigma). -// - -class Jacobian - { -public: - // basic meta-info - patch_system& my_patch_system(); - -public: - // access a matrix element - virtual const fp& operator()(int II, int JJ) const; // rvalue - virtual fp& operator()(int II, int JJ) const; // lvalue - - // zero the whole matrix or a single row - virtual void zero(); - virtual void zero_row(int II); - - // add to an element: Jacobian(II,JJ) += delta - virtual void add_to_element(int II, int JJ, fp delta); - }; diff --git a/src/gr/cg.hh b/src/gr/cg.hh index 3095a5a..2688a93 100644 --- a/src/gr/cg.hh +++ b/src/gr/cg.hh @@ -3,7 +3,9 @@ // // This header file defines the "virtual machine" used by machine-generated -// code. +// code. It is "dangerous" in that it #defines macros for all the gridfns, +// which will likely produce syntax errors if they are used outside their +// intentended environment. // // prerequisites // gfn.hh diff --git a/src/gr/driver.cc b/src/gr/driver.cc index 3982be2..0c04dfe 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -3,7 +3,7 @@ // // <<<prototypes for functions local to this file>>> // AHFinderDirect_driver - top-level driver -/// setup_ellipsoid_initial_guess - set up ellipsoid in h gridfn +/// setup_Kerr_KerrSchild_initial_guess - set up Kerr/Kerr-Schild initial guess // #include <stdio.h> @@ -32,6 +32,7 @@ using jtutil::error_exit; #include "../util/patch_interp.hh" #include "../util/ghost_zone.hh" #include "../util/patch_system.hh" +#include "../util/Jacobian.hh" #include "gfn.hh" #include "AHFinderDirect.hh" @@ -43,10 +44,9 @@ using jtutil::error_exit; // namespace { -void setup_ellipsoid_initial_guess - (patch_system& ps, - fp global_center_x, fp global_center_y, fp global_center_z, - fp radius_x, fp radius_y, fp radius_z); +void setup_Kerr_KerrSchild_horizon(patch_system& ps, + fp mass, fp spin, + fp xposn, fp yposn, fp zposn); }; //****************************************************************************** @@ -113,42 +113,22 @@ cgi.coord_delta[2] = cctk_delta_space[2]; cgi.gridfn_dims[0] = cctk_lsh[0]; cgi.gridfn_dims[1] = cctk_lsh[1]; cgi.gridfn_dims[2] = cctk_lsh[2]; -cgi.g_dd_11_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gxx") - ); -cgi.g_dd_12_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gxy") - ); -cgi.g_dd_13_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gxz") - ); -cgi.g_dd_22_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gyy") - ); -cgi.g_dd_23_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gyz") - ); -cgi.g_dd_33_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gzz") - ); -cgi.K_dd_11_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kxx") - ); -cgi.K_dd_12_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kxy") - ); -cgi.K_dd_13_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kxz") - ); -cgi.K_dd_22_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kyy") - ); -cgi.K_dd_23_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kyz") - ); -cgi.K_dd_33_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kzz") - ); +// n.b. The cgi.[gK]_dd_??_data are actually const fp * pointers, +// since we won't modify the 3-D gridfn data! But static_cast<...> +// won't change const modifiers, so we just cast to fp* and let +// the assignment take care of the const part... +cgi.g_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxx")); +cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxy")); +cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxz")); +cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyy")); +cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyz")); +cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gzz")); +cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxx")); +cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxy")); +cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxz")); +cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyy")); +cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyz")); +cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kzz")); // @@ -175,15 +155,14 @@ if (STRING_EQUAL(initial_guess_method, "read from file")) initial_guess__read_from_file__file_name, false); // no ghost zones } -else if (STRING_EQUAL(initial_guess_method, "ellipsoid")) +else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild")) then { - setup_ellipsoid_initial_guess(ps, - initial_guess__ellipsoid__center_global_x, - initial_guess__ellipsoid__center_global_y, - initial_guess__ellipsoid__center_global_z, - initial_guess__ellipsoid__radius_x, - initial_guess__ellipsoid__radius_y, - initial_guess__ellipsoid__radius_z); + setup_Kerr_KerrSchild_horizon(ps, + initial_guess__Kerr_KerrSchild__mass, + initial_guess__Kerr_KerrSchild__spin, + initial_guess__Kerr_KerrSchild__xposn, + initial_guess__Kerr_KerrSchild__yposn, + initial_guess__Kerr_KerrSchild__zposn); ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, true, ghosted_gfns::gfn__h, "h.dat", @@ -204,6 +183,12 @@ if (STRING_EQUAL(method, "horizon")) true, ghosted_gfns::gfn__h, "H.dat"); } +if (STRING_EQUAL(method, "horizon Jacobian")) + then { + Jacobian& Jac = create_Jacobian(ps, Jacobian_type); + horizon_function(ps, cgi, gii, true); + horizon_Jacobian(ps, Jac); + } else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown method=\"%s\"!", method); /*NOTREACHED*/ @@ -212,38 +197,30 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, //****************************************************************************** // -// This function sets up an ellipsoid in the gridfn h, using the -// formulas in "ellipsoid.maple" and the Maple-generated C code in -// "ellipsoid.c": -// -// ellipsoid has center (A,B,C), radius (a,b,c) -// angular coordinate system has center (U,V,W) -// -// direction cosines wrt angular coordinate center are (xcos,ycos,zcos) -// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) +// This function sets up the horizon of a Kerr black hole in Kerr-Schild +// coordinates, on the nominal grid, in the h gridfn. // -// then the equation of the ellipsoid is -// (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 -// ----------------- + ---------------- + ----------------- = 1 -// a^2 b^2 c^2 +// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903, +// and the horizon is worked out on page 13 of my AHFinderDirect notes. // -// to solve this, we introduce intermediate variables -// AU = A - U -// BV = B - V -// CW = C - W +// Arguments: +// (mass,spin) = Describe the Kerr black hole. +// [xyz]_posn = The position of the Kerr black hole. // namespace { -void setup_ellipsoid_initial_guess - (patch_system& ps, - fp global_center_x, fp global_center_y, fp global_center_z, - fp radius_x, fp radius_y, fp radius_z) +void setup_Kerr_KerrSchild_horizon(patch_system& ps, + fp mass, fp spin, + fp xposn, fp yposn, fp zposn) { CCTK_VInfo(CCTK_THORNSTRING, - " h = ellipsoid: global_center=(%g,%g,%g)", - global_center_x, global_center_y, global_center_z); + " setting up Kerr/Kerr-Schild horizon"); CCTK_VInfo(CCTK_THORNSTRING, - " radius=(%g,%g,%g)", - radius_x, radius_y, radius_z); + " mass=%g, spin=%g, posn=(%g,%g,%g)", + double(mass), double(spin), + double(xposn), double(yposn), double(zposn)); + +// horizon in Kerr-Schild is coordinate sphere $r = (1 + \sqrt{1-a^2}) m$ +const fp r = (1.0 + sqrt(1.0 - spin*spin)) * mass; for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { @@ -255,44 +232,19 @@ CCTK_VInfo(CCTK_THORNSTRING, isigma <= p.max_isigma() ; ++isigma) { - const fp rho = p.rho_of_irho(irho); + const fp rho = p.rho_of_irho(irho); const fp sigma = p.sigma_of_isigma(isigma); - fp xcos, ycos, zcos; - p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos); - - // set up variables used by Maple-generated code - const fp AU = global_center_x - ps.origin_x(); - const fp BV = global_center_y - ps.origin_y(); - const fp CW = global_center_z - ps.origin_z(); - const fp a = radius_x; - const fp b = radius_y; - const fp c = radius_z; - - // compute the solutions r_plus and r_minus - fp r_plus, r_minus; - #include "ellipsoid.c" - // exactly one of the solutions (call it r) should be positive - fp r; - if ((r_plus > 0.0) && (r_minus < 0.0)) - then r = r_plus; - else if ((r_plus < 0.0) && (r_minus > 0.0)) - then r = r_minus; - else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "\n" -" expected exactly one r>0 solution, got 0 or 2!\n" -" %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n" -" direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n" -" ==> r_plus=%g r_minus=%g\n" - , - p.name(), irho, isigma, - double(rho), double(sigma), - double(xcos), double(ycos), double(zcos), - double(r_plus), double(r_minus)); - /*NOTREACHED*/ + fp theta, phi; + p.theta_phi_of_rho_sigma(rho,sigma, theta,phi); + fp Kerr_x, Kerr_y, Kerr_z; + p.xyz_of_r_rho_sigma(r,rho,sigma, Kerr_x,Kerr_y,Kerr_z); - // r = horizon radius at this grid point - p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r; + const fp KS_x = Kerr_x - spin*sin(theta)*sin(phi); + const fp KS_y = Kerr_y + spin*sin(theta)*cos(phi); + const fp KS_z = Kerr_z; + p.ghosted_gridfn(ghosted_gfns::gfn__h,irho,isigma) + = jtutil::hypot3(KS_x, KS_y, KS_z); } } } diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index 23390ab..c7d397e 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -2,7 +2,9 @@ // $Id$ // // <<<prototypes for functions local to this file>>> -// horizon_Jacobian - top-level driver +// create_Jacobian - create a Jacobian matrix of the appropriate type +// horizon_Jacobian - top-level driver to compute the Jacobian +/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian // #include <stdio.h> @@ -30,9 +32,9 @@ using jtutil::error_exit; #include "../util/patch_interp.hh" #include "../util/ghost_zone.hh" #include "../util/patch_system.hh" +#include "../util/Jacobian.hh" #include "gfn.hh" -#include "Jacobian.hh" #include "AHFinderDirect.hh" //****************************************************************************** @@ -42,14 +44,29 @@ using jtutil::error_exit; // namespace { -void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp, - int x_irho, int x_isigma, - int xm_irho, int xm_isigma, - fp mol); +void add_ghost_zone_Jacobian(const patch_system& ps, + const patch& xp, const ghost_zone& xmgz, + int x_II, + int xm_irho, int xm_isigma, + fp mol, Jacobian& Jac); } //****************************************************************************** -//****************************************************************************** + +// +// This function decodes the Jacobian_type parameter and creates the +// appropriate type of Jacobian matrix object. +// +Jacobian& create_Jacobian(const patch_system& ps, + const char Jacobian_type[]) +{ +if (STRING_EQUAL(Jacobian_type, "dense matrix")) + then return *new dense_Jacobian(ps); +else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + "unknown Jacobian_type=\"%s\"!", + Jacobian_type); /*NOTREACHED*/ +} + //****************************************************************************** // @@ -65,8 +82,8 @@ void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp, // // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. -void horizon_Jacobian(const patch_system& ps, - Jacobian& Jac); +// +void horizon_Jacobian(patch_system& ps, Jacobian& Jac) { CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); @@ -80,8 +97,8 @@ Jac.zero(); for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho) { for (int x_isigma = xp.min_isigma() ; - x_isigma <= xp.max_isigma() ; - ++x_isigma) + x_isigma <= xp.max_isigma() ; + ++x_isigma) { // // compute the Jacobian row for this grid point, i.e. @@ -96,17 +113,25 @@ Jac.zero(); // Newton's method really wants the latter // + // Jacobian row index + const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma); + // Jacobian coefficients for this point const fp Jacobian_coeff_rho - = xp.gridfn(gfn__partial_H_wrt_partial_d_h_1, x_irho,x_isigma); + = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_1, + x_irho, x_isigma); const fp Jacobian_coeff_sigma - = xp.gridfn(gfn__partial_H_wrt_partial_d_h_2, x_irho,x_isigma); + = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_2, + x_irho, x_isigma); const fp Jacobian_coeff_rho_rho - = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_11, x_irho,x_isigma); + = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_11, + x_irho, x_isigma); const fp Jacobian_coeff_rho_sigma - = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_12, x_irho,x_isigma); + = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_12, + x_irho, x_isigma); const fp Jacobian_coeff_sigma_sigma - = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_22, x_irho,x_isigma); + = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, + x_irho, x_isigma); // partial_rho, partial_rho_rho for (int m_irho = xp.molecule_min_m() ; @@ -115,13 +140,16 @@ Jac.zero(); { const int xm_irho = x_irho + m_irho; const fp Jac_rho = Jacobian_coeff_rho - * xp.partial_rho_coeff(m_rho); + * xp.partial_rho_coeff(m_irho); const fp Jac_rho_rho = Jacobian_coeff_rho_rho - * xp.partial_rho_rho_coeff(m_rho); - add_molecule_point_to_Jacobian(ps, xp, - x_irho, x_isigma, - xm_irho, x_isigma, - Jac_rho + Jac_rho_rho); + * xp.partial_rho_rho_coeff(m_irho); + const fp Jac_sum = Jac_rho + Jac_rho_rho; + if (xp.is_in_nominal_grid(xm_irho, x_isigma)) + then Jac(II, xp,xm_irho,x_isigma) += Jac_sum; + else add_ghost_zone_Jacobian + (ps, xp, xp.minmax_rho_ghost_zone(m_irho < 0), + II, xm_irho, x_isigma, + Jac_sum, Jac); } // partial_sigma, partial_sigma_sigma @@ -131,13 +159,16 @@ Jac.zero(); { const int xm_isigma = x_isigma + m_isigma; const fp Jac_sigma = Jacobian_coeff_sigma - * xp.partial_sigma_coeff(m_sigma); + * xp.partial_sigma_coeff(m_isigma); const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma - * xp.partial_sigma_sigma_coeff(m_sigma); - add_molecule_point_to_Jacobian(ps, xp, - x_irho, x_isigma, - x_irho, xm_isigma, - Jac_sigma + Jac_sigma_sigma); + * xp.partial_sigma_sigma_coeff(m_isigma); + const fp Jac_sum = Jac_sigma + Jac_sigma_sigma; + if (xp.is_in_nominal_grid(x_irho, xm_isigma)) + then Jac(II, xp,x_irho,xm_isigma) += Jac_sum; + else add_ghost_zone_Jacobian + (ps, xp, xp.minmax_sigma_ghost_zone(m_isigma < 0), + II, x_irho, xm_isigma, + Jac_sum, Jac); } // partial_rho_sigma @@ -152,12 +183,16 @@ Jac.zero(); const int xm_irho = x_irho + m_irho; const int xm_isigma = x_isigma + m_isigma; const fp Jac_rho_sigma - = Jacobian_coeff_rho_sigma - * xp.partial_rho_sigma_coeff(m_rho, m_sigma); - add_molecule_point_to_Jacobian(ps, xp, - x_irho, x_isigma, - xm_irho, xm_isigma, - Jac_rho_sigma); + = Jacobian_coeff_rho_sigma + * xp.partial_rho_sigma_coeff(m_irho, m_isigma); + if (xp.is_in_nominal_grid(xm_irho, xm_isigma)) + then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma; + else add_ghost_zone_Jacobian + (ps, xp, + xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0, + xm_irho, xm_isigma), + II, xm_irho, xm_isigma, + Jac_rho_sigma, Jac); } } @@ -169,59 +204,52 @@ Jac.zero(); //****************************************************************************** // -// This function adds all the elements for a single molecule point to -// a Jacobian matrix, including any contributions from other patches if -// the specified point lies in a ghost zone. +// This function adds the ghost-zone Jacobian dependency contributions +// for a single ghost-zone point, to a Jacobian matrix. // // Arguments: // ps = The patch system. // xp = The patch containing the center point of the molecule. -// x_(irho,isigma) = The coordinates of the center point of the molecule. -// xm_(irho,isigma) = The coordinates of the x+m point of the molecule. +// xmgz = If the x+m point is in a ghost zone, this must be that ghost zone. +// If the x+m point is not in a ghost zone, this argument is ignored. +// x_II = The Jacobian row of the x point. +// xm_(irho,isigma) = The coordinates (in xp) of the x+m point of the molecule. // mol = The molecule coefficient. +// Jac = The Jacobian matrix. // namespace { -void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp, - int x_irho, int x_isigma, - int xm_irho, int xm_isigma, - fp mol) +void add_ghost_zone_Jacobian(const patch_system& ps, + const patch& xp, const ghost_zone& xmgz, + int x_II, + int xm_irho, int xm_isigma, + fp mol, Jacobian& Jac) { -// Jacobian row -const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma); - -if (xp.is_in_nominal_grid(xm_irho, xm_isigma)) - then { - const int JJ = ps.gpn_of_patch_irho_isigma(xp, xm_irho, xm_isigma); - Jac.add_to_element(II,JJ, Jac_coeff); - } - else { - const ghost_zone& xmgz - = xp.ghost_zone_containing_point(xm_irho, xm_isigma); - 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); - - // 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); - for (int ym_ipar_m = min_ym_ipar_m ; - ym_ipar_m <= max_ym_ipar_m ; - ++ym_ipar_m) - { - const int ym_ipar = ym_ipar_posn + ym_ipar_m; - 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.add_to_element(II,JJ, mol*sync_Jac); - } +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); + +// 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); + +// add the Jacobian contributions from the ym points + for (int ym_ipar_m = min_ym_ipar_m ; + ym_ipar_m <= max_ym_ipar_m ; + ++ym_ipar_m) + { + const int ym_ipar = ym_ipar_posn + ym_ipar_m; + 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; } } } diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 50e9254..36d5746 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -33,6 +33,7 @@ using jtutil::error_exit; #include "../util/patch_interp.hh" #include "../util/ghost_zone.hh" #include "../util/patch_system.hh" +#include "../util/Jacobian.hh" #include "gfn.hh" #include "AHFinderDirect.hh" @@ -98,8 +99,8 @@ void horizon_function(patch_system& ps, { CCTK_VInfo(CCTK_THORNSTRING, " horizon function"); -// fill in values of ghosted gridfns in ghost zones -ps.synchronize(ghosted_gfns::gfn__h, ghosted_gfns::gfn__h); +// fill in values of all ghosted gridfns in ghost zones +ps.synchronize(); // set up xyz positions of grid points setup_xyz_posns(ps); |