diff options
Diffstat (limited to 'src/gr')
| -rw-r--r-- | src/gr/Schwarzschild_EF.cc | 3 | ||||
| -rw-r--r-- | src/gr/cg.hh | 3 | ||||
| -rw-r--r-- | src/gr/expansion.cc | 491 | ||||
| -rw-r--r-- | src/gr/expansion_Jacobian.cc | 415 | ||||
| -rw-r--r-- | src/gr/gfns.hh | 17 | ||||
| -rw-r--r-- | src/gr/gr.hh | 75 | ||||
| -rw-r--r-- | src/gr/gr_gfas.minc | 6 | ||||
| -rw-r--r-- | src/gr/horizon.maple | 34 | ||||
| -rw-r--r-- | src/gr/make.code.defn | 9 | ||||
| -rw-r--r-- | src/gr/maple.log | 265 | ||||
| -rw-r--r-- | src/gr/misc-gr.cc | 1 | ||||
| -rw-r--r-- | src/gr/setup_gr_gfas.maple | 7 |
12 files changed, 1109 insertions, 217 deletions
diff --git a/src/gr/Schwarzschild_EF.cc b/src/gr/Schwarzschild_EF.cc index b27d5a4..63bc01f 100644 --- a/src/gr/Schwarzschild_EF.cc +++ b/src/gr/Schwarzschild_EF.cc @@ -16,6 +16,7 @@ #include "util_Table.h" #include "cctk.h" +#include "cctk_Arguments.h" #include "config.h" #include "stdc.h" @@ -23,6 +24,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "../patch/coords.hh" #include "../patch/grid.hh" @@ -41,7 +43,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** diff --git a/src/gr/cg.hh b/src/gr/cg.hh index 50d1468..da331ac 100644 --- a/src/gr/cg.hh +++ b/src/gr/cg.hh @@ -106,7 +106,8 @@ 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) +#define old_Theta p.gridfn(gfns::gfn__old_Theta, irho,isigma) +#define Delta_h p.gridfn(gfns::gfn__Delta_h, irho,isigma) //****************************************************************************** diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc index 5a03783..cdbd57c 100644 --- a/src/gr/expansion.cc +++ b/src/gr/expansion.cc @@ -37,9 +37,11 @@ #include <stdio.h> #include <assert.h> #include <math.h> +#include <string.h> #include "util_Table.h" #include "cctk.h" +#include "cctk_Arguments.h" #include "config.h" #include "stdc.h" @@ -47,6 +49,10 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; +using jtutil::pow2; +using jtutil::pow3; +using jtutil::pow4; #include "../patch/coords.hh" #include "../patch/grid.hh" @@ -65,9 +71,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; -using jtutil::pow2; -using jtutil::pow4; //****************************************************************************** //****************************************************************************** @@ -99,8 +102,13 @@ bool geometry_is_finite(patch_system& ps, const struct error_info& error_info, bool initial_flag, bool print_msg_flag); -bool compute_Theta(patch_system& ps, fp add_to_expansion, - bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, +bool compute_Theta(patch_system& ps, const struct what_to_compute& comput_info, + bool Jacobian_flag, + jtutil::norm<fp>* Theta_norms_ptr, + jtutil::norm<fp>* expansion_Theta_norms_ptr, + jtutil::norm<fp>* inner_expansion_Theta_norms_ptr, + jtutil::norm<fp>* product_expansion_Theta_norms_ptr, + jtutil::norm<fp>* mean_curvature_Theta_norms_ptr, const struct error_info& error_info, bool initial_flag, bool print_msg_flag); } @@ -172,13 +180,18 @@ bool compute_Theta(patch_system& ps, fp add_to_expansion, // succeeded or failed, and if the latter, what caused the failure. // enum expansion_status - expansion(patch_system* ps_ptr, fp add_to_expansion, + expansion(patch_system* ps_ptr, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct error_info& error_info, bool initial_flag, bool Jacobian_flag /* = false */, bool print_msg_flag /* = false */, - jtutil::norm<fp>* Theta_norms_ptr /* = NULL */) + jtutil::norm<fp>* Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* expansion_Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* inner_expansion_Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* product_expansion_Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* mean_curvature_Theta_norms_ptr /* = NULL */) { const bool active_flag = (ps_ptr != NULL); if (print_msg_flag) @@ -230,6 +243,8 @@ if (gi.hardwire_Schwarzschild_EF_geometry) if (active_flag) then { + + if (gi.check_that_geometry_is_finite && !geometry_is_finite(*ps_ptr, error_info, initial_flag, @@ -237,15 +252,134 @@ if (active_flag) then return expansion_failure__geometry_nonfinite; // *** ERROR RETURN *** + // Ensure that there is a norm object + const bool want_norms = Theta_norms_ptr; + jtutil::norm<fp> norms; + if (compute_info.surface_selection != selection_definition) + then if (! Theta_norms_ptr) Theta_norms_ptr = &norms; + + // compute remaining gridfns --> $\Theta$ // and optionally also the Jacobian coefficients // by algebraic ops and angular finite differencing - if (!compute_Theta(*ps_ptr, add_to_expansion, + what_to_compute this_compute_info (compute_info); + this_compute_info.surface_selection = selection_definition; + if (!compute_Theta(*ps_ptr, this_compute_info, Jacobian_flag, Theta_norms_ptr, + expansion_Theta_norms_ptr, + inner_expansion_Theta_norms_ptr, + product_expansion_Theta_norms_ptr, + mean_curvature_Theta_norms_ptr, error_info, initial_flag, print_msg_flag)) then return expansion_failure__gij_not_positive_definite; // *** ERROR RETURN *** + + if (compute_info.surface_selection != selection_definition) { + // + // Apply correction to find a surface by its areal radius + // + // get mean expansion + fp mean_expansion; + fp areal_radius; + switch (compute_info.surface_selection) { + case selection_mean_coordinate_radius: { + const int np = ps_ptr->N_grid_points(); + fp sum_expansion = 0; + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + sum_expansion += p.gridfn(gfns::gfn__Theta, irho,isigma); + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_expansion = sum_expansion / np; + areal_radius = sum_radius / np; + break; + } + case selection_areal_radius: { + // get surface area + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + mean_expansion = Theta_norms_ptr->mean(); + areal_radius = sqrt(area / (4.0*PI)); + break; + } + case selection_expansion_mean_coordinate_radius: { + const int np = ps_ptr->N_grid_points(); + fp sum_expansion = 0; + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + sum_expansion += p.gridfn(gfns::gfn__Theta, irho,isigma); + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_expansion = sum_expansion / np; + areal_radius = mean_expansion * sum_radius / np; + break; + } + case selection_expansion_areal_radius: { + // get surface area + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + mean_expansion = Theta_norms_ptr->mean(); + areal_radius = mean_expansion * sqrt(area / (4.0*PI)); + break; + } + default: + assert (0); + } // switch areal_radius_definition + + if (! ps_ptr->N_additional_points()) { + // calculate correction + const fp correction + = (- mean_expansion + + areal_radius - compute_info.desired_value); + // apply correction + ps_ptr->add_to_gridfn(-correction, gfns::gfn__Theta); + } else { + const int np = ps_ptr->N_grid_points(); + const int gnp = ps_ptr->ghosted_N_grid_points(); + // apply correction + const fp correction + = ps_ptr->ghosted_gridfn_data(gfns::gfn__h)[gnp]; + ps_ptr->add_to_gridfn(-correction, gfns::gfn__Theta); + ps_ptr->gridfn_data(gfns::gfn__Theta)[np] + = (mean_expansion - correction + - areal_radius + compute_info.desired_value); + } + if (want_norms) { + // recalculate norms + ps_ptr->gridfn_norms (gfns::gfn__Theta, *Theta_norms_ptr); + } // if want_norms + } else { + // + // do not apply correction + // + if (ps_ptr->N_additional_points()) { + const int np = ps_ptr->N_grid_points(); + const int gnp = ps_ptr->ghosted_N_grid_points(); + ps_ptr->gridfn_data(gfns::gfn__Theta)[np] + = ps_ptr->ghosted_gridfn_data(gfns::gfn__h)[gnp]; + } + } // if use-areal-radius } return expansion_success; // *** NORMAL RETURN *** @@ -256,25 +390,19 @@ return expansion_success; // *** NORMAL RETURN *** //****************************************************************************** // -// This function sets up the xyz-position gridfns: -// * global_[xyz] (used by interplate_geometry() ) -// * global_{xx,xy,xz,yy,yz,zz} (used by higher-level code to compute -// quadrupole moments of horizons) -// -// Bugs: -// * We initialize the global_{xx,xy,xz,yy,yz,zz} gridfns every time -// this function is called, i.e. at each horizon-finding Newton -// iteration, even though they're only needed at the end of the -// horizon-finding process. In practice the extra cost is small, -// though, so it's probably not worth fixing this... +// This function sets up the global xyz positions of the grid points +// in the gridfns global_[xyz]. These will be used by interplate_geometry(). // namespace { -void setup_xyz_posns(patch_system& ps, bool print_msg_flag) +void setup_xyz_posns(patch_system& ps, const bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " xyz positions and derivative coefficients"); +#ifdef GEOMETRY_INTERP_DEBUG2 + printf ("AH exp\n"); +#endif for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -285,12 +413,15 @@ if (print_msg_flag) isigma <= p.max_isigma() ; ++isigma) { - const fp r = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + const fp r = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); const fp rho = p.rho_of_irho(irho); const fp sigma = p.sigma_of_isigma(isigma); - fp local_x, local_y, local_z; p.xyz_of_r_rho_sigma(r,rho,sigma, local_x,local_y,local_z); +#ifdef GEOMETRY_INTERP_DEBUG2 + printf (" pn=%d irho=%d isigma=%d x=%g y=%g z=%g\n", + pn, irho, isigma, local_x, local_y, local_z); +#endif const fp global_x = ps.origin_x() + local_x; const fp global_y = ps.origin_y() + local_y; @@ -353,12 +484,15 @@ if (print_msg_flag) // global_[xyz] # xyz positions of grid points // // Inputs (Cactus 3-D gridfns): +// ahmask # excision mask // gxx,gxy,gxz,gyy,gyz,gzz # 3-metric $g_{ij}$ // # (may be either physical or conformal) // kxx,kxy,kxz,kyy,kyz,kzz # extrinsic curvature $K_{ij}$ // psi # optional conformal factor $\psi$ // // Outputs (angular gridfns, all on the nominal grid): +// mask # excision mask +// partial_d_mask_[123] # derivatives of the mask // g_dd_{11,12,13,22,23,33} # $\stored{g}_{ij}$ // K_dd_{11,12,13,22,23,33} # $K_{ij}$ // partial_d_g_dd_[123]{11,12,13,22,23,33} # $\partial_k \stored{g}_{ij}$ @@ -428,6 +562,7 @@ const void* const interp_coords[N_GRID_DIMS] const CCTK_INT input_array_variable_indices[] = { + cgi.mask_varindex, cgi.g_dd_11_varindex, cgi.g_dd_12_varindex, cgi.g_dd_13_varindex, cgi.g_dd_22_varindex, cgi.g_dd_23_varindex, cgi.g_dd_33_varindex, @@ -436,19 +571,6 @@ const CCTK_INT input_array_variable_indices[] cgi.K_dd_33_varindex, cgi.psi_varindex, }; -const CCTK_INT input_array_type_codes[] - = { - // 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, - // psi - CCTK_VARIABLE_REAL, - }; const int N_input_arrays_for_psi = 1; const int N_input_arrays_dim = sizeof(input_array_variable_indices) / sizeof(input_array_variable_indices[0]); @@ -463,6 +585,8 @@ const int N_input_arrays_use const CCTK_INT output_array_type_codes[] = { + // mask $\partial_x$ mask $\partial_y$ mask $\partial_z$ mask + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, // $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, @@ -480,18 +604,20 @@ const CCTK_INT output_array_type_codes[] const CCTK_INT operand_indices[] = { - 0, 0, 0, 0, // g_dd_11, partial_[xyz] g_dd_11 - 1, 1, 1, 1, // g_dd_12, partial_[xyz] g_dd_12 - 2, 2, 2, 2, // g_dd_13, partial_[xyz] g_dd_13 - 3, 3, 3, 3, // g_dd_22, partial_[xyz] g_dd_22 - 4, 4, 4, 4, // g_dd_23, partial_[xyz] g_dd_23 - 5, 5, 5, 5, // g_dd_33, partial_[xyz] g_dd_33 - 6, 7, 8, 9, 10, 11, // K_dd_{11,12,13,22,23,33} - 12, 12, 12, 12, // psi, partial_[xyz] psi + 0, 0, 0, 0, // mask, partial_[xyz] mask + 1, 1, 1, 1, // g_dd_11, partial_[xyz] g_dd_11 + 2, 2, 2, 2, // g_dd_12, partial_[xyz] g_dd_12 + 3, 3, 3, 3, // g_dd_13, partial_[xyz] g_dd_13 + 4, 4, 4, 4, // g_dd_22, partial_[xyz] g_dd_22 + 5, 5, 5, 5, // g_dd_23, partial_[xyz] g_dd_23 + 6, 6, 6, 6, // g_dd_33, partial_[xyz] g_dd_33 + 7, 8, 9, 10, 11, 12, // K_dd_{11,12,13,22,23,33} + 13, 13, 13, 13, // psi, partial_[xyz] psi }; #define DERIV(x) x const CCTK_INT operation_codes[] = { + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // mask, partial_[xyz] mask DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_11, partial_[xyz] g_dd_11 DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_12, partial_[xyz] g_dd_12 DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_13, partial_[xyz] g_dd_13 @@ -505,6 +631,10 @@ const CCTK_INT operation_codes[] void* const output_arrays[] = { + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__mask)), + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_mask_1)), + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_mask_2)), + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_mask_3)), CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_11)), CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_111)), CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_211)), @@ -625,6 +755,21 @@ if (print_msg_flag) " calling geometry interpolator (%s%d points)", (active_flag ? "" : "dummy: "), N_interp_points); +#ifdef GEOMETRY_INTERP_DEBUG2 + { + printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() coordinates are:\n", + int(CCTK_MyProc(cgi.GH))); + for (int pt = 0 ; pt < N_interp_points ; ++ pt) + { + printf(" pt=%d [x,y,z]=[%g,%g,%g]\n", + pt, + double(((const CCTK_REAL*)interp_coords[0])[pt]), + double(((const CCTK_REAL*)interp_coords[1])[pt]), + double(((const CCTK_REAL*)interp_coords[2])[pt])); + } + } +#endif /* GEOMETRY_INTERP_DEBUG2 */ + #ifdef GEOMETRY_INTERP_DEBUG printf("AHFinderDirect:: proc %d: initializing interpolator outputs to 999.999\n", int(CCTK_MyProc(cgi.GH))); @@ -676,8 +821,11 @@ fflush(stdout); { for (int pt = 0 ; pt < N_interp_points ; pt = 2*pt + (pt == 0)) { - printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() results for pt=%d:\n", - int(CCTK_MyProc(cgi.GH)), pt); + printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() results for pt=%d at [x,y,z]=[%g,%g,%g]:\n", + int(CCTK_MyProc(cgi.GH)), pt, + double(((const CCTK_REAL*)interp_coords[0])[pt]), + double(((const CCTK_REAL*)interp_coords[1])[pt]), + double(((const CCTK_REAL*)interp_coords[2])[pt])); for (int out = 0 ; out < N_output_arrays_use ; ++out) { const CCTK_REAL* const out_ptr @@ -743,8 +891,8 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE) CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, "interpolate_geometry():\n" -" one or more points on the trial horizon surface point" -" is/are outside the grid (or too close to the grid boundary)" +" one or more points on the trial horizon surface point\n" +" is/are outside the grid (or too close to the grid boundary)\n" " (in a single-processor run, this may also mean that\n" " driver::ghost_size is too small for this geometry interpolator)\n"); } @@ -763,6 +911,51 @@ else if (status < 0) "***** interpolate_geometry(): error return %d from interpolator!\n", status); /*NOTREACHED*/ +// +// ***** check the interpolated excision mask ***** +// +{ +bool did_use_excised_gridpoint = false; +if (active_flag) + then { + for (int pn = 0 ; pn < ps_ptr->N_patches() ; ++pn) + { + patch& p = ps_ptr->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 m = p.gridfn(gfns::gfn__mask, irho,isigma); + const fp m1 = p.gridfn(gfns::gfn__partial_d_mask_1, irho,isigma); + const fp m2 = p.gridfn(gfns::gfn__partial_d_mask_2, irho,isigma); + const fp m3 = p.gridfn(gfns::gfn__partial_d_mask_3, irho,isigma); + if (fabs(m) > 1.0e-12 + || fabs(m1) > 1.0e-12 || fabs(m2) > 1.0e-12 || fabs(m3) > 1.0e-12) + then did_use_excised_gridpoint = true; + } + } + } +if (did_use_excised_gridpoint) + then { + if (print_msg_flag) + then { + // see if we can get further info + const int warn_level + = initial_flag + ? error_info.warn_level__point_outside__initial + : error_info.warn_level__point_outside__subsequent; + + + CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, +"interpolate_geometry():\n" +" one or more points on the trial horizon surface point\n" +" is/are in an excised region (or too close to the excision boundary)\n"); + } + + return expansion_failure__surface_in_excised_region; // *** ERROR RETURN *** + } +} + return expansion_success; // *** NORMAL RETURN *** } } @@ -1166,14 +1359,59 @@ return true; // *** no check possible *** // g_ij isn't positive definite). // namespace { -bool compute_Theta(patch_system& ps, fp add_to_expansion, - bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, +bool compute_Theta(patch_system& ps, const struct what_to_compute& compute_info, + bool Jacobian_flag, + jtutil::norm<fp>* Theta_norms_ptr, + jtutil::norm<fp>* expansion_Theta_norms_ptr, + jtutil::norm<fp>* inner_expansion_Theta_norms_ptr, + jtutil::norm<fp>* product_expansion_Theta_norms_ptr, + jtutil::norm<fp>* mean_curvature_Theta_norms_ptr, const struct error_info& error_info, bool initial_flag, bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " computing Theta(h)"); +#if 0 + fp mean_radius, areal_radius; + switch (compute_info.surface_modification) { + case modification_none: + case modification_radius: + case modification_radius2: + // do nothing + break; + case modification_mean_radius: { + // get average coordinate radius + const int np = ps.N_grid_points(); + fp sum_radius = 0; + 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) { + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_radius = sum_radius / np; + break; + } + case modification_areal_radius: { + // get surface area + const fp area = ps.integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + areal_radius = sqrt(area / (4.0*PI)); + break; + } + default: + assert (0); + } +#endif + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -1284,24 +1522,169 @@ CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, return false; // *** ERROR RETURN *** } + assert (compute_info.surface_selection == selection_definition); + // compute H via equation (14) of my 1996 horizon finding paper 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 - + add_to_expansion; + const fp Theta_X = + Theta_A/(Theta_D*sqrt_Theta_D) + + Theta_B/sqrt_Theta_D; + const fp Theta_Y = + Theta_C/Theta_D + - K; + +#define mean_curvature p.gridfn(gfns::gfn__mean_curvature, irho,isigma) + mean_curvature = Theta_X; +#undef mean_curvature + + switch (compute_info.surface_definition) { + case definition_expansion: + Theta = + Theta_X + Theta_Y; + break; + case definition_inner_expansion: + Theta = - Theta_X + Theta_Y; + break; + case definition_mean_curvature: + Theta = + Theta_X; + break; + case definition_expansion_product: + Theta = (+ Theta_X + Theta_Y) * (- Theta_X + Theta_Y); + break; + default: + assert (0); + } + + switch (compute_info.surface_modification) { + case modification_none: + // do nothing + break; + case modification_radius: + // multiply by radius + Theta *= r; + break; + case modification_radius2: + // multiply by radius^2 + Theta *= pow2(r); + break; +#if 0 + case modification_mean_radius: + // multiply by average coordinate radius + Theta *= mean_radius; + break; + case modification_areal_radius: + // multiply by areal radius + Theta *= areal_radius; + break; +#endif + default: + assert (0); + } + + Theta -= compute_info.desired_value; + // update running norms of Theta(h) function if (Theta_norms_ptr != NULL) then Theta_norms_ptr->data(Theta); + if (expansion_Theta_norms_ptr != NULL) + then expansion_Theta_norms_ptr->data(+ Theta_X + Theta_Y); + + if (inner_expansion_Theta_norms_ptr != NULL) + then inner_expansion_Theta_norms_ptr->data(- Theta_X + Theta_Y); + + if (product_expansion_Theta_norms_ptr != NULL) + then product_expansion_Theta_norms_ptr->data((+ Theta_X + Theta_Y) * (- Theta_X + Theta_Y)); + + if (mean_curvature_Theta_norms_ptr != NULL) + then mean_curvature_Theta_norms_ptr->data(+ Theta_X); + + fp partial_Theta_X_wrt_partial_d_h_1; + fp partial_Theta_X_wrt_partial_d_h_2; + fp partial_Theta_X_wrt_partial_dd_h_11; + fp partial_Theta_X_wrt_partial_dd_h_12; + fp partial_Theta_X_wrt_partial_dd_h_22; + fp partial_Theta_Y_wrt_partial_d_h_1; + fp partial_Theta_Y_wrt_partial_d_h_2; + fp partial_Theta_Y_wrt_partial_dd_h_11; + fp partial_Theta_Y_wrt_partial_dd_h_12; + fp partial_Theta_Y_wrt_partial_dd_h_22; + if (Jacobian_flag) then { // partial_Theta_wrt_partial_d_h, // partial_Theta_wrt_partial_dd_h #include "../gr.cg/expansion_Jacobian.c" } + + if (Jacobian_flag) { + switch (compute_info.surface_definition) { + + case definition_expansion: + partial_Theta_wrt_partial_d_h_1 + = (+ partial_Theta_X_wrt_partial_d_h_1 + + partial_Theta_Y_wrt_partial_d_h_1); + partial_Theta_wrt_partial_d_h_2 + = (+ partial_Theta_X_wrt_partial_d_h_2 + + partial_Theta_Y_wrt_partial_d_h_2); + partial_Theta_wrt_partial_dd_h_11 + = (+ partial_Theta_X_wrt_partial_dd_h_11 + + partial_Theta_Y_wrt_partial_dd_h_11); + partial_Theta_wrt_partial_dd_h_12 + = (+ partial_Theta_X_wrt_partial_dd_h_12 + + partial_Theta_Y_wrt_partial_dd_h_12); + partial_Theta_wrt_partial_dd_h_22 + = (+ partial_Theta_X_wrt_partial_dd_h_22 + + partial_Theta_Y_wrt_partial_dd_h_22); + break; + + case definition_inner_expansion: + partial_Theta_wrt_partial_d_h_1 + = (- partial_Theta_X_wrt_partial_d_h_1 + + partial_Theta_Y_wrt_partial_d_h_1); + partial_Theta_wrt_partial_d_h_2 + = (- partial_Theta_X_wrt_partial_d_h_2 + + partial_Theta_Y_wrt_partial_d_h_2); + partial_Theta_wrt_partial_dd_h_11 + = (- partial_Theta_X_wrt_partial_dd_h_11 + + partial_Theta_Y_wrt_partial_dd_h_11); + partial_Theta_wrt_partial_dd_h_12 + = (- partial_Theta_X_wrt_partial_dd_h_12 + + partial_Theta_Y_wrt_partial_dd_h_12); + partial_Theta_wrt_partial_dd_h_22 + = (- partial_Theta_X_wrt_partial_dd_h_22 + + partial_Theta_Y_wrt_partial_dd_h_22); + break; + + case definition_mean_curvature: + partial_Theta_wrt_partial_d_h_1 + = + partial_Theta_X_wrt_partial_d_h_1; + partial_Theta_wrt_partial_d_h_2 + = + partial_Theta_X_wrt_partial_d_h_2; + partial_Theta_wrt_partial_dd_h_11 + = + partial_Theta_X_wrt_partial_dd_h_11; + partial_Theta_wrt_partial_dd_h_12 + = + partial_Theta_X_wrt_partial_dd_h_12; + partial_Theta_wrt_partial_dd_h_22 + = + partial_Theta_X_wrt_partial_dd_h_22; + break; + + case definition_expansion_product: { +#define f(x,y,dx,dy) (- x*x + y*y) +#define df(x,y,dx,dy) (- 2*x*dx + 2*y*dy) + partial_Theta_wrt_partial_d_h_1 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_d_h_1 , partial_Theta_Y_wrt_partial_d_h_1 ); + partial_Theta_wrt_partial_d_h_2 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_d_h_2 , partial_Theta_Y_wrt_partial_d_h_2 ); + partial_Theta_wrt_partial_dd_h_11 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_dd_h_11, partial_Theta_Y_wrt_partial_dd_h_11); + partial_Theta_wrt_partial_dd_h_12 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_dd_h_12, partial_Theta_Y_wrt_partial_dd_h_12); + partial_Theta_wrt_partial_dd_h_22 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_dd_h_22, partial_Theta_Y_wrt_partial_dd_h_22); +#undef f +#undef df + break; + } + + default: + assert (0); + } + } + } } } diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc index 372d114..04e00b5 100644 --- a/src/gr/expansion_Jacobian.cc +++ b/src/gr/expansion_Jacobian.cc @@ -17,6 +17,7 @@ #include "util_Table.h" #include "cctk.h" +#include "cctk_Arguments.h" #include "config.h" #include "stdc.h" @@ -24,6 +25,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "../patch/coords.hh" #include "../patch/grid.hh" @@ -42,7 +44,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** @@ -53,7 +54,8 @@ using jtutil::error_exit; namespace { enum expansion_status expansion_Jacobian_NP - (patch_system& ps, Jacobian& Jac, fp add_to_expansion, + (patch_system& ps, Jacobian& Jac, + const struct what_to_compute& comput_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -73,7 +75,8 @@ void add_ghost_zone_Jacobian(const patch_system& ps, int xm_irho, int xm_isigma); enum expansion_status expansion_Jacobian_dr_FD - (patch_system* ps_ptr, Jacobian* Jac_ptr, fp add_to_expansion, + (patch_system* ps_ptr, Jacobian* Jac_ptr, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -109,7 +112,7 @@ enum expansion_status // enum expansion_status expansion_Jacobian(patch_system* ps_ptr, Jacobian* Jac_ptr, - fp add_to_expansion, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -125,7 +128,7 @@ case Jacobian__numerical_perturbation: if (active_flag) then { status = expansion_Jacobian_NP(*ps_ptr, *Jac_ptr, - add_to_expansion, + compute_info, cgi, gi, Jacobian_info, error_info, initial_flag, print_msg_flag); @@ -153,7 +156,8 @@ case Jacobian__symbolic_diff_with_FD_dr: // this function looks at ps_ptr and Jac_ptr (non-NULL vs NULL) // to choose a normal vs dummy computation { - status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr, add_to_expansion, + status = expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr, + compute_info, cgi, gi, Jacobian_info, error_info, initial_flag, print_msg_flag); @@ -171,6 +175,327 @@ default: int(Jacobian_info.Jacobian_compute_method)); /*NOTREACHED*/ } + if (active_flag) { + + switch (compute_info.surface_modification) { + + case modification_none: + // do nothing + break; + + case modification_radius: { + // multiply with the coordinate radius + // H_{(r)i} = H_i h_i + // J_{(r)ij} = J_{ij} h_i + H_i \delta_{ij} + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + const fp radius = p.ghosted_gridfn(gfns::gfn__h, irho, isigma); + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * radius); + } + } + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / radius; + Jac_ptr->sum_into_element (i, i, Theta); + } + } + } + break; + } + + case modification_radius2: { + // multiply with the square of the coordinate radius + // H_{(r2)i} = H_i h_i^2 + // J_{(r2)ij} = J_{ij} h_i^2 + 2 H_i h_i \delta_{ij} + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + const fp radius = p.ghosted_gridfn(gfns::gfn__h, irho, isigma); + const fp radius2 = radius * radius; + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * radius2); + } + } + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / radius2; + Jac_ptr->sum_into_element (i, i, 2 * Theta * radius); + } + } + } + break; + } + +#if 0 + case modification_mean_radius: { + // multiply with the average coordinate radius + // H_{(\bar r)i} = H_i \bar r + // J_{(\bar r)ij} = J_{ij} \bar r + H_i / N + // calculate average coordinate radius + const int np = ps_ptr->N_grid_points(); + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_radius = sum_radius / np; + // correct Jacobian + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * mean_radius); + } + } +#error "unfinished" + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / areal_radius; + const fp dRdh = 0.5 * areal_radius; + Jac_ptr->sum_into_element (i, i, Theta * dRdh); + } + } + } + break; + } + + case modification_areal_radius: { + // multiply with the areal radius + // H_{(R)i} = H_i R + // J_{(R)ij} = J_{ij} R + H_i dR/dh_j + // get surface area + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + const fp areal_radius = sqrt(area / (4.0*PI)); + // correct Jacobian + const int np = ps_ptr->N_grid_points(); + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); + for (int j=0; j<np; ++j) { + if (Jac_ptr->is_explicitly_stored (i, j)) { + const fp val = Jac_ptr->element (i, j); + Jac_ptr->set_element (i, j, val * areal_radius); + } + } + const fp Theta = (p.gridfn(gfns::gfn__Theta, irho, isigma) + + compute_info.desired_value) / areal_radius; + const fp dRdh = 0.5 * areal_radius; + Jac_ptr->sum_into_element (i, i, Theta * dRdh); + } + } + } + break; + } +#endif + + default: + assert (0); + } // switch surface_modification + + if (ps_ptr->N_additional_points()) { + switch (compute_info.surface_selection) { + + case selection_definition: { + // we want nothing special + const int np = ps_ptr->N_grid_points(); + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, 0.0); + } + for (int j=0; j<np; ++j) { + Jac_ptr->set_element (np, j, 0.0); + } + Jac_ptr->set_element (np, np, 1.0); + break; + } + + case selection_mean_coordinate_radius: { + // Jac_ptr->set_element (II, JJ, x) == dTheta(II)/dh(JJ) + // \frac{\partial R}{\partial h_j} = 1 / N + const int np = ps_ptr->N_grid_points(); + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, -1.0); + } + for (int j=0; j<np; ++j) { + fp val = 0; + for (int k=0; k<np; ++k) { + val += Jac_ptr->element (k, j) / np; + } + val -= 1.0 / np; + Jac_ptr->set_element (np, j, val); + } + Jac_ptr->set_element (np, np, -1.0); + break; + } + + case selection_areal_radius: { + // \frac{\partial R_a}{\partial h_j} + // = \sqrt{1 / 16 \pi A} \sum_k \sqrt{q_k} dS_k + // The "trapezoid" method is faster +// const enum patch::integration_method method +// = patch::integration_method__automatic_choice; + const enum patch::integration_method method + = patch::integration_method__trapezoid; + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method); + const int np = ps_ptr->N_grid_points(); + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, -1.0); + } + for (int j=0; j<np; ++j) { + fp val = 0; + for (int k=0; k<np; ++k) { + val += Jac_ptr->element (k, j) / np; + } + Jac_ptr->set_element (np, j, val); + } + for (int jpn = 0; jpn < ps_ptr->N_patches(); ++jpn) { + patch& jp = ps_ptr->ith_patch(jpn); + for (int jrho = jp.min_irho(); jrho <= jp.max_irho(); ++jrho) { + for (int jsigma = jp.min_isigma(); jsigma <= jp.max_isigma(); ++jsigma) { + const int j = ps_ptr->gpn_of_patch_irho_isigma(jp, jrho,jsigma); + // const fp radius = jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma); + const fp epsilon = Jacobian_info.perturbation_amplitude; + fp val1, val2; +#if 0 + // Re-calculate all points + // (this is slow, but it works) + val1 = area; + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) += epsilon; + val2 = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method); + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) -= epsilon; +#else + // Re-calculate all points with non-zero Jacobian entries + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) -= epsilon/2; + val1 = 0; + for (int ipn = 0; ipn < ps_ptr->N_patches(); ++ipn) { + patch& ip = ps_ptr->ith_patch(ipn); + for (int irho = ip.min_irho(); irho <= ip.max_irho(); ++irho) { + for (int isigma = ip.min_isigma(); isigma <= ip.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(ip, irho,isigma); + if (Jac_ptr->is_explicitly_stored (i, j)) { + val1 += ps_ptr->integrate_gridpoint + (gfns::gfn__one, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method, + ipn, irho, isigma); + } + } + } + } + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) += epsilon; + val2 = 0; + for (int ipn = 0; ipn < ps_ptr->N_patches(); ++ipn) { + patch& ip = ps_ptr->ith_patch(ipn); + for (int irho = ip.min_irho(); irho <= ip.max_irho(); ++irho) { + for (int isigma = ip.min_isigma(); isigma <= ip.max_isigma(); ++isigma) { + const int i = ps_ptr->gpn_of_patch_irho_isigma(ip, irho,isigma); + if (Jac_ptr->is_explicitly_stored (i, j)) { + val2 += ps_ptr->integrate_gridpoint + (gfns::gfn__one, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + method, + ipn, irho, isigma); + } + } + } + } + jp.ghosted_gridfn(gfns::gfn__h, jrho, jsigma) -= epsilon/2; +#endif + const fp val = 1 / sqrt(16*PI*area) * ps_ptr->integrate_correction(true, true, true) * (val2 - val1) / epsilon; + Jac_ptr->sum_into_element (np, j, -val); + } + } + } + Jac_ptr->set_element (np, np, -1.0); + break; + } + + case selection_expansion_mean_coordinate_radius: { + // Jac_ptr->set_element (II, JJ, x) == dTheta(II)/dh(JJ) + const int np = ps_ptr->N_grid_points(); + fp sum_expansion = 0; + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->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) { + sum_expansion += p.gridfn(gfns::gfn__Theta, irho,isigma); + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + for (int i=0; i<np; ++i) { + Jac_ptr->set_element (i, np, -1.0); + } + for (int j=0; j<np; ++j) { + fp val = 0; + for (int k=0; k<np; ++k) { + val += (Jac_ptr->element (k, j) / np) * (1.0 - sum_radius / np); + } + val -= (sum_expansion / np) / np; + Jac_ptr->set_element (np, j, val); + } + Jac_ptr->set_element (np, np, -1.0); + break; + } + + case selection_expansion_areal_radius: { + CCTK_WARN (0, "selection_expansion_areal_radius not implemented"); + break; + } + + default: + assert (0); + } // switch surface_selection + } else { + assert (compute_info.surface_selection == selection_definition); + } + + } // if active + return expansion_success; // *** NORMAL RETURN *** } @@ -214,7 +539,8 @@ return expansion_success; // *** NORMAL RETURN *** namespace { enum expansion_status expansion_Jacobian_NP - (patch_system& ps, Jacobian& Jac, fp add_to_expansion, + (patch_system& ps, Jacobian& Jac, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -227,6 +553,7 @@ if (print_msg_flag) const fp epsilon = Jacobian_info.perturbation_amplitude; ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); +ps.gridfn_copy(gfns::gfn__mean_curvature, gfns::gfn__save_mean_curvature); for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) { @@ -247,7 +574,8 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); 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; const - enum expansion_status status = expansion(&ps, add_to_expansion, + enum expansion_status status = expansion(&ps, + compute_info, cgi, gi, error_info, initial_flag); if (status != expansion_success) @@ -275,12 +603,19 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); } } + if (ps.N_additional_points()) + { + const int np = ps.N_grid_points(); + Jac.set_element(np,JJ, 0.0); // insert dummy value + } + yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y; } } } ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta); +ps.gridfn_copy(gfns::gfn__save_mean_curvature, gfns::gfn__mean_curvature); return expansion_success; // *** NORMAL RETURN *** } } @@ -440,6 +775,12 @@ ps.compute_synchronize_Jacobian(); } } + if (ps.N_additional_points()) + { + const int np = ps.N_grid_points(); + Jac.set_element(II,np, 0.0); // insert dummy value + } + } } } @@ -519,11 +860,12 @@ patch& yp = ye.my_patch(); // It's illegal for one but not both of ps_ptr and Jac_ptr to be NULL. // // The basic algorithm is that -// Jac += diag[ (Theta(h+epsilon) - Theta(h)) / epsilon ] +// Jac += diag[ (Theta(h+epsilon/2) - Theta(h-epsilon/2)) / epsilon ] // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface // Theta # Theta(h) assumed to already be computed +// # (saved and restored, but not used) // // Outputs: // Jac += d/dr terms @@ -535,7 +877,8 @@ patch& yp = ye.my_patch(); namespace { enum expansion_status expansion_Jacobian_dr_FD - (patch_system* ps_ptr, Jacobian* Jac_ptr, fp add_to_expansion, + (patch_system* ps_ptr, Jacobian* Jac_ptr, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, @@ -550,18 +893,52 @@ if (print_msg_flag) const fp epsilon = Jacobian_info.perturbation_amplitude; -// compute Theta(h+epsilon) +what_to_compute this_compute_info (compute_info); +this_compute_info.surface_modification = modification_none; +this_compute_info.surface_selection = selection_definition; +this_compute_info.desired_value = 0.0; + +fp additional_save_Theta; + +// compute Theta(h-epsilon/2) if (active_flag) then { ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta); - ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h); + ps_ptr->gridfn_copy(gfns::gfn__mean_curvature, gfns::gfn__save_mean_curvature); + if (ps_ptr->N_additional_points()) + then { + const int np = ps_ptr->N_grid_points(); + additional_save_Theta = ps_ptr->gridfn_data(gfns::gfn__Theta)[np]; + } + ps_ptr->add_to_ghosted_gridfn(-epsilon/2, gfns::gfn__h); } const - enum expansion_status status = expansion(ps_ptr, add_to_expansion, + enum expansion_status status = expansion(ps_ptr, + this_compute_info, cgi, gi, error_info, initial_flag); if (status != expansion_success) - then return status; // *** ERROR RETURN *** + then { + expansion(NULL, + this_compute_info, + cgi, gi, + error_info, false); + return status; // *** ERROR RETURN *** + } + +// compute Theta(h+epsilon/2) +if (active_flag) + then { + ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__old_Theta); + ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h); + } +const + enum expansion_status status2 = expansion(ps_ptr, + this_compute_info, + cgi, gi, + error_info, initial_flag); +if (status2 != expansion_success) + then return status2; // *** ERROR RETURN *** if (active_flag) then { @@ -575,7 +952,7 @@ if (active_flag) ++isigma) { const int II = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma); - const fp old_Theta = p.gridfn(gfns::gfn__save_Theta, + const fp old_Theta = p.gridfn(gfns::gfn__old_Theta, irho,isigma); const fp new_Theta = p.gridfn(gfns::gfn__Theta, irho,isigma); @@ -586,8 +963,14 @@ if (active_flag) } // restore h and Theta - ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h); + ps_ptr->add_to_ghosted_gridfn(-epsilon/2, gfns::gfn__h); ps_ptr->gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta); + ps_ptr->gridfn_copy(gfns::gfn__save_mean_curvature, gfns::gfn__mean_curvature); + if (ps_ptr->N_additional_points()) + then { + const int np = ps_ptr->N_grid_points(); + ps_ptr->gridfn_data(gfns::gfn__Theta)[np] = additional_save_Theta; + } } return expansion_success; // *** NORMAL RETURN *** diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh index ee80c90..c223d03 100644 --- a/src/gr/gfns.hh +++ b/src/gr/gfns.hh @@ -16,10 +16,11 @@ namespace gfns // ghosted gridfns enum { - ghosted_min_gfn = -1, // must set this by hand so + ghosted_min_gfn = -3, // must set this by hand so // ghosted_max_gfn is still < 0 gfn__h = ghosted_min_gfn, - ghosted_max_gfn = gfn__h + gfn__save_h, + ghosted_max_gfn = gfn__save_h }; // nominal gridfns @@ -38,7 +39,7 @@ enum { gfn__global_x = nominal_min_gfn, // no access macro gfn__global_y, // no access macro gfn__global_z, // no access macro - + gfn__global_xx, // no access macro gfn__global_xy, // no access macro gfn__global_xz, // no access macro @@ -46,6 +47,11 @@ enum { gfn__global_yz, // no access macro gfn__global_zz, // no access macro + gfn__mask, // no access macro + gfn__partial_d_mask_1, // no access macro + gfn__partial_d_mask_2, // no access macro + gfn__partial_d_mask_3, // no access macro + gfn__g_dd_11, gfn__g_dd_12, gfn__g_dd_13, @@ -82,6 +88,9 @@ enum { gfn__partial_d_psi_2, // no access macro gfn__partial_d_psi_3, // no access macro + gfn__mean_curvature, // no access macro + gfn__save_mean_curvature, // no access macro + gfn__Theta, gfn__partial_Theta_wrt_partial_d_h_1, gfn__partial_Theta_wrt_partial_d_h_2, @@ -90,6 +99,8 @@ enum { gfn__partial_Theta_wrt_partial_dd_h_22, gfn__Delta_h, gfn__save_Theta, + gfn__old_Theta, + gfn__zero, gfn__one, nominal_max_gfn = gfn__one // no comma }; diff --git a/src/gr/gr.hh b/src/gr/gr.hh index 00f106c..3564438 100644 --- a/src/gr/gr.hh +++ b/src/gr/gr.hh @@ -1,3 +1,4 @@ + // gr.hh -- header file for general relativity code // $Header$ @@ -14,6 +15,68 @@ namespace AHFinderDirect enum { N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 }; // +// this enum specifies what kind of surface we want +// +enum a_surface_definition + { + definition_error, // this value is illegal + definition_expansion, // apparent horizon + definition_inner_expansion, // expansion Theta_(l), ingoing null normal + definition_mean_curvature, // mean curvature + definition_expansion_product // product of Theta_(n) and Theta_(l) + }; + +// +// this enum specifies how the surface definition is modified +// +enum a_surface_modification + { + modification_error, // this value is illegal + modification_none, // no modification + modification_radius, // times coordinate radius + modification_radius2 // times coordinate radius^2 +#if 0 + modification_mean_radius, // times mean coordinate radius + modification_areal_radius // times areal radius +#endif + }; + +// +// this enum specifies how we select the surface +// +enum a_surface_selection + { + selection_error, // this value is illegal + selection_definition, // use the surface's definition + selection_mean_coordinate_radius, // mean coordinate radius (cheap) + selection_areal_radius, // areal radius + selection_expansion_mean_coordinate_radius, // expansion times mean coordinate radius + selection_expansion_areal_radius // expansion times areal radius + }; + +// +// this struct specifies what to calculate +// +struct what_to_compute + { + // how Theta is calculated + a_surface_definition surface_definition; + // how Theta is modified + a_surface_modification surface_modification; + // what is solved for + a_surface_selection surface_selection; + // the desired value (expansion, areal radius, etc.) + fp desired_value; + + what_to_compute () + : surface_definition (definition_error), + surface_modification (modification_error), + surface_selection (selection_error), + desired_value (0.0) + { } + }; + +// // this enum holds the (a) decoded Jacobian_compute_method parameter, // i.e. it specifies how we compute the (a) Jacobian matrix // @@ -76,6 +139,7 @@ struct cactus_grid_info bool use_Cactus_conformal_metric; // Cactus variable indices of geometry variables + int mask_varindex; int g_dd_11_varindex, g_dd_12_varindex, g_dd_13_varindex, g_dd_22_varindex, g_dd_23_varindex, g_dd_33_varindex; @@ -169,18 +233,23 @@ struct error_info // expansion.cc enum expansion_status - expansion(patch_system* ps_ptr, fp add_to_expansion, + expansion(patch_system* ps_ptr, + const struct what_to_compute& comput_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct error_info& error_info, bool initial_flag, bool Jacobian_flag = false, bool print_msg_flag = false, - jtutil::norm<fp>* H_norms_ptr = NULL); + jtutil::norm<fp>* H_norms_ptr = NULL, + jtutil::norm<fp>* expansion_H_norms_ptr = NULL, + jtutil::norm<fp>* inner_expansion_H_norms_ptr = NULL, + jtutil::norm<fp>* product_expansion_H_norms_ptr = NULL, + jtutil::norm<fp>* mean_curvature_H_norms_ptr = NULL); // expansion_Jacobian.cc enum expansion_status expansion_Jacobian(patch_system* ps_ptr, Jacobian* Jac_ptr, - fp add_to_expansion, + const struct what_to_compute& comput_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct Jacobian_info& Jacobian_info, diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc index 34da085..b72531c 100644 --- a/src/gr/gr_gfas.minc +++ b/src/gr/gr_gfas.minc @@ -40,5 +40,7 @@ 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 +partial_Theta_X_wrt_partial_d_h, partial_Theta_X_wrt_partial_d_h__fnd, +partial_Theta_Y_wrt_partial_d_h, partial_Theta_Y_wrt_partial_d_h__fnd, +partial_Theta_X_wrt_partial_dd_h, partial_Theta_X_wrt_partial_dd_h__fnd, +partial_Theta_Y_wrt_partial_dd_h, partial_Theta_Y_wrt_partial_dd_h__fnd # no comma diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple index 7284253..423642f 100644 --- a/src/gr/horizon.maple +++ b/src/gr/horizon.maple @@ -258,7 +258,7 @@ global @include "../maple/gfa.minc", @include "../gr/gr_gfas.minc"; local u,v, - temp; + temp1,temp2; printf("%a...\n", procname); @@ -284,14 +284,16 @@ assert_fnd_exists(Theta_D); := frontend('diff', [Theta_D__fnd, Diff(h,y_rs[u])]); # equation (A1a) in my 1996 apparent horizon finding paper - 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] + temp1 := + (3/2)*Theta_A/Theta_D^(5/2) + + (1/2)*Theta_B/Theta_D^(3/2); + partial_Theta_X_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; + - partial_Theta_D_wrt_partial_d_h__fnd[u] * temp1; + temp2 := + Theta_C/Theta_D^2; + partial_Theta_Y_wrt_partial_d_h__fnd[u] + := + partial_Theta_C_wrt_partial_d_h__fnd[u] / Theta_D + - partial_Theta_D_wrt_partial_d_h__fnd[u] * temp2; end do; # Jacobian coefficients of Theta_[AB] and Theta wrt Diff(h,y_rs[u],y_rs[v]) @@ -305,17 +307,23 @@ assert_fnd_exists(Theta_D); := frontend('diff', [Theta_B__fnd, Diff(h,y_rs[u],y_rs[v])]); # equation (A1b) in my 1996 apparent horizon finding paper - partial_Theta_wrt_partial_dd_h__fnd[u,v] + partial_Theta_X_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) + + partial_Theta_B_wrt_partial_dd_h__fnd[u,v] / Theta_D^(1/2); + partial_Theta_Y_wrt_partial_dd_h__fnd[u,v] + := 0; end do; end do; 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'], + then codegen2([partial_Theta_X_wrt_partial_d_h__fnd, + partial_Theta_Y_wrt_partial_d_h__fnd, + partial_Theta_X_wrt_partial_dd_h__fnd, + partial_Theta_Y_wrt_partial_dd_h__fnd], + ['partial_Theta_X_wrt_partial_d_h', + 'partial_Theta_Y_wrt_partial_d_h', + 'partial_Theta_X_wrt_partial_dd_h', + 'partial_Theta_Y_wrt_partial_dd_h'], "../gr.cg/expansion_Jacobian.c"); fi; diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn index 37e9db6..71bf658 100644 --- a/src/gr/make.code.defn +++ b/src/gr/make.code.defn @@ -10,16 +10,9 @@ SRCS = expansion.cc \ # Subdirectories containing source files SUBDIRS = -# disable automatic template instantiation on DEC Alphas cxx compiler +# disable automatic template instantiation on DEC Alphas ifeq ($(shell uname), OSF1) ifeq ($(CXX), cxx) CXXFLAGS += -nopt endif endif - -# disable automagic template instantiation on SGI Irix CC compiler -ifneq (,$(findstring IRIX,$(shell uname))) - ifeq ($(notdir $(CXX)), CC) - CXXFLAGS += -no_auto_include - endif -endif diff --git a/src/gr/maple.log b/src/gr/maple.log index 54c762f..1a91a6c 100644 --- a/src/gr/maple.log +++ b/src/gr/maple.log @@ -1,10 +1,11 @@ +RedHat v <9 or other Linux present, starting standard mode... |\^/| Maple 7 (IBM INTEL LINUX) ._|\| |/|_. Copyright (c) 2001 by Waterloo Maple Inc. \ MAPLE / All rights reserved. Maple is a registered trademark of <____ ____> Waterloo Maple Inc. | Type ? for help. # top-level Maple file to read/run all code in this directory -# $Header: /numrelcvs/AEIDevelopment/AHFinderDirect/src/gr/doit.maple,v 1.5 2002/09/13 14:12:18 jthorn Exp $ +# $Header: /numrelcvs/AEIThorns/AHFinderDirect/src/gr/doit.maple,v 1.5 2002/09/13 14:12:18 jthorn Exp $ > > read "../maple/setup.mm"; msum := proc(fn::algebraic) @@ -205,19 +206,22 @@ 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_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; option remember; var_list := [args[2 .. nargs]]; if type(operand, indexed) and op(0, operand) = 'X_ud' and @@ -597,9 +601,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_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); @@ -635,9 +641,13 @@ partial_Theta_wrt_partial_dd_h__fnd; [1 .. N_ang, 1 .. N_ang], symmetric); make_gfa('partial_Theta_B_wrt_partial_dd_h', {inert, fnd}, [1 .. N_ang, 1 .. N_ang], symmetric); - make_gfa('partial_Theta_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], + make_gfa('partial_Theta_X_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], + none); + make_gfa('partial_Theta_Y_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none); - make_gfa('partial_Theta_wrt_partial_dd_h', {inert, fnd}, + make_gfa('partial_Theta_X_wrt_partial_dd_h', {inert, fnd}, + [1 .. N_ang, 1 .. N_ang], symmetric); + make_gfa('partial_Theta_Y_wrt_partial_dd_h', {inert, fnd}, [1 .. N_ang, 1 .. N_ang], symmetric); NULL end proc @@ -660,9 +670,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; option remember; var_list := [args[2 .. nargs]]; if type(operand, indexed) and op(0, operand) = 'g_dd' and @@ -696,9 +708,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(g_dd); assert_fnd_exists(g_uu, fnd); @@ -727,9 +741,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(g_uu); assert_fnd_exists(K_dd); @@ -771,9 +787,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(g_dd); assert_fnd_exists(g_uu); @@ -809,9 +827,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(g_dd); assert_fnd_exists(g_uu); @@ -852,9 +872,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(h); assert_fnd_exists(X_ud); @@ -883,9 +905,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(h); assert_fnd_exists(X_ud); @@ -923,9 +947,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(g_uu); assert_fnd_exists(K_uu); @@ -958,7 +984,7 @@ partial_Theta_wrt_partial_dd_h__fnd; end proc expansion_Jacobian := proc(cg_flag::boolean) -local u, v, temp; +local u, v, temp1, temp2; 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, @@ -975,9 +1001,11 @@ 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; +partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_X_wrt_partial_d_h, +partial_Theta_X_wrt_partial_d_h__fnd, partial_Theta_Y_wrt_partial_d_h, +partial_Theta_Y_wrt_partial_d_h__fnd, partial_Theta_X_wrt_partial_dd_h, +partial_Theta_X_wrt_partial_dd_h__fnd, partial_Theta_Y_wrt_partial_dd_h, +partial_Theta_Y_wrt_partial_dd_h__fnd; printf("%a...\n", procname); assert_fnd_exists(g_uu); assert_fnd_exists(K_uu); @@ -996,13 +1024,15 @@ partial_Theta_wrt_partial_dd_h__fnd; 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] := + temp1 := 3/2*Theta_A/Theta_D^(5/2) + 1/2*Theta_B/Theta_D^(3/2); + partial_Theta_X_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 + - partial_Theta_D_wrt_partial_d_h__fnd[u]*temp1; + temp2 := Theta_C/Theta_D^2; + partial_Theta_Y_wrt_partial_d_h__fnd[u] := + partial_Theta_C_wrt_partial_d_h__fnd[u]/Theta_D + - partial_Theta_D_wrt_partial_d_h__fnd[u]*temp2 end do; for u to N_ang do for v from u to N_ang do partial_Theta_A_wrt_partial_dd_h__fnd[u, v] := @@ -1011,16 +1041,22 @@ partial_Theta_wrt_partial_dd_h__fnd; 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_X_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) + partial_Theta_B_wrt_partial_dd_h__fnd[u, v]/Theta_D^(1/2); + partial_Theta_Y_wrt_partial_dd_h__fnd[u, v] := 0 end do end do; - 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") + if cg_flag then codegen2([partial_Theta_X_wrt_partial_d_h__fnd, + partial_Theta_Y_wrt_partial_d_h__fnd, + partial_Theta_X_wrt_partial_dd_h__fnd, + partial_Theta_Y_wrt_partial_dd_h__fnd], [ + 'partial_Theta_X_wrt_partial_d_h', + 'partial_Theta_Y_wrt_partial_d_h', + 'partial_Theta_X_wrt_partial_dd_h', + 'partial_Theta_Y_wrt_partial_dd_h'], + "../gr.cg/expansion_Jacobian.c") end if; NULL end proc @@ -1040,7 +1076,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=1000088, alloc=917336, time=0.12 +bytes used=1000776, alloc=917336, time=0.13 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc --> `codegen2/unindex` @@ -1053,7 +1089,7 @@ codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c" convert --> equation list --> `codegen2/eqnlist` optimizing computation sequence -bytes used=2000352, alloc=1441528, time=0.17 +bytes used=2001072, alloc=1441528, time=0.20 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` @@ -1061,39 +1097,39 @@ bytes used=2000352, alloc=1441528, time=0.17 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc --> `codegen2/unindex` -bytes used=3000584, alloc=1638100, time=0.22 +bytes used=3001364, alloc=1638100, time=0.26 convert p/q --> RATIONAL(p/q) --> `codegen2/fix_rationals` writing C code > curvature(true); inverse_metric_gradient... -bytes used=4000860, alloc=1703624, time=0.28 +bytes used=4001680, alloc=1703624, time=0.34 codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c" --> `codegen2/input` convert --> equation list --> `codegen2/eqnlist` optimizing computation sequence -bytes used=5001020, alloc=1703624, time=0.38 +bytes used=5001836, alloc=1703624, time=0.42 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc -bytes used=6001324, alloc=1769148, time=0.44 +bytes used=6002212, alloc=1769148, time=0.49 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc --> `codegen2/unindex` -bytes used=7001528, alloc=1769148, time=0.51 +bytes used=7002736, alloc=1769148, time=0.55 convert p/q --> RATIONAL(p/q) --> `codegen2/fix_rationals` writing C code -bytes used=8001896, alloc=1769148, time=0.57 +bytes used=8003012, alloc=1769148, time=0.62 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=9002080, alloc=1769148, time=0.64 +bytes used=9003256, alloc=1769148, time=0.75 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` @@ -1108,29 +1144,29 @@ codegen/C/expression: Unknown function: RATIONAL will be left as is. > horizon(true); non_unit_normal... non_unit_normal_deriv... -bytes used=10002316, alloc=1769148, time=0.69 +bytes used=10005164, alloc=1769148, time=0.82 expansion... -bytes used=11002748, alloc=1834672, time=0.77 +bytes used=11005388, alloc=1834672, time=0.89 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=12003104, alloc=1834672, time=0.83 -bytes used=13003356, alloc=2031244, time=0.95 -bytes used=14003520, alloc=2031244, time=1.05 +bytes used=12005892, alloc=1834672, time=0.96 +bytes used=13006112, alloc=2031244, time=1.09 +bytes used=14006296, alloc=2031244, time=1.20 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc -bytes used=15003896, alloc=2096768, time=1.11 +bytes used=15006844, alloc=2096768, time=1.27 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc -bytes used=16004136, alloc=2096768, time=1.18 +bytes used=16007112, alloc=2096768, time=1.34 --> `codegen2/unindex` -bytes used=17004388, alloc=2096768, time=1.25 +bytes used=17007616, alloc=2096768, time=1.40 convert p/q --> RATIONAL(p/q) -bytes used=18004580, alloc=2096768, time=1.31 +bytes used=18007784, alloc=2096768, time=1.46 --> `codegen2/fix_rationals` writing C code codegen/C/expression: Unknown function: PARTIAL_RHO will be left as is. @@ -1141,57 +1177,58 @@ 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=19004732, alloc=2096768, time=1.37 -bytes used=20004972, alloc=2096768, time=1.47 +bytes used=19008004, alloc=2096768, time=1.51 +bytes used=20008500, alloc=2096768, time=1.62 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" +bytes used=21009180, alloc=2096768, time=1.70 +codegen2([partial_Theta_X_wrt_partial_d_h, partial_Theta_Y_wrt_partial_d_h, partial_Theta_X_wrt_partial_dd_h, partial_Theta_Y_wrt_partial_dd_h]) --> "../gr.cg/expansion_Jacobian.c" --> `codegen2/input` convert --> equation list -bytes used=22005440, alloc=2096768, time=1.63 +bytes used=22009512, alloc=2096768, time=1.77 --> `codegen2/eqnlist` optimizing computation sequence -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 +bytes used=23009696, alloc=2162292, time=1.84 +bytes used=24009880, alloc=2293340, time=1.91 +bytes used=25010716, alloc=2293340, time=1.97 +bytes used=26010940, alloc=2293340, time=2.07 +bytes used=27011284, alloc=2293340, time=2.21 +bytes used=28011604, alloc=2293340, time=2.37 +bytes used=29011784, alloc=2293340, time=2.47 +bytes used=30012024, alloc=2424388, time=2.60 +bytes used=31012336, alloc=2424388, time=2.69 +bytes used=32012604, alloc=2489912, time=2.80 --> `codegen2/optimize` find temporary variables +bytes used=33012892, alloc=2620960, time=2.87 --> `codegen2/temps` convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc -bytes used=33010028, alloc=2555436, time=2.77 -bytes used=34010432, alloc=2555436, time=2.84 -bytes used=35010740, alloc=2555436, time=2.91 +bytes used=34013116, alloc=2620960, time=2.94 +bytes used=35013316, alloc=2620960, time=3.01 +bytes used=36013600, alloc=2620960, time=3.08 --> `codegen2/fix_Diff` -bytes used=36021172, alloc=2555436, time=2.99 convert R_dd[2,3] --> R_dd_23 etc -bytes used=37021516, alloc=2555436, time=3.05 -bytes used=38021792, alloc=2555436, time=3.12 -bytes used=39022208, alloc=2555436, time=3.20 +bytes used=37014080, alloc=2620960, time=3.15 +bytes used=38014560, alloc=2620960, time=3.22 +bytes used=39014928, alloc=2620960, time=3.29 --> `codegen2/unindex` -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 +bytes used=40015272, alloc=2620960, time=3.36 +bytes used=41015476, alloc=2620960, time=3.43 +bytes used=42015660, alloc=2620960, time=3.49 +bytes used=43016244, alloc=2620960, time=3.56 +bytes used=44016600, alloc=2620960, time=3.63 convert p/q --> RATIONAL(p/q) -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 +bytes used=45016764, alloc=2620960, time=3.68 +bytes used=46016940, alloc=2620960, time=3.74 +bytes used=47017228, alloc=2620960, time=3.81 +bytes used=48017676, alloc=2620960, time=3.87 --> `codegen2/fix_rationals` writing C code -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 +bytes used=49017944, alloc=2620960, time=3.94 +bytes used=50018104, alloc=2620960, time=3.99 +bytes used=51018308, alloc=2620960, time=4.08 +bytes used=52018568, alloc=2620960, time=4.20 +bytes used=53018900, alloc=2620960, time=4.35 +bytes used=54019192, alloc=2620960, time=4.53 +bytes used=55019432, alloc=2620960, time=4.68 > quit -bytes used=54895808, alloc=2555436, time=4.56 +bytes used=55315640, alloc=2620960, time=4.73 diff --git a/src/gr/misc-gr.cc b/src/gr/misc-gr.cc index be9ccd0..eaaaf9a 100644 --- a/src/gr/misc-gr.cc +++ b/src/gr/misc-gr.cc @@ -11,6 +11,7 @@ #include "util_Table.h" #include "cctk.h" +#include "cctk_Arguments.h" #include "config.h" #include "stdc.h" diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple index 6af4aed..b2fed25 100644 --- a/src/gr/setup_gr_gfas.maple +++ b/src/gr/setup_gr_gfas.maple @@ -74,8 +74,11 @@ make_gfa('partial_Theta_B_wrt_partial_dd_h', {inert,fnd}, [1..N_ang, 1..N_ang], symmetric); # 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}, +make_gfa('partial_Theta_X_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +make_gfa('partial_Theta_Y_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +make_gfa('partial_Theta_X_wrt_partial_dd_h', {inert,fnd}, + [1..N_ang, 1..N_ang], symmetric); +make_gfa('partial_Theta_Y_wrt_partial_dd_h', {inert,fnd}, [1..N_ang, 1..N_ang], symmetric); NULL; |