aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/elliptic/Jacobian.cc24
-rw-r--r--src/elliptic/Jacobian.hh9
-rw-r--r--src/gr/AHFinderDirect.hh80
-rw-r--r--src/gr/Newton.cc19
-rw-r--r--src/gr/README4
-rw-r--r--src/gr/Schwarzschild_EF.cc366
-rw-r--r--src/gr/cg.hh96
-rw-r--r--src/gr/driver.cc180
-rw-r--r--src/gr/horizon_Jacobian.cc261
-rw-r--r--src/gr/horizon_function.cc162
-rw-r--r--src/gr/make.code.defn1
-rw-r--r--src/jtutil/miscfp.cc5
-rw-r--r--src/make.configuration.defn2
-rw-r--r--src/patch/coords.cc12
-rw-r--r--src/patch/coords.hh1
15 files changed, 866 insertions, 356 deletions
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 <stdio.h>
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<fp>* H_norms_ptr,
- bool msg_flag = true);
+ const struct geometry_info& gi,
+ bool Jacobian_flag = false,
+ bool msg_flag = false,
+ jtutil::norm<fp>* 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<fp> 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<fp> 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$
+//
+// <<<prototypes for functions local to this file>>>
+// 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 <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#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 @@
//
// <<<prototypes for functions local to this file>>>
// 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<fp*>(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<fp> 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<fp> 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,
@@ -286,6 +304,31 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
}
//******************************************************************************
+
+//
+// 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$
//
// <<<prototypes for functions local to this file>>>
-// 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 <stdio.h>
#include <assert.h>
@@ -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,
@@ -56,26 +68,117 @@ void add_ghost_zone_Jacobian(const patch_system& ps,
}
//******************************************************************************
+
+//
+// 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 @@
// <<<prototypes for functions local to this file>>>
// 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 <stdio.h>
#include <assert.h>
@@ -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<fp>* 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<fp>* H_norms_ptr,
- bool msg_flag = true)
+ const struct geometry_info& gi,
+ bool Jacobian_flag = false,
+ bool msg_flag = false,
+ jtutil::norm<fp>* 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<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_x)),
- static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_y)),
- static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_z)),
+ static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_x)),
+ static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_y)),
+ static_cast<const void*>(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<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_11)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_111)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_211)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_311)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_12)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_112)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_212)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_312)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_13)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_113)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_213)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_313)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_22)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_122)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_222)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_322)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_23)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_123)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_223)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_323)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_33)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_133)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_233)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_333)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_11)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_12)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_13)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_22)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_23)),
- static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_33)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_11)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_211)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_311)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_12)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_112)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_212)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_312)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_13)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_113)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_213)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_313)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_22)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_122)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_222)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_322)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_23)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_123)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_223)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_323)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_33)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_133)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_233)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_333)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_11)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_12)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_13)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_22)),
+ static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_23)),
+ static_cast<void*>(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
@@ -324,6 +324,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)
{
return arctan_xy(z , hypot(x,y));
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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-29 19:10:57 +0000