diff options
-rw-r--r-- | src/gr/AHFinderDirect.hh | 156 | ||||
-rw-r--r-- | src/gr/Newton.cc | 180 | ||||
-rw-r--r-- | src/gr/README | 16 | ||||
-rw-r--r-- | src/gr/Schwarzschild_EF.cc | 2 | ||||
-rw-r--r-- | src/gr/driver.cc | 623 | ||||
-rw-r--r-- | src/gr/ellipsoid.c | 28 | ||||
-rw-r--r-- | src/gr/ellipsoid.maple | 37 | ||||
-rw-r--r-- | src/gr/ellipsoid.out | 38 | ||||
-rw-r--r-- | src/gr/horizon_Jacobian.cc | 422 | ||||
-rw-r--r-- | src/gr/horizon_function.cc | 77 | ||||
-rw-r--r-- | src/gr/make.code.defn | 7 | ||||
-rw-r--r-- | src/gr/makefile | 7 |
12 files changed, 57 insertions, 1536 deletions
diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh deleted file mode 100644 index 009f5f0..0000000 --- a/src/gr/AHFinderDirect.hh +++ /dev/null @@ -1,156 +0,0 @@ -// AHFinderDirect.hh -- misc global-within-this-thorn stuff -// $Header$ - -//****************************************************************************** - -// -// 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 enum holds the (a) decoded Jacobian_method parameter, -// i.e. it specifies how we compute the (a) Jacobian matrix -// -enum Jacobian_method - { - Jacobian__numerical_perturb, - Jacobian__symbolic_diff_with_FD_dr, - Jacobian__symbolic_diff // no comma - }; - -// -// this enum holds the (a) decoded geometry_method parameter, -// i.e. it specifies how we compute the (a) geometry -// -enum geometry_method - { - geometry__interpolate_from_Cactus_grid, - geometry__Schwarzschild_EF // 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 - { - enum geometry_method geometry_method; - - // parameters for geometry_method = interpolate_from_Cactus_grid - int operator_handle; // Cactus handle to interpolation op - int param_table_handle; // Cactus handle to parameter table - // for the interpolator - - // parameters for geometry_method = Schwarzschild_EF - fp geometry__Schwarzschild_EF__x_posn; // x posn of Schwarzschild BH - fp geometry__Schwarzschild_EF__y_posn; // y posn of Schwarzschild BH - fp geometry__Schwarzschild_EF__z_posn; // z posn of Schwarzschild BH - fp geometry__Schwarzschild_EF__mass; // mass of Schwarzschild BH - fp geometry__Schwarzschild_EF__epsilon; // use z axis limits if - // (x^2+y^2)/r^2 <= this - fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD - // computation of partial_k g_ij - }; - -// -// This structure holds all the information we need about the Cactus grid -// and gridfns outside the top-level driver. At present we use the -// CCTK_InterpLocalUniform() local interpolator to access the Cactus -// gridfns. Much of the information in this structure is specific to -// that API, and (alas) will need changing when we eventually switch to -// a "global" multiprocessor/distributed interpolator. -// -struct cactus_grid_info - { - cGH *GH; // --> Cactus grid hierarchy - - // Cactus coordinate system - fp coord_origin[N_GRID_DIMS]; // (x,y,z) of grid posn (0,0,0) - fp coord_delta[N_GRID_DIMS]; // (x,y,z) grid spacing - - // dimensions of gridfn data, viewed as a 3-D array - // n.b. storage ordering is Fortran, - // i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj - CCTK_INT gridfn_dims[N_GRID_DIMS]; - - // --> Cactus gridfn data for grid posn (0,0,0) - const fp* g_dd_11_data; - const fp* g_dd_12_data; - const fp* g_dd_13_data; - const fp* g_dd_22_data; - const fp* g_dd_23_data; - const fp* g_dd_33_data; - const fp* K_dd_11_data; - const fp* K_dd_12_data; - const fp* K_dd_13_data; - const fp* K_dd_22_data; - const fp* K_dd_23_data; - const fp* K_dd_33_data; - }; - -// -// This struct holds parameters for computing the Jacobian matrix. -// -struct Jacobian_info - { - enum Jacobian_method Jacobian_method; - fp perturbation_amplitude; - }; - -// -// This struct holds parameters for solving the H(h) = 0 equations. -// -struct solver_info - { - // stuff for Newton_solve() - int max_Newton_iterations; - bool output_h_and_H_at_each_Newton_iteration; - const char *h_file_name; - const char *H_of_h_file_name; - fp H_norm_for_convergence; - fp Delta_h_norm_for_convergence; - bool final_H_update_if_exit_x_H_small; - }; - -//****************************************************************************** - -// -// prototypes for functions visible outside their source files -// - -// driver.cc -extern "C" - void AHFinderDirect_driver(CCTK_ARGUMENTS); - -// horizon_function.cc -void horizon_function(patch_system& ps, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - bool Jacobian_flag = false, - bool msg_flag = false, - jtutil::norm<fp>* H_norms_ptr = NULL); - -// horizon_Jacobian.cc -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); - -// 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_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 deleted file mode 100644 index 79f1bab..0000000 --- a/src/gr/Newton.cc +++ /dev/null @@ -1,180 +0,0 @@ -// Newton.cc -- solve H(h) = 0 via Newton's method -// $Id$ -// -// <<<prototypes for functions local to this file>>> -// Newton_solve - driver to solve H(h) = 0 via Newton's method -// - -#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 { - } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function solves H(h) = 0 for h via Newton's method. -// -bool Newton_solve(patch_system& ps, - const struct cactus_grid_info& cgi, - const struct geometry_info& gi, - const struct Jacobian_info& Jacobian_info, - const struct solver_info& solver_info, - Jacobian& Jac) -{ - for (int iteration = 1 ; - iteration <= solver_info.max_Newton_iterations ; - ++iteration) - { - CCTK_VInfo(CCTK_THORNSTRING, - "Newton iteration %d", iteration); - - jtutil::norm<fp> 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()); - - if (solver_info.output_h_and_H_at_each_Newton_iteration) - then { - const char N_buffer = 100; - char file_name_buffer[100]; - - snprintf(file_name_buffer, N_buffer, - "%s.it%d", - solver_info.h_file_name, iteration); - CCTK_VInfo(CCTK_THORNSTRING, - " writing h to \"%s\"", - file_name_buffer); - ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h, - true, gfns::gfn__h, - file_name_buffer, - false); // no ghost zones - - snprintf(file_name_buffer, N_buffer, - "%s.it%d", - solver_info.H_of_h_file_name, iteration); - CCTK_VInfo(CCTK_THORNSTRING, - " writing H(h) to \"%s\"", - file_name_buffer); - ps.print_gridfn_with_xyz(gfns::gfn__H, - true, gfns::gfn__h, - file_name_buffer); - } - - if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence) - then { - CCTK_VInfo(CCTK_THORNSTRING, - "==> finished (||H|| < %g)", - double(solver_info.H_norm_for_convergence)); - return true; // *** NORMAL RETURN *** - } - - // compute the Newton step - horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac); - CCTK_VInfo(CCTK_THORNSTRING, - "solving linear system for Delta_h (%d unknowns)", - Jac.NN()); - const fp rcond = Jac.solve_linear_system(gfns::gfn__H, - gfns::gfn__Delta_h); - if (rcond == 0.0) - then { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Newton_solve: Jacobian matrix is numerically singular!"); - return false; // *** ERROR RETURN *** - } - if (rcond > 0.0) - then CCTK_VInfo(CCTK_THORNSTRING, - "done solving linear system (condition number %.2e)", - double(rcond)); - else CCTK_VInfo(CCTK_THORNSTRING, - "done solving linear system"); - - // take the Newton step - jtutil::norm<fp> Delta_h_norms; - 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 Delta_h = p.gridfn(gfns::gfn__Delta_h, irho,isigma); - Delta_h_norms.data(Delta_h); - - p.ghosted_gridfn(gfns::gfn__h, irho,isigma) -= Delta_h; - } - } - } - CCTK_VInfo(CCTK_THORNSTRING, - "h += Delta_h (rms-norm=%.2e, infinity-norm=%.2e)", - Delta_h_norms.rms_norm(), Delta_h_norms.infinity_norm()); - - if (Delta_h_norms.infinity_norm() <= solver_info - .Delta_h_norm_for_convergence) - then { - CCTK_VInfo(CCTK_THORNSTRING, - "==> finished (||Delta_h|| < %g)", - double(solver_info.Delta_h_norm_for_convergence)); - if (solver_info.final_H_update_if_exit_x_H_small) - then { - CCTK_VInfo(CCTK_THORNSTRING, - " doing final H(h) evaluation"); - horizon_function(ps, cgi, gi, false, true); - } - return true; // *** NORMAL RETURN *** - } - } - - -CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Newton_solve: no convergence in %d iterations!", - solver_info.max_Newton_iterations); -if (solver_info.final_H_update_if_exit_x_H_small) - then { - CCTK_VInfo(CCTK_THORNSTRING, - " doing final H(h) evaluation"); - horizon_function(ps, cgi, gi, false, true); - } -return false; // *** ERROR RETURN *** -} diff --git a/src/gr/README b/src/gr/README index 32b834f..fdd3450 100644 --- a/src/gr/README +++ b/src/gr/README @@ -3,16 +3,10 @@ general relativity. The main source files in this directory are as follows: -driver.cc - Top-level driver for the whole thorn. - horizon_function.cc Computes the LHS of the AH equation, $H = H(h)$. -horizon_Jacobian.cc - Computes the Jacobian of the $H(h)$ function, $J[H(h)]$. - -AHFinderDirect.hh +gr.hh Overall header file for all the external routines in this directory. gfns.hh @@ -28,6 +22,9 @@ cg.hh doit.maple Top-level driver for all the Maple code. +maple.log + This is a log file (transcript) of the Maple computation. + setup_gr_gfas.maple gr_gfas.minc These files define the relativity gridfns for the Maple code, @@ -39,5 +36,6 @@ horizon.maple These files compute the $H(h)$ function and the $J[H(h)]$ coefficients in Maple. -maple.log - This is a log file (transcript) of the Maple computation. +makefile + This drives the Maple code-generation process; the resulting + (Maple-generated) C code lives in ../gr.cg/ diff --git a/src/gr/Schwarzschild_EF.cc b/src/gr/Schwarzschild_EF.cc index 1bce005..628e641 100644 --- a/src/gr/Schwarzschild_EF.cc +++ b/src/gr/Schwarzschild_EF.cc @@ -39,7 +39,7 @@ using jtutil::error_exit; #include "../elliptic/Jacobian.hh" #include "gfns.hh" -#include "AHFinderDirect.hh" +#include "gr.hh" //****************************************************************************** diff --git a/src/gr/driver.cc b/src/gr/driver.cc deleted file mode 100644 index 5da7b3f..0000000 --- a/src/gr/driver.cc +++ /dev/null @@ -1,623 +0,0 @@ -// driver.cc -- top level driver for finding apparent horizons -// $Id$ -// -// <<<prototypes for functions local to this file>>> -// AHFinderDirect_driver - top-level driver -/// decode_Jacobian_method - decode the Jacobian_method parameter -/// decode_geometry_method - decode the geometry_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 -/// -/// print_Jacobians - print a pair of Jacobians -/// - -#include <stdio.h> -#include <assert.h> -#include <math.h> -#include <vector> - -#include "util_Table.h" -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.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 { -enum Jacobian_method - decode_Jacobian_method(const char Jacobian_method_string[]); -enum geometry_method - decode_geometry_method(const char geometry_method_string[]); - -void setup_Kerr_horizon(patch_system& ps, - fp x_posn, fp y_posn, fp z_posn, - fp m, fp a, - bool Kerr_Schild_flag); -void setup_ellipsoid(patch_system& ps, - fp x_center, fp y_center, fp z_center, - fp x_radius, fp y_radius, fp z_radius); - -void print_Jacobians(const patch_system& ps, - const Jacobian* Jac_NP, - bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr, - const char file_name[]); - }; - -//****************************************************************************** - -// -// This function is the Cactus interface for the test driver. -// -extern "C" - void AHFinderDirect_driver(CCTK_ARGUMENTS) -{ -DECLARE_CCTK_ARGUMENTS -DECLARE_CCTK_PARAMETERS - -CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures"); - - -// -// set up the geometry parameters and the geometry interpolator -// -struct geometry_info gi; -gi.geometry_method = decode_geometry_method(geometry_method); - -// parameters for geometry_method = "interpolate from Cactus grid" -CCTK_VInfo(CCTK_THORNSTRING, " setting up geometry interpolator"); -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*/ -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*/ - -// parameters for geometry_method = "Schwarzschild/EF" -gi.geometry__Schwarzschild_EF__mass = geometry__Schwarzschild_EF__mass; -gi.geometry__Schwarzschild_EF__x_posn = geometry__Schwarzschild_EF__x_posn; -gi.geometry__Schwarzschild_EF__y_posn = geometry__Schwarzschild_EF__y_posn; -gi.geometry__Schwarzschild_EF__z_posn = geometry__Schwarzschild_EF__z_posn; -gi.geometry__Schwarzschild_EF__epsilon = geometry__Schwarzschild_EF__epsilon; -gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz; - - -// -// set up the interpatch interpolator -// -CCTK_VInfo(CCTK_THORNSTRING, " setting up interpatch interpolator"); -const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name); -if (interp_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "couldn't find interpolator \"%s\"!", - interpatch_interpolator_name); /*NOTREACHED*/ -const int interp_param_table_handle - = Util_TableCreateFromString(interpatch_interpolator_pars); -if (interp_param_table_handle < 0) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "bad interpatch-interpolator parameter(s) \"%s\"!", - interpatch_interpolator_pars); /*NOTREACHED*/ - - -// -// set up the Cactus grid info -// -CCTK_VInfo(CCTK_THORNSTRING, " setting up Cactus grid info"); -struct cactus_grid_info cgi; -cgi.GH = cctkGH; -cgi.coord_origin[0] = cctk_origin_space[0]; -cgi.coord_origin[1] = cctk_origin_space[1]; -cgi.coord_origin[2] = cctk_origin_space[2]; -cgi.coord_delta[0] = cctk_delta_space[0]; -cgi.coord_delta[1] = cctk_delta_space[1]; -cgi.coord_delta[2] = cctk_delta_space[2]; -cgi.gridfn_dims[0] = cctk_lsh[0]; -cgi.gridfn_dims[1] = cctk_lsh[1]; -cgi.gridfn_dims[2] = cctk_lsh[2]; -// n.b. The cgi.[gK]_dd_??_data are actually const fp * pointers, -// since we won't modify the 3-D gridfn data! But static_cast<...> -// won't change const modifiers, so we just cast to fp* and let -// the assignment take care of the const part... -cgi.g_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxx")); -cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxy")); -cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxz")); -cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyy")); -cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyz")); -cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gzz")); -cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxx")); -cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxy")); -cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxz")); -cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyy")); -cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyz")); -cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz")); - - -// -// create the patch system and initialize the xyz derivative coefficients -// -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, - gfns::nominal_min_gfn, gfns::nominal_max_gfn, - gfns::ghosted_min_gfn, gfns::ghosted_max_gfn, - interp_handle, interp_param_table_handle); - - -// -// set up the initial guess for the apparent horizon shape -// -if (STRING_EQUAL(initial_guess_method, "read from file")) - then { - CCTK_VInfo(CCTK_THORNSTRING, - "reading initial guess h from \"%s\"", - h_file_name); - ps.read_ghosted_gridfn(gfns::gfn__h, - h_file_name, - false); // no ghost zones - } - else { - if (STRING_EQUAL(initial_guess_method, "ellipsoid")) - then setup_ellipsoid(ps, - initial_guess__ellipsoid__x_center, - initial_guess__ellipsoid__y_center, - initial_guess__ellipsoid__z_center, - initial_guess__ellipsoid__x_radius, - initial_guess__ellipsoid__y_radius, - initial_guess__ellipsoid__z_radius); - else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr")) - then setup_Kerr_horizon(ps, - 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__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\"!", - initial_guess_method); /*NOTREACHED*/ - - // write initial guess back to this file - CCTK_VInfo(CCTK_THORNSTRING, - "writing initial guess h to \"%s\"", - h_file_name); - ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h, - true, gfns::gfn__h, - h_file_name, - false); // no ghost zones - } - - -// -// decode/copy parameters into structures -// - -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.output_h_and_H_at_each_Newton_iteration - = (output_h_and_H_at_each_Newton_iteration != 0); -solver_info.h_file_name = h_file_name; -solver_info.H_of_h_file_name = H_of_h_file_name; -solver_info.H_norm_for_convergence = H_norm_for_convergence; -solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence; -solver_info.final_H_update_if_exit_x_H_small - = (final_H_update_if_exit_x_H_small != 0); - -// -// 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, gi, false, true, &H_norms); - CCTK_TimerStopI(timer_handle); - printf("timer stats for evaluating H(h):\n"); - CCTK_TimerPrintDataI(timer_handle, -1); - - CCTK_VInfo(CCTK_THORNSTRING, - " H(h) rms-norm %.2e, infinity-norm %.2e", - H_norms.rms_norm(), H_norms.infinity_norm()); - CCTK_VInfo(CCTK_THORNSTRING, - "writing H(h) to \"%s\"", - H_of_h_file_name); - ps.print_gridfn_with_xyz(gfns::gfn__H, - true, gfns::gfn__h, - H_of_h_file_name); - } -else if ( STRING_EQUAL(method, "Jacobian test") - || STRING_EQUAL(method, "Jacobian test (NP only)") ) - then { - const bool NP_only = STRING_EQUAL(method, "Jacobian test (NP only)"); - const bool do_non_NP_Jac = ! NP_only; - - // numerical perturbation - Jacobian* const Jac_NP = & new_Jacobian(ps, Jacobian_storage_method); - horizon_function(ps, cgi, gi, true); - Jacobian_info.Jacobian_method = Jacobian__numerical_perturb; - horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_NP); - - // symbolic differentiation with finite diff d/dr - Jacobian* Jac_SD_FDdr = NULL; - if (do_non_NP_Jac) - then { - Jac_SD_FDdr = & new_Jacobian(ps, Jacobian_storage_method); - horizon_function(ps, cgi, gi, true); - Jacobian_info.Jacobian_method - = Jacobian__symbolic_diff_with_FD_dr; - horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_SD_FDdr); - } - - CCTK_VInfo(CCTK_THORNSTRING, - "writing Jacobian%s to \"%s\"", - (NP_only ? "" : "s"), Jacobian_file_name); - print_Jacobians(ps, - Jac_NP, - do_non_NP_Jac, Jac_SD_FDdr, - Jacobian_file_name); - } -else if (STRING_EQUAL(method, "Newton solve")) - then { - Jacobian& Jac = new_Jacobian(ps, Jacobian_storage_method); - - const int timer_handle = CCTK_TimerCreateI(); - CCTK_TimerStartI(timer_handle); - 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); - - CCTK_VInfo(CCTK_THORNSTRING, - "writing final horizon shape h to \"%s\"", - h_file_name); - 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(gfns::gfn__H, - true, gfns::gfn__h, - H_of_h_file_name); - } -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "unknown method=\"%s\"!", - 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 Jacobian__numerical_perturb; -else if (STRING_EQUAL(Jacobian_method_string, - "symbolic differentiation with finite diff d/dr")) - then return Jacobian__symbolic_diff_with_FD_dr; -else if (STRING_EQUAL(Jacobian_method_string, - "symbolic differentiation")) - then return Jacobian__symbolic_diff; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "unknown Jacobian_method_string=\"%s\"!", - Jacobian_method_string); /*NOTREACHED*/ -} - } - -//****************************************************************************** - -// -// This function decodes the geometry_method parameter (string) into -// an internal enum for future use. -// -namespace { -enum geometry_method - decode_geometry_method(const char geometry_method_string[]) -{ -if (STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid")) - then return geometry__interpolate_from_Cactus_grid; -else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF")) - then return geometry__Schwarzschild_EF; -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "unknown geometry_method_string=\"%s\"!", - geometry_method_string); /*NOTREACHED*/ -} - } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function sets up the horizon of a Kerr black hole in Kerr or -// Kerr-Schild coordinates, on the nominal grid, in the h gridfn. -// -// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903, -// and the horizon is worked out on page 13.2 of my AHFinderDirect notes. -// -// Arguments: -// [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, -// true to use Kerr-Schild coordinates -// -namespace { -void setup_Kerr_horizon(patch_system& ps, - fp x_posn, fp y_posn, fp z_posn, - fp m, fp a, - bool Kerr_Schild_flag) -{ -const char* const name = Kerr_Schild_flag ? "Kerr-Schild" : "Kerr"; - -CCTK_VInfo(CCTK_THORNSTRING, - "setting up horizon for Kerr in %s coords", - name); -CCTK_VInfo(CCTK_THORNSTRING, - " 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)); - -// horizon in Kerr-Schild coordinates is coordinate ellipsoid -const fp z_radius = r; -const fp xy_radius = Kerr_Schild_flag ? r * sqrt(1.0 + a*a*m*m/(r*r)) : r; - -CCTK_VInfo(CCTK_THORNSTRING, - " horizon is coordinate %s", - Kerr_Schild_flag ? "ellipsoid" : "sphere"); -setup_ellipsoid(ps, - x_posn, y_posn, z_posn, - xy_radius, xy_radius, z_radius); -} - } - -//****************************************************************************** - -// -// This function sets up an ellipsoid in the gridfn h, using the -// formulas in "ellipsoid.maple" and the Maple-generated C code in -// "ellipsoid.c": -// -// ellipsoid has center (A,B,C), radius (a,b,c) -// angular coordinate system has center (U,V,W) -// -// direction cosines wrt angular coordinate center are (xcos,ycos,zcos) -// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) -// -// then the equation of the ellipsoid is -// (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 -// ----------------- + ---------------- + ----------------- = 1 -// a^2 b^2 c^2 -// -// to solve this, we introduce intermediate variables -// AU = A - U -// BV = B - V -// CW = C - W -// -namespace { -void setup_ellipsoid(patch_system& ps, - fp x_center, fp y_center, fp z_center, - fp x_radius, fp y_radius, fp z_radius) -{ -CCTK_VInfo(CCTK_THORNSTRING, - "setting h = ellipsoid: center=(%g,%g,%g)", - double(x_center), double(y_center), double(z_center)); -CCTK_VInfo(CCTK_THORNSTRING, - " radius=(%g,%g,%g)", - double(x_radius), double(y_radius), double(z_radius)); - - 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 rho = p.rho_of_irho(irho); - const fp sigma = p.sigma_of_isigma(isigma); - fp xcos, ycos, zcos; - p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos); - - // set up variables used by Maple-generated code - const fp AU = x_center - ps.origin_x(); - const fp BV = y_center - ps.origin_y(); - const fp CW = z_center - ps.origin_z(); - const fp a = x_radius; - const fp b = y_radius; - const fp c = z_radius; - - // compute the solutions r_plus and r_minus - fp r_plus, r_minus; - #include "ellipsoid.c" - - // exactly one of the solutions (call it r) should be positive - fp r; - if ((r_plus > 0.0) && (r_minus < 0.0)) - then r = r_plus; - else if ((r_plus < 0.0) && (r_minus > 0.0)) - then r = r_minus; - else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "\n" -" expected exactly one r>0 solution to quadratic, got 0 or 2!\n" -" %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n" -" direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n" -" ==> r_plus=%g r_minus=%g\n" - , - p.name(), irho, isigma, - double(rho), double(sigma), - double(xcos), double(ycos), double(zcos), - double(r_plus), double(r_minus)); - /*NOTREACHED*/ - - // r = horizon radius at this grid point - p.ghosted_gridfn(gfns::gfn__h, irho,isigma) = r; - } - } - } -} - } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function prints one or two Jacobian matrices (and their difference -// in the latter case) to a named output file. -// -// Arguments: -// do_non_NP_Jac = true if we should also print *Jac_SD_FDdr -// false to skip this (and print only *Jac_NP) -// -namespace { -void print_Jacobians(const patch_system& ps, - const Jacobian* Jac_NP, - bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr, - const char file_name[]) -{ -FILE *fileptr = fopen(file_name, "w"); -if (fileptr == NULL) - then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "print_Jacobians(): can't open output file \"%s\"!", - file_name); /*NOTREACHED*/ - -fprintf(fileptr, "# column 1 = x patch name\n"); -fprintf(fileptr, "# column 2 = x patch number\n"); -fprintf(fileptr, "# column 3 = x irho\n"); -fprintf(fileptr, "# column 4 = x isigma\n"); -fprintf(fileptr, "# column 5 = x II\n"); -fprintf(fileptr, "# column 6 = y patch name\n"); -fprintf(fileptr, "# column 7 = y patch number\n"); -fprintf(fileptr, "# column 8 = y irho\n"); -fprintf(fileptr, "# column 9 = y isigma\n"); -fprintf(fileptr, "# column 10 = y JJ\n"); -fprintf(fileptr, "# column 11 = Jac_NP(II,JJ)\n"); -if (do_non_NP_Jac) - then { - fprintf(fileptr, "# column 12 = Jac_SD_FDdr(II,JJ)\n"); - fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n"); - fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n"); - } - - 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); - - for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) - { - patch& yp = ps.ith_patch(ypn); - - 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); - - if (! Jac_NP->is_explicitly_stored(II,JJ)) - then continue; // skip sparse points - - const fp NP = (*Jac_NP)(II,JJ); - const fp SD_FDdr = do_non_NP_Jac ? (*Jac_SD_FDdr)(II,JJ) : 0.0; - - if ((NP == 0.0) && (SD_FDdr == 0.0)) - then continue; // skip zero values (if == ) - - fprintf(fileptr, - "%s %d %d %d %d\t%s %d %d %d %d\t%.10g", - xp.name(), xpn, x_irho, x_isigma, II, - yp.name(), ypn, y_irho, y_isigma, JJ, - double(NP)); - - if (do_non_NP_Jac) - then { - 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, - "\t%.10g\t%e\t%e", - double(SD_FDdr), - double(SD_FDdr_abs_error), double(SD_FDdr_rel_error)); - } - - fprintf(fileptr, "\n"); - } - } - } - } - } - } - -fclose(fileptr); -} - } diff --git a/src/gr/ellipsoid.c b/src/gr/ellipsoid.c deleted file mode 100644 index 49a6457..0000000 --- a/src/gr/ellipsoid.c +++ /dev/null @@ -1,28 +0,0 @@ -fp t1, t2, t3, t5, t6, t7, t9, t10, t12, t28; -fp t30, t33, t35, t36, t40, t42, t43, t48, t49, t52; -fp t55; - t1 = a*a; - t2 = b*b; - t3 = t1*t2; - t5 = t3*zcos*CW; - t6 = c*c; - t7 = t1*t6; - t9 = t7*ycos*BV; - t10 = t2*t6; - t12 = t10*xcos*AU; - t28 = xcos*xcos; - t30 = CW*CW; - t33 = BV*BV; - t35 = t10*t28; - t36 = ycos*ycos; - t40 = AU*AU; - t42 = t7*t36; - t43 = zcos*zcos; - t48 = t3*t43; - t49 = -2.0*t1*zcos*CW*ycos*BV-2.0*t2*zcos*CW*xcos*AU-2.0*t6*ycos*BV*xcos* -AU+t2*t28*t30+t6*t28*t33-t35+t1*t36*t30+t6*t36*t40-t42+t1*t43*t33+t2*t43*t40- -t48; - t52 = sqrt(-t3*t6*t49); - t55 = 1/(t35+t42+t48); - r_plus = (t5+t9+t12+t52)*t55; - r_minus = (t5+t9+t12-t52)*t55; diff --git a/src/gr/ellipsoid.maple b/src/gr/ellipsoid.maple deleted file mode 100644 index 349cc8d..0000000 --- a/src/gr/ellipsoid.maple +++ /dev/null @@ -1,37 +0,0 @@ -# ellipsoid.maple -- compute equations for offset ellipsoid setup -# $Header$ - -# -# This program finds the intersection of an ellipsoid with an offset ray. -# -# ellipsoid has center (A,B,C), radius (a,b,c) -# angular coordinate system has center (U,V,W) -# -# direction cosines wrt angular coordinate center are (alpha,beta,gamma) -# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos) -# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) -# -# then the equation of the ellipsoid is -# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 -# ----------------- + ---------------- + ----------------- = 1 -# a^2 b^2 c^2 -# -# to solve this, we introduce intermediate variables -# AU = A - U -# BV = B - V -# CW = C - W -# -eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1; - -read "../maple/util.mm"; -read "../maple/codegen2.mm"; - -[solve(eqn, r)]; -map(simplify, %); -[r_plus = %[1], r_minus = %[2]]; -solnlist := [codegen[optimize](%)]; - -ftruncate("ellipsoid.c"); -print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]), - "fp", "ellipsoid.c"); -codegen[C](solnlist, filename="ellipsoid.c"); diff --git a/src/gr/ellipsoid.out b/src/gr/ellipsoid.out deleted file mode 100644 index 4a0ee01..0000000 --- a/src/gr/ellipsoid.out +++ /dev/null @@ -1,38 +0,0 @@ - |\^/| Maple 7 (IBM INTEL LINUX) -._|\| |/|_. Copyright (c) 2001 by Waterloo Maple Inc. - \ MAPLE / All rights reserved. Maple is a registered trademark of - <____ ____> Waterloo Maple Inc. - | Type ? for help. -# ellipsoid.maple -- compute equations for offset ellipsoid setup -# $Header: /numrelcvs/ThornburgCVS/AHFinderDirect/src/gr/ellipsoid.maple,v 1.2 2002/04/22 15:10:32 jthorn Exp $ -> -# -# ellipsoid has center (A,B,C), radius (a,b,c) -# angular coordinate system has center (U,V,W) -# -# direction cosines wrt angular coordinate center are (alpha,beta,gamma) -# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos) -# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) -# -# then the equation of the ellipsoid is -# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 -# ----------------- + ---------------- + ----------------- = 1 -# a^2 b^2 c^2 -# -# to solve this, we introduce intermediate variables -# AU = A - U -# BV = B - V -# CW = C - W -# -> eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1; - 2 2 2 - (xcos r - AU) (ycos r - BV) (zcos r - CW) - eqn := -------------- + -------------- + -------------- = 1 - 2 2 2 - a b c - -> -> read "../maple/util.mm"; -Error, unable to read `../maple/util.mm` -> quit -bytes used=129844, alloc=196572, time=0.03 diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc deleted file mode 100644 index 24a59a7..0000000 --- a/src/gr/horizon_Jacobian.cc +++ /dev/null @@ -1,422 +0,0 @@ -// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h) -// $Id$ -// -// <<<prototypes for functions local to this file>>> -// -// 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 -/// - -#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 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, - int xm_irho, int xm_isigma, - 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 Jacobian__numerical_perturb: - horizon_Jacobian_NP(ps, cgi, gi, Jacobian_info, Jac); - break; -case Jacobian__symbolic_diff_with_FD_dr: - horizon_Jacobian_SD(ps, cgi, gi, Jacobian_info, Jac); - break; -case Jacobian__symbolic_diff: - // 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 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 -// -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) -{ -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); -} - } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function computes the Jacobian matrix of the LHS function H(h) -// by symbolic differentiation from the Jacobian coefficient (angular) -// gridfns. The computation is done a Jacobian row at a time, using -// equation (25) of my 1996 apparent horizon finding paper, except that -// the d/dr term is done by numerical perturbation. -// -// The overall algorithm is -// save_H = H -// h += perturbation_amplitude -// evaluate H(h) (silently) -// Jac = diagonal[ (H(h) - save_H)/perturbation_amplitude ] -// h -= perturbation_amplitude -// H = save_H -// Jac += Jacobian coeffs * molecule coeffs -// -// Inputs (angular gridfns, on ghosted grid): -// h # shape of trial surface -// H # H(h) assumed to already be computed -// partial_H_wrt_partial_d_h # Jacobian coefficients -// partial_H_wrt_partial_dd_h # (also assumed to already be computed) -// -// Outputs: -// The Jacobian matrix is stored in the Jacobian object 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 (symbolic differentiation)"); - -ps.compute_synchronize_Jacobian(); -Jac.zero_matrix(); - - -// -// compute d/dr term in Jacobian by numerical perturbation -// - -CCTK_VInfo(CCTK_THORNSTRING, - " computing d/dr terms by numerical perturbation"); -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); - for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) - { - for (int isigma = p.min_isigma() ; - isigma <= p.max_isigma() ; - ++isigma) - { - const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma); - const fp old_H = p.gridfn(gfns::gfn__save_H, irho,isigma); - const fp new_H = p.gridfn(gfns::gfn__H, irho,isigma); - Jac(II,II) += (new_H - old_H) / epsilon; - } - } - } -ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h); -ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H); - - -// -// now the main symbolic-differentiation Jacobian computation -// - -CCTK_VInfo(CCTK_THORNSTRING, - " computing main terms by symbolic differentiation"); - 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) - { - // - // compute the main Jacobian terms for this grid point, i.e. - // partial H(this point x, Jacobian row II) - // --------------------------------------------- - // partial h(other points y, Jacobian column JJ) - // - - // Jacobian row index - const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma); - - // Jacobian coefficients for this point - const fp Jacobian_coeff_rho - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1, - x_irho, x_isigma); - const fp Jacobian_coeff_sigma - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2, - x_irho, x_isigma); - const fp Jacobian_coeff_rho_rho - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11, - x_irho, x_isigma); - const fp Jacobian_coeff_rho_sigma - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12, - x_irho, x_isigma); - const fp Jacobian_coeff_sigma_sigma - = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22, - x_irho, x_isigma); - - // partial_rho, partial_rho_rho - for (int m_irho = xp.molecule_min_m() ; - m_irho <= xp.molecule_max_m() ; - ++m_irho) - { - const int xm_irho = x_irho + m_irho; - const fp Jac_rho = Jacobian_coeff_rho - * xp.partial_rho_coeff(m_irho); - const fp Jac_rho_rho = Jacobian_coeff_rho_rho - * xp.partial_rho_rho_coeff(m_irho); - const fp Jac_sum = Jac_rho + Jac_rho_rho; - if (xp.is_in_nominal_grid(xm_irho, x_isigma)) - then Jac(II, xp,xm_irho,x_isigma) += Jac_sum; - else add_ghost_zone_Jacobian - (ps, xp, xp.minmax_rho_ghost_zone(m_irho < 0), - II, xm_irho, x_isigma, - Jac_sum, Jac); - } - - // partial_sigma, partial_sigma_sigma - for (int m_isigma = xp.molecule_min_m() ; - m_isigma <= xp.molecule_max_m() ; - ++m_isigma) - { - const int xm_isigma = x_isigma + m_isigma; - const fp Jac_sigma = Jacobian_coeff_sigma - * xp.partial_sigma_coeff(m_isigma); - const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma - * xp.partial_sigma_sigma_coeff(m_isigma); - const fp Jac_sum = Jac_sigma + Jac_sigma_sigma; - if (xp.is_in_nominal_grid(x_irho, xm_isigma)) - then Jac(II, xp,x_irho,xm_isigma) += Jac_sum; - else add_ghost_zone_Jacobian - (ps, xp, xp.minmax_sigma_ghost_zone(m_isigma < 0), - II, x_irho, xm_isigma, - Jac_sum, Jac); - } - - // partial_rho_sigma - for (int m_irho = xp.molecule_min_m() ; - m_irho <= xp.molecule_max_m() ; - ++m_irho) - { - for (int m_isigma = xp.molecule_min_m() ; - m_isigma <= xp.molecule_max_m() ; - ++m_isigma) - { - const int xm_irho = x_irho + m_irho; - const int xm_isigma = x_isigma + m_isigma; - const fp Jac_rho_sigma - = Jacobian_coeff_rho_sigma - * xp.partial_rho_sigma_coeff(m_irho, m_isigma); - if (xp.is_in_nominal_grid(xm_irho, xm_isigma)) - then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma; - else add_ghost_zone_Jacobian - (ps, xp, - xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0, - xm_irho, xm_isigma), - II, xm_irho, xm_isigma, - Jac_rho_sigma, Jac); - } - } - - } - } - } -} - } - -//****************************************************************************** - -// -// This function adds the ghost-zone Jacobian dependency contributions -// for a single ghost-zone point, to a Jacobian matrix. -// -// Arguments: -// ps = The patch system. -// xp = The patch containing the center point of the molecule. -// xmgz = If the x+m point is in a ghost zone, this must be that ghost zone. -// If the x+m point is not in a ghost zone, this argument is ignored. -// x_II = The Jacobian row of the x point. -// xm_(irho,isigma) = The coordinates (in xp) of the x+m point of the molecule. -// mol = The molecule coefficient. -// Jac = The Jacobian matrix. -// -namespace { -void add_ghost_zone_Jacobian(const patch_system& ps, - const patch& xp, const ghost_zone& xmgz, - int x_II, - int xm_irho, int xm_isigma, - fp mol, Jacobian& Jac) -{ -const patch_edge& xme = xmgz.my_edge(); -const int xm_iperp = xme.iperp_of_irho_isigma(xm_irho, xm_isigma); -const int xm_ipar = xme. ipar_of_irho_isigma(xm_irho, xm_isigma); - -// FIXME: this won't change from one call to another -// ==> it would be more efficient to reuse the same buffer -// across multiple calls on this function -int global_min_ym, global_max_ym; -ps.synchronize_Jacobian_global_minmax_ym(global_min_ym, global_max_ym); -jtutil::array1d<fp> Jacobian_buffer(global_min_ym, global_max_ym); - -// on what other points y does this molecule point xm depend -// via the patch_system::synchronize() operation? -int y_iperp; -int y_posn, min_ym, max_ym; -const patch_edge& ye = ps.synchronize_Jacobian(xmgz, - xm_iperp, xm_ipar, - y_iperp, - y_posn, min_ym, max_ym, - Jacobian_buffer); -const patch& yp = ye.my_patch(); - -// add the Jacobian contributions from the ym points - for (int ym = min_ym ; ym <= max_ym ; ++ym) - { - const int y_ipar = y_posn + ym; - const int y_irho = ye. irho_of_iperp_ipar(y_iperp,y_ipar); - const int y_isigma = ye.isigma_of_iperp_ipar(y_iperp,y_ipar); - const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma); - Jac(x_II,JJ) += mol*Jacobian_buffer(ym); - } -} - } diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 41eaf23..3749ea0 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -39,7 +39,7 @@ using jtutil::pow4; #include "../elliptic/Jacobian.hh" #include "gfns.hh" -#include "AHFinderDirect.hh" +#include "gr.hh" //****************************************************************************** @@ -48,19 +48,19 @@ using jtutil::pow4; // namespace { -void setup_xyz_posns(patch_system& ps, bool msg_flag); -void interpolate_geometry(patch_system& ps, +void setup_xyz_posns(patch_system& ps, bool print_msg_flag); +bool interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, - bool msg_flag); + bool print_msg_flag); void Schwarzschild_EF_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, - bool msg_flag); + bool print_msg_flag); void compute_H(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* H_norms_ptr, - bool msg_flag); + bool print_msg_flag); } //****************************************************************************** @@ -100,37 +100,43 @@ 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. +// print_msg_flag = true to print status messages, +// false to skip this. // H_norms_ptr = (out) If this pointer is non-NULL, the norm object it // points to is updated with all the H values in the // grid. This norm object can then be used to compute // various (gridwise) norms of H. // -void horizon_function(patch_system& ps, +// Results: +// This function returns an error flag: normally this is false, but if +// the geometry interpolation would need data outside the Cactus grid, +// or data from an excised region, then this function returns true. +// +bool horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, bool Jacobian_flag = false, - bool msg_flag = false, + bool print_msg_flag = false, jtutil::norm<fp>* H_norms_ptr = NULL) { -if (msg_flag) +if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " horizon function"); // fill in values of all ghosted gridfns in ghost zones ps.synchronize(); // set up xyz positions of grid points -setup_xyz_posns(ps, msg_flag); +setup_xyz_posns(ps, print_msg_flag); // compute the "geometry" g_ij, K_ij, and partial_k g_ij switch (gi.geometry_method) { case geometry__interpolate_from_Cactus_grid: - interpolate_geometry(ps, cgi, gi, msg_flag); + if (interpolate_geometry(ps, cgi, gi, print_msg_flag)) + then return true; // *** ERROR RETURN *** break; case geometry__Schwarzschild_EF: - Schwarzschild_EF_geometry(ps, gi, msg_flag); + Schwarzschild_EF_geometry(ps, gi, print_msg_flag); break; default: error_exit(PANIC_EXIT, @@ -142,7 +148,9 @@ default: // compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing -compute_H(ps, Jacobian_flag, H_norms_ptr, msg_flag); +compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag); + +return false; // *** NORMAL RETURN *** } //****************************************************************************** @@ -152,9 +160,9 @@ compute_H(ps, Jacobian_flag, H_norms_ptr, msg_flag); // in the gridfns global_[xyz]. These will be used by interplate_geometry(). // namespace { -void setup_xyz_posns(patch_system& ps, bool msg_flag) +void setup_xyz_posns(patch_system& ps, bool print_msg_flag) { -if (msg_flag) +if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " xyz positions and derivative coefficients"); @@ -213,13 +221,18 @@ if (msg_flag) // On the first call this function also modifies the interpolator // parameter table. // +// Results: +// This function returns an error flag: normally this is false, but if +// the interpolation would need data outside the Cactus grid, or data +// from an excised region, then this function returns true. +// namespace { -void interpolate_geometry(patch_system& ps, +bool interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_info& gi, - bool msg_flag) + bool print_msg_flag) { -if (msg_flag) +if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " interpolating g_ij and Kij from Cactus 3-D grid"); @@ -343,7 +356,7 @@ if (first_time) first_time = false; // store derivative info in interpolator parameter table - if (msg_flag) + if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " setting up interpolator derivative info"); @@ -370,7 +383,7 @@ if (first_time) status); /*NOTREACHED*/ } -if (msg_flag) +if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " calling interpolator (%d points)", N_interp_points); @@ -425,13 +438,13 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE) assert((out_of_range_end == -1) || (out_of_range_end == +1)); const char end = (out_of_range_end == -1) ? '-' : '+'; - error_exit(ERROR_EXIT, -"***** interpolate_geometry():\n" -" the point (%g,%g,%g) on the trial horizon surface\n" -" is outside the grid in the %c%c direction!\n" -, - global_x, global_y, global_z, - end, axis); /*NOTREACHED*/ + CCTK_VInfo(CCTK_THORNSTRING, + "*** the trial-horizon-surface point (%g,%g,%g)", + global_x, global_y, global_z); + CCTK_VInfo(CCTK_THORNSTRING, + "*** is outside the grid in the %c%c direction!", + end, axis); + return true; // *** ERROR RETURN *** } if (status < 0) then error_exit(ERROR_EXIT, @@ -439,6 +452,8 @@ if (status < 0) " CCTK_InterpLocalUniform() status=%d\n" , status); /*NOTREACHED*/ + +return false; // *** NORMAL RETURN *** } } @@ -458,9 +473,9 @@ namespace { void compute_H(patch_system& ps, bool Jacobian_flag, jtutil::norm<fp>* H_norms_ptr, - bool msg_flag) + bool print_msg_flag) { -if (msg_flag) +if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); for (int pn = 0 ; pn < ps.N_patches() ; ++pn) diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn index ef5576b..0fcbabb 100644 --- a/src/gr/make.code.defn +++ b/src/gr/make.code.defn @@ -2,11 +2,8 @@ # $Header$ # Source files in this directory -SRCS = driver.cc \ - horizon_function.cc \ - Schwarzschild_EF.cc \ - horizon_Jacobian.cc \ - Newton.cc +SRCS = horizon_function.cc \ + Schwarzschild_EF.cc # Subdirectories containing source files SUBDIRS = diff --git a/src/gr/makefile b/src/gr/makefile index 0b980fb..d3618f9 100644 --- a/src/gr/makefile +++ b/src/gr/makefile @@ -1,5 +1,5 @@ # Makefile for GR code -# $Id: makefile,v 1.6 2002-04-22 11:02:27 jthorn Exp $ +# $Id: makefile,v 1.7 2002-09-03 15:14:03 jthorn Exp $ # # Environment Variables: @@ -11,7 +11,6 @@ # cg ==> run all the *.mm files to generate C code # mm ==> preprocess all *.maple files to produce *.mm files # clean ==> delete *.mm (mpp Maple preprocessor output) -# ellipsoid ==> compute ellipsoid equations # ifneq ($(MAPLE_VERSION),) @@ -33,7 +32,3 @@ mm : $(patsubst %.maple, %.mm, $(wildcard *.maple)) .PHONY : clean clean : -rm *.mm - -.PHONY : ellipsoid -ellipsoid: - maple <ellipsoid.maple >ellipsoid.out |