diff options
| -rw-r--r-- | param.ccl | 66 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Jacobian-small.par | 7 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset-small-ahf.par | 10 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset-small.par | 5 | ||||
| -rw-r--r-- | run/test-ahfinderdirect/Kerr-offset.par | 5 | ||||
| -rw-r--r-- | src/elliptic/Jacobian.cc | 146 | ||||
| -rw-r--r-- | src/elliptic/Jacobian.hh | 18 | ||||
| -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 | ||||
| -rw-r--r-- | src/jtutil/miscfp.cc | 81 | ||||
| -rw-r--r-- | src/jtutil/util.hh | 10 | ||||
| -rw-r--r-- | src/make.configuration.defn | 4 | ||||
| -rw-r--r-- | src/patch/patch_system.cc | 51 | ||||
| -rw-r--r-- | src/patch/patch_system.hh | 12 |
19 files changed, 440 insertions, 406 deletions
@@ -11,7 +11,7 @@ keyword method "top-level method used to find the apparent horizon" "horizon function" :: "evaluate the LHS function H(h)" "Jacobian test" :: "test the J[H(h)] Jacobian matrix (SD vs NP methods)" "Newton solve" :: "find the horizon via Newton's method" -} "horizon" +} "horizon function" keyword Jacobian_type \ "what type of storage scheme do we use for the Jacobian matrix?" @@ -24,10 +24,40 @@ int max_Newton_iterations "maximum number of Newton iterations before giving up" (0:* :: "any positive integer" } 10 +# +# we declare convergence if *either* of the following two criteria are met +# real H_norm_for_convergence "declare convergence if ||H||_inf <= this" { (0.0:* :: "any positive real number" -} 1.0e-6 +} 1.0e-10 + +real Delta_h_norm_for_convergence \ + "declare convergence after any Newton step with ||Delta_h||_inf <= this" +{ +(0.0:* :: "any positive real number" +} 1.0e-10 + +################################################################################ + +# +# input/output file name parameters +# + +string h_file_name "file name from/to which to read/write horizon shape h" +{ +.* :: "any string" +} "h.dat" + +string H_of_h_file_name "file name to which to write final H(h) function" +{ +.* :: "any string" +} "H.dat" + +string Jacobian_file_name "name of the Jacobian output file" +{ +.+ :: "any valid file name" +} "Jacobian.dat" ################################################################################ @@ -46,22 +76,6 @@ real NP_Jacobian__perturbation_amplitude \ (0.0:* :: "any real number > 0" } 1.0e-6 -# -# true ==> gives a more thorough test of the Jacobian, -# but makes the test run much slower -# false ==> gives a slightly less thorough test, but runs faster -# -boolean NP_Jacobian__perturb_all_y_patch_points \ - "should we perturb at *all* points in the y patch, or just those with the \ - iperp which is (supposedly) involved in the interpatch interpolation?" -{ -} true; - -string Jacobian_file_name "name of the Jacobian output file" -{ -.+ :: "any valid file name" -} "Jacobian.dat" - ################################################################################ # @@ -158,13 +172,6 @@ keyword initial_guess_method \ "Kerr/Kerr-Schild" :: "set to horizon of Kerr spacetime in Kerr-Schild coords" } "read from file" -# parameters for initial_guess_method = "read from file" -string initial_guess__read_from_file__file_name \ - "file name from which to read initial guess for apparent horizon shape" -{ -.* :: "any string" -} "h.dat" - # parameters for initial_guess_method = "ellipsoid" real initial_guess__ellipsoid__x_center "x coordinate of ellipsoid center" { @@ -235,4 +242,13 @@ int which_derivs "bit flags to specify which derivatives to test" 0:63 :: "any set of bit flags" } 63 +# true ==> gives a more thorough test of the Jacobian, +# but makes the test run much slower +# false ==> gives a slightly less thorough test, but runs faster +boolean NP_Jacobian__perturb_all_y_patch_points \ + "should we perturb at *all* points in the y patch, or just those with the \ + iperp which is (supposedly) involved in the interpatch interpolation?" +{ +} true; + ################################################################################ diff --git a/run/test-ahfinderdirect/Jacobian-small.par b/run/test-ahfinderdirect/Jacobian-small.par index 39fb556..7f0e72d 100644 --- a/run/test-ahfinderdirect/Jacobian-small.par +++ b/run/test-ahfinderdirect/Jacobian-small.par @@ -26,6 +26,10 @@ Exact::exact_model = "Kerr/Kerr-Schild" Exact::Kerr_KerrSchild__mass = 1.0 Exact::Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::method = "Jacobian test" +AHFinderDirect::Jacobian_type = "dense matrix" +AHFinderDirect::Jacobian_file_name = "Jacobian-small.dat" + AHFinderDirect::origin_x = 0.5 AHFinderDirect::origin_y = 0.7 AHFinderDirect::origin_z = 0.9 @@ -42,6 +46,3 @@ AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 -AHFinderDirect::method = "Jacobian (SD+NP)" -AHFinderDirect::NP_Jacobian__perturbation_amplitude = 1.0e-6 -AHFinderDirect::Jacobian_file_name = "Jacobian-small.dat" diff --git a/run/test-ahfinderdirect/Kerr-offset-small-ahf.par b/run/test-ahfinderdirect/Kerr-offset-small-ahf.par index 9550fbd..fcfd442 100644 --- a/run/test-ahfinderdirect/Kerr-offset-small-ahf.par +++ b/run/test-ahfinderdirect/Kerr-offset-small-ahf.par @@ -26,6 +26,10 @@ Exact::exact_model = "Kerr/Kerr-Schild" Exact::Kerr_KerrSchild__mass = 1.0 Exact::Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::method = "Newton solve" +AHFinderDirect::h_file_name = "Kerr-offset-small-ahf.h.dat" +AHFinderDirect::H_of_h_file_name = "Kerr-offset-small-ahf.H.dat" + AHFinderDirect::origin_x = 0.5 AHFinderDirect::origin_y = 0.7 AHFinderDirect::origin_z = 0.9 @@ -43,9 +47,3 @@ AHFinderDirect::initial_guess_method = "ellipsoid" AHFinderDirect::initial_guess__ellipsoid__x_radius = 1.8 AHFinderDirect::initial_guess__ellipsoid__y_radius = 1.8 AHFinderDirect::initial_guess__ellipsoid__z_radius = 1.8 - -AHFinderDirect::method = "Newton (NP Jacobian)" -AHFinderDirect::Jacobian_type = "dense matrix" -AHFinderDirect::max_Newton_iterations = 10 -AHFinderDirect::H_norm_for_convergence = 1.0e-10 -AHFinderDirect::NP_Jacobian__perturbation_amplitude = 1.0e-6 diff --git a/run/test-ahfinderdirect/Kerr-offset-small.par b/run/test-ahfinderdirect/Kerr-offset-small.par index c83c33e..0a873b6 100644 --- a/run/test-ahfinderdirect/Kerr-offset-small.par +++ b/run/test-ahfinderdirect/Kerr-offset-small.par @@ -26,6 +26,10 @@ Exact::exact_model = "Kerr/Kerr-Schild" Exact::Kerr_KerrSchild__mass = 1.0 Exact::Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::method = "horizon function" +AHFinderDirect::h_file_name = "Kerr-offset-small.h.dat" +AHFinderDirect::H_of_h_file_name = "Kerr-offset-small.H.dat" + AHFinderDirect::origin_x = 0.5 AHFinderDirect::origin_y = 0.7 AHFinderDirect::origin_z = 0.9 @@ -42,4 +46,3 @@ AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 -AHFinderDirect::method = "horizon" diff --git a/run/test-ahfinderdirect/Kerr-offset.par b/run/test-ahfinderdirect/Kerr-offset.par index 61b562b..301e84f 100644 --- a/run/test-ahfinderdirect/Kerr-offset.par +++ b/run/test-ahfinderdirect/Kerr-offset.par @@ -26,6 +26,10 @@ Exact::exact_model = "Kerr/Kerr-Schild" Exact::Kerr_KerrSchild__mass = 1.0 Exact::Kerr_KerrSchild__spin = 0.6 +AHFinderDirect::method = "horizon function" +AHFinderDirect::h_file_name = "Kerr-offset.h.dat" +AHFinderDirect::H_of_h_file_name = "Kerr-offset.H.dat" + AHFinderDirect::origin_x = 0.5 AHFinderDirect::origin_y = 0.7 AHFinderDirect::origin_z = 0.9 @@ -42,4 +46,3 @@ AHFinderDirect::interpatch_interpolator_pars = "order=3" AHFinderDirect::initial_guess_method = "Kerr/Kerr-Schild" AHFinderDirect::initial_guess__Kerr_KerrSchild__mass = 1.0 AHFinderDirect::initial_guess__Kerr_KerrSchild__spin = 0.6 -AHFinderDirect::method = "horizon" diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 4d3e411..58009ec 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -4,9 +4,6 @@ // // Jacobian::Jacobian // -// print_Jacobian -// print_Jacobians -// // dense_Jacobian::dense_Jacobian // dense_Jacobian::~dense_Jacobian // dense_Jacobian::zero @@ -77,12 +74,6 @@ void CCTK_FCALL #endif //************************************** -namespace { -void print_Jacobians_internal(const char file_name[], - const Jacobian& SD_Jac, const Jacobian& NP_Jac, - bool pair_flag); - }; - //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -100,130 +91,6 @@ Jacobian::Jacobian(patch_system& ps) //****************************************************************************** // -// This function prints a Jacobian matrix to a named output file. -// -void print_Jacobian(const char file_name[], const Jacobian& Jac) -{ -print_Jacobians_internal(file_name, Jac, Jac, false); -} - -//****************************************************************************** - -// -// This function prints a pair of Jacobian matrices (and their difference) -// to a named output file. -// -void print_Jacobians(const char file_name[], - const Jacobian& SD_Jac, const Jacobian& NP_Jac) -{ -print_Jacobians_internal(file_name, SD_Jac, NP_Jac, true); -} - -//****************************************************************************** - -// -// If pair_flag = false, this prints SD_Jac. -// If pair_flag = true, this prints both Jacobians, and the error in SD_Jac. -// -namespace { -void print_Jacobians_internal(const char file_name[], - const Jacobian& SD_Jac, const Jacobian& NP_Jac, - bool pair_flag) -{ -const patch_system& ps = SD_Jac.my_patch_system(); - -FILE *fileptr = fopen(file_name, "w"); -if (fileptr == NULL) - then error_exit(ERROR_EXIT, -"***** dense_Jacobian::print(): 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"); -if (pair_flag) - then { - fprintf(fileptr, "# column 9 = SD_Jac(II,JJ)\n"); - fprintf(fileptr, "# column 10 = NP_Jac(II,JJ)\n"); - fprintf(fileptr, "# column 11 = abs error\n"); - fprintf(fileptr, "# column 12 = rel error\n"); - } - else fprintf(fileptr, "# column 9 = Jac(II,JJ)\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 (! SD_Jac.is_explicitly_stored(II,JJ)) - then continue; // skip sparse points - - const fp SD = SD_Jac(II,JJ); - const fp NP = NP_Jac(II,JJ); - const fp abs_error = SD - NP; - - if (pair_flag ? ((SD == 0.0) && (NP == 0.0)) - : (SD == 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", - II, xpn, x_irho, x_isigma, - JJ, ypn, y_irho, y_isigma); - if (pair_flag) - then fprintf(fileptr, - "%.10g\t%.10g\t%e\t%e\n", - double(SD), double(NP), - double(abs_error), double(rel_error)); - else fprintf(fileptr, - "%.10g\n", - double(SD)); - } - } - } - } - } - } - -fclose(fileptr); -} - }; - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// // This function constructs a dense_Jacobian object. // dense_Jacobian::dense_Jacobian(patch_system& ps) @@ -249,7 +116,13 @@ delete[] pivot_; // void dense_Jacobian::zero() { -jtutil::array_zero(matrix_.N_array(), matrix_.data_array()); + for (int JJ = 0 ; JJ < NN() ; ++JJ) + { + for (int II = 0 ; II < NN() ; ++II) + { + matrix_(JJ,II) = 0.0; + } + } } //****************************************************************************** @@ -282,7 +155,10 @@ fp *x = my_patch_system().gridfn_data(x_gfn); // [sd]gesv() use an "in out" design, where the same argument is used for // both rhs and x ==> first copy rhs to x so we can pass that to [sd]gesv() // -jtutil::array_copy(NN(), rhs, x); + for (int II = 0 ; II < NN() ; ++II) + { + x[II] = rhs[II]; + } integer N = NN(); integer NRHS = 1; diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh index d948f07..6d20e63 100644 --- a/src/elliptic/Jacobian.hh +++ b/src/elliptic/Jacobian.hh @@ -3,10 +3,6 @@ // // Jacobian -- abstract base class to describe a Jacobian matrix -// -// print_Jacobian - print a Jacobian -// print_Jacobians - print two Jacobians -// // dense_Jacobian - Jacobian stored as a dense matrix // @@ -29,8 +25,6 @@ // //****************************************************************************** -//****************************************************************************** -//****************************************************************************** // // A Jacobian object stores the (a) Jacobian matrix for a patch system. @@ -101,18 +95,6 @@ protected: //****************************************************************************** // -// These functions print 1 or 2 Jacobians to a named output file; -// in the latter case they also print the error in the SD Jacobian. -// -void print_Jacobian(const char file_name[], const Jacobian& Jac); -void print_Jacobians(const char file_name[], - const Jacobian& SD_Jac, const Jacobian& NP_Jac); - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// // This class stores the Jacobian as a dense matrix in Fortran (column) // order. // 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 diff --git a/src/jtutil/miscfp.cc b/src/jtutil/miscfp.cc index 4263d2c..20060a2 100644 --- a/src/jtutil/miscfp.cc +++ b/src/jtutil/miscfp.cc @@ -6,12 +6,6 @@ // jtutil::arctan_xy - 4-quadrant arc tangent // jtutil::modulo_reduce - reduce an angle modulo 2*pi radians (360 degrees) // -// jtutil::array_zero - zero an array -// jtutil::array_copy - copy one array to another -// jtutil::array_scale - scale an array by a scalar -// -// ***** template instantiations ***** -// #include <math.h> #include "stdc.h" @@ -121,78 +115,3 @@ if (! (fuzzy<double>::GE(xx, xmin) && fuzzy<double>::LE(xx, xmax)) ) return xx; } } // namespace jtutil:: - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function zeros the size-N array dst[]. -// -namespace jtutil - { -template <typename fp> - void array_zero(int N, fp dst[]) -{ - for (int i = 0 ; i < N ; ++i) - { - dst[i] = 0.0; - } -} - } // namespace jtutil:: - -//****************************************************************************** - -// -// This function copies the size-N array src[] to the size-N array dst[], -// analogously to the BLAS routines xCOPY(). -// -namespace jtutil - { -template <typename fp> - void array_copy(int N, const fp src[], fp dst[]) -{ - for (int i = 0 ; i < N ; ++i) - { - dst[i] = src[i]; - } -} - } // namespace jtutil:: - -//****************************************************************************** - -// -// This function multiplies the size-N array dst[] by the scalar alpha, -// analogously to the BLAS routines xSCAL(). -// -namespace jtutil - { -template <typename fp> - void array_scale(int N, fp alpha, fp dst[]) -{ - for (int i = 0 ; i < N ; ++i) - { - dst[i] *= alpha; - } -} - } // namespace jtutil:: - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// ***** template instantiations ***** -// - -namespace jtutil - { -template void array_zero<float >(int, float []); -template void array_zero<double>(int, double[]); - -template void array_copy<float >(int, const float [], float []); -template void array_copy<double>(int, const double[], double[]); - -template void array_scale<float >(int, float , float []); -template void array_scale<double>(int, double, double[]); - } // namespace jtutil:: diff --git a/src/jtutil/util.hh b/src/jtutil/util.hh index 89d6bde..9273ebe 100644 --- a/src/jtutil/util.hh +++ b/src/jtutil/util.hh @@ -62,16 +62,6 @@ double arctan_xy(double x, double y); double modulo_reduce(double x, double xmod, double xmin, double xmax); // -// simple BLAS-like array arithmetic routines -// -template <typename fp> - void array_zero(int N, fp dst[]); // dst[] = 0 -template <typename fp> - void array_copy(int N, const fp src[], fp dst[]); // dst[] = src[] -template <typename fp> - void array_scale(int N, fp alpha, fp dst[]); // dst[] *= alpha - -// // more misc math stuff, valid only if <math.h> has been #included // #ifdef PI // from "jt/stdc.h" diff --git a/src/make.configuration.defn b/src/make.configuration.defn index 06830d3..ecc4f59 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 -LIBS += lapack blas +##LIBDIRS += /usr/local/lib +LIBS := lapack blas $(LIBS) diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index fa0d28d..2fb1f17 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -22,6 +22,9 @@ // patch_system::type_of_name // patch_system::patch_number_of_name // +// patch_system::gridfn_copy +// patch_system::add_to_ghosted_gridfn +// // patch_system::synchronize // patch_system::compute_synchronize_Jacobian // patch_system::synchronize_Jacobian_y_ipar_posn @@ -1060,6 +1063,54 @@ error_exit(ERROR_EXIT, //****************************************************************************** // +// This function copies one (nominal-grid) gridfn to another. +// +void patch_system::gridfn_copy(int src_gfn, int dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = 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) + { + p.gridfn(dst_gfn, irho,isigma) = p.gridfn(src_gfn, irho,isigma); + } + } + } +} + +//****************************************************************************** + +// +// This function adds a scalar to a ghosted gridfn. +// +void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; + ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) += delta; + } + } + } +} + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// // This function "synchronizes" all ghost zones of all patches, i.e. it // update the ghost-zone values of the specified gridfns via the appropriate // sequence of symmetry operations and interpatch interpolations. This diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh index a1b4227..bc1cdbb 100644 --- a/src/patch/patch_system.hh +++ b/src/patch/patch_system.hh @@ -128,6 +128,7 @@ public: patch& ghosted_patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma) const; + // // ***** meta-info about gridfns ***** // @@ -142,6 +143,17 @@ public: // + // ***** gridfn operations ***** + // + + // dst = src + void gridfn_copy(int src_gfn, int dst_gfn); + + // dst += delta + void add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn); + + + // // ***** access to gridfns as 1-D arrays[gpn] ***** // ... n.b. this interface implicitly assumes that gridfn data // arrays are contiguous across patches; this is ensured by |