diff options
Diffstat (limited to 'src/gr/horizon_Jacobian.cc')
| -rw-r--r-- | src/gr/horizon_Jacobian.cc | 106 |
1 files changed, 73 insertions, 33 deletions
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index 99d8d44..f16d9c5 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -8,6 +8,7 @@ /// 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> @@ -85,25 +86,75 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // 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, except that -// the d/dr term is done by numerical finite differencing. +// the d/dr term is done by numerical perturbation. +// +// The overall algorithm is +// save_H = H +// h += perturbation_amplitude +// evaluate H(h) (silently) +// Jac = diagonal[ (H(h) - save_H)/perturbation_amplitude ] +// h -= perturbation_amplitude +// H = save_H +// Jac += Jacobian coeffs * molecule coeffs // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface // H # H(h) assumed to already be computed // partial_H_wrt_partial_d_h # Jacobian coefficients -// partial_H_wrt_partial_dd_h +// partial_H_wrt_partial_dd_h # (also assumed to already be computed) // // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. // -void horizon_Jacobian_SD(patch_system& ps, - Jacobian& Jac) +void horizon_Jacobian(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_interpolator_info& gii, + const struct Jacobian_info& Jac_info, + Jacobian& Jac) { -CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_SD"); +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); ps.compute_synchronize_Jacobian(); Jac.zero(); + +// +// compute d/dr term in Jacobian by numerical perturbation +// + +CCTK_VInfo(CCTK_THORNSTRING, + " computing d/dr terms by numerical perturbation"); +const fp epsilon = Jac_info.perturbation_amplitude; +ps.gridfn_copy(nominal_gfns::gfn__H, nominal_gfns::gfn__save_H); +ps.add_to_ghosted_gridfn(epsilon, ghosted_gfns::gfn__h); +horizon_function(ps, cgi, gii, false, NULL, false); + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + { + patch& p = ps.ith_patch(pn); + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma); + const fp old_H = p.gridfn(nominal_gfns::gfn__save_H, + irho,isigma); + const fp new_H = p.gridfn(nominal_gfns::gfn__H, irho,isigma); + Jac(II,II) += (new_H - old_H) / epsilon; + } + } + } +ps.add_to_ghosted_gridfn(-epsilon, ghosted_gfns::gfn__h); +ps.gridfn_copy(nominal_gfns::gfn__save_H, nominal_gfns::gfn__H); + + +// +// now the main symbolic-differentiation Jacobian computation +// + +CCTK_VInfo(CCTK_THORNSTRING, + " computing main terms by symbolic differentiation"); for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) { patch& xp = ps.ith_patch(xpn); @@ -213,10 +264,6 @@ Jac.zero(); } } } - -// compute d/dr term by numerical finite differencing -perturb_h(Jacobian_info.perturbation_amplitude); -horizon_function(ps, cgi, gii, false, NULL, false); } //****************************************************************************** @@ -281,36 +328,31 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn // 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 +// we assume that H = H(h) has already been evaluated +// save_H = H // for each point (y,JJ) // { // const fp save_h_y = h at y; // h at y += perturbation_amplitude; -// evaluate H(h) +// evaluate H(h) (silently) // 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) +// H = save_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) + const struct Jacobian_info& Jac_info, + Jacobian& Jac) { CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_NP"); +const fp epsilon = Jac_info.perturbation_amplitude; -// 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); +ps.gridfn_copy(nominal_gfns::gfn__H, nominal_gfns::gfn__save_H); for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) { @@ -327,10 +369,8 @@ jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H); 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); + yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) += epsilon; + horizon_function(ps, cgi, ii, false, NULL, false); for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) { patch& xp = ps.ith_patch(xpn); @@ -344,18 +384,18 @@ jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H); ++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; + const fp old_H = xp.gridfn(nominal_gfns::gfn__save_H, + x_irho,x_isigma); + const fp new_H = xp.gridfn(nominal_gfns::gfn__H, + x_irho,x_isigma); + Jac(II,JJ) = (new_H - old_H) / epsilon; } } } - yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) - = save_h_y; + 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; +ps.gridfn_copy(nominal_gfns::gfn__save_H, nominal_gfns::gfn__H); } |