diff options
Diffstat (limited to 'src/gr/horizon_Jacobian.cc')
| -rw-r--r-- | src/gr/horizon_Jacobian.cc | 115 |
1 files changed, 107 insertions, 8 deletions
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index 7718796..bcaec78 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -3,9 +3,11 @@ // // <<<prototypes for functions local to this file>>> // create_Jacobian - create a Jacobian matrix of the appropriate type -// horizon_Jacobian - top-level driver to compute the Jacobian +// +// horizon_Jacobian_SD - compute the Jacobian by symbolic differentiation /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian // +// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation #include <stdio.h> #include <assert.h> @@ -53,12 +55,18 @@ void add_ghost_zone_Jacobian(const patch_system& ps, } //****************************************************************************** +//****************************************************************************** +//****************************************************************************** // // This function decodes the Jacobian_type parameter and creates the // appropriate type of Jacobian matrix object. // -Jacobian& create_Jacobian(const patch_system& ps, +// FIXME: the patch system shouldn't really have to be non-const, but +// the Jacobian constructors all require this to allow the +// linear solvers to directly update gridfns +// +Jacobian& create_Jacobian(patch_system& ps, const char Jacobian_type[]) { if (STRING_EQUAL(Jacobian_type, "dense matrix")) @@ -69,12 +77,14 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, } //****************************************************************************** +//****************************************************************************** +//****************************************************************************** // // This function computes the Jacobian matrix of the LHS function H(h) -// from the Jacobian coefficient (angular) gridfns. The computation is -// done a Jacobian row at a time, using equation (25) of my 1996 apparent -// horizon finding paper. +// by symbolic differentiation from the Jacobian coefficient (angular) +// gridfns. The computation is done a Jacobian row at a time, using +// equation (25) of my 1996 apparent horizon finding paper. // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface @@ -84,9 +94,10 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. // -void horizon_Jacobian(patch_system& ps, Jacobian& Jac) +void horizon_Jacobian_SD(patch_system& ps, + Jacobian& Jac) { -CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_SD"); ps.compute_synchronize_Jacobian(); Jac.zero(); @@ -189,7 +200,7 @@ Jac.zero(); 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, + (ps, xp, xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0, xm_irho, xm_isigma), II, xm_irho, xm_isigma, @@ -254,3 +265,91 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn } } } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function computes the Jacobian matrix of the LHS function H(h) +// by numerical perturbation of the H(h) function. The algorithm is as +// follows: +// +// evaluate H(h) +// array-copy H(h) to scratch array save_H +// for each point (y,JJ) +// { +// const fp save_h_y = h at y; +// h at y += perturbation_amplitude; +// evaluate H(h) +// for each point (x,II) +// { +// Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude; +// } +// h at y = save_h_y; +// } +// array-copy save_H back to H(h) +// +void horizon_Jacobian_NP(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_interpolator_info& ii, + Jacobian& Jac, + fp perturbation_amplitude) +{ +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_NP"); + +// evaluate H(h) +jtutil::norm<fp> H_norms; +horizon_function(ps, cgi, ii, false, H_norms); + +// array-copy H(h) to scratch array save_H +fp *save_H = new fp[Jac.NN()]; +jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H); + + for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) + { + patch& yp = ps.ith_patch(ypn); + CCTK_VInfo(CCTK_THORNSTRING, " perturbing in %s patch", yp.name()); + + for (int y_irho = yp.min_irho() ; y_irho <= yp.max_irho() ; ++y_irho) + { + for (int y_isigma = yp.min_isigma() ; + y_isigma <= yp.max_isigma() ; + ++y_isigma) + { + const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); + + const fp save_h_y = yp.ghosted_gridfn(ghosted_gfns::gfn__h, + y_irho,y_isigma); + yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) + += perturbation_amplitude; + H_norms.reset(); + horizon_function(ps, cgi, ii, false, H_norms, false); + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) + { + patch& xp = ps.ith_patch(xpn); + + 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) + { + const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); + const fp new_xH = xp.gridfn(nominal_gfns::gfn__H, + x_irho,x_isigma); + Jac(II,JJ) = (new_xH - save_H[II]) / perturbation_amplitude; + } + } + } + yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) + = save_h_y; + } + } + } + +jtutil::array_copy(Jac.NN(), save_H, ps.gridfn_data(nominal_gfns::gfn__H)); +delete[] save_H; +} |