aboutsummaryrefslogtreecommitdiff
path: root/src/gr
diff options
context:
space:
mode:
Diffstat (limited to 'src/gr')
-rw-r--r--src/gr/cg.hh36
-rw-r--r--src/gr/expansion.cc (renamed from src/gr/horizon_function.cc)121
-rw-r--r--src/gr/expansion_Jacobian.cc (renamed from src/gr/horizon_Jacobian.cc)209
-rw-r--r--src/gr/gfns.hh14
-rw-r--r--src/gr/gr.hh30
-rw-r--r--src/gr/gr_gfas.minc42
-rw-r--r--src/gr/horizon.maple195
-rw-r--r--src/gr/make.code.defn4
-rw-r--r--src/gr/maple.log529
-rw-r--r--src/gr/setup_gr_gfas.maple55
10 files changed, 635 insertions, 600 deletions
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index 02bf370..affd54b 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -89,20 +89,20 @@
#define partial_d_g_dd_323 p.gridfn(gfns::gfn__partial_d_g_dd_323, irho,isigma)
#define partial_d_g_dd_333 p.gridfn(gfns::gfn__partial_d_g_dd_333, irho,isigma)
-#define H p.gridfn(gfns::gfn__H, irho,isigma)
-
-#define partial_H_wrt_partial_d_h_1 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1, irho,isigma)
-#define partial_H_wrt_partial_d_h_2 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2, irho,isigma)
-#define partial_H_wrt_partial_dd_h_11 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11, irho,isigma)
-#define partial_H_wrt_partial_dd_h_12 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12, irho,isigma)
-#define partial_H_wrt_partial_dd_h_22 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma)
-
-#define save_H p.gridfn(gfns::gfn__save_H, irho,isigma)
+#define Theta p.gridfn(gfns::gfn__Theta, irho,isigma)
+
+#define partial_Theta_wrt_partial_d_h_1 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_d_h_1, irho,isigma)
+#define partial_Theta_wrt_partial_d_h_2 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_d_h_2, irho,isigma)
+#define partial_Theta_wrt_partial_dd_h_11 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_11, irho,isigma)
+#define partial_Theta_wrt_partial_dd_h_12 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_12, irho,isigma)
+#define partial_Theta_wrt_partial_dd_h_22 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_22, irho,isigma)
+
+#define save_Theta p.gridfn(gfns::gfn__save_Theta, irho,isigma)
#define Delta_h p.gridfn(gfns::gfn__Delta_h, irho,isigma)
//******************************************************************************
@@ -148,7 +148,7 @@ fp partial_d_g_uu_322;
fp partial_d_g_uu_323;
fp partial_d_g_uu_333;
-fp HA;
-fp HB;
-fp HC;
-fp HD;
+fp Theta_A;
+fp Theta_B;
+fp Theta_C;
+fp Theta_D;
diff --git a/src/gr/horizon_function.cc b/src/gr/expansion.cc
index 0821eaf..9cafc2b 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/expansion.cc
@@ -1,8 +1,8 @@
-// horizon_function.cc -- evaluate LHS function H(h)
+// expansion.cc -- evaluate expansion Theta of trial horizon surface
// $Header$
//
// <<<prototypes for functions local to this file>>>
-// horizon_function - top-level driver
+// expansion - top-level driver
// decode_geometry_method - decode the geometry_method parameter
//
/// setup_xyz_posns - setup global xyz posns of grid points
@@ -12,7 +12,7 @@
/// h_is_finite - does function h contain any NaNs/infinities?
/// geometry_is_finite - do geometry vars contain NaN/infty?
///
-/// compute_H - compute H(h) given earlier setup
+/// compute_Theta - compute Theta(h) given earlier setup
///
#include <stdio.h>
@@ -73,10 +73,10 @@ bool h_is_finite(patch_system& ps, bool print_msg_flag);
bool geometry_is_finite(patch_system& ps,
bool print_msg_flag);
-bool compute_H(patch_system& ps,
- bool Jacobian_flag,
- jtutil::norm<fp>* H_norms_ptr,
- bool print_msg_flag);
+bool compute_Theta(patch_system& ps,
+ bool Jacobian_flag,
+ jtutil::norm<fp>* Theta_norms_ptr,
+ bool print_msg_flag);
}
//******************************************************************************
@@ -84,7 +84,7 @@ bool compute_H(patch_system& ps,
//******************************************************************************
//
-// This function computes the LHS function H(h), and optionally also
+// This function computes the LHS function Theta(h), and optionally also
// its Jacobian coefficients (from which the Jacobian matrix may be
// computed later).
//
@@ -113,17 +113,18 @@ bool compute_H(patch_system& ps,
// 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}$
-// H # $H = H(h)$
+// ## computed at the nominal grid points
+// Theta # $\Theta = \Theta(h)$
//
// Arguments:
// Jacobian_flag = true to compute the Jacobian coefficients,
// false to skip this.
// print_msg_flag = true to print status messages,
// false to skip this.
-// H_norms_ptr = (out) If this pointer is non-NULL, the norm object it
-// points to is updated with all the H values in the
-// grid. This norm object can then be used to compute
-// various (gridwise) norms of H.
+// Theta_norms_ptr = (out) If this pointer is non-NULL, the norm object it
+// points to is updated with all the Theta values
+// in the grid. This norm object can then be used
+// to compute various (gridwise) norms of Theta.
//
// Results:
// This function returns true for a successful computation, or false
@@ -132,17 +133,17 @@ bool compute_H(patch_system& ps,
// or from an excised region
// FIXME: excision isn't implemented yet :(
// - NaNs are found in h or in the interpolate geometry variables
-// - HD <= 0 (this means the interpolated g_ij aren't positive definite)
+// - Theta_D <= 0 (this means the interpolated g_ij aren't positive definite)
//
-bool horizon_function(patch_system& ps,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- bool Jacobian_flag /* = false */,
- bool print_msg_flag /* = false */,
- jtutil::norm<fp>* H_norms_ptr /* = NULL */)
+bool expansion(patch_system& ps,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ bool Jacobian_flag /* = false */,
+ bool print_msg_flag /* = false */,
+ jtutil::norm<fp>* Theta_norms_ptr /* = NULL */)
{
if (print_msg_flag)
- then CCTK_VInfo(CCTK_THORNSTRING, " horizon function");
+ then CCTK_VInfo(CCTK_THORNSTRING, " expansion");
// fill in values of all ghosted gridfns in ghost zones
ps.synchronize();
@@ -168,7 +169,7 @@ case geometry__Schwarzschild_EF:
break;
default:
error_exit(PANIC_EXIT,
-"***** horizon_function(): unknown gi.geometry_method=(int)%d!\n"
+"***** expansion(): unknown gi.geometry_method=(int)%d!\n"
" (this should never happen!)\n"
,
int(gi.geometry_method)); /*NOTREACHED*/
@@ -178,9 +179,10 @@ if (gi.check_that_geometry_is_finite
&& !geometry_is_finite(ps, print_msg_flag))
then return false; // *** ERROR RETURN ***
-// compute remaining gridfns --> $H$ and optionally Jacobian coefficients
+// compute remaining gridfns --> $\Theta$
+// and optionally also the Jacobian coefficients
// by algebraic ops and angular finite differencing
-if (!compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag))
+if (!compute_Theta(ps, Jacobian_flag, Theta_norms_ptr, print_msg_flag))
then return false; // *** ERROR RETURN ***
return true; // *** NORMAL RETURN ***
@@ -586,7 +588,7 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
if (print_msg_flag)
then {
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" the trial-horizon-surface point"
@@ -743,7 +745,7 @@ if (print_msg_flag)
// or more NaNs or infinities.
//
// Results:
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
// This function returns true if all the h values are finite, false
// otherwise (i.e. if it contains any NaNs or infinities).
#else
@@ -757,7 +759,7 @@ bool h_is_finite(patch_system& ps, bool print_msg_flag)
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING, " checking that h is finite");
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
patch& p = ps.ith_patch(pn);
@@ -775,7 +777,7 @@ if (print_msg_flag)
const fp sigma = p.sigma_of_isigma(isigma);
const fp drho = jtutil::degrees_of_radians(rho);
const fp dsigma = jtutil::degrees_of_radians(sigma);
- CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+ CCTK_VWarn(expansion_warning_levels::nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" h=%g isn't finite!\n"
@@ -791,7 +793,7 @@ if (print_msg_flag)
}
return true; // *** all values finite ***
#else
-CCTK_VWarn(horizon_function_warning_levels::skip_nonfinite,
+CCTK_VWarn(expansion_warning_levels::skip_nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
" no finite() fn ==> skipping is-h-finite check");
return true; // *** no check possible ***
@@ -810,7 +812,7 @@ return true; // *** no check possible ***
// or more NaNs or infinities.
//
// Results:
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
// This function returns true if all the geometry variables are finite,
// false otherwise (i.e. if they contain any NaNs or infinities).
#else
@@ -824,7 +826,7 @@ bool geometry_is_finite(patch_system& ps, bool print_msg_flag)
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING, " checking that geometry is finite");
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
patch& p = ps.ith_patch(pn);
@@ -921,7 +923,7 @@ if (print_msg_flag)
const fp global_x = ps.global_x_of_local_x(local_x);
const fp global_y = ps.global_y_of_local_y(local_y);
const fp global_z = ps.global_z_of_local_z(local_z);
- CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+ CCTK_VWarn(expansion_warning_levels::nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" geometry isn't finite at %s patch\n"
@@ -973,7 +975,7 @@ if (print_msg_flag)
}
return true; // *** no NaNs found ***
#else
-CCTK_VWarn(horizon_function_warning_levels::skip_nonfinite,
+CCTK_VWarn(expansion_warning_levels::skip_nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
" no finite() ==> skipping is-geometry-finite check");
return true; // *** no check possible ***
@@ -986,10 +988,11 @@ return true; // *** no check possible ***
//******************************************************************************
//
-// This function computes H(h), and optionally its Jacobian coefficients,
-// (from which the Jacobian matrix may be computed later). This function
-// uses a mixture of algebraic operations and (rho,sigma) finite differencing.
-// The computation is done (entirely) on the nominal angular grid.
+// This function computes the expansion Theta(h), and optionally also
+// its Jacobian coefficients, (from which the Jacobian matrix may be
+// computed later). This function uses a mixture of algebraic operations
+// and (rho,sigma) finite differencing. The computation is done entirely
+// on the nominal angular grid.
//
// N.b. This function #includes "cg.hh", which defines "dangerous" macros
// which will stay in effect for the rest of this compilation unit!
@@ -1000,17 +1003,17 @@ return true; // *** no check possible ***
//
// Results:
// This function returns true for a successful computation, or false
-// if the computation failed because HD <= 0 (this means the interpolated
+// if the computation failed because Theta_D <= 0 (this means the interpolated
// g_ij isn't positive definite).
//
namespace {
-bool compute_H(patch_system& ps,
- bool Jacobian_flag,
- jtutil::norm<fp>* H_norms_ptr,
- bool print_msg_flag)
+bool compute_Theta(patch_system& ps,
+ bool Jacobian_flag,
+ jtutil::norm<fp>* Theta_norms_ptr,
+ bool print_msg_flag)
{
if (print_msg_flag)
- then CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)");
+ then CCTK_VInfo(CCTK_THORNSTRING, " computing Theta(h)");
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
@@ -1084,37 +1087,41 @@ if (print_msg_flag)
}
{
- // HA, HB, HC, HD
- #include "../gr.cg/horizon_function.c"
+ // Theta_A, Theta_B, Theta_C, Theta_D
+ #include "../gr.cg/expansion.c"
}
- if (HD <= 0)
+ if (Theta_D <= 0)
then {
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
- "HD <= 0 at %s patch rho=%g sigma=%g!",
+ "Theta_D <= 0 at %s patch rho=%g sigma=%g!",
p.name(), double(rho), double(sigma));
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
"(this probably means the interpolated g_ij");
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
" isn't positive definite)");
return false; // *** ERROR RETURN ***
}
// compute H via equation (14) of my 1996 horizon finding paper
- const fp sqrt_HD = sqrt(HD);
- H = HA/(HD*sqrt_HD) + HB/sqrt_HD + HC/HD - K;
+ const fp sqrt_Theta_D = sqrt(Theta_D);
+ Theta = + Theta_A/(Theta_D*sqrt_Theta_D)
+ + Theta_B/sqrt_Theta_D
+ + Theta_C/Theta_D
+ - K;
- // update running norms of H(h) function
- if (H_norms_ptr != NULL)
- then H_norms_ptr->data(H);
+ // update running norms of Theta(h) function
+ if (Theta_norms_ptr != NULL)
+ then Theta_norms_ptr->data(Theta);
if (Jacobian_flag)
then {
- // partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
- #include "../gr.cg/horizon_Jacobian.c"
+ // partial_Theta_wrt_partial_d_h,
+ // partial_Theta_wrt_partial_dd_h
+ #include "../gr.cg/expansion_Jacobian.c"
}
}
}
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index 2d4e8e3..c9cb8bc 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -1,14 +1,14 @@
-// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
+// expansion_Jacobian.cc -- evaluate Jacobian matrix of LHS function Theta(h)
// $Header$
//
// <<<prototypes for functions local to this file>>>
//
-// horizon_Jacobian - top-level driver to compute the Jacobian
+// expansion_Jacobian - top-level driver to compute the Jacobian
// decode_geometry_method - decode the geometry_method parameter
///
-/// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
+/// expansion_Jacobian_NP - compute the Jacobian by numerical perturbation
///
-/// horizon_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff
+/// expansion_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff
/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
///
@@ -51,58 +51,58 @@ using jtutil::error_exit;
//
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);
+bool expansion_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 expansion_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);
+bool expansion_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.
+// J[Theta(h)] of the expansion Theta(h). It just decodes the
+// Jacobian_method parameter and calls the appropriate subfunction.
//
-// We assume that H(h) has already been computed.
+// We assume that Theta(h) has already been computed.
//
// Results:
// This function returns true for a successful computation, or false
-// for failure. See horizon_function() (in "horizon_function.cc")
+// for failure. See expansion() (in "expansion.cc")
// for details of the 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)
+bool expansion_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,
+ if (! expansion_Jacobian_NP(ps, Jac,
cgi, gi, Jacobian_info,
print_msg_flag))
then return false; // *** ERROR RETURN ***
@@ -110,22 +110,22 @@ case Jacobian_method__numerical_perturb:
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\"");
+" expansion_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))
+ expansion_Jacobian_partial_SD(ps, Jac,
+ cgi, gi, Jacobian_info,
+ print_msg_flag);
+ if (! expansion_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"
+"***** expansion_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n"
" (this should never happen!)\n"
,
int(Jacobian_info.Jacobian_method)); /*NOTREACHED*/
@@ -161,51 +161,52 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
//******************************************************************************
//
-// This function computes the Jacobian matrix of the LHS function H(h)
-// by numerical perturbation of the H(h) function. The algorithm is as
+// This function computes the Jacobian matrix of the expansion Theta(h)
+// by numerical perturbation of the Theta(h) function. The algorithm is as
// follows:
//
-// we assume that H = H(h) has already been evaluated
-// save_H = H
+// we assume that Theta = Theta(h) has already been evaluated
+// save_Theta = Theta
// for each point (y,JJ)
// {
// const fp save_h_y = h at y;
// h at y += perturbation_amplitude;
-// evaluate H(h) (silently)
+// evaluate Theta(h) (silently)
// for each point (x,II)
// {
-// Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude;
+// Jac(II,JJ) = (Theta(II) - save_Theta(II))
+// / perturbation_amplitude;
// }
// h at y = save_h_y;
// }
-// H = save_H
+// Theta = save_Theta
//
// Inputs (angular gridfns, on ghosted grid):
-// h # shape of trial surface
-// H # H(h) assumed to already be computed
+// h # shape of trial surface
+// Theta # Theta(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_function() (in "horizon_function.cc")
+// for failure. See expansion() (in "expansion.cc")
// for details of the 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)
+bool expansion_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);
+ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
{
@@ -225,7 +226,7 @@ ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
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))
+ if (! expansion(ps, cgi, ii))
then return false; // *** ERROR RETURN ***
for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
{
@@ -240,9 +241,11 @@ ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
++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;
+ const fp old_Theta = xp.gridfn(gfns::gfn__save_Theta,
+ x_irho,x_isigma);
+ const fp new_Theta = xp.gridfn(gfns::gfn__Theta,
+ x_irho,x_isigma);
+ Jac(II,JJ) = (new_Theta - old_Theta) / epsilon;
}
}
}
@@ -251,7 +254,7 @@ ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
}
}
-ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
+ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
return true; // *** NORMAL RETURN ***
}
}
@@ -262,27 +265,27 @@ 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
+// matrix of the expansion Theta(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)
+// h # shape of trial surface
+// Theta # Theta(h) assumed to already be computed
+// partial_Theta_wrt_partial_d_h # Jacobian coefficients
+// partial_Theta_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)
+void expansion_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,
@@ -303,7 +306,7 @@ ps.compute_synchronize_Jacobian();
{
//
// compute the main Jacobian terms for this grid point, i.e.
- // partial H(this point x, Jacobian row II)
+ // partial Theta(this point x, Jacobian row II)
// ---------------------------------------------
// partial h(other points y, Jacobian column JJ)
//
@@ -313,19 +316,19 @@ ps.compute_synchronize_Jacobian();
// Jacobian coefficients for this point
const fp Jacobian_coeff_rho
- = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1,
+ = xp.gridfn(gfns::gfn__partial_Theta_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,
+ = xp.gridfn(gfns::gfn__partial_Theta_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,
+ = xp.gridfn(gfns::gfn__partial_Theta_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,
+ = xp.gridfn(gfns::gfn__partial_Theta_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,
+ = xp.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_22,
x_irho, x_isigma);
// partial_rho, partial_rho_rho
@@ -471,14 +474,14 @@ const patch& yp = ye.my_patch();
//
// This function sums the d/dr terms into the Jacobian matrix of the
-// LHS function H(h), computing those terms by finite differencing.
+// expansion Theta(h), computing those terms by finite differencing.
//
// The basic algorithm is that
-// Jac += diag[ (H(h+epsilon) - H(h)) / epsilon ]
+// Jac += diag[ (Theta(h+epsilon) - Theta(h)) / epsilon ]
//
// Inputs (angular gridfns, on ghosted grid):
-// h # shape of trial surface
-// H # H(h) assumed to already be computed
+// h # shape of trial surface
+// Theta # Theta(h) assumed to already be computed
//
// Outputs:
// Jac += d/dr terms
@@ -490,12 +493,12 @@ const patch& yp = ye.my_patch();
// 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)
+bool expansion_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,
@@ -503,10 +506,10 @@ if (print_msg_flag)
const fp epsilon = Jacobian_info.perturbation_amplitude;
-// compute H(h+epsilon)
-ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
+// compute Theta(h+epsilon)
+ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
-if (! horizon_function(ps, cgi, gi))
+if (! expansion(ps, cgi, gi))
then return false; // *** ERROR RETURN ***
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
@@ -519,16 +522,18 @@ if (! horizon_function(ps, cgi, gi))
++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;
+ const fp old_Theta = p.gridfn(gfns::gfn__save_Theta,
+ irho,isigma);
+ const fp new_Theta = p.gridfn(gfns::gfn__Theta,
+ irho,isigma);
+ Jac(II,II) += (new_Theta - old_Theta) / epsilon;
}
}
}
-// restore h and H
+// restore h and Theta
ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
-ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
+ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
return true; // *** NORMAL RETURN ***
}
diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh
index 5c3f6d3..f64fa56 100644
--- a/src/gr/gfns.hh
+++ b/src/gr/gfns.hh
@@ -60,14 +60,14 @@ enum {
gfn__partial_d_psi_2, // no access macro
gfn__partial_d_psi_3, // no access macro
- gfn__H,
- gfn__partial_H_wrt_partial_d_h_1,
- gfn__partial_H_wrt_partial_d_h_2,
- gfn__partial_H_wrt_partial_dd_h_11,
- gfn__partial_H_wrt_partial_dd_h_12,
- gfn__partial_H_wrt_partial_dd_h_22,
+ gfn__Theta,
+ gfn__partial_Theta_wrt_partial_d_h_1,
+ gfn__partial_Theta_wrt_partial_d_h_2,
+ gfn__partial_Theta_wrt_partial_dd_h_11,
+ gfn__partial_Theta_wrt_partial_dd_h_12,
+ gfn__partial_Theta_wrt_partial_dd_h_22,
gfn__Delta_h,
- gfn__save_H,
+ gfn__save_Theta,
gfn__one,
nominal_max_gfn = gfn__one // no comma
};
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 45771ef..8ba54f8 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -129,12 +129,12 @@ struct Jacobian_info
// prototypes for functions visible outside their source files
//
-// horizon_function.cc
+// expansion.cc
// ... returns true for successful computation,
// ... returns false for failure (eg horizon outside grid);
// if (print_msg_flag) then also reports CCTK_VWarn() at the
// appropriate warning level
-namespace horizon_function_warning_levels
+namespace expansion_warning_levels
{
enum {
failure = 2, // point outside grid etc
@@ -143,23 +143,23 @@ namespace horizon_function_warning_levels
// ==> skipping nonfinite check
};
}
-bool horizon_function(patch_system& ps,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- bool Jacobian_flag = false,
- bool print_msg_flag = false,
- jtutil::norm<fp>* H_norms_ptr = NULL);
+bool expansion(patch_system& ps,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ bool Jacobian_flag = false,
+ bool print_msg_flag = false,
+ jtutil::norm<fp>* H_norms_ptr = NULL);
enum geometry_method
decode_geometry_method(const char geometry_method_string[]);
-// horizon_Jacobian.cc
+// expansion_Jacobian.cc
// ... returns true for successful computation, false for failure
-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);
+bool expansion_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);
enum Jacobian_method
decode_Jacobian_method(const char Jacobian_method_string[]);
diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc
index 7579e4d..34da085 100644
--- a/src/gr/gr_gfas.minc
+++ b/src/gr/gr_gfas.minc
@@ -21,24 +21,24 @@ partial_d_s_d, partial_d_s_d__fnd,
n_u, n_u__fnd,
-HA, HA__fnd,
-HB, HB__fnd,
-HC, HC__fnd,
-HD, HD__fnd,
-H, H__fnd,
-
-partial_d_HA, partial_d_HA__fnd,
-partial_d_HB, partial_d_HB__fnd,
-partial_d_HC, partial_d_HC__fnd,
-partial_d_HD, partial_d_HD__fnd,
-partial_d_H, partial_d_H__fnd,
-
-partial_HA_wrt_partial_d_h, partial_HA_wrt_partial_d_h__fnd,
-partial_HB_wrt_partial_d_h, partial_HB_wrt_partial_d_h__fnd,
-partial_HC_wrt_partial_d_h, partial_HC_wrt_partial_d_h__fnd,
-partial_HD_wrt_partial_d_h, partial_HD_wrt_partial_d_h__fnd,
-partial_HA_wrt_partial_dd_h, partial_HA_wrt_partial_dd_h__fnd,
-partial_HB_wrt_partial_dd_h, partial_HB_wrt_partial_dd_h__fnd,
-
-partial_H_wrt_partial_d_h, partial_H_wrt_partial_d_h__fnd,
-partial_H_wrt_partial_dd_h, partial_H_wrt_partial_dd_h__fnd # no comma
+Theta_A, Theta_A__fnd,
+Theta_B, Theta_B__fnd,
+Theta_C, Theta_C__fnd,
+Theta_D, Theta_D__fnd,
+Theta, Theta__fnd,
+
+partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd,
+partial_d_Theta_C, partial_d_Theta_C__fnd,
+partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd,
+
+partial_Theta_A_wrt_partial_d_h, partial_Theta_A_wrt_partial_d_h__fnd,
+partial_Theta_B_wrt_partial_d_h, partial_Theta_B_wrt_partial_d_h__fnd,
+partial_Theta_C_wrt_partial_d_h, partial_Theta_C_wrt_partial_d_h__fnd,
+partial_Theta_D_wrt_partial_d_h, partial_Theta_D_wrt_partial_d_h__fnd,
+partial_Theta_A_wrt_partial_dd_h, partial_Theta_A_wrt_partial_dd_h__fnd,
+partial_Theta_B_wrt_partial_dd_h, partial_Theta_B_wrt_partial_dd_h__fnd,
+
+partial_Theta_wrt_partial_d_h, partial_Theta_wrt_partial_d_h__fnd,
+partial_Theta_wrt_partial_dd_h, partial_Theta_wrt_partial_dd_h__fnd # no comma
diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple
index 433c723..7284253 100644
--- a/src/gr/horizon.maple
+++ b/src/gr/horizon.maple
@@ -4,8 +4,8 @@
# horizon - overall driver
## non_unit_normal - compute s_d
## non_unit_normal_deriv - compute partial_d_s_d
-## horizon_function - compute H[ABCD],H = fn(s_d__fnd) = fn(h)
-## horizon_Jacobian - compute Jacobian coeffs for H[ABCD], H
+## expansion - compute Theta_[ABCD],Theta = fn(s_d__fnd) = fn(h)
+## expansion_Jacobian - compute Jacobian coeffs for Theta_[ABCD], Theta
#
################################################################################
@@ -20,15 +20,17 @@
# h
#
# Outputs:
-# H = H__fnd(h, ...)
+# Theta = Theta__fnd(h, ...)
+# partial_Theta_wrt_partial_d_h__fnd = partial_Theta_wrt_partial_d_h(h, ...)
+# partial_Theta_wrt_partial_dd_h__fnd = partial_Theta_wrt_partial_dd_h(h, ...)
#
horizon :=
proc(cg_flag::boolean)
non_unit_normal();
non_unit_normal_deriv();
-horizon_function(cg_flag);
-horizon_Jacobian(cg_flag);
+expansion(cg_flag);
+expansion_Jacobian(cg_flag);
NULL;
end proc;
@@ -124,22 +126,21 @@ end proc;
################################################################################
#
-# This function computes the LHS of the apparent horizon equation as a
-# function of 1st and 2nd angular partial derivatives of the apparent
-# horizon radius $h$, using equations (14) and (15[abcd]) in my 1996
-# apparent horizon finding paper, but using partial_d_g_uu[k,i,j]
-# instead of Diff(g_uu[i,j], x_xyz[k]).
+# This function computes the expansion $\Theta$ of a trial horizon surface
+# as a function of 1st and 2nd angular partial derivatives of the surface
+# radius $h$, using equations (14) and (15[abcd]) in my 1996 AH-finding
+# paper, but using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]).
#
-# These equations give H[ABCD] and H as functions of s_d, but here we
-# use s_d__fnd, which we assume is already given in terms of angular
-# derivatives of h. The result is that we compute H[ABCD] and H directly
-# in terms of 1st and 2nd angular derivatives of h, without s_d ever
-# appearing in our final results.
+# These equations give Theta_[ABCD] and Theta as functions of s_d, but
+# here we use s_d__fnd, which we assume is already given in terms of angular
+# derivatives of h. The result is that we compute Theta_[ABCD] and Theta
+# directly in terms of 1st and 2nd angular derivatives of h, without s_d
+# ever appearing in our final results.
#
# This function also optionally generates C code for the computation
-# of H[ABCD]. It does *not* compute C code for the computation of H
-# itself, since we may want to check that HD > 0 in the C code before
-# we compute H itself.
+# of Theta_[ABCD]. It does *not* compute C code for the computation of
+# Theta itself, since we may want to check that Theta_D > 0 in the C
+# code before we compute Theta itself.
#
# Inputs:
# s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) )
@@ -147,17 +148,17 @@ end proc;
# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
#
# Outputs:
-# HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HB = HB__fnd( x_xyz, X_ud, X_udd,
-# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
-# HD = HD__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
+# Theta_A = Theta_A__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_B = Theta_B__fnd( x_xyz, X_ud, X_udd,
+# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_C = Theta_C__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
+# Theta_D = Theta_D__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
# ---------------------------------------------------------------
-# H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd )
+# Theta = Theta__fnd(Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd)
#
-horizon_function :=
+expansion :=
proc(cg_flag::boolean)
global
@include "../maple/coords.minc",
@@ -180,33 +181,38 @@ assert_fnd_exists(partial_d_s_d, fnd);
#
# (15a) in my paper
-HA__fnd := - msum('g_uu[i,k]*s_d__fnd[k] * g_uu[j,l]*s_d__fnd[l]
- * partial_d_s_d__fnd[i,j]',
- 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N)
- - (1/2)*msum('g_uu[i,j]*s_d__fnd[j]
- * partial_d_g_uu[i,k,l] * s_d__fnd[k]*s_d__fnd[l]',
- 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N);
+Theta_A__fnd
+ := - msum('g_uu[i,k]*s_d__fnd[k] * g_uu[j,l]*s_d__fnd[l]
+ * partial_d_s_d__fnd[i,j]',
+ 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N)
+ - (1/2)*msum('g_uu[i,j]*s_d__fnd[j]
+ * partial_d_g_uu[i,k,l] * s_d__fnd[k]*s_d__fnd[l]',
+ 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N);
# (15b) in my paper
-HB__fnd := + msum('partial_d_g_uu[i,i,j]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N)
- + msum('g_uu[i,j]*partial_d_s_d__fnd[i,j]', 'i'=1..N, 'j'=1..N)
- + msum('g_uu[i,j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]',
- 'i'=1..N, 'j'=1..N);
+Theta_B__fnd
+ := + msum('partial_d_g_uu[i,i,j]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N)
+ + msum('g_uu[i,j]*partial_d_s_d__fnd[i,j]', 'i'=1..N, 'j'=1..N)
+ + msum('g_uu[i,j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]',
+ 'i'=1..N, 'j'=1..N);
# (15c) in my paper
-HC__fnd := msum('K_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
+Theta_C__fnd := msum('K_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
# (15d) in my paper
-HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
+Theta_D__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
# n.b. no simplify() here, it would try to put things over a common
-# denominator, which would make the equation much more complicated
-H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K;
+# denominator, which would make the equation *much* more complicated
+Theta__fnd := + Theta_A__fnd/Theta_D__fnd^(3/2)
+ + Theta_B__fnd/Theta_D__fnd^(1/2)
+ + Theta_C__fnd/Theta_D__fnd
+ - K;
if (cg_flag)
- then codegen2([ HA__fnd, HB__fnd, HC__fnd, HD__fnd],
- ['HA', 'HB', 'HC', 'HD' ],
- "../gr.cg/horizon_function.c");
+ then codegen2([ Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd],
+ ['Theta_A', 'Theta_B', 'Theta_C', 'Theta_D' ],
+ "../gr.cg/expansion.c");
fi;
NULL;
@@ -215,12 +221,12 @@ end proc;
################################################################################
#
-# This function computes the Jacobian coefficients for the LHS of the
-# apparent horizon equation with respect to 1st and 2nd angular partial
-# derivatives of the apparent horizon radius $h$. These coefficients
+# This function computes the Jacobian coefficients for the expansion
+# $\Theta$ of a trial horizon surface with respect to 1st and 2nd angular
+# partial derivatives of the surface radius $h$. These coefficients
# are (or at least should be :) the same as those in equation (A1) in
# my 1996 apparent horizon finding paper, bue here we compute them in
-# a different manner: we have Maple directly differentiate H__fnd with
+# a different manner: we have Maple directly differentiate Theta__fnd with
# respect to Diff(h,y_rs[u]) and Diff(h,y_rs[u],y_rs[v]). We use the
# Maple frontend() function to do this.
#
@@ -228,24 +234,24 @@ end proc;
# coefficients.
#
# Inputs:
-# HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HB = HB__fnd( x_xyz, X_ud, X_udd,
-# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
-# HD = HD__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
-# H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd )
+# Theta_A = Theta_A__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_B = Theta_B__fnd( x_xyz, X_ud, X_udd,
+# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_C = Theta_C__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
+# Theta_D = Theta_D__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
+# Theta = Theta__fnd(Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd)
#
# Outputs:
-# partial_HA_wrt_partial_d_h partial_HA_wrt_partial_dd_h
-# partial_HB_wrt_partial_d_h partial_HB_wrt_partial_dd_h
-# partial_HC_wrt_partial_d_h
-# partial_HD_wrt_partial_d_h
+# partial_Theta_A_wrt_partial_d_h partial_Theta_A_wrt_partial_dd_h
+# partial_Theta_B_wrt_partial_d_h partial_Theta_B_wrt_partial_dd_h
+# partial_Theta_C_wrt_partial_d_h
+# partial_Theta_D_wrt_partial_d_h
# ---------------------------------------------------------------
-# partial_H_wrt_partial_d_h partial_H_wrt_partial_dd_h
+# partial_Theta_wrt_partial_d_h partial_Theta_wrt_partial_dd_h
#
-horizon_Jacobian :=
+expansion_Jacobian :=
proc(cg_flag::boolean)
global
@include "../maple/coords.minc",
@@ -260,54 +266,57 @@ assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
assert_fnd_exists(partial_d_ln_sqrt_g);
assert_fnd_exists(partial_d_g_uu);
-assert_fnd_exists(HA);
-assert_fnd_exists(HB);
-assert_fnd_exists(HC);
-assert_fnd_exists(HD);
+assert_fnd_exists(Theta_A);
+assert_fnd_exists(Theta_B);
+assert_fnd_exists(Theta_C);
+assert_fnd_exists(Theta_D);
-# Jacobian coefficients of H[ABCD] and H wrt Diff(h,y_rs[u])
+# Jacobian coefficients of Theta_[ABCD] and Theta wrt Diff(h,y_rs[u])
for u from 1 to N_ang
do
- partial_HA_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HA__fnd, Diff(h,y_rs[u])]);
- partial_HB_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HB__fnd, Diff(h,y_rs[u])]);
- partial_HC_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HC__fnd, Diff(h,y_rs[u])]);
- partial_HD_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HD__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_A_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_A__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_B_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_B__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_C_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_C__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_D_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_D__fnd, Diff(h,y_rs[u])]);
# equation (A1a) in my 1996 apparent horizon finding paper
- temp := (3/2)*HA/HD^(5/2) + (1/2)*HB/HD^(3/2) + HC/HD^2;
- partial_H_wrt_partial_d_h__fnd[u]
- := + partial_HA_wrt_partial_d_h__fnd[u] / HD^(3/2)
- + partial_HB_wrt_partial_d_h__fnd[u] / HD^(1/2)
- + partial_HC_wrt_partial_d_h__fnd[u] / HD
- - partial_HD_wrt_partial_d_h__fnd[u] * temp;
+ temp := + (3/2)*Theta_A/Theta_D^(5/2)
+ + (1/2)*Theta_B/Theta_D^(3/2)
+ + Theta_C/Theta_D^2;
+ partial_Theta_wrt_partial_d_h__fnd[u]
+ := + partial_Theta_A_wrt_partial_d_h__fnd[u] / Theta_D^(3/2)
+ + partial_Theta_B_wrt_partial_d_h__fnd[u] / Theta_D^(1/2)
+ + partial_Theta_C_wrt_partial_d_h__fnd[u] / Theta_D
+ - partial_Theta_D_wrt_partial_d_h__fnd[u] * temp;
end do;
-# Jacobian coefficients of H[AB] and H wrt Diff(h,y_rs[u],y_rs[v])
+# Jacobian coefficients of Theta_[AB] and Theta wrt Diff(h,y_rs[u],y_rs[v])
for u from 1 to N_ang
do
for v from u to N_ang
do
- partial_HA_wrt_partial_dd_h__fnd[u,v]
- := frontend('diff', [HA__fnd, Diff(h,y_rs[u],y_rs[v])]);
- partial_HB_wrt_partial_dd_h__fnd[u,v]
- := frontend('diff', [HB__fnd, Diff(h,y_rs[u],y_rs[v])]);
+ partial_Theta_A_wrt_partial_dd_h__fnd[u,v]
+ := frontend('diff', [Theta_A__fnd, Diff(h,y_rs[u],y_rs[v])]);
+ partial_Theta_B_wrt_partial_dd_h__fnd[u,v]
+ := frontend('diff', [Theta_B__fnd, Diff(h,y_rs[u],y_rs[v])]);
# equation (A1b) in my 1996 apparent horizon finding paper
- partial_H_wrt_partial_dd_h__fnd[u,v]
- := + partial_HA_wrt_partial_dd_h__fnd[u,v] / HD^(3/2)
- + partial_HB_wrt_partial_dd_h__fnd[u,v] / HD^(1/2)
+ partial_Theta_wrt_partial_dd_h__fnd[u,v]
+ := + partial_Theta_A_wrt_partial_dd_h__fnd[u,v] / Theta_D^(3/2)
+ + partial_Theta_B_wrt_partial_dd_h__fnd[u,v] / Theta_D^(1/2)
end do;
end do;
if (cg_flag)
- then codegen2([partial_H_wrt_partial_d_h__fnd,
- partial_H_wrt_partial_dd_h__fnd],
- ['partial_H_wrt_partial_d_h', 'partial_H_wrt_partial_dd_h'],
- "../gr.cg/horizon_Jacobian.c");
+ then codegen2([partial_Theta_wrt_partial_d_h__fnd,
+ partial_Theta_wrt_partial_dd_h__fnd],
+ ['partial_Theta_wrt_partial_d_h',
+ 'partial_Theta_wrt_partial_dd_h'],
+ "../gr.cg/expansion_Jacobian.c");
fi;
NULL;
diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn
index a371c50..e93414f 100644
--- a/src/gr/make.code.defn
+++ b/src/gr/make.code.defn
@@ -2,8 +2,8 @@
# $Header$
# Source files in this directory
-SRCS = horizon_function.cc \
- horizon_Jacobian.cc \
+SRCS = expansion.cc \
+ expansion_Jacobian.cc \
Schwarzschild_EF.cc
# Subdirectories containing source files
diff --git a/src/gr/maple.log b/src/gr/maple.log
index 0a28bb5..54c762f 100644
--- a/src/gr/maple.log
+++ b/src/gr/maple.log
@@ -586,19 +586,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
make_gfa('g_dd', {inert}, [1 .. N, 1 .. N], symmetric);
make_gfa('K_dd', {inert}, [1 .. N, 1 .. N], symmetric);
make_gfa('g_uu', {inert, fnd}, [1 .. N, 1 .. N], symmetric);
@@ -612,31 +613,31 @@ partial_H_wrt_partial_dd_h__fnd;
make_gfa('h', {inert, fnd}, [], none);
make_gfa('s_d', {inert, fnd}, [1 .. N], none);
make_gfa('partial_d_s_d', {inert, fnd}, [1 .. N, 1 .. N], none);
- make_gfa('HA', {inert, fnd}, [], none);
- make_gfa('HB', {inert, fnd}, [], none);
- make_gfa('HC', {inert, fnd}, [], none);
- make_gfa('HD', {inert, fnd}, [], none);
- make_gfa('H', {inert, fnd}, [], none);
- make_gfa('partial_d_HA', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_HB', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_HC', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_HD', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_H', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_HA_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HB_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HC_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HD_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HA_wrt_partial_dd_h', {inert, fnd},
+ make_gfa('Theta_A', {inert, fnd}, [], none);
+ make_gfa('Theta_B', {inert, fnd}, [], none);
+ make_gfa('Theta_C', {inert, fnd}, [], none);
+ make_gfa('Theta_D', {inert, fnd}, [], none);
+ make_gfa('Theta', {inert, fnd}, [], none);
+ make_gfa('partial_d_Theta_A', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta_B', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta_C', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta_D', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_Theta_A_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_B_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_C_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_D_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_A_wrt_partial_dd_h', {inert, fnd},
[1 .. N_ang, 1 .. N_ang], symmetric);
- make_gfa('partial_HB_wrt_partial_dd_h', {inert, fnd},
+ make_gfa('partial_Theta_B_wrt_partial_dd_h', {inert, fnd},
[1 .. N_ang, 1 .. N_ang], symmetric);
- make_gfa('partial_H_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_H_wrt_partial_dd_h', {inert, fnd},
+ make_gfa('partial_Theta_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_wrt_partial_dd_h', {inert, fnd},
[1 .. N_ang, 1 .. N_ang], symmetric);
NULL
end proc
@@ -648,19 +649,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
option remember;
var_list := [args[2 .. nargs]];
if type(operand, indexed) and op(0, operand) = 'g_dd' and
@@ -683,19 +685,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_dd);
assert_fnd_exists(g_uu, fnd);
@@ -713,19 +716,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_uu);
assert_fnd_exists(K_dd);
@@ -756,19 +760,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_dd);
assert_fnd_exists(g_uu);
@@ -793,19 +798,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_dd);
assert_fnd_exists(g_uu);
@@ -823,8 +829,8 @@ end proc
horizon := proc(cg_flag::boolean)
non_unit_normal();
non_unit_normal_deriv();
- horizon_function(cg_flag);
- horizon_Jacobian(cg_flag);
+ expansion(cg_flag);
+ expansion_Jacobian(cg_flag);
NULL
end proc
@@ -835,19 +841,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(h);
assert_fnd_exists(X_ud);
@@ -865,19 +872,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(h);
assert_fnd_exists(X_ud);
@@ -897,26 +905,27 @@ partial_H_wrt_partial_dd_h__fnd;
NULL
end proc
-horizon_function := proc(cg_flag::boolean)
+expansion := proc(cg_flag::boolean)
local i, j, k, l;
global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
@@ -924,88 +933,94 @@ partial_H_wrt_partial_dd_h__fnd;
assert_fnd_exists(partial_d_g_uu);
assert_fnd_exists(s_d, fnd);
assert_fnd_exists(partial_d_s_d, fnd);
- HA__fnd := -msum('g_uu[i, k]*s_d__fnd[k]*g_uu[j, l]*s_d__fnd[l]*
+ Theta_A__fnd := -msum('g_uu[i, k]*s_d__fnd[k]*g_uu[j, l]*s_d__fnd[l]*
partial_d_s_d__fnd[i, j]', 'i' = 1 .. N, 'j' = 1 .. N, 'k' = 1 .. N,
'l' = 1 .. N) - 1/2*msum('g_uu[i, j]*s_d__fnd[j]*
partial_d_g_uu[i, k, l]*s_d__fnd[k]*s_d__fnd[l]', 'i' = 1 .. N,
'j' = 1 .. N, 'k' = 1 .. N, 'l' = 1 .. N);
- HB__fnd := msum('partial_d_g_uu[i, i, j]*s_d__fnd[j]', 'i' = 1 .. N,
- 'j' = 1 .. N) + msum('g_uu[i, j]*partial_d_s_d__fnd[i, j]',
+ Theta_B__fnd := msum('partial_d_g_uu[i, i, j]*s_d__fnd[j]',
'i' = 1 .. N, 'j' = 1 .. N) + msum(
- 'g_uu[i, j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]', 'i' = 1 .. N,
+ 'g_uu[i, j]*partial_d_s_d__fnd[i, j]', 'i' = 1 .. N, 'j' = 1 .. N)
+ + msum('g_uu[i, j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]',
+ 'i' = 1 .. N, 'j' = 1 .. N);
+ Theta_C__fnd := msum('K_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
'j' = 1 .. N);
- HC__fnd := msum('K_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
+ Theta_D__fnd := msum('g_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
'j' = 1 .. N);
- HD__fnd := msum('g_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
- 'j' = 1 .. N);
- H__fnd :=
- HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K
- ;
- if cg_flag then codegen2([HA__fnd, HB__fnd, HC__fnd, HD__fnd],
- ['HA', 'HB', 'HC', 'HD'], "../gr.cg/horizon_function.c")
+ Theta__fnd := Theta_A__fnd/Theta_D__fnd^(3/2)
+ + Theta_B__fnd/Theta_D__fnd^(1/2) + Theta_C__fnd/Theta_D__fnd - K;
+ if cg_flag then codegen2(
+ [Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd],
+ ['Theta_A', 'Theta_B', 'Theta_C', 'Theta_D'],
+ "../gr.cg/expansion.c")
end if;
NULL
end proc
-horizon_Jacobian := proc(cg_flag::boolean)
+expansion_Jacobian := proc(cg_flag::boolean)
local u, v, temp;
global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
assert_fnd_exists(partial_d_ln_sqrt_g);
assert_fnd_exists(partial_d_g_uu);
- assert_fnd_exists(HA);
- assert_fnd_exists(HB);
- assert_fnd_exists(HC);
- assert_fnd_exists(HD);
+ assert_fnd_exists(Theta_A);
+ assert_fnd_exists(Theta_B);
+ assert_fnd_exists(Theta_C);
+ assert_fnd_exists(Theta_D);
for u to N_ang do
- partial_HA_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HA__fnd, Diff(h, y_rs[u])]);
- partial_HB_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HB__fnd, Diff(h, y_rs[u])]);
- partial_HC_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HC__fnd, Diff(h, y_rs[u])]);
- partial_HD_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HD__fnd, Diff(h, y_rs[u])]);
- temp := 3/2*HA/HD^(5/2) + 1/2*HB/HD^(3/2) + HC/HD^2;
- partial_H_wrt_partial_d_h__fnd[u] :=
- partial_HA_wrt_partial_d_h__fnd[u]/HD^(3/2)
- + partial_HB_wrt_partial_d_h__fnd[u]/HD^(1/2)
- + partial_HC_wrt_partial_d_h__fnd[u]/HD
- - partial_HD_wrt_partial_d_h__fnd[u]*temp
+ partial_Theta_A_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_A__fnd, Diff(h, y_rs[u])]);
+ partial_Theta_B_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_B__fnd, Diff(h, y_rs[u])]);
+ partial_Theta_C_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_C__fnd, Diff(h, y_rs[u])]);
+ partial_Theta_D_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_D__fnd, Diff(h, y_rs[u])]);
+ temp := 3/2*Theta_A/Theta_D^(5/2) + 1/2*Theta_B/Theta_D^(3/2)
+ + Theta_C/Theta_D^2;
+ partial_Theta_wrt_partial_d_h__fnd[u] :=
+ partial_Theta_A_wrt_partial_d_h__fnd[u]/Theta_D^(3/2)
+ + partial_Theta_B_wrt_partial_d_h__fnd[u]/Theta_D^(1/2)
+ + partial_Theta_C_wrt_partial_d_h__fnd[u]/Theta_D
+ - partial_Theta_D_wrt_partial_d_h__fnd[u]*temp
end do;
for u to N_ang do for v from u to N_ang do
- partial_HA_wrt_partial_dd_h__fnd[u, v] :=
- frontend('diff', [HA__fnd, Diff(h, y_rs[u], y_rs[v])]);
- partial_HB_wrt_partial_dd_h__fnd[u, v] :=
- frontend('diff', [HB__fnd, Diff(h, y_rs[u], y_rs[v])]);
- partial_H_wrt_partial_dd_h__fnd[u, v] :=
- partial_HA_wrt_partial_dd_h__fnd[u, v]/HD^(3/2)
- + partial_HB_wrt_partial_dd_h__fnd[u, v]/HD^(1/2)
+ partial_Theta_A_wrt_partial_dd_h__fnd[u, v] :=
+ frontend('diff', [Theta_A__fnd, Diff(h, y_rs[u], y_rs[v])])
+ ;
+ partial_Theta_B_wrt_partial_dd_h__fnd[u, v] :=
+ frontend('diff', [Theta_B__fnd, Diff(h, y_rs[u], y_rs[v])])
+ ;
+ partial_Theta_wrt_partial_dd_h__fnd[u, v] :=
+ partial_Theta_A_wrt_partial_dd_h__fnd[u, v]/Theta_D^(3/2)
+ +
+ partial_Theta_B_wrt_partial_dd_h__fnd[u, v]/Theta_D^(1/2)
end do
end do;
- if cg_flag then codegen2(
- [partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h__fnd],
- ['partial_H_wrt_partial_d_h', 'partial_H_wrt_partial_dd_h'],
- "../gr.cg/horizon_Jacobian.c")
+ if cg_flag then codegen2([partial_Theta_wrt_partial_d_h__fnd,
+ partial_Theta_wrt_partial_dd_h__fnd],
+ ['partial_Theta_wrt_partial_d_h', 'partial_Theta_wrt_partial_dd_h']
+ , "../gr.cg/expansion_Jacobian.c")
end if;
NULL
end proc
@@ -1025,7 +1040,7 @@ codegen2(g_uu) --> "../gr.cg/inverse_metric.c"
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=1006876, alloc=917336, time=0.09
+bytes used=1000088, alloc=917336, time=0.12
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
--> `codegen2/unindex`
@@ -1038,7 +1053,7 @@ codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c"
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=2007252, alloc=1441528, time=0.17
+bytes used=2000352, alloc=1441528, time=0.17
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
@@ -1046,39 +1061,39 @@ bytes used=2007252, alloc=1441528, time=0.17
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
--> `codegen2/unindex`
-bytes used=3007472, alloc=1638100, time=0.25
+bytes used=3000584, alloc=1638100, time=0.22
convert p/q --> RATIONAL(p/q)
--> `codegen2/fix_rationals`
writing C code
> curvature(true);
inverse_metric_gradient...
-bytes used=4007696, alloc=1638100, time=0.33
+bytes used=4000860, alloc=1703624, time=0.28
codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c"
--> `codegen2/input`
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=5007880, alloc=1703624, time=0.43
+bytes used=5001020, alloc=1703624, time=0.38
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=6008088, alloc=1769148, time=0.49
+bytes used=6001324, alloc=1769148, time=0.44
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
--> `codegen2/unindex`
-bytes used=7008488, alloc=1769148, time=0.56
+bytes used=7001528, alloc=1769148, time=0.51
convert p/q --> RATIONAL(p/q)
--> `codegen2/fix_rationals`
writing C code
-bytes used=8008672, alloc=1769148, time=0.62
+bytes used=8001896, alloc=1769148, time=0.57
metric_det_gradient...
codegen2(partial_d_ln_sqrt_g) --> "../gr.cg/metric_det_gradient.c"
--> `codegen2/input`
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=9008904, alloc=1769148, time=0.75
+bytes used=9002080, alloc=1769148, time=0.64
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
@@ -1093,90 +1108,90 @@ codegen/C/expression: Unknown function: RATIONAL will be left as is.
> horizon(true);
non_unit_normal...
non_unit_normal_deriv...
-bytes used=10009180, alloc=1769148, time=0.83
-horizon_function...
-bytes used=11009748, alloc=1834672, time=0.91
-codegen2([HA, HB, HC, HD]) --> "../gr.cg/horizon_function.c"
+bytes used=10002316, alloc=1769148, time=0.69
+expansion...
+bytes used=11002748, alloc=1834672, time=0.77
+codegen2([Theta_A, Theta_B, Theta_C, Theta_D]) --> "../gr.cg/expansion.c"
--> `codegen2/input`
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=12010468, alloc=1834672, time=0.98
-bytes used=13010628, alloc=2031244, time=1.08
-bytes used=14010940, alloc=2031244, time=1.16
+bytes used=12003104, alloc=1834672, time=0.83
+bytes used=13003356, alloc=2031244, time=0.95
+bytes used=14003520, alloc=2031244, time=1.05
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=15011240, alloc=2096768, time=1.23
+bytes used=15003896, alloc=2096768, time=1.11
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
-bytes used=16011692, alloc=2096768, time=1.31
+bytes used=16004136, alloc=2096768, time=1.18
--> `codegen2/unindex`
-bytes used=17012040, alloc=2096768, time=1.37
+bytes used=17004388, alloc=2096768, time=1.25
convert p/q --> RATIONAL(p/q)
-bytes used=18012780, alloc=2096768, time=1.44
+bytes used=18004580, alloc=2096768, time=1.31
--> `codegen2/fix_rationals`
writing C code
codegen/C/expression: Unknown function: PARTIAL_RHO will be left as is.
codegen/C/expression: Unknown function: PARTIAL_SIGMA will be left as is.
-bytes used=19012960, alloc=2096768, time=1.52
codegen/C/expression: Unknown function: PARTIAL_RHO_RHO
will be left as is.
codegen/C/expression: Unknown function: PARTIAL_RHO_SIGMA
will be left as is.
codegen/C/expression: Unknown function: PARTIAL_SIGMA_SIGMA
will be left as is.
-bytes used=20013320, alloc=2096768, time=1.66
-horizon_Jacobian...
-bytes used=21013548, alloc=2096768, time=1.73
-codegen2([partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h]) --> "../gr.cg/horizon_Jacobian.c"
+bytes used=19004732, alloc=2096768, time=1.37
+bytes used=20004972, alloc=2096768, time=1.47
+expansion_Jacobian...
+bytes used=21005184, alloc=2096768, time=1.56
+codegen2([partial_Theta_wrt_partial_d_h, partial_Theta_wrt_partial_dd_h]) --> "../gr.cg/expansion_Jacobian.c"
--> `codegen2/input`
convert --> equation list
-bytes used=22013808, alloc=2096768, time=1.79
+bytes used=22005440, alloc=2096768, time=1.63
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=23014336, alloc=2096768, time=1.87
-bytes used=24014548, alloc=2293340, time=1.92
-bytes used=25014936, alloc=2293340, time=1.98
-bytes used=26015132, alloc=2293340, time=2.08
-bytes used=27015308, alloc=2293340, time=2.20
-bytes used=28015956, alloc=2293340, time=2.36
-bytes used=29016184, alloc=2293340, time=2.48
-bytes used=30016340, alloc=2424388, time=2.62
-bytes used=31016540, alloc=2424388, time=2.72
-bytes used=32016792, alloc=2555436, time=2.82
+bytes used=23005788, alloc=2162292, time=1.71
+bytes used=24006060, alloc=2293340, time=1.78
+bytes used=25006304, alloc=2293340, time=1.84
+bytes used=26006456, alloc=2293340, time=1.94
+bytes used=27007216, alloc=2293340, time=2.10
+bytes used=28007500, alloc=2293340, time=2.26
+bytes used=29007664, alloc=2293340, time=2.39
+bytes used=30007868, alloc=2424388, time=2.50
+bytes used=31008068, alloc=2424388, time=2.60
+bytes used=32009692, alloc=2555436, time=2.71
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=33016952, alloc=2555436, time=2.89
-bytes used=34017300, alloc=2555436, time=2.96
-bytes used=35017536, alloc=2555436, time=3.03
+bytes used=33010028, alloc=2555436, time=2.77
+bytes used=34010432, alloc=2555436, time=2.84
+bytes used=35010740, alloc=2555436, time=2.91
--> `codegen2/fix_Diff`
+bytes used=36021172, alloc=2555436, time=2.99
convert R_dd[2,3] --> R_dd_23 etc
-bytes used=36018544, alloc=2555436, time=3.10
-bytes used=37018944, alloc=2555436, time=3.17
-bytes used=38019240, alloc=2555436, time=3.24
-bytes used=39019452, alloc=2555436, time=3.31
+bytes used=37021516, alloc=2555436, time=3.05
+bytes used=38021792, alloc=2555436, time=3.12
+bytes used=39022208, alloc=2555436, time=3.20
--> `codegen2/unindex`
-bytes used=40019836, alloc=2555436, time=3.38
-bytes used=41019988, alloc=2555436, time=3.45
-bytes used=42020192, alloc=2555436, time=3.51
-bytes used=43020452, alloc=2555436, time=3.58
+bytes used=40022676, alloc=2555436, time=3.27
+bytes used=41022924, alloc=2555436, time=3.33
+bytes used=42023096, alloc=2555436, time=3.40
+bytes used=43023256, alloc=2555436, time=3.47
+bytes used=44023688, alloc=2555436, time=3.53
convert p/q --> RATIONAL(p/q)
-bytes used=44021232, alloc=2555436, time=3.65
-bytes used=45021924, alloc=2555436, time=3.71
-bytes used=46022240, alloc=2555436, time=3.77
-bytes used=47023068, alloc=2555436, time=3.84
-bytes used=48023224, alloc=2555436, time=3.92
+bytes used=45023996, alloc=2555436, time=3.60
+bytes used=46024268, alloc=2555436, time=3.67
+bytes used=47024688, alloc=2555436, time=3.73
+bytes used=48025108, alloc=2555436, time=3.80
--> `codegen2/fix_rationals`
writing C code
-bytes used=49023460, alloc=2555436, time=3.99
-bytes used=50023664, alloc=2555436, time=4.04
-bytes used=51023844, alloc=2555436, time=4.20
-bytes used=52024240, alloc=2555436, time=4.32
-bytes used=53024400, alloc=2555436, time=4.48
-bytes used=54024568, alloc=2555436, time=4.58
+bytes used=49025260, alloc=2555436, time=3.86
+bytes used=50025744, alloc=2555436, time=3.91
+bytes used=51026008, alloc=2555436, time=4.03
+bytes used=52026168, alloc=2555436, time=4.18
+bytes used=53026328, alloc=2555436, time=4.32
+bytes used=54026584, alloc=2555436, time=4.44
> quit
-bytes used=54791456, alloc=2555436, time=4.67
+bytes used=54895808, alloc=2555436, time=4.56
diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple
index ff39681..6af4aed 100644
--- a/src/gr/setup_gr_gfas.maple
+++ b/src/gr/setup_gr_gfas.maple
@@ -44,39 +44,38 @@ make_gfa('h', {inert,fnd}, [], none);
make_gfa('s_d', {inert,fnd}, [1..N], none);
make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none);
-# subexpressions for computing LHS of apparent horizon equation
-# ... these are defined by (15) in my 1996 apparent horizon finding paper
-make_gfa('HA', {inert,fnd}, [], none);
-make_gfa('HB', {inert,fnd}, [], none);
-make_gfa('HC', {inert,fnd}, [], none);
-make_gfa('HD', {inert,fnd}, [], none);
-
-# LHS of apparent horizon equation
-make_gfa('H', {inert,fnd}, [], none);
-
-# 1st derivatives of H[ABCD] and of H
-make_gfa('partial_d_HA', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_HB', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_HC', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_HD', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_H', {inert,fnd}, [1..N], none);
-
-# Jacobian coefficients for H[ABCD]
-# these are the partial derivatives of H[ABCD]
+# expansion of horizon and subexpressions for computing it
+# ... these are defined by (14) and (15) in my 1996 AH-finding paper
+# except I now use Theta instead of H for the LHS
+make_gfa('Theta_A', {inert,fnd}, [], none);
+make_gfa('Theta_B', {inert,fnd}, [], none);
+make_gfa('Theta_C', {inert,fnd}, [], none);
+make_gfa('Theta_D', {inert,fnd}, [], none);
+make_gfa('Theta', {inert,fnd}, [], none);
+
+# 1st derivatives of Theta[ABCD] and of Theta
+make_gfa('partial_d_Theta_A', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta_B', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta_C', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta_D', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta', {inert,fnd}, [1..N], none);
+
+# Jacobian coefficients for Theta[ABCD]
+# these are the partial derivatives of Theta[ABCD]
# ... wrt Diff(h,x_rs[x])
-make_gfa('partial_HA_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_HB_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_HC_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_HD_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_A_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_B_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_C_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_D_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
# ... wrt Diff(h,x_rs[x],x_rs[y])
-make_gfa('partial_HA_wrt_partial_dd_h', {inert,fnd},
+make_gfa('partial_Theta_A_wrt_partial_dd_h', {inert,fnd},
[1..N_ang, 1..N_ang], symmetric);
-make_gfa('partial_HB_wrt_partial_dd_h', {inert,fnd},
+make_gfa('partial_Theta_B_wrt_partial_dd_h', {inert,fnd},
[1..N_ang, 1..N_ang], symmetric);
-# Jacobian coefficients for H itself
-make_gfa('partial_H_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_H_wrt_partial_dd_h', {inert,fnd},
+# Jacobian coefficients for Theta itself
+make_gfa('partial_Theta_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_wrt_partial_dd_h', {inert,fnd},
[1..N_ang, 1..N_ang], symmetric);
NULL;
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 00:15:22 +0000