From 720ce53a79ef5bdff53b3db18bc45c46101aa0fa Mon Sep 17 00:00:00 2001 From: jthorn Date: Fri, 26 Jul 2002 14:57:48 +0000 Subject: re-sync changes from laptop - parameter to control how Jacobian is computed - can hardwire geometry to Schwarzschild/EF git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@661 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/Jacobian.cc | 24 +++ src/elliptic/Jacobian.hh | 9 +- src/gr/AHFinderDirect.hh | 80 +++++++--- src/gr/Newton.cc | 19 ++- src/gr/README | 4 +- src/gr/Schwarzschild_EF.cc | 366 ++++++++++++++++++++++++++++++++++++++++++++ src/gr/cg.hh | 96 ++++++------ src/gr/driver.cc | 180 ++++++++++++++-------- src/gr/horizon_Jacobian.cc | 261 +++++++++++++++++-------------- src/gr/horizon_function.cc | 162 ++++++++++---------- src/gr/make.code.defn | 1 + src/jtutil/miscfp.cc | 5 +- src/make.configuration.defn | 2 +- src/patch/coords.cc | 12 ++ src/patch/coords.hh | 1 + 15 files changed, 866 insertions(+), 356 deletions(-) create mode 100644 src/gr/Schwarzschild_EF.cc diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 58009ec..603abac 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -10,6 +10,8 @@ // dense_Jacobian::zero_row // dense_Jacobian::solve_linear_system // +// new_Jacobian +// #include using std::fopen; @@ -180,3 +182,25 @@ if (info != 0) rhs_gfn, x_gfn, int(info)); /*NOTREACHED*/ } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function is an "object factory" for Jacobians: it constructs +// and returns a new-allocated Jacobian object of a specified type. +// +// FIXME: the patch system shouldn't really have to be non-const, but +// the Jacobian constructors all require this to allow the +// linear solvers to directly update gridfns +// +Jacobian& new_Jacobian(patch_system& ps, + const char Jacobian_type[]) +{ +if (STRING_EQUAL(Jacobian_type, "dense matrix")) + then return *new dense_Jacobian(ps); +else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + "unknown Jacobian_type=\"%s\"!", + Jacobian_type); /*NOTREACHED*/ +} diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh index 6d20e63..7bd37db 100644 --- a/src/elliptic/Jacobian.hh +++ b/src/elliptic/Jacobian.hh @@ -3,7 +3,8 @@ // // Jacobian -- abstract base class to describe a Jacobian matrix -// dense_Jacobian - Jacobian stored as a dense matrix +// dense_Jacobian -- Jacobian stored as a dense matrix +// new_Jacobian - factory method // // @@ -132,3 +133,9 @@ private: // pivot vector for LAPACK routines integer *pivot_; }; + +//****************************************************************************** + +// construct and return new-allocated Jacobian object of specified type +Jacobian& new_Jacobian(patch_system& ps, + const char Jacobian_type[]); diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index ac79f14..403ea39 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -3,20 +3,57 @@ //****************************************************************************** +// // number of spatial dimensions in the main Cactus grid // and in our trial-horizon-surface grid +// enum { N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 }; // -// This structure holds various information about the interpolator -// used to interpolate the "geometry" Cactus 3-D gridfns ($g_{ij}$ and -// $K_{ij}$) to the apparent horizon position. +// this enum holds the (a) decoded Jacobian_method parameter, +// i.e. it specifies how we compute the (a) Jacobian matrix // -struct geometry_interpolator_info +enum Jacobian_method { + numerical_perturbation, + symbolic_differentiation_with_FD_dr, + symbolic_differentiation // no comma + }; + + +// +// This structure holds information for computing the spacetime geometry. +// This is normally done by interpolating $g_{ij}$ and $K_{ij}$ from the +// usual Cactus grid, but can optionally instead by done by hard-wiring +// the Schwarzschild geometry in Eddington-Finkelstein coordinates. +// +struct geometry_info + { + // + // parameters for hard-wiring Schwarzschild/EF geometry + // + bool hardwire_Schwarzschild_EF; // should we hard-wire the + // Schwarzschild/EF geometry? + fp hardwire_Schwarzschild_EF__x_posn; // x posn of Schwarzschild BH + fp hardwire_Schwarzschild_EF__y_posn; // y posn of Schwarzschild BH + fp hardwire_Schwarzschild_EF__z_posn; // z posn of Schwarzschild BH + fp hardwire_Schwarzschild_EF__mass; // mass of Schwarzschild BH + fp hardwire_Schwarzschild_EF__epsilon; // threshold for sin^2 theta + // = (x^2+y^2)/r^2 below which + // we use z axis limits + fp Delta_xyz; // pseudo-grid spacing for finite differencing + // computation of $\partial_k g_{ij}$ + + // + // parameters for normal interpolation from Cactus grid + // int operator_handle; // Cactus handle to interpolation op int param_table_handle; // Cactus handle to parameter table // for the interpolator + + // this doesn't really belong in this structure (it doesn't + // have any logical connection to the geometry calculations), + // but it's convenient to put it here anyway... }; // @@ -56,17 +93,16 @@ struct cactus_grid_info }; // -// This struct holds the various parameters used in computing the Jacobian -// matrix. +// This struct holds parameters for computing the Jacobian matrix. // struct Jacobian_info { + enum Jacobian_method Jacobian_method; fp perturbation_amplitude; }; // -// This struct holds the various parameters used in solving the H(h) = 0 -// equations. +// This struct holds parameters for solving the H(h) = 0 equations. // struct solver_info { @@ -89,30 +125,28 @@ extern "C" // horizon_function.cc void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - bool Jacobian_flag, - jtutil::norm* H_norms_ptr, - bool msg_flag = true); + const struct geometry_info& gi, + bool Jacobian_flag = false, + bool msg_flag = false, + jtutil::norm* H_norms_ptr = NULL); // horizon_Jacobian.cc -Jacobian& create_Jacobian(patch_system& ps, - const char Jacobian_type[]); void horizon_Jacobian(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - const struct Jacobian_info& Jac_info, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, Jacobian& Jac); -void horizon_Jacobian_NP(patch_system& ps, - const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - const struct Jacobian_info& Jac_info, - Jacobian& Jac); + +// Schwarzschild_EF.cc +void Schwarzschild_EF_geometry(patch_system& ps, + const struct geometry_info& gi, + bool msg_flag); // Newton.cc // return true for success, false for failure to converge bool Newton_solve(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - const struct Jacobian_info& Jac_info, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, const struct solver_info& solver_info, Jacobian& Jac); diff --git a/src/gr/Newton.cc b/src/gr/Newton.cc index 4c2eb04..d82c752 100644 --- a/src/gr/Newton.cc +++ b/src/gr/Newton.cc @@ -33,7 +33,7 @@ using jtutil::error_exit; #include "../elliptic/Jacobian.hh" -#include "gfn.hh" +#include "gfns.hh" #include "AHFinderDirect.hh" //****************************************************************************** @@ -54,8 +54,8 @@ namespace { // bool Newton_solve(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - const struct Jacobian_info& Jac_info, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, const struct solver_info& solver_info, Jacobian& Jac) { @@ -67,7 +67,7 @@ bool Newton_solve(patch_system& ps, "Newton iteration %d", iteration); jtutil::norm H_norms; - horizon_function(ps, cgi, gii, true, &H_norms); + horizon_function(ps, cgi, gi, true, true, &H_norms); CCTK_VInfo(CCTK_THORNSTRING, " H rms-norm=%.2e, infinity-norm=%.2e", H_norms.rms_norm(), H_norms.infinity_norm()); @@ -81,11 +81,11 @@ bool Newton_solve(patch_system& ps, } // take a Newton step - horizon_Jacobian(ps, cgi, gii, Jac_info, Jac); + horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac); CCTK_VInfo(CCTK_THORNSTRING, "solving linear system for %d unknowns", Jac.NN()); - Jac.solve_linear_system(nominal_gfns::gfn__H, // rhs gridfn - nominal_gfns::gfn__Delta_h); // soln gridfn + Jac.solve_linear_system(gfns::gfn__H, // rhs gridfn + gfns::gfn__Delta_h); // soln gridfn CCTK_VInfo(CCTK_THORNSTRING, "done solving linear system"); jtutil::norm Delta_h_norms; @@ -99,11 +99,10 @@ bool Newton_solve(patch_system& ps, isigma <= p.max_isigma() ; ++isigma) { - const fp Delta_h = p.gridfn(nominal_gfns::gfn__Delta_h, - irho,isigma); + const fp Delta_h = p.gridfn(gfns::gfn__Delta_h, irho,isigma); Delta_h_norms.data(Delta_h); - p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) -= Delta_h; + p.ghosted_gridfn(gfns::gfn__h, irho,isigma) -= Delta_h; } } } diff --git a/src/gr/README b/src/gr/README index e44caab..5944ef3 100644 --- a/src/gr/README +++ b/src/gr/README @@ -15,8 +15,8 @@ horizon_Jacobian.cc AHFinderDirect.hh Overall header file for all the external routines in this directory. -gfn.hh - Defines the gfns of all the gridfns. +gfns.hh + Defines the gfns of all the gridfns as enums in namespace gfns::. cg.hh This is a "dangerous" header file which defines the "virtual machine" diff --git a/src/gr/Schwarzschild_EF.cc b/src/gr/Schwarzschild_EF.cc new file mode 100644 index 0000000..f0e18d4 --- /dev/null +++ b/src/gr/Schwarzschild_EF.cc @@ -0,0 +1,366 @@ +// Schwarzschild_EF.cc -- set up Schwarzschild/EF geometry variables +// $Id$ +// +// <<>> +// Schwarzschild_EF_geometry - top-level driver +/// Schwarzschild_EF_gij_xyz - (x,y,z) g_ij +/// Schwarzschild_EF_Kij_xyz - (x,y,z) g_ij +/// Schwarzschild_EF_gij_rthetaphi - (r,theta,phi) g_ij +/// Schwarzschild_EF_Kij_rthetaphi - (r,theta,phi) g_ij +/// tensor_xform_rthetaphi_to_xyz - tensor xform diag T_dd +/// + +#include +#include +#include +#include + +#include "util_Table.h" +#include "cctk.h" +#include "cctk_Arguments.h" + +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" +using jtutil::error_exit; + +#include "../util/coords.hh" +#include "../util/grid.hh" +#include "../util/fd_grid.hh" +#include "../util/patch.hh" +#include "../util/patch_edge.hh" +#include "../util/patch_interp.hh" +#include "../util/ghost_zone.hh" +#include "../util/patch_system.hh" + +#include "../elliptic/Jacobian.hh" + +#include "gfns.hh" +#include "AHFinderDirect.hh" + +//****************************************************************************** + +// +// ***** prototypes for functions local to this file ***** +// + +namespace { +void Schwarzschild_EF_gij_xyz(fp m, fp epsilon, + fp x, fp y, fp z, + fp& g_xx, fp& g_xy, fp& g_xz, + fp& g_yy, fp& g_yz, + fp& g_zz); +void Schwarzschild_EF_Kij_xyz(fp m, fp epsilon, + fp x, fp y, fp z, + fp& K_xx, fp& K_xy, fp& K_xz, + fp& K_yy, fp& K_yz, + fp& K_zz); + +void Schwarzschild_EF_gij_rthetaphi(fp m, + fp r, fp theta, fp phi, + fp& g_rr, fp& g_theta_theta); +void Schwarzschild_EF_Kij_rthetaphi(fp m, + fp r, fp theta, fp phi, + fp& K_rr, fp& K_theta_theta); + +void xform_from_rthetaphi_to_xyz(fp epsilon, + fp x, fp y, fp z, + fp T_rr, fp T_theta_theta, + fp& T_xx, fp& T_xy, fp& T_xz, + fp& T_yy, fp& T_yz, + fp& T_zz); + }; + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function sets the geometry gridfns to their analytic values +// for a unit-mass Schwarzschild spacetime in Eddington-Finkelstein +// coordinates. +// +// The (r,theta,phi) $g_{ij}$ and $K_{ij}$ are as given in appendix A2 +// of my (Jonathan Thornburg's) Ph.D thesis, available on my web page at +// http://www.aei.mpg.de/~jthorn/phd/html/phd.html +// The (r,theta,phi) --> (x,y,z) tensor transformation is as worked out +// on pages 16-17 of my AHFinderDirect working notes. +// The $\partial_k g_{kl}$ is done by numerical finite differencing. +// +// Inputs (angular gridfns, on ghosted grid): +// ... defined on ghosted grid +// ... only values on nominal grid are actually used as input +// h # shape of trial surface +// +// Inputs (angular gridfns, all on the nominal grid): +// global_[xyz] # xyz positions of grid points +// +// Outputs (angular gridfns, all on the nominal grid): +// g_dd_{11,12,13,22,23,33} # $g_{ij}$ +// K_dd_{11,12,13,22,23,33} # $K_{ij}$ +// partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k g_{ij}$ +// +void Schwarzschild_EF_geometry(patch_system& ps, + const struct geometry_info& gi, + bool msg_flag) +{ +const fp mass = gi.hardwire_Schwarzschild_EF__mass; +const fp epsilon = gi.hardwire_Schwarzschild_EF__epsilon; +const fp x_posn = gi.hardwire_Schwarzschild_EF__x_posn; +const fp y_posn = gi.hardwire_Schwarzschild_EF__y_posn; +const fp z_posn = gi.hardwire_Schwarzschild_EF__z_posn; + +if (msg_flag) + then { + CCTK_VInfo(CCTK_THORNSTRING, + " setting up Schwarzschild/EF geometry"); + CCTK_VInfo(CCTK_THORNSTRING, + " posn=(%g,%g,%g) mass=%g", + double(x_posn),double(y_posn),double(z_posn), double(mass)); + } + + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + { + patch& p = ps.ith_patch(pn); + + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + const fp r = p.ghosted_gridfn(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); + + const fp x_rel = ps.global_x_of_local_x(local_x) - x_posn; + const fp y_rel = ps.global_y_of_local_y(local_y) - y_posn; + const fp z_rel = ps.global_z_of_local_z(local_z) - z_posn; + + // compute g_ij and K_ij + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel, y_rel, z_rel, + p.gridfn(gfns::gfn__g_dd_11, irho,isigma), + p.gridfn(gfns::gfn__g_dd_12, irho,isigma), + p.gridfn(gfns::gfn__g_dd_13, irho,isigma), + p.gridfn(gfns::gfn__g_dd_22, irho,isigma), + p.gridfn(gfns::gfn__g_dd_23, irho,isigma), + p.gridfn(gfns::gfn__g_dd_33, irho,isigma)); + Schwarzschild_EF_Kij_xyz(mass, epsilon, + x_rel, y_rel, z_rel, + p.gridfn(gfns::gfn__K_dd_11, irho,isigma), + p.gridfn(gfns::gfn__K_dd_12, irho,isigma), + p.gridfn(gfns::gfn__K_dd_13, irho,isigma), + p.gridfn(gfns::gfn__K_dd_22, irho,isigma), + p.gridfn(gfns::gfn__K_dd_23, irho,isigma), + p.gridfn(gfns::gfn__K_dd_33, irho,isigma)); + + fp g11p, g12p, g13p, g22p, g23p, g33p; + fp g11m, g12m, g13m, g22m, g23m, g33m; + + // compute partial_x g_ij by finite differencing in xyz + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel+gi.Delta_xyz, y_rel, z_rel, + g11p, g12p, g13p, + g22p, g23p, + g33p); + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel-gi.Delta_xyz, y_rel, z_rel, + g11m, g12m, g13m, + g22m, g23m, + g33m); + const fp fx = 1.0 / (2.0*gi.Delta_xyz); + p.gridfn(gfns::gfn__partial_d_g_dd_111,irho,isigma) = fx * (g11p-g11m); + p.gridfn(gfns::gfn__partial_d_g_dd_112,irho,isigma) = fx * (g12p-g12m); + p.gridfn(gfns::gfn__partial_d_g_dd_113,irho,isigma) = fx * (g13p-g13m); + p.gridfn(gfns::gfn__partial_d_g_dd_122,irho,isigma) = fx * (g22p-g22m); + p.gridfn(gfns::gfn__partial_d_g_dd_123,irho,isigma) = fx * (g23p-g23m); + p.gridfn(gfns::gfn__partial_d_g_dd_133,irho,isigma) = fx * (g33p-g33m); + + // compute partial_y g_ij by finite differencing in xyz + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel, y_rel+gi.Delta_xyz, z_rel, + g11p, g12p, g13p, + g22p, g23p, + g33p); + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel, y_rel-gi.Delta_xyz, z_rel, + g11m, g12m, g13m, + g22m, g23m, + g33m); + const fp fy = 1.0 / (2.0*gi.Delta_xyz); + p.gridfn(gfns::gfn__partial_d_g_dd_211,irho,isigma) = fy * (g11p-g11m); + p.gridfn(gfns::gfn__partial_d_g_dd_212,irho,isigma) = fy * (g12p-g12m); + p.gridfn(gfns::gfn__partial_d_g_dd_213,irho,isigma) = fy * (g13p-g13m); + p.gridfn(gfns::gfn__partial_d_g_dd_222,irho,isigma) = fy * (g22p-g22m); + p.gridfn(gfns::gfn__partial_d_g_dd_223,irho,isigma) = fy * (g23p-g23m); + p.gridfn(gfns::gfn__partial_d_g_dd_233,irho,isigma) = fy * (g33p-g33m); + + // compute partial_x g_ij by finite differencing in xyz + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel, y_rel, z_rel+gi.Delta_xyz, + g11p, g12p, g13p, + g22p, g23p, + g33p); + Schwarzschild_EF_gij_xyz(mass, epsilon, + x_rel, y_rel, z_rel-gi.Delta_xyz, + g11m, g12m, g13m, + g22m, g23m, + g33m); + const fp fz = 1.0 / (2.0*gi.Delta_xyz); + p.gridfn(gfns::gfn__partial_d_g_dd_311,irho,isigma) = fz * (g11p-g11m); + p.gridfn(gfns::gfn__partial_d_g_dd_312,irho,isigma) = fz * (g12p-g12m); + p.gridfn(gfns::gfn__partial_d_g_dd_313,irho,isigma) = fz * (g13p-g13m); + p.gridfn(gfns::gfn__partial_d_g_dd_322,irho,isigma) = fz * (g22p-g22m); + p.gridfn(gfns::gfn__partial_d_g_dd_323,irho,isigma) = fz * (g23p-g23m); + p.gridfn(gfns::gfn__partial_d_g_dd_333,irho,isigma) = fz * (g33p-g33m); + } + } + } +} + +//****************************************************************************** + +// +// This function computes the Schwarzschild/Eddington-Finkelstein +// 3-metric $g_{ij}$ in $(x,y,z)$ spatial coordinates. +// +namespace { +void Schwarzschild_EF_gij_xyz(fp m, fp epsilon, + fp x, fp y, fp z, + fp& g_xx, fp& g_xy, fp& g_xz, + fp& g_yy, fp& g_yz, + fp& g_zz) +{ +fp r, theta, phi; +local_coords::r_theta_phi_of_xyz(x,y,z, r,theta,phi); + +fp g_rr, g_theta_theta, g_phi_phi; +Schwarzschild_EF_gij_rthetaphi(m, + r, theta, phi, + g_rr, g_theta_theta); +xform_from_rthetaphi_to_xyz(epsilon, + x, y, z, + g_rr, g_theta_theta, + g_xx, g_xy, g_xz, + g_yy, g_yz, + g_zz); +} + } + +//****************************************************************************** + +// +// This function computes the Schwarzschild/Eddington-Finkelstein +// 3-extrinsic curvature $K_{ij}$ in $(x,y,z)$ spatial coordinates. +// +namespace { +void Schwarzschild_EF_Kij_xyz(fp m, fp epsilon, + fp x, fp y, fp z, + fp& K_xx, fp& K_xy, fp& K_xz, + fp& K_yy, fp& K_yz, + fp& K_zz) +{ +fp r, theta, phi; +local_coords::r_theta_phi_of_xyz(x,y,z, r,theta,phi); + +fp K_rr, K_theta_theta, K_phi_phi; +Schwarzschild_EF_Kij_rthetaphi(m, + r, theta, phi, + K_rr, K_theta_theta); +xform_from_rthetaphi_to_xyz(epsilon, + x, y, z, + K_rr, K_theta_theta, + K_xx, K_xy, K_xz, + K_yy, K_yz, + K_zz); +} + } + +//****************************************************************************** + +// +// This function computes the two independent components of the 3-metric +// $g_{ij}$ for Schwarzschild spacetime in Eddington-Finkelstein coordinates +// $(r,theta,phi)$, as per equation (A2.2) of my Ph.D thesis, rescaled in +// the obvious way for non-unit mass. +// +// By the spherical symmetry, $g_{\phi\phi} = g_{\theta\theta} \sin^2\theta$. +// +namespace { +void Schwarzschild_EF_gij_rthetaphi(fp m, + fp r, fp theta, fp phi, + fp& g_rr, fp& g_theta_theta) +{ +g_rr = 1.0 + 2.0*m/r; +g_theta_theta = r*r; +} + } + +//****************************************************************************** + +// +// This function computes the two independent components of the 3-extrinsic +// curvature $K_{ij}$ for Schwarzschild spacetime in Eddington-Finkelstein +// coordinates $(r,theta,phi)$, as per equation (A2.2) of my Ph.D thesis, +// rescaled in the obvious way for non-unit mass. +// +// By the spherical symmetry, $K_{\phi\phi} = K_{\theta\theta} \sin^2\theta$. +// +// +namespace { +void Schwarzschild_EF_Kij_rthetaphi(fp m, + fp r, fp theta, fp phi, + fp& K_rr, fp& K_theta_theta) +{ +fp temp = 1.0 / sqrt(1.0 + 2.0*m/r); +K_rr = -(2.0/(r*r)) * (1.0 + m/r) * temp; +K_theta_theta = 2.0 * temp; +} + } + +//****************************************************************************** + +// +// This function transforms a diagonal T_dd (rank 2 covariant) tensor +// in spherical symmetry, +// T_dd = diag[ T_rr T_theta_theta T_theta_theta*sin(theta)**2 ] +// from (r,theta,phi) to (x,y,z) coordinates, as per pages 16-17 of my +// AHFInderDirect working notes. +// +// Arguments: +// epsilon = Tolerance parameter for deciding when to switch from +// generic expressions to z-axis limits: we switch iff +// sin^2 theta = (x^2+y^2)/r^2 is < epsilon +// +namespace { +void xform_from_rthetaphi_to_xyz(fp epsilon, + fp x, fp y, fp z, + fp T_rr, fp T_theta_theta, + fp& T_xx, fp& T_xy, fp& T_xz, + fp& T_yy, fp& T_yz, + fp& T_zz) +{ +const fp x2 = x*x; +const fp y2 = y*y; +const fp x2py2 = x2 + y2; +const fp z2 = z*z; +const fp r2 = x2 + y2 + z2; const fp r4 = r2*r2; +const bool z_flag = x2py2/r2 < epsilon; + +T_xx = (x2/r2)*T_rr + (z_flag ? 1.0/r2 : (x2*z2/r2 + y2) / (r2*x2py2)) + *T_theta_theta; +T_yy = (y2/r2)*T_rr + (z_flag ? 1.0/r2 : (y2*z2/r2 + x2) / (r2*x2py2)) + *T_theta_theta; +T_zz = (z2/r2)*T_rr + (x2py2/r4)*T_theta_theta; +T_xy = (x*y/r2)*T_rr + (z_flag ? 0.0 : (z2/r2 - 1.0) * x*y / (r2*x2py2)) + *T_theta_theta; +T_xz = (x*z/r2)*T_rr - (x*z/r4)*T_theta_theta; +T_yz = (y*z/r2)*T_rr - (y*z/r4)*T_theta_theta; +} + } diff --git a/src/gr/cg.hh b/src/gr/cg.hh index 6535753..b63ff2a 100644 --- a/src/gr/cg.hh +++ b/src/gr/cg.hh @@ -10,7 +10,7 @@ // amount of other code affected by this file). // // prerequisites -// gfn.hh +// gfns.hh // //****************************************************************************** @@ -33,15 +33,15 @@ // partial derivatives // #define PARTIAL_RHO(ghosted_gridfn_name) \ - p.partial_rho(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) + p.partial_rho(gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) #define PARTIAL_SIGMA(ghosted_gridfn_name) \ - p.partial_sigma(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) + p.partial_sigma(gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) #define PARTIAL_RHO_RHO(ghosted_gridfn_name) \ - p.partial_rho_rho(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) + p.partial_rho_rho(gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) #define PARTIAL_RHO_SIGMA(ghosted_gridfn_name) \ - p.partial_rho_sigma(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) -#define PARTIAL_SIGMA_SIGMA(ghosted_gridfn_name) \ - p.partial_sigma_sigma(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) + p.partial_rho_sigma(gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) +#define PARTIAL_SIGMA_SIGMA(ghosted_gridfn_name)\ + p.partial_sigma_sigma(gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) //****************************************************************************** @@ -49,7 +49,7 @@ // ghosted-grid gridfns // n.b. since we're always evaluating on the surface, r == h // -#define h p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) +#define h p.ghosted_gridfn(gfns::gfn__h, irho,isigma) #define r h //****************************************************************************** @@ -57,53 +57,53 @@ // // nominal-grid gridfns // -#define g_dd_11 p.gridfn(nominal_gfns::gfn__g_dd_11, irho,isigma) -#define g_dd_12 p.gridfn(nominal_gfns::gfn__g_dd_12, irho,isigma) -#define g_dd_13 p.gridfn(nominal_gfns::gfn__g_dd_13, irho,isigma) -#define g_dd_22 p.gridfn(nominal_gfns::gfn__g_dd_22, irho,isigma) -#define g_dd_23 p.gridfn(nominal_gfns::gfn__g_dd_23, irho,isigma) -#define g_dd_33 p.gridfn(nominal_gfns::gfn__g_dd_33, irho,isigma) -#define K_dd_11 p.gridfn(nominal_gfns::gfn__K_dd_11, irho,isigma) -#define K_dd_12 p.gridfn(nominal_gfns::gfn__K_dd_12, irho,isigma) -#define K_dd_13 p.gridfn(nominal_gfns::gfn__K_dd_13, irho,isigma) -#define K_dd_22 p.gridfn(nominal_gfns::gfn__K_dd_22, irho,isigma) -#define K_dd_23 p.gridfn(nominal_gfns::gfn__K_dd_23, irho,isigma) -#define K_dd_33 p.gridfn(nominal_gfns::gfn__K_dd_33, irho,isigma) - -#define partial_d_g_dd_111 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_111, irho,isigma) -#define partial_d_g_dd_112 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_112, irho,isigma) -#define partial_d_g_dd_113 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_113, irho,isigma) -#define partial_d_g_dd_122 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_122, irho,isigma) -#define partial_d_g_dd_123 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_123, irho,isigma) -#define partial_d_g_dd_133 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_133, irho,isigma) -#define partial_d_g_dd_211 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_211, irho,isigma) -#define partial_d_g_dd_212 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_212, irho,isigma) -#define partial_d_g_dd_213 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_213, irho,isigma) -#define partial_d_g_dd_222 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_222, irho,isigma) -#define partial_d_g_dd_223 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_223, irho,isigma) -#define partial_d_g_dd_233 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_233, irho,isigma) -#define partial_d_g_dd_311 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_311, irho,isigma) -#define partial_d_g_dd_312 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_312, irho,isigma) -#define partial_d_g_dd_313 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_313, irho,isigma) -#define partial_d_g_dd_322 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_322, irho,isigma) -#define partial_d_g_dd_323 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_323, irho,isigma) -#define partial_d_g_dd_333 p.gridfn(nominal_gfns::gfn__partial_d_g_dd_333, irho,isigma) - -#define H p.gridfn(nominal_gfns::gfn__H, irho,isigma) +#define g_dd_11 p.gridfn(gfns::gfn__g_dd_11, irho,isigma) +#define g_dd_12 p.gridfn(gfns::gfn__g_dd_12, irho,isigma) +#define g_dd_13 p.gridfn(gfns::gfn__g_dd_13, irho,isigma) +#define g_dd_22 p.gridfn(gfns::gfn__g_dd_22, irho,isigma) +#define g_dd_23 p.gridfn(gfns::gfn__g_dd_23, irho,isigma) +#define g_dd_33 p.gridfn(gfns::gfn__g_dd_33, irho,isigma) +#define K_dd_11 p.gridfn(gfns::gfn__K_dd_11, irho,isigma) +#define K_dd_12 p.gridfn(gfns::gfn__K_dd_12, irho,isigma) +#define K_dd_13 p.gridfn(gfns::gfn__K_dd_13, irho,isigma) +#define K_dd_22 p.gridfn(gfns::gfn__K_dd_22, irho,isigma) +#define K_dd_23 p.gridfn(gfns::gfn__K_dd_23, irho,isigma) +#define K_dd_33 p.gridfn(gfns::gfn__K_dd_33, irho,isigma) + +#define partial_d_g_dd_111 p.gridfn(gfns::gfn__partial_d_g_dd_111, irho,isigma) +#define partial_d_g_dd_112 p.gridfn(gfns::gfn__partial_d_g_dd_112, irho,isigma) +#define partial_d_g_dd_113 p.gridfn(gfns::gfn__partial_d_g_dd_113, irho,isigma) +#define partial_d_g_dd_122 p.gridfn(gfns::gfn__partial_d_g_dd_122, irho,isigma) +#define partial_d_g_dd_123 p.gridfn(gfns::gfn__partial_d_g_dd_123, irho,isigma) +#define partial_d_g_dd_133 p.gridfn(gfns::gfn__partial_d_g_dd_133, irho,isigma) +#define partial_d_g_dd_211 p.gridfn(gfns::gfn__partial_d_g_dd_211, irho,isigma) +#define partial_d_g_dd_212 p.gridfn(gfns::gfn__partial_d_g_dd_212, irho,isigma) +#define partial_d_g_dd_213 p.gridfn(gfns::gfn__partial_d_g_dd_213, irho,isigma) +#define partial_d_g_dd_222 p.gridfn(gfns::gfn__partial_d_g_dd_222, irho,isigma) +#define partial_d_g_dd_223 p.gridfn(gfns::gfn__partial_d_g_dd_223, irho,isigma) +#define partial_d_g_dd_233 p.gridfn(gfns::gfn__partial_d_g_dd_233, irho,isigma) +#define partial_d_g_dd_311 p.gridfn(gfns::gfn__partial_d_g_dd_311, irho,isigma) +#define partial_d_g_dd_312 p.gridfn(gfns::gfn__partial_d_g_dd_312, irho,isigma) +#define partial_d_g_dd_313 p.gridfn(gfns::gfn__partial_d_g_dd_313, irho,isigma) +#define partial_d_g_dd_322 p.gridfn(gfns::gfn__partial_d_g_dd_322, irho,isigma) +#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(nominal_gfns::gfn__partial_H_wrt_partial_d_h_1, irho,isigma) + p.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1, irho,isigma) #define partial_H_wrt_partial_d_h_2 \ - p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_2, irho,isigma) + p.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2, irho,isigma) #define partial_H_wrt_partial_dd_h_11 \ - p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_11, irho,isigma) + p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11, irho,isigma) #define partial_H_wrt_partial_dd_h_12 \ - p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_12, irho,isigma) + p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12, irho,isigma) #define partial_H_wrt_partial_dd_h_22 \ - p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma) + p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma) -#define save_H p.gridfn(nominal_gfns::gfn__save_H, irho,isigma) -#define Delta_h p.gridfn(nominal_gfns::gfn__Delta_h, irho,isigma) +#define save_H p.gridfn(gfns::gfn__save_H, irho,isigma) +#define Delta_h p.gridfn(gfns::gfn__Delta_h, irho,isigma) //****************************************************************************** diff --git a/src/gr/driver.cc b/src/gr/driver.cc index 00695e9..a9b9a13 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -3,6 +3,7 @@ // // <<>> // AHFinderDirect_driver - top-level driver +/// decode_Jacobian_method - decode the Jacobian_method parameter /// /// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords) /// setup_ellipsoid - setup up a coordiante ellipsoid in h @@ -39,7 +40,7 @@ using jtutil::error_exit; #include "../elliptic/Jacobian.hh" -#include "gfn.hh" +#include "gfns.hh" #include "AHFinderDirect.hh" //****************************************************************************** @@ -49,8 +50,10 @@ using jtutil::error_exit; // namespace { +enum Jacobian_method + decode_Jacobian_method(const char Jacobian_method_string[]); void setup_Kerr_horizon(patch_system& ps, - fp x_center, fp y_center, fp z_center, + fp x_posn, fp y_posn, fp z_posn, fp m, fp a, bool Kerr_Schild_flag); void setup_ellipsoid(patch_system& ps, @@ -58,7 +61,7 @@ void setup_ellipsoid(patch_system& ps, fp x_radius, fp y_radius, fp z_radius); void print_Jacobians(const patch_system& ps, - const Jacobian& Jac, const Jacobian& NP_Jac, + const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr, const char file_name[]); }; @@ -77,18 +80,25 @@ CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures"); // -// set up the geometry interpolator +// set up the geometry parameters and the geometry interpolator // -struct geometry_interpolator_info gii; +struct geometry_info gi; +gi.hardwire_Schwarzschild_EF = (hardwire_Schwarzschild_EF != 0); +gi.hardwire_Schwarzschild_EF__x_posn = hardwire_Schwarzschild_EF__x_posn; +gi.hardwire_Schwarzschild_EF__y_posn = hardwire_Schwarzschild_EF__y_posn; +gi.hardwire_Schwarzschild_EF__z_posn = hardwire_Schwarzschild_EF__z_posn; +gi.hardwire_Schwarzschild_EF__mass = hardwire_Schwarzschild_EF__mass; +gi.hardwire_Schwarzschild_EF__epsilon = hardwire_Schwarzschild_EF__epsilon; +gi.Delta_xyz = hardwire_Schwarzschild_EF__Delta_xyz; + CCTK_VInfo(CCTK_THORNSTRING, " setting up geometry interpolator"); -gii.operator_handle = CCTK_InterpHandle(geometry_interpolator_name); -if (gii.operator_handle < 0) +gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name); +if (gi.operator_handle < 0) then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "couldn't find interpolator \"%s\"!", geometry_interpolator_name); /*NOTREACHED*/ - -gii.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars); -if (gii.param_table_handle < 0) +gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars); +if (gi.param_table_handle < 0) then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "bad geometry-interpolator parameter(s) \"%s\"!", geometry_interpolator_pars); /*NOTREACHED*/ @@ -150,8 +160,8 @@ cgi.K_dd_33_data = static_cast(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz")); patch_system ps(origin_x, origin_y, origin_z, patch_system::type_of_name(patch_system_type), N_ghost_points, N_overlap_points, delta_drho_dsigma, - nominal_gfns::min_gfn, nominal_gfns::max_gfn, - ghosted_gfns::min_gfn, ghosted_gfns::max_gfn, + gfns::nominal_min_gfn, gfns::nominal_max_gfn, + gfns::ghosted_min_gfn, gfns::ghosted_max_gfn, interp_handle, interp_param_table_handle); @@ -163,7 +173,7 @@ if (STRING_EQUAL(initial_guess_method, "read from file")) CCTK_VInfo(CCTK_THORNSTRING, "reading initial guess h from \"%s\"", h_file_name); - ps.read_ghosted_gridfn(ghosted_gfns::gfn__h, + ps.read_ghosted_gridfn(gfns::gfn__h, h_file_name, false); // no ghost zones } @@ -178,19 +188,19 @@ if (STRING_EQUAL(initial_guess_method, "read from file")) initial_guess__ellipsoid__z_radius); else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr")) then setup_Kerr_horizon(ps, - initial_guess__Kerr_KerrSchild__x_center, - initial_guess__Kerr_KerrSchild__y_center, - initial_guess__Kerr_KerrSchild__z_center, - initial_guess__Kerr_KerrSchild__mass, - initial_guess__Kerr_KerrSchild__spin, + initial_guess__Kerr__x_posn, + initial_guess__Kerr__y_posn, + initial_guess__Kerr__z_posn, + initial_guess__Kerr__mass, + initial_guess__Kerr__spin, false); // use Kerr coords else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild")) then setup_Kerr_horizon(ps, - initial_guess__Kerr_KerrSchild__x_center, - initial_guess__Kerr_KerrSchild__y_center, - initial_guess__Kerr_KerrSchild__z_center, - initial_guess__Kerr_KerrSchild__mass, - initial_guess__Kerr_KerrSchild__spin, + initial_guess__Kerr__x_posn, + initial_guess__Kerr__y_posn, + initial_guess__Kerr__z_posn, + initial_guess__Kerr__mass, + initial_guess__Kerr__spin, true); // use Kerr-Schild coords else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown initial_guess_method=\"%s\"!", @@ -200,31 +210,36 @@ if (STRING_EQUAL(initial_guess_method, "read from file")) CCTK_VInfo(CCTK_THORNSTRING, "writing initial guess h to \"%s\"", h_file_name); - ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, - true, ghosted_gfns::gfn__h, + ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h, + true, gfns::gfn__h, h_file_name, false); // no ghost zones } // -// find the apparent horizon +// decode/copy parameters into structures // -struct Jacobian_info Jac_info; -Jac_info.perturbation_amplitude = NP_Jacobian__perturbation_amplitude; + +struct Jacobian_info Jacobian_info; +Jacobian_info.Jacobian_method = decode_Jacobian_method(Jacobian_method); +Jacobian_info.perturbation_amplitude = Jacobian_perturbation_amplitude; struct solver_info solver_info; solver_info.max_Newton_iterations = max_Newton_iterations; solver_info.H_norm_for_convergence = H_norm_for_convergence; solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence; +// +// find the apparent horizon +// if (STRING_EQUAL(method, "horizon function")) then { jtutil::norm H_norms; const int timer_handle = CCTK_TimerCreateI(); CCTK_TimerStartI(timer_handle); - horizon_function(ps, cgi, gii, false, &H_norms); + horizon_function(ps, cgi, gi, false, true, &H_norms); CCTK_TimerStopI(timer_handle); printf("timer stats for evaluating H(h):\n"); CCTK_TimerPrintDataI(timer_handle, -1); @@ -235,33 +250,36 @@ if (STRING_EQUAL(method, "horizon function")) CCTK_VInfo(CCTK_THORNSTRING, "writing H(h) to \"%s\"", H_of_h_file_name); - ps.print_gridfn_with_xyz(nominal_gfns::gfn__H, - true, ghosted_gfns::gfn__h, + ps.print_gridfn_with_xyz(gfns::gfn__H, + true, gfns::gfn__h, H_of_h_file_name); } else if (STRING_EQUAL(method, "Jacobian test")) then { - jtutil::norm H_norms; - Jacobian& Jac = create_Jacobian(ps, Jacobian_type); - horizon_function(ps, cgi, gii, true, &H_norms); - horizon_Jacobian(ps, cgi, gii, Jac_info, Jac); - - Jacobian& NP_Jac = create_Jacobian(ps, Jacobian_type); - horizon_function(ps, cgi, gii, true, &H_norms); - horizon_Jacobian_NP(ps, cgi, gii, Jac_info, NP_Jac); + // symbolic differentiation with finite diff d/dr + Jacobian& Jac_SD_FDdr = new_Jacobian(ps, Jacobian_storage_method); + horizon_function(ps, cgi, gi); + Jacobian_info.Jacobian_method = symbolic_differentiation_with_FD_dr; + horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_SD_FDdr); + + // numerical perturbation + Jacobian& Jac_NP = new_Jacobian(ps, Jacobian_storage_method); + horizon_function(ps, cgi, gi); + Jacobian_info.Jacobian_method = numerical_perturbation; + horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_NP); CCTK_VInfo(CCTK_THORNSTRING, "writing Jacobians to \"%s\"", Jacobian_file_name); - print_Jacobians(ps, Jac, NP_Jac, Jacobian_file_name); + print_Jacobians(ps, Jac_NP, Jac_SD_FDdr, Jacobian_file_name); } else if (STRING_EQUAL(method, "Newton solve")) then { - Jacobian& Jac = create_Jacobian(ps, Jacobian_type); + Jacobian& Jac = new_Jacobian(ps, Jacobian_storage_method); const int timer_handle = CCTK_TimerCreateI(); CCTK_TimerStartI(timer_handle); - Newton_solve(ps, cgi, gii, Jac_info, solver_info, Jac); + Newton_solve(ps, cgi, gi, Jacobian_info, solver_info, Jac); CCTK_TimerStopI(timer_handle); printf("timer stats for finding the apparent horizon:\n"); CCTK_TimerPrintDataI(timer_handle, -1); @@ -269,15 +287,15 @@ else if (STRING_EQUAL(method, "Newton solve")) CCTK_VInfo(CCTK_THORNSTRING, "writing final horizon shape h to \"%s\"", h_file_name); - ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, - true, ghosted_gfns::gfn__h, + ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h, + true, gfns::gfn__h, h_file_name, false); // no ghost zones CCTK_VInfo(CCTK_THORNSTRING, "writing H(h) to \"%s\"", H_of_h_file_name); - ps.print_gridfn_with_xyz(nominal_gfns::gfn__H, - true, ghosted_gfns::gfn__h, + ps.print_gridfn_with_xyz(gfns::gfn__H, + true, gfns::gfn__h, H_of_h_file_name); } else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, @@ -285,6 +303,31 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, method); /*NOTREACHED*/ } +//****************************************************************************** + +// +// This function decodes the Jacobian_method parameter (string) into +// an internal enum for future use. +// +namespace { +enum Jacobian_method + decode_Jacobian_method(const char Jacobian_method_string[]) +{ +if (STRING_EQUAL(Jacobian_method_string, + "numerical perturbation")) + then return numerical_perturbation; +else if (STRING_EQUAL(Jacobian_method_string, + "symbolic differentiation with finite diff d/dr")) + then return symbolic_differentiation_with_FD_dr; +else if (STRING_EQUAL(Jacobian_method_string, + "symbolic differentiation")) + then return symbolic_differentiation; +else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + "unknown Jacobian_method_string=\"%s\"!", + Jacobian_method_string); /*NOTREACHED*/ +} + } + //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -297,7 +340,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // and the horizon is worked out on page 13.2 of my AHFinderDirect notes. // // Arguments: -// [xyz]_center = The position of the Kerr black hole. +// [xyz]_posn = The position of the Kerr black hole. // (m,a) = Describe the Kerr black hole. Note that my convention has // a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a). // Kerr_Schild_flag = false to use Kerr coordinates, @@ -305,7 +348,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // namespace { void setup_Kerr_horizon(patch_system& ps, - fp x_center, fp y_center, fp z_center, + fp x_posn, fp y_posn, fp z_posn, fp m, fp a, bool Kerr_Schild_flag) { @@ -315,9 +358,9 @@ CCTK_VInfo(CCTK_THORNSTRING, "setting up horizon for Kerr in %s coords", name); CCTK_VInfo(CCTK_THORNSTRING, - " mass=%g, spin=J/m^2=%g, posn=(%g,%g,%g)", - double(m), double(a), - double(x_center), double(y_center), double(z_center)); + " posn=(%g,%g,%g) mass=%g spin=J/m^2=%g", + double(x_posn), double(y_posn), double(z_posn), + double(m), double(a)); // horizon in Kerr coordinates is coordinate sphere const fp r = m * (1.0 + sqrt(1.0 - a*a)); @@ -330,7 +373,7 @@ CCTK_VInfo(CCTK_THORNSTRING, " horizon is coordinate %s", Kerr_Schild_flag ? "ellipsoid" : "sphere"); setup_ellipsoid(ps, - x_center, y_center, z_center, + x_posn, y_posn, z_posn, xy_radius, xy_radius, z_radius); } } @@ -417,7 +460,7 @@ CCTK_VInfo(CCTK_THORNSTRING, /*NOTREACHED*/ // r = horizon radius at this grid point - p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r; + p.ghosted_gridfn(gfns::gfn__h, irho,isigma) = r; } } } @@ -434,7 +477,7 @@ CCTK_VInfo(CCTK_THORNSTRING, // namespace { void print_Jacobians(const patch_system& ps, - const Jacobian& Jac, const Jacobian& NP_Jac, + const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr, const char file_name[]) { FILE *fileptr = fopen(file_name, "w"); @@ -451,10 +494,10 @@ fprintf(fileptr, "# column 5 = y JJ\n"); fprintf(fileptr, "# column 6 = y patch number\n"); fprintf(fileptr, "# column 7 = y irho\n"); fprintf(fileptr, "# column 8 = y isigma\n"); -fprintf(fileptr, "# column 9 = Jac(II,JJ)\n"); -fprintf(fileptr, "# column 10 = NP_Jac(II,JJ)\n"); -fprintf(fileptr, "# column 11 = abs error = Jac - NP_Jac\n"); -fprintf(fileptr, "# column 12 = rel error = abs error / max(|Jac|,|NP_Jac|)\n"); +fprintf(fileptr, "# column 9 = Jac_NP(II,JJ)\n"); +fprintf(fileptr, "# column 10 = Jac_SD_FDdr(II,JJ)\n"); +fprintf(fileptr, "# column 11 = abs error in Jac_SD_FDdr\n"); +fprintf(fileptr, "# column 12 = rel error in Jac_SD_FDdr\n"); for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) { @@ -482,26 +525,27 @@ fprintf(fileptr, "# column 12 = rel error = abs error / max(|Jac|,|NP_Jac|)\n"); { const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); - if (! Jac.is_explicitly_stored(II,JJ)) + if (! Jac_SD_FDdr.is_explicitly_stored(II,JJ)) then continue; // skip sparse points - const fp SD = Jac(II,JJ); - const fp NP = NP_Jac(II,JJ); - const fp abs_error = SD - NP; + const fp NP = Jac_NP (II,JJ); + const fp SD_FDdr = Jac_SD_FDdr(II,JJ); - if ((SD == 0.0) && (NP == 0.0)) + if ((NP == 0.0) && (SD_FDdr == 0.0)) then continue; // skip zero values (if == ) - const fp abs_SD = jtutil::abs(SD); - const fp abs_NP = jtutil::abs(NP); - const fp rel_error = abs_error / jtutil::max(abs_SD, abs_NP); + const fp abs_NP = jtutil::abs(NP ); + const fp abs_SD_FDdr = jtutil::abs(SD_FDdr); + const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP); + const fp SD_FDdr_abs_error = SD_FDdr - NP; + const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs; fprintf(fileptr, "%d %d %d %d\t%d %d %d %d\t%.10g\t%.10g\t%e\t%e\n", II, xpn, x_irho, x_isigma, JJ, ypn, y_irho, y_isigma, - double(SD), double(NP), - double(abs_error), double(rel_error)); + double(NP), double(SD_FDdr), + double(SD_FDdr_abs_error), double(SD_FDdr_rel_error)); } } } diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index f16d9c5..0b9e0e0 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -2,13 +2,14 @@ // $Id$ // // <<>> -// create_Jacobian - create a Jacobian matrix of the appropriate type // -// horizon_Jacobian - compute the Jacobian (symbolic differentiation) +// horizon_Jacobian - top-level driver to compute the Jacobian +/// +/// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation +/// +/// horizon_Jacobian_SD - compute the Jacobian (symbolic differentiation) /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian -// -// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation -// +/// #include #include @@ -38,7 +39,7 @@ using jtutil::error_exit; #include "../elliptic/Jacobian.hh" -#include "gfn.hh" +#include "gfns.hh" #include "AHFinderDirect.hh" //****************************************************************************** @@ -48,6 +49,17 @@ using jtutil::error_exit; // namespace { +void horizon_Jacobian_NP(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& ii, + const struct Jacobian_info& Jacobian_info, + Jacobian& Jac); + +void horizon_Jacobian_SD(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + Jacobian& Jac); void add_ghost_zone_Jacobian(const patch_system& ps, const patch& xp, const ghost_zone& xmgz, int x_II, @@ -55,27 +67,118 @@ void add_ghost_zone_Jacobian(const patch_system& ps, fp mol, Jacobian& Jac); } +//****************************************************************************** + +// +// 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. +void horizon_Jacobian(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + Jacobian& Jac) +{ +switch (Jacobian_info.Jacobian_method) + { +case numerical_perturbation: + horizon_Jacobian_NP(ps, cgi, gi, Jacobian_info, Jac); + break; +case symbolic_differentiation_with_FD_dr: + horizon_Jacobian_SD(ps, cgi, gi, Jacobian_info, Jac); + break; +case symbolic_differentiation: + // FIXME FIXME: need true symbolic diff for d/dr terms + horizon_Jacobian_SD(ps, cgi, gi, Jacobian_info, Jac); + break; +default: + error_exit(PANIC_EXIT, +"***** horizon_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n" +" (this should never happen!)\n" +, + int(Jacobian_info.Jacobian_method)); /*NOTREACHED*/ + } +} + //****************************************************************************** //****************************************************************************** //****************************************************************************** // -// This function decodes the Jacobian_type parameter and creates the -// appropriate type of Jacobian matrix object. +// This function computes the Jacobian matrix of the LHS function H(h) +// by numerical perturbation of the H(h) function. The algorithm is as +// follows: // -// FIXME: the patch system shouldn't really have to be non-const, but -// the Jacobian constructors all require this to allow the -// linear solvers to directly update gridfns +// we assume that H = H(h) has already been evaluated +// save_H = H +// for each point (y,JJ) +// { +// const fp save_h_y = h at y; +// h at y += perturbation_amplitude; +// evaluate H(h) (silently) +// for each point (x,II) +// { +// Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude; +// } +// h at y = save_h_y; +// } +// H = save_H // -Jacobian& create_Jacobian(patch_system& ps, - const char Jacobian_type[]) +namespace { +void horizon_Jacobian_NP(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& ii, + const struct Jacobian_info& Jacobian_info, + Jacobian& Jac) { -if (STRING_EQUAL(Jacobian_type, "dense matrix")) - then return *new dense_Jacobian(ps); -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "unknown Jacobian_type=\"%s\"!", - Jacobian_type); /*NOTREACHED*/ +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); + + for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) + { + patch& yp = ps.ith_patch(ypn); + CCTK_VInfo(CCTK_THORNSTRING, " perturbing in %s patch", yp.name()); + + for (int y_irho = yp.min_irho() ; y_irho <= yp.max_irho() ; ++y_irho) + { + for (int y_isigma = yp.min_isigma() ; + y_isigma <= yp.max_isigma() ; + ++y_isigma) + { + const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); + + 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; + horizon_function(ps, cgi, ii); + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) + { + patch& xp = ps.ith_patch(xpn); + + for (int x_irho = xp.min_irho() ; + x_irho <= xp.max_irho() ; + ++x_irho) + { + for (int x_isigma = xp.min_isigma() ; + x_isigma <= xp.max_isigma() ; + ++x_isigma) + { + 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; + } + } + } + yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y; + } + } + } + +ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H); } + } //****************************************************************************** //****************************************************************************** @@ -106,13 +209,14 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. // -void horizon_Jacobian(patch_system& ps, - const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - const struct Jacobian_info& Jac_info, - Jacobian& Jac) +namespace { +void horizon_Jacobian_SD(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + const struct Jacobian_info& Jacobian_info, + Jacobian& Jac) { -CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian (symbolic differentiation)"); ps.compute_synchronize_Jacobian(); Jac.zero(); @@ -124,10 +228,10 @@ Jac.zero(); CCTK_VInfo(CCTK_THORNSTRING, " computing d/dr terms by numerical perturbation"); -const fp epsilon = Jac_info.perturbation_amplitude; -ps.gridfn_copy(nominal_gfns::gfn__H, nominal_gfns::gfn__save_H); -ps.add_to_ghosted_gridfn(epsilon, ghosted_gfns::gfn__h); -horizon_function(ps, cgi, gii, false, NULL, false); +const fp epsilon = Jacobian_info.perturbation_amplitude; +ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H); +ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h); +horizon_function(ps, cgi, gi); for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -138,15 +242,14 @@ horizon_function(ps, cgi, gii, false, NULL, false); ++isigma) { const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma); - const fp old_H = p.gridfn(nominal_gfns::gfn__save_H, - irho,isigma); - const fp new_H = p.gridfn(nominal_gfns::gfn__H, 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; } } } -ps.add_to_ghosted_gridfn(-epsilon, ghosted_gfns::gfn__h); -ps.gridfn_copy(nominal_gfns::gfn__save_H, nominal_gfns::gfn__H); +ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h); +ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H); // @@ -183,19 +286,19 @@ CCTK_VInfo(CCTK_THORNSTRING, // Jacobian coefficients for this point const fp Jacobian_coeff_rho - = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_1, + = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1, x_irho, x_isigma); const fp Jacobian_coeff_sigma - = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_2, + = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2, x_irho, x_isigma); const fp Jacobian_coeff_rho_rho - = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_11, + = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11, x_irho, x_isigma); const fp Jacobian_coeff_rho_sigma - = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_12, + = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12, x_irho, x_isigma); const fp Jacobian_coeff_sigma_sigma - = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, + = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22, x_irho, x_isigma); // partial_rho, partial_rho_rho @@ -265,6 +368,7 @@ CCTK_VInfo(CCTK_THORNSTRING, } } } + } //****************************************************************************** @@ -318,84 +422,3 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn } } } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function computes the Jacobian matrix of the LHS function H(h) -// by numerical perturbation of the H(h) function. The algorithm is as -// follows: -// -// we assume that H = H(h) has already been evaluated -// save_H = H -// for each point (y,JJ) -// { -// const fp save_h_y = h at y; -// h at y += perturbation_amplitude; -// evaluate H(h) (silently) -// for each point (x,II) -// { -// Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude; -// } -// h at y = save_h_y; -// } -// H = save_H -// -void horizon_Jacobian_NP(patch_system& ps, - const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& ii, - const struct Jacobian_info& Jac_info, - Jacobian& Jac) -{ -CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_NP"); -const fp epsilon = Jac_info.perturbation_amplitude; - -ps.gridfn_copy(nominal_gfns::gfn__H, nominal_gfns::gfn__save_H); - - for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) - { - patch& yp = ps.ith_patch(ypn); - CCTK_VInfo(CCTK_THORNSTRING, " perturbing in %s patch", yp.name()); - - for (int y_irho = yp.min_irho() ; y_irho <= yp.max_irho() ; ++y_irho) - { - for (int y_isigma = yp.min_isigma() ; - y_isigma <= yp.max_isigma() ; - ++y_isigma) - { - const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); - - const fp save_h_y = yp.ghosted_gridfn(ghosted_gfns::gfn__h, - y_irho,y_isigma); - yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) += epsilon; - horizon_function(ps, cgi, ii, false, NULL, false); - for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) - { - patch& xp = ps.ith_patch(xpn); - - for (int x_irho = xp.min_irho() ; - x_irho <= xp.max_irho() ; - ++x_irho) - { - for (int x_isigma = xp.min_isigma() ; - x_isigma <= xp.max_isigma() ; - ++x_isigma) - { - const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); - const fp old_H = xp.gridfn(nominal_gfns::gfn__save_H, - x_irho,x_isigma); - const fp new_H = xp.gridfn(nominal_gfns::gfn__H, - x_irho,x_isigma); - Jac(II,JJ) = (new_H - old_H) / epsilon; - } - } - } - yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) = save_h_y; - } - } - } - -ps.gridfn_copy(nominal_gfns::gfn__save_H, nominal_gfns::gfn__H); -} diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 9f05698..02310ce 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -4,9 +4,9 @@ // <<>> // horizon_function - top-level driver /// setup_xyz_posns - setup global xyz posns of grid points -/// interpolate_geometry - interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid +/// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid /// compute_H - compute H(h) given earlier setup -// +/// #include #include @@ -24,6 +24,8 @@ #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" using jtutil::error_exit; +using jtutil::pow2; +using jtutil::pow4; #include "../util/coords.hh" #include "../util/grid.hh" @@ -36,7 +38,7 @@ using jtutil::error_exit; #include "../elliptic/Jacobian.hh" -#include "gfn.hh" +#include "gfns.hh" #include "AHFinderDirect.hh" //****************************************************************************** @@ -49,8 +51,12 @@ namespace { void setup_xyz_posns(patch_system& ps, bool msg_flag); void interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& ii, + const struct geometry_info& gi, bool msg_flag); +void Schwarzschild_EF_geometry(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_info& gi, + bool msg_flag); void compute_H(patch_system& ps, bool Jacobian_flag, jtutil::norm* H_norms_ptr, @@ -58,8 +64,6 @@ void compute_H(patch_system& ps, } //****************************************************************************** -//****************************************************************************** -//****************************************************************************** // // This function computes the LHS function H(h), and optionally also @@ -77,7 +81,7 @@ void compute_H(patch_system& ps, // // Outputs (temporaries computed at each grid point) // ## computed by hand-written code -// xx,yy,zz # xyz positions of grid points +// global_[xyz] # xyz positions of grid points // X_ud_*, X_udd_* # xyz derivative coefficients // ## computed by Maple-generated code // g_uu_{11,12,13,22,23,33} # $g^{ij}$ @@ -96,19 +100,19 @@ void compute_H(patch_system& ps, // Arguments: // Jacobian_flag = true to compute the Jacobian coefficients, // false to skip this. +// 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. -// msg_flag = true to print status messages, -// false to skip this. // void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, - bool Jacobian_flag, - jtutil::norm* H_norms_ptr, - bool msg_flag = true) + const struct geometry_info& gi, + bool Jacobian_flag = false, + bool msg_flag = false, + jtutil::norm* H_norms_ptr = NULL) { if (msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " horizon function"); @@ -119,8 +123,10 @@ ps.synchronize(); // set up xyz positions of grid points setup_xyz_posns(ps, msg_flag); -// interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$ -interpolate_geometry(ps, cgi, gii, msg_flag); +// compute g_ij, K_ij, and partial_k g_ij +if (gi.hardwire_Schwarzschild_EF) + then Schwarzschild_EF_geometry(ps, gi, msg_flag); + else interpolate_geometry(ps, cgi, gi, msg_flag); // compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing @@ -150,8 +156,7 @@ if (msg_flag) isigma <= p.max_isigma() ; ++isigma) { - const fp r = p.ghosted_gridfn(ghosted_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; @@ -160,9 +165,9 @@ if (msg_flag) const fp global_y = ps.global_y_of_local_y(local_y); const fp global_z = ps.global_z_of_local_z(local_z); - p.gridfn(nominal_gfns::gfn__global_x, irho,isigma) = global_x; - p.gridfn(nominal_gfns::gfn__global_y, irho,isigma) = global_y; - p.gridfn(nominal_gfns::gfn__global_z, irho,isigma) = global_z; + p.gridfn(gfns::gfn__global_x, irho,isigma) = global_x; + p.gridfn(gfns::gfn__global_y, irho,isigma) = global_y; + p.gridfn(gfns::gfn__global_z, irho,isigma) = global_z; } } } @@ -172,9 +177,9 @@ if (msg_flag) //****************************************************************************** // -// This function interpolates the Cactus $g_{ij}$ and $K_{ij}$, and -// the spatial derivatives $\partial_k g_{ij}$, to the trial horizon -// surface position given by h. +// This function interpolates the Cactus gridfns $g_{ij}$ and $K_{ij}$ +// to determine $g_{ij}$, $K_{ij}$, and the spatial derivatives +// $\partial_k g_{ij}$, on the trial horizon surface position given by h. // // Inputs (angular gridfns, on ghosted grid): // ... defined on ghosted grid @@ -182,16 +187,12 @@ if (msg_flag) // h # shape of trial surface // // Inputs (angular gridfns, all on the nominal grid): -// xx,yy,zz # xyz positions of grid points +// global_[xyz] # xyz positions of grid points // // Inputs (Cactus 3-D gridfns): // gxx,gxy,gxz,gyy,gyz,gzz # 3-metric $g_{ij}$ // kxx,kxy,kxz,kyy,kyz,kzz # extrinsic curvature $K_{ij}$ // -// Inputs (angular gridfns, all on the nominal grid): -// ## computed directly from h -// xx,yy,zz # xyz positions -// // Outputs (angular gridfns, all on the nominal grid): // g_dd_{11,12,13,22,23,33} # $g_{ij}$ // K_dd_{11,12,13,22,23,33} # $K_{ij}$ @@ -203,12 +204,12 @@ if (msg_flag) namespace { void interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii, + const struct geometry_info& gi, bool msg_flag) { if (msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid"); + " interpolating g_ij and Kij from Cactus 3-D grid"); int status; @@ -216,19 +217,19 @@ const int N_interp_points = ps.N_grid_points(); const int interp_coords_type_code = CCTK_VARIABLE_REAL; const void* const interp_coords[N_GRID_DIMS] = { - static_cast(ps.gridfn_data(nominal_gfns::gfn__global_x)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__global_y)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__global_z)), + static_cast(ps.gridfn_data(gfns::gfn__global_x)), + static_cast(ps.gridfn_data(gfns::gfn__global_y)), + static_cast(ps.gridfn_data(gfns::gfn__global_z)), }; const int N_INPUT_ARRAYS = 12; const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS] = { - // $g_{ij}$ + // g_ij CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, - // $K_{ij}$ + // K_ij CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, @@ -253,7 +254,7 @@ const int N_OUTPUT_ARRAYS = 30; const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS] = { - // $g_{ij}$ $\partial_x g_{ij}$ $\partial_y g_{ij}$ $\partial_z g_{ij}$ + // g_ij partial_x g_ij partial_y g_ij partial_z g_ij CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, @@ -292,36 +293,36 @@ const CCTK_INT operation_codes[N_OUTPUT_ARRAYS] }; void* const output_arrays[N_OUTPUT_ARRAYS] = { - static_cast(ps.gridfn_data(nominal_gfns::gfn__g_dd_11)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_111)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_211)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_311)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__g_dd_12)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_112)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_212)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_312)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__g_dd_13)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_113)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_213)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_313)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__g_dd_22)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_122)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_222)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_322)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__g_dd_23)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_123)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_223)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_323)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__g_dd_33)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_133)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_233)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_333)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__K_dd_11)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__K_dd_12)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__K_dd_13)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__K_dd_22)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__K_dd_23)), - static_cast(ps.gridfn_data(nominal_gfns::gfn__K_dd_33)), + static_cast(ps.gridfn_data(gfns::gfn__g_dd_11)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_211)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_311)), + static_cast(ps.gridfn_data(gfns::gfn__g_dd_12)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_112)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_212)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_312)), + static_cast(ps.gridfn_data(gfns::gfn__g_dd_13)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_113)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_213)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_313)), + static_cast(ps.gridfn_data(gfns::gfn__g_dd_22)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_122)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_222)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_322)), + static_cast(ps.gridfn_data(gfns::gfn__g_dd_23)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_123)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_223)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_323)), + static_cast(ps.gridfn_data(gfns::gfn__g_dd_33)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_133)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_233)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_333)), + static_cast(ps.gridfn_data(gfns::gfn__K_dd_11)), + static_cast(ps.gridfn_data(gfns::gfn__K_dd_12)), + static_cast(ps.gridfn_data(gfns::gfn__K_dd_13)), + static_cast(ps.gridfn_data(gfns::gfn__K_dd_22)), + static_cast(ps.gridfn_data(gfns::gfn__K_dd_23)), + static_cast(ps.gridfn_data(gfns::gfn__K_dd_33)), }; bool first_time = true; @@ -332,9 +333,9 @@ if (first_time) // store derivative info in interpolator parameter table if (msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, - " setting up derivative info for interpolator"); + " setting up interpolator derivative info"); - status = Util_TableSetIntArray(gii.param_table_handle, + status = Util_TableSetIntArray(gi.param_table_handle, N_OUTPUT_ARRAYS, operand_indices, "operand_indices"); if (status < 0) @@ -345,7 +346,7 @@ if (first_time) , status); /*NOTREACHED*/ - status = Util_TableSetIntArray(gii.param_table_handle, + status = Util_TableSetIntArray(gi.param_table_handle, N_OUTPUT_ARRAYS, operation_codes, "operation_codes"); if (status < 0) @@ -362,7 +363,7 @@ if (msg_flag) " calling interpolator (%d points)", N_interp_points); status = CCTK_InterpLocalUniform(N_GRID_DIMS, - gii.operator_handle, gii.param_table_handle, + gi.operator_handle, gi.param_table_handle, cgi.coord_origin, cgi.coord_delta, N_interp_points, interp_coords_type_code, @@ -379,13 +380,13 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE) // look in interpolator output table entries // to see *which* point is out-of-range CCTK_INT out_of_range_pt, out_of_range_axis, out_of_range_end; - if ( (Util_TableGetInt(gii.param_table_handle, + if ( (Util_TableGetInt(gi.param_table_handle, &out_of_range_pt, "out_of_range_pt") < 0) - || (Util_TableGetInt(gii.param_table_handle, + || (Util_TableGetInt(gi.param_table_handle, &out_of_range_axis, "out_of_range_axis") < 0) - || (Util_TableGetInt(gii.param_table_handle, + || (Util_TableGetInt(gi.param_table_handle, &out_of_range_end, "out_of_range_end") < 0) ) then error_exit(ERROR_EXIT, @@ -398,12 +399,12 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE) assert(out_of_range_pt >= 0); assert(out_of_range_pt < ps.N_grid_points()); - const double global_x = ps.gridfn_data(nominal_gfns::gfn__global_x) - [out_of_range_pt]; - const double global_y = ps.gridfn_data(nominal_gfns::gfn__global_y) - [out_of_range_pt]; - const double global_z = ps.gridfn_data(nominal_gfns::gfn__global_z) - [out_of_range_pt]; + const double global_x + = ps.gridfn_data(gfns::gfn__global_x)[out_of_range_pt]; + const double global_y + = ps.gridfn_data(gfns::gfn__global_y)[out_of_range_pt]; + const double global_z + = ps.gridfn_data(gfns::gfn__global_z)[out_of_range_pt]; assert(out_of_range_axis >= 0); assert(out_of_range_axis < N_GRID_DIMS); @@ -464,8 +465,7 @@ if (msg_flag) // compute the X_ud and X_udd derivative coefficients // ... n.b. this uses the *local* (x,y,z) coordinates // - const fp r = p.ghosted_gridfn(ghosted_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 xx, yy, zz; diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn index 8e9745a..ef5576b 100644 --- a/src/gr/make.code.defn +++ b/src/gr/make.code.defn @@ -4,6 +4,7 @@ # Source files in this directory SRCS = driver.cc \ horizon_function.cc \ + Schwarzschild_EF.cc \ horizon_Jacobian.cc \ Newton.cc diff --git a/src/jtutil/miscfp.cc b/src/jtutil/miscfp.cc index 20060a2..3d19276 100644 --- a/src/jtutil/miscfp.cc +++ b/src/jtutil/miscfp.cc @@ -73,9 +73,8 @@ namespace jtutil { double arctan_xy(double x, double y) { -return ((x == 0.0) && (y == 0.0)) - ? - : std::atan2(y, x); // note reversed argument order (y,x)!! +// note reversed argument order (y,x) in std::atan2() function +return ((x == 0.0) && (y == 0.0)) ? 0.0 : std::atan2(y,x); } } // namespace jtutil:: diff --git a/src/make.configuration.defn b/src/make.configuration.defn index ecc4f59..24804b8 100644 --- a/src/make.configuration.defn +++ b/src/make.configuration.defn @@ -1,5 +1,5 @@ # additional configuration information for Cactus to build this thorn # $Header$ -##LIBDIRS += /usr/local/lib +LIBDIRS += /usr/local/lib LIBS := lapack blas $(LIBS) diff --git a/src/patch/coords.cc b/src/patch/coords.cc index a39ba2a..5b02c5b 100644 --- a/src/patch/coords.cc +++ b/src/patch/coords.cc @@ -322,6 +322,18 @@ y = r * sin(theta) * sin(phi); //************************************** +namespace local_coords + { +void r_theta_phi_of_xyz(fp x, fp y, fp z, fp& r, fp& theta, fp& phi) +{ +r = r_of_xyz(x, y, z); +theta = theta_of_xyz(x, y, z); +phi = phi_of_xy(x, y); +} + } + +//************************************** + namespace local_coords { fp theta_of_xyz(fp x, fp y, fp z) diff --git a/src/patch/coords.hh b/src/patch/coords.hh index 7ebc49e..54889da 100644 --- a/src/patch/coords.hh +++ b/src/patch/coords.hh @@ -167,6 +167,7 @@ fp partial2_phi_wrt_yy(fp x, fp y); // usual polar spherical (r,theta,phi) <--> (x,y,z) void xyz_of_r_theta_phi(fp r, fp theta, fp phi, fp& x, fp& y, fp& z); +void r_theta_phi_of_xyz(fp x, fp y, fp z, fp& r, fp& theta, fp& phi); // ... already have r_of_xyz() // ... already have phi_of_xy() fp theta_of_xyz(fp x, fp y, fp z); -- cgit v1.2.3