// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
// $Id$
//
// <<<prototypes for functions local to this file>>>
// horizon_Jacobian - top-level driver
//

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>

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

#include "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
using jtutil::error_exit;

#include "../config.hh"

#include "../util/coords.hh"
#include "../util/grid.hh"
#include "../util/fd_grid.hh"
#include "../util/patch.hh"
#include "../util/patch_edge.hh"
#include "../util/patch_interp.hh"
#include "../util/ghost_zone.hh"
#include "../util/patch_system.hh"

#include "gfn.hh"
#include "Jacobian.hh"
#include "AHFinderDirect.hh"

//******************************************************************************

//
// ***** prototypes for functions local to this file *****
//

namespace {
void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
				    int x_irho, int x_isigma,
				    int xm_irho, int xm_isigma,
				    fp mol);
	  }

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// 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.
//
// Inputs (angular gridfns, on ghosted grid):
//	h				# shape of trial surface
//	partial_H_wrt_partial_d_h,	# Jacobian coefficients
//	partial_H_wrt_partial_dd_h,
//
// Outputs:
//	The Jacobian matrix is stored in the Jacobian object Jac.
void horizon_Jacobian(const patch_system& ps,
		      Jacobian& Jac);
{
CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian");

ps.compute_synchronize_Jacobian();
Jac.zero();

    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 Jacobian row for this grid point, i.e.
	//	  partial H(this point x, Jacobian row II)
	//	---------------------------------------------
	//	partial h(other points y, Jacobian column JJ)
	//
	// FIXME FIXME: we still have to take into account the
	//		position-dependence of the coefficients,
	//		cf the difference between J[3H(h)] and J[2H(h)];
	//		here we're sort of computing the former, but
	//		Newton's method really wants the latter
	//

	// Jacobian coefficients for this point
	const fp Jacobian_coeff_rho
	   = xp.gridfn(gfn__partial_H_wrt_partial_d_h_1, x_irho,x_isigma);
	const fp Jacobian_coeff_sigma
	   = xp.gridfn(gfn__partial_H_wrt_partial_d_h_2, x_irho,x_isigma);
	const fp Jacobian_coeff_rho_rho
	   = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_11, x_irho,x_isigma);
	const fp Jacobian_coeff_rho_sigma
	   = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_12, x_irho,x_isigma);
	const fp Jacobian_coeff_sigma_sigma
	   = xp.gridfn(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_rho);
	    const fp Jac_rho_rho = Jacobian_coeff_rho_rho
				   * xp.partial_rho_rho_coeff(m_rho);
	    add_molecule_point_to_Jacobian(ps, xp,
					   x_irho, x_isigma,
					   xm_irho, x_isigma,
					   Jac_rho + Jac_rho_rho);
	    }

	// 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_sigma);
	    const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma
				       * xp.partial_sigma_sigma_coeff(m_sigma);
	    add_molecule_point_to_Jacobian(ps, xp,
					   x_irho, x_isigma,
					   x_irho, xm_isigma,
					   Jac_sigma + Jac_sigma_sigma);
	    }

	// 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_rho, m_sigma);
	    add_molecule_point_to_Jacobian(ps, xp,
					   x_irho, x_isigma,
					   xm_irho, xm_isigma,
					   Jac_rho_sigma);
	    }
	    }

	}
	}
    }
}

//******************************************************************************

//
// This function adds all the elements for a single molecule point to
// a Jacobian matrix, including any contributions from other patches if
// the specified point lies in a ghost zone.
//
// Arguments:
// ps = The patch system.
// xp = The patch containing the center point of the molecule.
// x_(irho,isigma) = The coordinates of the center point of the molecule.
// xm_(irho,isigma) = The coordinates of the x+m point of the molecule.
// mol = The molecule coefficient.
//
namespace {
void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
				    int x_irho, int x_isigma,
				    int xm_irho, int xm_isigma,
				    fp mol)
{
// Jacobian row
const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);

if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
   then {
	const int JJ = ps.gpn_of_patch_irho_isigma(xp, xm_irho, xm_isigma);
	Jac.add_to_element(II,JJ, Jac_coeff);
	}
   else {
	const ghost_zone& xmgz
		= xp.ghost_zone_containing_point(xm_irho, 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);

	// on what other points ym does this molecule point xm depend
	// via the patch_system::synchronize() operation?
	const patch&      ymp = ps.synchronize_Jacobian_y_patch(xmgz);
	const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz);
	const int min_ym_ipar_m = ps.synchronize_Jacobian_min_y_ipar_m(xmgz);
	const int max_ym_ipar_m = ps.synchronize_Jacobian_max_y_ipar_m(xmgz);
	const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp);
	const int ym_ipar_posn
	   = ps.synchronize_Jacobian_y_ipar_posn(xmgz, xm_iperp, xm_ipar);
	    for (int ym_ipar_m = min_ym_ipar_m ;
		 ym_ipar_m <= max_ym_ipar_m ;
		 ++ym_ipar_m)
	    {
	    const int ym_ipar = ym_ipar_posn + ym_ipar_m;
	    const int ym_irho   = yme.  irho_of_iperp_ipar(ym_iperp,ym_ipar);
	    const int ym_isigma = yme.isigma_of_iperp_ipar(ym_iperp,ym_ipar);
	    const int JJ = ps.gpn_of_patch_irho_isigma(ymp, ym_irho, ym_isigma);
	    const fp sync_Jac = ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar,
							ym_ipar_m);
	    Jac.add_to_element(II,JJ, mol*sync_Jac);
	    }
	}
}
	  }