path: root/src/gr/horizon_function.cc

// horizon_function.cc -- evaluate LHS function H(h)
// $Id$
//
// <<<prototypes>>>
// horizon_function - top-level driver
/// xyz_posns_and_deriv_coeffs - init xyz posns of grid points, xyz deriv coeffs
/// interpolate_geometry - interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid
//

#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 "AHFinderDirect.hh"

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

//
// ***** static data structures *****
//

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

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

namespace {
void xyz_posns_and_deriv_coeffs(patch_system& ps);
void interpolate_geometry(patch_system& ps,
			  const struct cactus_grid_info& cgi,
			  const struct geometry_interpolator_info& ii);
void compute_H(patch_system& ps);
	  };

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

//
// This function computes the LHS function H(h).
//
// Inputs (angular gridfns, on ghosted grid):
// ... defined on ghosted grid
// ... only values on nominal grid are actually used as input
//	h				# shape of trial surface
//
// Inputs (Cactus 3-D gridfns):
//	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
//	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
//
// Outputs (angular gridfns, all on the nominal grid):
//	## computed directly from h
//	xx,yy,zz			# xyz positions of grid points
//	X_ud_*, X_udd_*			# xyz derivative coefficients
//	## interpolated from 3-D Cactus grid
//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
//	## computed by Maple-generated code
//	g_uu_{11,12,13,22,23,33}	# $g^{ij}$
//	K				# $K$
//	K_dd_{11,12,13,22,23,33}	# $K^{ij}$
//	partial_d_ln_sqrt_g_d		# $\partial_i \ln \sqrt{g}$
//	partial_d_g_uu_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g^{ij}$
//	H				# $H = H(h)$
//
void horizon_function(patch_system& ps,
		      const struct cactus_grid_info& cgi,
		      const struct geometry_interpolator_info& ii)
{
CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");

// fill in values of ghosted gridfns in ghost zones
ps.synchronize_ghost_zones(ghosted_gfns::gfn__h, ghosted_gfns::gfn__h);

// compute xyz positions of grid points and xyz derivative coeffs
xyz_posns_and_deriv_coeffs(ps);

// interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$
interpolate_geometry(ps, cgi, ii);

// compute remaining gridfns --> $H$
// by algebraic ops and angular finite differencing
compute_H(ps);
}

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

//
// This function computes
// - the xyz positions of the grid points (gridfns xx, yy, and zz)
// - the xyz derivative coefficients (gridfns X_ud and X_udd).
// on the nominal grid.
//
namespace {
void xyz_posns_and_deriv_coeffs(patch_system& ps)
{
CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");

	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 fp r
			= p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma);
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);
		fp x, y, z;
		p.xyz_of_r_rho_sigma(r,rho,sigma, x,y,z);

		// xyz positions of grid points
		p.gridfn(nominal_gfns::gfn__xx, irho,isigma) = x;
		p.gridfn(nominal_gfns::gfn__yy, irho,isigma) = y;
		p.gridfn(nominal_gfns::gfn__zz, irho,isigma) = z;

		// 1st derivative coefficient gridfns X_ud
		p.gridfn(nominal_gfns::gfn__X_ud_11, irho,isigma)
			= p.partial_rho_wrt_x(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_ud_12, irho,isigma)
			= p.partial_rho_wrt_y(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_ud_13, irho,isigma)
			= p.partial_rho_wrt_z(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_ud_21, irho,isigma)
			= p.partial_sigma_wrt_x(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_ud_22, irho,isigma)
			= p.partial_sigma_wrt_y(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_ud_23, irho,isigma)
			= p.partial_sigma_wrt_z(x,y,z);

		// 2nd derivative coefficient gridfns X_udd
		p.gridfn(nominal_gfns::gfn__X_udd_111, irho,isigma)
			= p.partial2_rho_wrt_xx(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_112, irho,isigma)
			= p.partial2_rho_wrt_xy(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_113, irho,isigma)
			= p.partial2_rho_wrt_xz(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_122, irho,isigma)
			= p.partial2_rho_wrt_yy(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_123, irho,isigma)
			= p.partial2_rho_wrt_yz(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_133, irho,isigma)
			= p.partial2_rho_wrt_zz(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_211, irho,isigma)
			= p.partial2_sigma_wrt_xx(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_212, irho,isigma)
			= p.partial2_sigma_wrt_xy(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_213, irho,isigma)
			= p.partial2_sigma_wrt_xz(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_222, irho,isigma)
			= p.partial2_sigma_wrt_yy(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_223, irho,isigma)
			= p.partial2_sigma_wrt_yz(x,y,z);
		p.gridfn(nominal_gfns::gfn__X_udd_233, irho,isigma)
			= p.partial2_sigma_wrt_zz(x,y,z);
		}
		}
	}
}
	  }

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

//
// This function interpolates the Cactus $g_{ij}$ and $K_{ij}$, and
// the spatial derivatives $\partial_k g_{ij}$, to the trial horizon
// surface position given by h.
//
// Inputs (angular gridfns, on ghosted grid):
// ... defined on ghosted grid
// ... only values on nominal grid are actually used as input
//	h				# shape of trial surface
//
// Inputs (angular gridfns, all on the nominal grid):
//	xx,yy,zz			# xyz positions of grid points
//
// Inputs (Cactus 3-D gridfns):
//	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
//	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
//
// Inputs (angular gridfns, all on the nominal grid):
//	## computed directly from h
//	xx,yy,zz			# xyz positions
//
// Outputs (angular gridfns, all on the nominal grid):
//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
//
// On the first call this function also modifies the interpolator
// parameter table.
//
namespace {
void interpolate_geometry(patch_system& ps,
			  const struct cactus_grid_info& cgi,
			  const struct geometry_interpolator_info& gii)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "      interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid");

