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

1 files changed, 0 insertions, 512 deletions

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$
-//
-// <<<prototypes for functions local to this file>>>
-//
-// 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 <stdio.h>
-#include <assert.h>
-#include <math.h>
-#include <vector.h>
-
-#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<fp> 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 ***
-}
-	  }