diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-16 09:26:48 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-16 09:26:48 +0000 |
commit | c0fea6ff486fe45a24f26df8cbb29926a434b500 (patch) | |
tree | bdda0af1cb67253215502a77c3a22d0c753ad92f /src/gr/horizon_Jacobian.cc | |
parent | 35d8115e452c8ff4171aff23ca86856dd9208331 (diff) |
more code to compute the Jacobian matrix
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@622 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src/gr/horizon_Jacobian.cc')
-rw-r--r-- | src/gr/horizon_Jacobian.cc | 161 |
1 files changed, 114 insertions, 47 deletions
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index 40fada6..39db258 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -43,7 +43,11 @@ 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); + } //****************************************************************************** //****************************************************************************** @@ -67,7 +71,8 @@ void horizon_Jacobian(const patch_system& ps, { CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); -Jac.zero(II); +ps.compute_synchronize_Jacobian(); +Jac.zero(); for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) { @@ -92,9 +97,6 @@ Jac.zero(II); // Newton's method really wants the latter // - // Jacobian row index for this point - 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); @@ -112,50 +114,115 @@ Jac.zero(II); m_irho <= xp.molecule_max_m() ; ++m_irho) { - const fp Jac_rho = Jacobian_coeff_rho * xp.partial_rho_coeff(m_rho); + const int xm_irho = x_irho + m_irho; + const fp Jac_rho = Jacobian_coeff_rho + * xp.partial_rho_coeff(m_rho); const fp Jac_rho_rho = Jacobian_coeff_rho_rho * xp.partial_rho_rho_coeff(m_rho); - const fp Jac_sum = Jac_rho + Jac_rho_rho; - if (xp.is_in_nominal_grid(x_irho+m_irho, x_isigma)) - then { - const int JJ - = ps.gpn_of_patch_irho_isigma(xp, x_irho+m_irho, - x_isigma); - Jac.add_to_element(II,JJ, Jac_sum); - } - else { - const patch_edge& xe = (m_irho < 0) - ? xp.min_rho_patch_edge() - : xp.max_rho_patch_edge() - const ghost_zone& gz = xp.ghost_zone_on_edge(xe); - const int x_iperp = xe.iperp_of_irho_isigma(x_irho, - x_isigma); - const int x_ipar = xe.ipar_of_irho_isigma(x_irho, x_isigma); - gz.compute_Jacobian(gfn__h, gfn__h); - const patch& yp = gz.Jacobian_y_patch(); - const patch_edge& ye = gz.Jacobian_y_edge(); - const int min_y_ipar_m = gz.Jacobian_min_y_ipar_m(); - const int max_y_ipar_m = gz.Jacobian_max_y_ipar_m(); - const int y_iperp = gz.Jacobian_y_iperp(x_iperp); - const int y_ipar_posn - = gz.Jacobian_y_ipar_posn(x_iperp, x_ipar); - for (int y_ipar_m = min_y_ipar_m ; - y_ipar_m <= max_y_ipar_m ; - ++y_ipar_m) - { - const int y_ipar = y_ipar_posn + y_ipar_m; - const int y_irho - = ye.irho_of_iperp_ipar(y_iperp, y_ipar); - const int y_isigma - = ye.isigma_of_iperp_ipar(y_iperp, y_ipar); - const int JJ - = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma); - const fp gz_Jac = gz.Jacobian(x_iperp,x_ipar, y_ipar_m); - Jac.add_to_element(II,JJ, Jac_sum*gz_Jac); - } - } - } + add_molecule_point_to_Jacobian(ps, xp, + x_irho, x_isigma, + xm_irho, x_isigma, + Jac_rho + Jac_rho_rho); + } + + // partial_sigma, partial_sigma_sigma + for (int m_isigma = xp.molecule_min_m() ; + m_isigma <= xp.molecule_max_m() ; + ++m_isigma) + { + const int xm_isigma = x_isigma + m_isigma; + const fp Jac_sigma = Jacobian_coeff_sigma + * xp.partial_sigma_coeff(m_sigma); + 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); + } + + // partial_rho_sigma + for (int m_irho = xp.molecule_min_m() ; + m_irho <= xp.molecule_max_m() ; + ++m_irho) + { + for (int m_isigma = xp.molecule_min_m() ; + m_isigma <= xp.molecule_max_m() ; + ++m_isigma) + { + 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); + } + } + + } } + } +} + +//****************************************************************************** + +// +// 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. +// +// 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. +// mol = The molecule coefficient. +// +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) +{ +// 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); + } } - } } + } |