int status;

const int N_interp_points = ps.N_grid_points();
const int interp_coords_type_code = CCTK_VARIABLE_REAL;
const void* const interp_coords[N_GRID_DIMS]
  = {
    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__xx)),
    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__yy)),
    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__zz)),
    };

const int N_INPUT_ARRAYS = 12;
const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
	= {
	  // $g_{ij}$
	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
						  CCTK_VARIABLE_REAL,
	    // $K_{ij}$
	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
						  CCTK_VARIABLE_REAL,
	  };
const void* const input_arrays[N_INPUT_ARRAYS]
	= {
	  static_cast<const void*>(cgi.g_dd_11_data),
	  static_cast<const void*>(cgi.g_dd_12_data),
	  static_cast<const void*>(cgi.g_dd_13_data),
	  static_cast<const void*>(cgi.g_dd_22_data),
	  static_cast<const void*>(cgi.g_dd_23_data),
	  static_cast<const void*>(cgi.g_dd_33_data),
	  static_cast<const void*>(cgi.K_dd_11_data),
	  static_cast<const void*>(cgi.K_dd_12_data),
	  static_cast<const void*>(cgi.K_dd_13_data),
	  static_cast<const void*>(cgi.K_dd_22_data),
	  static_cast<const void*>(cgi.K_dd_23_data),
	  static_cast<const void*>(cgi.K_dd_33_data),
	  };

const int N_OUTPUT_ARRAYS = 30;

const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
	= {
 // $g_{ij}$         $\partial_x g_{ij}$ $\partial_y g_{ij}$ $\partial_z g_{ij}$
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 // $K_{ij}$
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
		     CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
					 CCTK_VARIABLE_REAL,
	  };
const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
	= {
	  0, 0, 0, 0,		// g_dd_11
	  1, 1, 1, 1,		// g_dd_12
	  2, 2, 2, 2,		// g_dd_13
	  3, 3, 3, 3,		// g_dd_22
	  4, 4, 4, 4,		// g_dd_23
	  5, 5, 5, 5,		// g_dd_33
	   6,  7,  8,		// K_dd_{11,12,13}
	       9, 10,		// K_dd_{22,23}
		  11,		// K_dd_33
	  };
#define DERIV(x)	x
const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
	= {
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_11
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_12
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_13
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_22
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_23
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_33
	  DERIV(0), DERIV(0), DERIV(0),			// K_dd_{11,12,13}
		    DERIV(0), DERIV(0),			// K_dd_{22,23}
			      DERIV(0),			// K_dd_{33}
	  };
void* const output_arrays[N_OUTPUT_ARRAYS]
  = {
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_11)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_111)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_211)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_311)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_12)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_112)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_212)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_312)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_13)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_113)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_213)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_313)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_22)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_122)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_222)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_322)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_23)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_123)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_223)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_323)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_33)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_133)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_233)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_333)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_11)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_12)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_13)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_22)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_23)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_33)),
    };

bool first_time = true;
if (first_time)
   then {
	// store derivative info in interpolator parameter table
	first_time = false;
	CCTK_VInfo(CCTK_THORNSTRING,
		   "         setting up derivative info for interpolator");

	status = Util_TableSetIntArray(gii.param_table_handle,
				       N_OUTPUT_ARRAYS, operand_indices,
				       "operand_indices");
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        unable to set operand_indices in interpolator parameter table!\n"
"        Util_TableSetIntArray() status=%d\n"
,
			   status);				/*NOTREACHED*/

	status = Util_TableSetIntArray(gii.param_table_handle,
				       N_OUTPUT_ARRAYS, operation_codes,
				       "operation_codes");
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        unable to set operation_codes in interpolator parameter table!\n"
"        Util_TableSetIntArray() status=%d\n"
,
			   status);				/*NOTREACHED*/
	}

CCTK_VInfo(CCTK_THORNSTRING,
	   "         calling interpolator (%d points)",
	   N_interp_points);
status = CCTK_InterpLocalUniform(N_GRID_DIMS,
				 gii.operator_handle, gii.param_table_handle,
				 cgi.coord_origin, cgi.coord_delta,
				 N_interp_points,
				    interp_coords_type_code,
				    interp_coords,
				 N_INPUT_ARRAYS,
				    cgi.gridfn_dims,
				    input_array_type_codes,
				    input_arrays,
				 N_OUTPUT_ARRAYS,
				    output_array_type_codes,
				    output_arrays);
if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
   then {
	// FIXME: should look at interpolator output table entries
	//	  and report *which* point is out-of-range
	error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        point out of range when interpolating geometry info from 3-D grid!\n"
"        ==> the trial horizon surface is (at least partially)\n"
"            outside the grid and/or in an excised region\n");	/*NOTREACHED*/
	}
if (status < 0)
   then error_exit(ERROR_EXIT,
"***** interpolate_geometry(): error return from interpolator!\n"
"        CCTK_InterpLocalUniform() status=%d\n"
,
		   status);					/*NOTREACHED*/
}
	  }

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

//
// This function computes H(h) by a mixture of algebraic operations
// and (rho,sigma) finite differencing, on the angular grid.
//
namespace {
void compute_H(patch_system& ps)
{
CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");

	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)
		{
		// "call" the Maple-generated code
		// ... each file has a separate set of temporary variables,
		//     and so must be inside its own set of { } braces
		#include "cg.hh"
		  {
		#include "cg/inverse_metric.c"
		  }
		  {
		#include "cg/extrinsic_curvature_trace_raise.c"
		  }
		  {
		#include "cg/inverse_metric_gradient.c"
		  }
		  {
		#include "cg/metric_det_gradient.c"
		  }
		  {
		#include "cg/horizon_function.c"
		  }
		}
		}
	}
}
	  }