diff options
Diffstat (limited to 'src/gr/horizon_function.cc')
| -rw-r--r-- | src/gr/horizon_function.cc | 58 |
1 files changed, 45 insertions, 13 deletions
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 6ea0a8d..745760b 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -57,7 +57,9 @@ void compute_H(patch_system& ps); //****************************************************************************** // -// This function computes the LHS function H(h). +// This function computes the LHS function H(h), and optionally also +// its Jacobian coefficients (from which the Jacobian matrix may be +// computed later). // // Inputs (angular gridfns, on ghosted grid): // ... defined on ghosted grid @@ -86,9 +88,14 @@ void compute_H(patch_system& ps); // partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k g_{ij}$ // H # $H = H(h)$ // +// Arguments: +// Jacobian_flag = true to compute the Jacobian coefficients, +// false to skip this. +// void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& ii) + const struct geometry_interpolator_info& ii, + bool Jacobian_flag) { CCTK_VInfo(CCTK_THORNSTRING, " horizon function"); @@ -101,9 +108,9 @@ setup_xyz_posns(ps); // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$ interpolate_geometry(ps, cgi, ii); -// compute remaining gridfns --> $H$ +// compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing -compute_H(ps); +compute_H(ps, Jacobian_flag); } //****************************************************************************** @@ -302,8 +309,9 @@ void* const output_arrays[N_OUTPUT_ARRAYS] bool first_time = true; if (first_time) then { - // store derivative info in interpolator parameter table first_time = false; + + // store derivative info in interpolator parameter table CCTK_VInfo(CCTK_THORNSTRING, " setting up derivative info for interpolator"); @@ -404,11 +412,17 @@ if (status < 0) //****************************************************************************** // -// This function computes H(h) by a mixture of algebraic operations -// and (rho,sigma) finite differencing, on the angular grid. +// 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. +// +// Arguments: +// Jacobian_flag = true to compute the Jacobian coefficients, +// false to skip this. // namespace { -void compute_H(patch_system& ps) +void compute_H(patch_system& ps, bool Jacobian_flag) { CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); @@ -460,22 +474,40 @@ CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); // ... each cg/*.c file has a separate set of temp variables, // and so must be inside its own set of { } braces // + + // gridfn #defins #include "cg.hh" + { - #include "cg/inverse_metric.c" + // g_uu + #include "../gr.cg/inverse_metric.c" } + { - #include "cg/extrinsic_curvature_trace_raise.c" + // K, K_uu + #include "../gr.cg/extrinsic_curvature_trace_raise.c" } + { - #include "cg/inverse_metric_gradient.c" + // partial_d_g_uu + #include "../gr.cg/inverse_metric_gradient.c" } + { - #include "cg/metric_det_gradient.c" + // partial_d_ln_sqrt_g + #include "../gr.cg/metric_det_gradient.c" } + { - #include "cg/horizon_function.c" + // HA, HB, HC, HD, H + #include "../gr.cg/horizon_function.c" } + + if (Jacobian_flag) + then { + // partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h + #include "../gr.cg/horizon_Jacobian.c" + } } } } |