path: root/src/driver/horizon_Jacobian.cc

// 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>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"

#include "stdc.h"
#include "config.hh"
#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 ***
}
	  }