diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-22 14:09:21 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-22 14:09:21 +0000 |
commit | f83ce9603ea6d88ff644c218a7585b0767b55a7b (patch) | |
tree | 7c1742ce1af9b72b6fb83324edfe9b8a336e3b51 /src | |
parent | 4e68d75167ec89786c3a8819f3bc851b8a6f3fcd (diff) |
add some code for the d/dr term in the Jacobian -- not finished yet
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@652 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src')
-rw-r--r-- | src/gr/AHFinderDirect.hh | 34 | ||||
-rw-r--r-- | src/gr/horizon_Jacobian.cc | 14 | ||||
-rw-r--r-- | src/gr/horizon_function.cc | 17 |
3 files changed, 46 insertions, 19 deletions
diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index 42b77c9..a8305d2 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -55,6 +55,26 @@ struct cactus_grid_info const fp* K_dd_33_data; }; +// +// This struct holds the various parameters used in computing the Jacobian +// matrix. +// +struct Jacobian_info + { + fp perturbation_amplitude; + }; + +// +// This struct holds the various parameters used in solving the H(h) = 0 +// equations. +// +struct solver_info + { + // stuff for Newton_solve() + int max_Newton_iterations; + fp H_norm_for_convergence; + }; + //****************************************************************************** // @@ -70,19 +90,17 @@ 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, + jtutil::norm<fp>* H_norms_ptr, bool msg_flag = true); // horizon_Jacobian.cc Jacobian& create_Jacobian(patch_system& ps, const char Jacobian_type[]); -void horizon_Jacobian_SD(patch_system& ps, - Jacobian& Jac); -void horizon_Jacobian_NP(patch_system& ps, - const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& ii, - Jacobian& Jac, - fp perturbation_amplitude); +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); // Newton.cc void Newton_solve(patch_system& ps, diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index bcaec78..99d8d44 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -4,7 +4,7 @@ // <<<prototypes for functions local to this file>>> // create_Jacobian - create a Jacobian matrix of the appropriate type // -// horizon_Jacobian_SD - compute the Jacobian by symbolic differentiation +// horizon_Jacobian - compute the Jacobian (symbolic differentiation) /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian // // horizon_Jacobian_NP - compute the Jacobian by numerical perturbation @@ -84,12 +84,14 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // 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. +// equation (25) of my 1996 apparent horizon finding paper, except that +// the d/dr term is done by numerical finite differencing. // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface -// partial_H_wrt_partial_d_h, # Jacobian coefficients -// partial_H_wrt_partial_dd_h, +// H # H(h) assumed to already be computed +// partial_H_wrt_partial_d_h # Jacobian coefficients +// partial_H_wrt_partial_dd_h // // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. @@ -211,6 +213,10 @@ Jac.zero(); } } } + +// compute d/dr term by numerical finite differencing +perturb_h(Jacobian_info.perturbation_amplitude); +horizon_function(ps, cgi, gii, false, NULL, false); } //****************************************************************************** diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 9c6cb17..95bb6b6 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -53,7 +53,7 @@ void interpolate_geometry(patch_system& ps, bool msg_flag); void compute_H(patch_system& ps, bool Jacobian_flag, - jtutil::norm<fp>& H_norms, + jtutil::norm<fp>& H_norms_ptr, bool msg_flag); } @@ -96,13 +96,16 @@ void compute_H(patch_system& ps, // Arguments: // Jacobian_flag = true to compute the Jacobian coefficients, // false to skip this. -// H_norms = (out) This is set to the gridwise norms of the H(h) function. +// 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, const struct cactus_grid_info& cgi, const struct geometry_interpolator_info& ii, bool Jacobian_flag, - jtutil::norm<fp>& H_norms, + jtutil::norm<fp>* H_norms_ptr, bool msg_flag = true) { if (msg_flag) @@ -119,7 +122,7 @@ interpolate_geometry(ps, cgi, ii, msg_flag); // compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing -compute_H(ps, Jacobian_flag, H_norms, msg_flag); +compute_H(ps, Jacobian_flag, H_norms_ptr, msg_flag); } //****************************************************************************** @@ -435,12 +438,11 @@ if (status < 0) // Arguments: // Jacobian_flag = true to compute the Jacobian coefficients, // false to skip this. -// H_norms = (out) This is set to the gridwise norms of the H(h) function. // namespace { void compute_H(patch_system& ps, bool Jacobian_flag, - jtutil::norm<fp>& H_norms, + jtutil::norm<fp>* H_norms_ptr, bool msg_flag) { if (msg_flag) @@ -524,7 +526,8 @@ if (msg_flag) } // update running norms of H(h) function - H_norms.data(H); + if (H_norms_ptr != NULL) + then H_norms_ptr->data(H); if (Jacobian_flag) then { |