diff options
Diffstat (limited to 'src/gr')
-rw-r--r-- | src/gr/AHFinderDirect.hh | 16 | ||||
-rw-r--r-- | src/gr/Newton.cc | 49 | ||||
-rw-r--r-- | src/gr/cg.hh | 1 | ||||
-rw-r--r-- | src/gr/driver.cc | 248 | ||||
-rw-r--r-- | src/gr/gfn.hh | 3 | ||||
-rw-r--r-- | src/gr/horizon_Jacobian.cc | 106 | ||||
-rw-r--r-- | src/gr/horizon_function.cc | 8 |
7 files changed, 307 insertions, 124 deletions
diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index a8305d2..ac79f14 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -73,6 +73,7 @@ struct solver_info // stuff for Newton_solve() int max_Newton_iterations; fp H_norm_for_convergence; + fp Delta_h_norm_for_convergence; }; //****************************************************************************** @@ -101,12 +102,17 @@ void horizon_Jacobian(patch_system& ps, const struct geometry_interpolator_info& gii, const struct Jacobian_info& Jac_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); // Newton.cc -void Newton_solve(patch_system& ps, +// 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 char Jacobian_type[], - fp perturbation_amplitude, - int max_Newton_iterations, - fp H_norm_for_convergence); + const struct Jacobian_info& Jac_info, + const struct solver_info& solver_info, + Jacobian& Jac); diff --git a/src/gr/Newton.cc b/src/gr/Newton.cc index 011f34c..4c2eb04 100644 --- a/src/gr/Newton.cc +++ b/src/gr/Newton.cc @@ -51,44 +51,39 @@ namespace { // // This function solves H(h) = 0 for h via Newton's method. -// At the moment the NP Jacobian is used. // -void Newton_solve(patch_system& ps, +bool Newton_solve(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_interpolator_info& gii, - const char Jacobian_type[], - fp perturbation_amplitude, - int max_Newton_iterations, - fp H_norm_for_convergence) + const struct Jacobian_info& Jac_info, + const struct solver_info& solver_info, + Jacobian& Jac) { -Jacobian& Jac = create_Jacobian(ps, Jacobian_type); - for (int iteration = 1 ; - iteration <= max_Newton_iterations ; + iteration <= solver_info.max_Newton_iterations ; ++iteration) { CCTK_VInfo(CCTK_THORNSTRING, "Newton iteration %d", iteration); jtutil::norm<fp> H_norms; - horizon_function(ps, cgi, gii, true, H_norms); + horizon_function(ps, cgi, gii, true, &H_norms); CCTK_VInfo(CCTK_THORNSTRING, " H rms-norm=%.2e, infinity-norm=%.2e", H_norms.rms_norm(), H_norms.infinity_norm()); - if (H_norms.infinity_norm() <= H_norm_for_convergence) + if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence) then { - CCTK_VInfo(CCTK_THORNSTRING, "==> finished!"); - return; // *** NORMAL RETURN *** + CCTK_VInfo(CCTK_THORNSTRING, + "==> finished (||H|| < %g)", + double(solver_info.H_norm_for_convergence)); + return true; // *** NORMAL RETURN *** } // take a Newton step - horizon_Jacobian_NP(ps, cgi, gii, Jac, perturbation_amplitude); - jtutil::array_scale(Jac.NN(), - -1.0, - ps.gridfn_data(nominal_gfns::gfn__H)); + horizon_Jacobian(ps, cgi, gii, Jac_info, Jac); CCTK_VInfo(CCTK_THORNSTRING, - "solving linear system"); + "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 CCTK_VInfo(CCTK_THORNSTRING, @@ -108,16 +103,26 @@ Jacobian& Jac = create_Jacobian(ps, Jacobian_type); irho,isigma); Delta_h_norms.data(Delta_h); - p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) += Delta_h; + p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) -= Delta_h; } } } CCTK_VInfo(CCTK_THORNSTRING, "moved h by Delta_h (rhs-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)); + return true; // *** NORMAL RETURN *** + } } -CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Newton_solve: failed to converge in %d iterations!", - max_Newton_iterations); +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Newton_solve: no convergence in %d iterations!\n", + solver_info.max_Newton_iterations); +return false; // *** ERROR RETURN *** } diff --git a/src/gr/cg.hh b/src/gr/cg.hh index cc0ab31..6535753 100644 --- a/src/gr/cg.hh +++ b/src/gr/cg.hh @@ -102,6 +102,7 @@ #define partial_H_wrt_partial_dd_h_22 \ p.gridfn(nominal_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) //****************************************************************************** diff --git a/src/gr/driver.cc b/src/gr/driver.cc index ca978ae..00695e9 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -7,6 +7,8 @@ /// 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> @@ -54,6 +56,10 @@ void setup_Kerr_horizon(patch_system& ps, 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, const Jacobian& NP_Jac, + const char file_name[]); }; //****************************************************************************** @@ -155,91 +161,124 @@ patch_system ps(origin_x, origin_y, origin_z, if (STRING_EQUAL(initial_guess_method, "read from file")) then { CCTK_VInfo(CCTK_THORNSTRING, - "reading initial guess from \"%s\"", - initial_guess__read_from_file__file_name); + "reading initial guess h from \"%s\"", + h_file_name); ps.read_ghosted_gridfn(ghosted_gfns::gfn__h, - initial_guess__read_from_file__file_name, + 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_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, - 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, - true); // use Kerr-Schild coords -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "unknown initial_guess_method=\"%s\"!", - initial_guess_method); /*NOTREACHED*/ -ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, - true, ghosted_gfns::gfn__h, - "h.dat", - 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_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, + 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, + 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(ghosted_gfns::gfn__h, + true, ghosted_gfns::gfn__h, + h_file_name, + false); // no ghost zones + } // // find the apparent horizon // -jtutil::norm<fp> H_norms; +struct Jacobian_info Jac_info; +Jac_info.perturbation_amplitude = NP_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; + if (STRING_EQUAL(method, "horizon function")) then { - horizon_function(ps, cgi, gii, false, H_norms); + jtutil::norm<fp> H_norms; + + const int timer_handle = CCTK_TimerCreateI(); + CCTK_TimerStartI(timer_handle); + horizon_function(ps, cgi, gii, false, &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\n", + " 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(nominal_gfns::gfn__H, true, ghosted_gfns::gfn__h, - "H.dat"); - } -else if (STRING_EQUAL(method, "Jacobian")) - then { - Jacobian& Jac = create_Jacobian(ps, Jacobian_type); - horizon_function(ps, cgi, gii, true, H_norms); - horizon_Jacobian_SD(ps, Jac); - print_Jacobian(Jacobian_file_name, Jac); + H_of_h_file_name); } else if (STRING_EQUAL(method, "Jacobian test")) then { - Jacobian& SD_Jac = create_Jacobian(ps, Jacobian_type); - horizon_function(ps, cgi, gii, true, H_norms); - horizon_Jacobian_SD(ps, SD_Jac); + 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, - NP_Jac, - NP_Jacobian__perturbation_amplitude); + horizon_function(ps, cgi, gii, true, &H_norms); + horizon_Jacobian_NP(ps, cgi, gii, Jac_info, NP_Jac); - print_Jacobians(Jacobian_file_name, SD_Jac, NP_Jac); + CCTK_VInfo(CCTK_THORNSTRING, + "writing Jacobians to \"%s\"", + Jacobian_file_name); + print_Jacobians(ps, Jac, NP_Jac, Jacobian_file_name); } -else if (STRING_EQUAL(method, "Newton (NP Jacobian)")) +else if (STRING_EQUAL(method, "Newton solve")) then { - Newton_solve(ps, cgi, gii, - Jacobian_type, - NP_Jacobian__perturbation_amplitude, - max_Newton_iterations, - H_norm_for_convergence); + Jacobian& Jac = create_Jacobian(ps, Jacobian_type); + + const int timer_handle = CCTK_TimerCreateI(); + CCTK_TimerStartI(timer_handle); + Newton_solve(ps, cgi, gii, Jac_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(ghosted_gfns::gfn__h, true, ghosted_gfns::gfn__h, - "h.dat", + 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, + H_of_h_file_name); } else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown method=\"%s\"!", @@ -384,3 +423,92 @@ CCTK_VInfo(CCTK_THORNSTRING, } } } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function prints a pair of Jacobian matrices (and their difference) +// to a named output file. +// +namespace { +void print_Jacobians(const patch_system& ps, + const Jacobian& Jac, const Jacobian& NP_Jac, + 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 II\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 = 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"); + + 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.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; + + if ((SD == 0.0) && (NP == 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); + + 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)); + } + } + } + } + } + } + +fclose(fileptr); +} + } diff --git a/src/gr/gfn.hh b/src/gr/gfn.hh index 9f54b5c..9bf5aac 100644 --- a/src/gr/gfn.hh +++ b/src/gr/gfn.hh @@ -4,7 +4,7 @@ namespace ghosted_gfns { enum { - min_gfn = -1, + min_gfn = -1, // must set this by hand so max_gfn is still < 0 gfn__h = min_gfn, max_gfn = gfn__h }; @@ -53,6 +53,7 @@ enum { gfn__partial_H_wrt_partial_dd_h_11, gfn__partial_H_wrt_partial_dd_h_12, gfn__partial_H_wrt_partial_dd_h_22, + gfn__save_H, gfn__Delta_h, max_gfn = gfn__Delta_h // no comma }; diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc index 99d8d44..f16d9c5 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -8,6 +8,7 @@ /// 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> @@ -85,25 +86,75 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // 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 finite differencing. +// 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 +// partial_H_wrt_partial_dd_h # (also assumed to already be computed) // // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. // -void horizon_Jacobian_SD(patch_system& ps, - Jacobian& 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) { -CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_SD"); +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); ps.compute_synchronize_Jacobian(); Jac.zero(); + +// +// compute d/dr term in Jacobian by numerical perturbation +// + +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); + 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(nominal_gfns::gfn__save_H, + irho,isigma); + const fp new_H = p.gridfn(nominal_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); + + +// +// 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); @@ -213,10 +264,6 @@ Jac.zero(); } } } - -// compute d/dr term by numerical finite differencing -perturb_h(Jacobian_info.perturbation_amplitude); -horizon_function(ps, cgi, gii, false, NULL, false); } //****************************************************************************** @@ -281,36 +328,31 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn // by numerical perturbation of the H(h) function. The algorithm is as // follows: // -// evaluate H(h) -// array-copy H(h) to scratch array save_H +// 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) +// 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; // } -// array-copy save_H back to H(h) +// H = save_H // void horizon_Jacobian_NP(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_interpolator_info& ii, - Jacobian& Jac, - fp perturbation_amplitude) + const struct Jacobian_info& Jac_info, + Jacobian& Jac) { CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_NP"); +const fp epsilon = Jac_info.perturbation_amplitude; -// evaluate H(h) -jtutil::norm<fp> H_norms; -horizon_function(ps, cgi, ii, false, H_norms); - -// array-copy H(h) to scratch array save_H -fp *save_H = new fp[Jac.NN()]; -jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H); +ps.gridfn_copy(nominal_gfns::gfn__H, nominal_gfns::gfn__save_H); for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) { @@ -327,10 +369,8 @@ jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H); 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) - += perturbation_amplitude; - H_norms.reset(); - horizon_function(ps, cgi, ii, false, H_norms, false); + 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); @@ -344,18 +384,18 @@ jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H); ++x_isigma) { const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); - const fp new_xH = xp.gridfn(nominal_gfns::gfn__H, - x_irho,x_isigma); - Jac(II,JJ) = (new_xH - save_H[II]) / perturbation_amplitude; + 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; + yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) = save_h_y; } } } -jtutil::array_copy(Jac.NN(), save_H, ps.gridfn_data(nominal_gfns::gfn__H)); -delete[] save_H; +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 95bb6b6..9f05698 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_ptr, + jtutil::norm<fp>* H_norms_ptr, bool msg_flag); } @@ -100,10 +100,12 @@ void compute_H(patch_system& ps, // 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& ii, + const struct geometry_interpolator_info& gii, bool Jacobian_flag, jtutil::norm<fp>* H_norms_ptr, bool msg_flag = true) @@ -118,7 +120,7 @@ ps.synchronize(); setup_xyz_posns(ps, msg_flag); // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$ -interpolate_geometry(ps, cgi, ii, msg_flag); +interpolate_geometry(ps, cgi, gii, msg_flag); // compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing |