From dd883a64eebd29a9a3baa0d724ee361424a704de Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 7 Oct 2002 15:05:21 +0000 Subject: Added Files: gr/horizon_Jacobian.cc Removed Files: driver/horizon_Jacobian.cc move this file from driver/ back to gr/ since it's not really a "driver", but rathera part of the GR computations (abeit one which doesn't know much about gr :) git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@800 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/driver/horizon_Jacobian.cc | 512 ----------------------------------------- 1 file changed, 512 deletions(-) delete mode 100644 src/driver/horizon_Jacobian.cc (limited to 'src/driver') diff --git a/src/driver/horizon_Jacobian.cc b/src/driver/horizon_Jacobian.cc deleted file mode 100644 index 09a6a4d..0000000 --- a/src/driver/horizon_Jacobian.cc +++ /dev/null @@ -1,512 +0,0 @@ -// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h) -// $Header$ -// -// <<>> -// -// horizon_Jacobian - top-level driver to compute the Jacobian -/// -/// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation -/// -/// horizon_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff -/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian -/// - -#include -#include -#include -#include - -#include "util_Table.h" -#include "cctk.h" -#include "cctk_Arguments.h" - -#include "config.h" -#include "stdc.h" -#include "../jtutil/util.hh" -#include "../jtutil/array.hh" -#include "../jtutil/cpm_map.hh" -#include "../jtutil/linear_map.hh" -using jtutil::error_exit; - -#include "../patch/coords.hh" -#include "../patch/grid.hh" -#include "../patch/fd_grid.hh" -#include "../patch/patch.hh" -#include "../patch/patch_edge.hh" -#include "../patch/patch_interp.hh" -#include "../patch/ghost_zone.hh" -#include "../patch/patch_system.hh" - -#include "../elliptic/Jacobian.hh" - -#include "../gr/gfns.hh" -#include "../gr/gr.hh" - -#include "driver.hh" - -//****************************************************************************** - -// -// ***** prototypes for functions local to this file ***** -// - -namespace { -bool horizon_Jacobian_NP(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& ii, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag); - -void horizon_Jacobian_partial_SD(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag); -void add_ghost_zone_Jacobian(const patch_system& ps, - Jacobian& Jac, - fp mol, - const patch& xp, const ghost_zone& xmgz, - int x_II, - int xm_irho, int xm_isigma); -bool horizon_Jacobian_dr_FD(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag); - } - -//****************************************************************************** - -// -// This function is a top-level driver to compute the Jacobian matrix -// J[H(h)] of the LHS function H(h). It just decodes the Jacobian_method -// parameter and calls the appropriate subfunction. -// -// We assume that H(h) has already been computed. -// -// Results: -// This function returns true for a successful computation, or false -// for failure. See horizon_Jacobian() (in ../gr/) for details on -// various ways the computation may fail. -// -bool horizon_Jacobian(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) -{ -switch (Jacobian_info.Jacobian_method) - { -case Jacobian_method__numerical_perturb: - if (! horizon_Jacobian_NP(ps, Jac, - cgi, gi, Jacobian_info, - print_msg_flag)) - then return false; // *** ERROR RETURN *** - break; -case Jacobian_method__symbolic_diff: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" horizon_Jacobian(): Jacobian_method == \"symbolic differentiation\"\n" -" isn't implemented (yet); reverting to\n" -" \"symbolic differentiation with finite diff d/dr\""); - // fall through -case Jacobian_method__symbolic_diff_with_FD_dr: - horizon_Jacobian_partial_SD(ps, Jac, - cgi, gi, Jacobian_info, - print_msg_flag); - if (! horizon_Jacobian_dr_FD(ps, Jac, - cgi, gi, Jacobian_info, - print_msg_flag)) - then return false; // *** ERROR RETURN *** - break; -default: - error_exit(PANIC_EXIT, -"***** horizon_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n" -" (this should never happen!)\n" -, - int(Jacobian_info.Jacobian_method)); /*NOTREACHED*/ - } -return true; // *** NORMAL RETURN *** -} - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// 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: -// -// 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) (silently) -// for each point (x,II) -// { -// Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude; -// } -// h at y = save_h_y; -// } -// H = save_H -// -// Inputs (angular gridfns, on ghosted grid): -// h # shape of trial surface -// H # H(h) assumed to already be computed -// -// Outputs: -// The Jacobian matrix is stored in the Jacobian object Jac. -// -// Results: -// This function returns true for a successful computation, or false -// for failure. See horizon_Jacobian() (in ../gr/) for details on -// various ways the computation may fail. -// -namespace { -bool horizon_Jacobian_NP(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& ii, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) -{ -if (print_msg_flag) - then CCTK_VInfo(CCTK_THORNSTRING, - " horizon Jacobian (numerical perturbation)"); -const fp epsilon = Jacobian_info.perturbation_amplitude; - -ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H); - - for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) - { - patch& yp = ps.ith_patch(ypn); - if (print_msg_flag) - then 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(gfns::gfn__h, y_irho,y_isigma); - yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon; - if (! horizon_function(ps, cgi, ii)) - then return false; // *** ERROR RETURN *** - 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 old_H = xp.gridfn(gfns::gfn__save_H, x_irho,x_isigma); - const fp new_H = xp.gridfn(gfns::gfn__H, x_irho,x_isigma); - Jac(II,JJ) = (new_H - old_H) / epsilon; - } - } - } - yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y; - } - } - } - -ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H); -return true; // *** NORMAL RETURN *** -} - } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function computes the partial derivative terms in the Jacobian -// matrix of the LHS function H(h), 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 -// H # H(h) assumed to already be computed -// partial_H_wrt_partial_d_h # Jacobian coefficients -// partial_H_wrt_partial_dd_h # (also assumed to already be computed) -// -// Outputs: -// The Jacobian matrix is stored in the Jacobian object Jac. -// -namespace { -void horizon_Jacobian_partial_SD(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) -{ -if (print_msg_flag) - then CCTK_VInfo(CCTK_THORNSTRING, - " horizon Jacobian: partial-deriv terms (symbolic diff)"); - -Jac.zero_matrix(); -ps.compute_synchronize_Jacobian(); - - 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) - { - // - // 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) - // - - // 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(gfns::gfn__partial_H_wrt_partial_d_h_1, - x_irho, x_isigma); - const fp Jacobian_coeff_sigma - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2, - x_irho, x_isigma); - const fp Jacobian_coeff_rho_rho - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11, - x_irho, x_isigma); - const fp Jacobian_coeff_rho_sigma - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12, - x_irho, x_isigma); - const fp Jacobian_coeff_sigma_sigma - = xp.gridfn(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() ; - m_irho <= xp.molecule_max_m() ; - ++m_irho) - { - const int xm_irho = x_irho + m_irho; - const fp Jac_rho = Jacobian_coeff_rho - * xp.partial_rho_coeff(m_irho); - const fp Jac_rho_rho = Jacobian_coeff_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, Jac, - Jac_sum, - xp, xp.minmax_rho_ghost_zone(m_irho < 0), - II, xm_irho, x_isigma); - } - } - - // 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_isigma); - const fp Jac_sigma_sigma = Jacobian_coeff_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, Jac, - Jac_sum, - xp, xp.minmax_sigma_ghost_zone(m_isigma < 0), - II, x_irho, xm_isigma); - } - } - - // 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_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 { - const ghost_zone& xmgz - = xp.corner_ghost_zone_containing_point - (m_irho < 0, m_isigma < 0, - xm_irho, xm_isigma); - add_ghost_zone_Jacobian(ps, Jac, - Jac_rho_sigma, - xp, xmgz, - II, xm_irho, xm_isigma); - } - } - } - } - - } - } - } -} - } - -//****************************************************************************** - -// -// This function adds the ghost-zone Jacobian dependency contributions -// for a single ghost-zone point, to a Jacobian matrix. -// -// Arguments: -// ps = The patch system. -// Jac = (out) The Jacobian matrix. -// mol = The molecule coefficient. -// xp = The patch containing the center 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. -// -namespace { -void add_ghost_zone_Jacobian(const patch_system& ps, - Jacobian& Jac, - fp mol, - const patch& xp, const ghost_zone& xmgz, - int x_II, - int xm_irho, int 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); - -// 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_ym, global_max_ym; -ps.synchronize_Jacobian_global_minmax_ym(global_min_ym, global_max_ym); -jtutil::array1d Jacobian_buffer(global_min_ym, global_max_ym); - -// on what other points y does this molecule point xm depend -// via the patch_system::synchronize() operation? -int y_iperp; -int y_posn, min_ym, max_ym; -const patch_edge& ye = ps.synchronize_Jacobian(xmgz, - xm_iperp, xm_ipar, - y_iperp, - y_posn, min_ym, max_ym, - Jacobian_buffer); -const patch& yp = ye.my_patch(); - -// add the Jacobian contributions from the ym points - for (int ym = min_ym ; ym <= max_ym ; ++ym) - { - const int y_ipar = y_posn + ym; - 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); - Jac(x_II,JJ) += mol*Jacobian_buffer(ym); - } -} - } - -//****************************************************************************** - -// -// This function sums the d/dr terms into the Jacobian matrix of the -// LHS function H(h), computing those terms by finite differencing. -// -// The basic algorithm is that -// Jac += diag[ (H(h+epsilon) - H(h)) / epsilon ] -// -// Inputs (angular gridfns, on ghosted grid): -// h # shape of trial surface -// H # H(h) assumed to already be computed -// -// Outputs: -// Jac += d/dr terms -// -// Results: -// This function returns true for a successful computation, or false -// if the computation failed because the geometry interpolation would -// need data outside the Cactus grid, or data from an excised region. -// FIXME: excision isn't implemented yet :( -// -namespace { -bool horizon_Jacobian_dr_FD(patch_system& ps, - Jacobian& Jac, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - bool print_msg_flag) -{ -if (print_msg_flag) - then CCTK_VInfo(CCTK_THORNSTRING, - " horizon Jacobian: d/dr terms (finite diff)"); - -const fp epsilon = Jacobian_info.perturbation_amplitude; - -// compute H(h+epsilon) -ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H); -ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h); -if (! horizon_function(ps, cgi, gi)) - then return false; // *** ERROR RETURN *** - - 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(gfns::gfn__save_H, irho,isigma); - const fp new_H = p.gridfn(gfns::gfn__H, irho,isigma); - Jac(II,II) += (new_H - old_H) / epsilon; - } - } - } - -// restore h and H -ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h); -ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H); - -return true; // *** NORMAL RETURN *** -} - } -- cgit v1.2.3