From b3c9f3f1269b213537b25a1ff0123871a5e07f9f Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 22 Jul 2002 12:29:42 +0000 Subject: add a whole bunch of changes from working at home --> AHFinderDirect now finds AH correctly for Kerr/offset!!! git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@648 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/Jacobian.cc | 238 +++++++++++++++++++++++++++++++++++++++++++-- src/elliptic/Jacobian.hh | 54 +++++++--- src/gr/AHFinderDirect.hh | 24 ++++- src/gr/cg.hh | 2 + src/gr/driver.cc | 192 +++++++++++++++++++++--------------- src/gr/gfn.hh | 3 +- src/gr/horizon_Jacobian.cc | 115 ++++++++++++++++++++-- src/gr/horizon_function.cc | 62 ++++++++---- src/gr/make.code.defn | 3 +- src/include/config.hh | 6 ++ src/jtutil/README | 15 +++ src/jtutil/miscfp.cc | 83 ++++++++++++++++ src/jtutil/util.hh | 24 +++-- src/patch/patch_system.cc | 45 +++++---- 14 files changed, 704 insertions(+), 162 deletions(-) (limited to 'src') diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 9b3e4d8..4d3e411 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -4,12 +4,22 @@ // // Jacobian::Jacobian // +// print_Jacobian +// print_Jacobians +// // dense_Jacobian::dense_Jacobian +// dense_Jacobian::~dense_Jacobian // dense_Jacobian::zero // dense_Jacobian::zero_row +// dense_Jacobian::solve_linear_system // #include +using std::fopen; +using std::printf; +using std::fprintf; +using std::fclose; +using std::FILE; #include #include #include @@ -36,6 +46,42 @@ using jtutil::error_exit; #include "../util/patch_system.hh" #include "Jacobian.hh" +//#include "lapack.h" +//************************************** +/* lapack.h -- C/C++ prototypes for (some) LAPACK routines */ +/* $Id$ */ + +/* + * prerequisites: + * "cctk.h" + * "config.hh" // for "integer" = Fortran integer + */ + +#ifdef __cplusplus +extern "C" { +#endif + +void CCTK_FCALL + CCTK_FNAME(sgesv)(const integer* N, const integer* NRHS, + float A[], const int* LDA, + integer IPIV[], + float B[], const integer* LDB, integer* info); +void CCTK_FCALL + CCTK_FNAME(dgesv)(const integer* N, const integer* NRHS, + double A[], const int* LDA, + integer IPIV[], + double B[], const integer* LDB, integer* info); + +#ifdef __cplusplus + }; /* extern "C" */ +#endif +//************************************** + +namespace { +void print_Jacobians_internal(const char file_name[], + const Jacobian& SD_Jac, const Jacobian& NP_Jac, + bool pair_flag); + }; //****************************************************************************** //****************************************************************************** @@ -44,31 +90,166 @@ using jtutil::error_exit; // // This function constructs a Jacobian object. // -Jacobian::Jacobian(const patch_system& ps) - : ps_(ps) +Jacobian::Jacobian(patch_system& ps) + : ps_(ps), + NN_(ps.N_grid_points()) { } //****************************************************************************** //****************************************************************************** //****************************************************************************** +// +// 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(const patch_system& ps) +dense_Jacobian::dense_Jacobian(patch_system& ps) : Jacobian(ps), - matrix_(0, ps.N_grid_points(), - 0, ps.N_grid_points()) + matrix_(0,NN()-1, 0,NN()-1), + pivot_(new integer[NN()]) { } //****************************************************************************** +// +// THis function destroys a dense_Jacobian object. +// +dense_Jacobian::~dense_Jacobian() +{ +delete[] pivot_; +} + +//****************************************************************************** + // // This function zeros a dense_Jacobian object. // void dense_Jacobian::zero() { -// FIXME FIXME +jtutil::array_zero(matrix_.N_array(), matrix_.data_array()); } //****************************************************************************** @@ -78,5 +259,48 @@ void dense_Jacobian::zero() // void dense_Jacobian::zero_row(int II) { -// FIXME FIXME + for (int JJ = 0 ; JJ < NN() ; ++JJ) + { + matrix_(JJ,II) = 0.0; + } +} + +//****************************************************************************** + +// +// This function solves the linear system J.x = rhs, with rhs and x +// being nominal-grid gridfns. The computation is done using the LAPACK +// [sd]gesv() subroutines; which modify the Jacobian matrix J in-place +// for the LU decomposition. +// +void dense_Jacobian::solve_linear_system(int rhs_gfn, int x_gfn) +{ +const fp *rhs = my_patch_system().gridfn_data(rhs_gfn); +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); + +integer N = NN(); +integer NRHS = 1; +integer info; + +#ifdef FP_IS_FLOAT + CCTK_FNAME(sgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info); +#endif +#ifdef FP_IS_DOUBLE + CCTK_FNAME(dgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info); +#endif + +if (info != 0) + then error_exit(ERROR_EXIT, +"\n" +"***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n" +" error return info=%d from [sd]gesv() LAPACK routine!\n" +, + rhs_gfn, x_gfn, + int(info)); /*NOTREACHED*/ } diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh index 2c64171..d948f07 100644 --- a/src/elliptic/Jacobian.hh +++ b/src/elliptic/Jacobian.hh @@ -4,6 +4,9 @@ // // 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 // @@ -46,10 +49,6 @@ class Jacobian { -public: - // basic meta-info - patch_system& my_patch_system(); - public: // access a matrix element via row/column indices virtual const fp& operator()(int II, int JJ) const = 0; // rvalue @@ -72,21 +71,43 @@ public: return operator()(II,JJ); } + // is a given element explicitly stored, or implicitly 0 via sparsity + virtual bool is_explicitly_stored(int II, int JJ) const = 0; + // zero the whole matrix or a single row virtual void zero() = 0; virtual void zero_row(int II) = 0; + // solve linear system J.x = rhs + // ... rhs and x are nominal-grid gridfns + // ... may modify Jacobian matrix (eg for LU decomposition) + virtual void solve_linear_system(int rhs_gfn, int x_gfn) = 0; + // basic meta-info - const patch_system& my_patch_system() const { return ps_; } + patch_system& my_patch_system() const { return ps_; } + int NN() const { return NN_; } // constructor, destructor - Jacobian(const patch_system& ps); + // ... for analyze-factor-solve sparse-matrix derived classes, + // construction may include the "analyze" operation + Jacobian(patch_system& ps); virtual ~Jacobian() { } -private: - const patch_system& ps_; +protected: + patch_system& ps_; + int NN_; // ps_.N_grid_points() }; +//****************************************************************************** + +// +// 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); + //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -105,16 +126,27 @@ public: fp& operator()(int II, int JJ) // lvalue { return matrix_(JJ,II); } + bool is_explicitly_stored(int II, int JJ) const + { return true; } + // zero the whole matrix or a single row void zero(); void zero_row(int II); + // solve linear system J.x = rhs + // ... rhs and x are nominal-grid gridfns + // ... may modify Jacobian matrix (eg for LU decomposition) + void solve_linear_system(int rhs_gfn, int x_gfn); + // constructor, destructor public: - dense_Jacobian(const patch_system& ps); - ~dense_Jacobian() { } + dense_Jacobian(patch_system& ps); + ~dense_Jacobian(); private: - // subscripts are (JJ,II) + // subscripts are (JJ,II) (both 0-origin) jtutil::array2d matrix_; + + // pivot vector for LAPACK routines + integer *pivot_; }; diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index 8e5fc8c..42b77c9 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -69,10 +69,26 @@ extern "C" void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_interpolator_info& gii, - bool Jacobian_flag); + bool Jacobian_flag, + jtutil::norm& H_norms, + bool msg_flag = true); // horizon_Jacobian.cc -Jacobian& create_Jacobian(const patch_system& ps, +Jacobian& create_Jacobian(patch_system& ps, const char Jacobian_type[]); -void horizon_Jacobian(patch_system& ps, - Jacobian& Jac); +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); + +// Newton.cc +void 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); diff --git a/src/gr/cg.hh b/src/gr/cg.hh index 9008159..cc0ab31 100644 --- a/src/gr/cg.hh +++ b/src/gr/cg.hh @@ -102,6 +102,8 @@ #define partial_H_wrt_partial_dd_h_22 \ p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, 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 546a488..ca978ae 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -3,9 +3,10 @@ // // <<>> // AHFinderDirect_driver - top-level driver +/// +/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords) /// setup_ellipsoid - setup up a coordiante ellipsoid in h -/// setup_Kerr_KerrSchild_horizon - set up Kerr/Kerr-Schild horizon in h -// +/// #include #include @@ -46,12 +47,13 @@ using jtutil::error_exit; // namespace { +void setup_Kerr_horizon(patch_system& ps, + fp x_center, fp y_center, fp z_center, + 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 setup_Kerr_KerrSchild_horizon(patch_system& ps, - fp x_center, fp y_center, fp z_center, - fp mass, fp spin); }; //****************************************************************************** @@ -139,7 +141,6 @@ cgi.K_dd_33_data = static_cast(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz")); // // create the patch system and initialize the xyz derivative coefficients // -CCTK_VInfo(CCTK_THORNSTRING, " creating patch system"); 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, @@ -154,7 +155,7 @@ 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\"", + "reading initial guess from \"%s\"", initial_guess__read_from_file__file_name); ps.read_ghosted_gridfn(ghosted_gfns::gfn__h, initial_guess__read_from_file__file_name, @@ -168,13 +169,22 @@ else if (STRING_EQUAL(initial_guess_method, "ellipsoid")) 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_KerrSchild_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); + 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*/ @@ -187,24 +197,105 @@ ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, // // find the apparent horizon // -if (STRING_EQUAL(method, "horizon")) +jtutil::norm H_norms; +if (STRING_EQUAL(method, "horizon function")) then { - horizon_function(ps, cgi, gii, false); + horizon_function(ps, cgi, gii, false, H_norms); + CCTK_VInfo(CCTK_THORNSTRING, + " H(h) rms-norm %.2e, infinity-norm %.2e\n", + H_norms.rms_norm(), H_norms.infinity_norm()); ps.print_gridfn_with_xyz(nominal_gfns::gfn__H, true, ghosted_gfns::gfn__h, "H.dat"); } -else if (STRING_EQUAL(method, "horizon Jacobian")) +else if (STRING_EQUAL(method, "Jacobian")) then { Jacobian& Jac = create_Jacobian(ps, Jacobian_type); - horizon_function(ps, cgi, gii, true); - horizon_Jacobian(ps, Jac); + horizon_function(ps, cgi, gii, true, H_norms); + horizon_Jacobian_SD(ps, Jac); + print_Jacobian(Jacobian_file_name, Jac); + } +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); + + 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); + + print_Jacobians(Jacobian_file_name, SD_Jac, NP_Jac); + } +else if (STRING_EQUAL(method, "Newton (NP Jacobian)")) + then { + Newton_solve(ps, cgi, gii, + Jacobian_type, + NP_Jacobian__perturbation_amplitude, + max_Newton_iterations, + H_norm_for_convergence); + ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, + true, ghosted_gfns::gfn__h, + "h.dat", + false); // no ghost zones } else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown method=\"%s\"!", method); /*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]_center = 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_center, fp y_center, fp z_center, + 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, + " mass=%g, spin=J/m^2=%g, posn=(%g,%g,%g)", + double(m), double(a), + double(x_center), double(y_center), double(z_center)); + +// 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_center, y_center, z_center, + xy_radius, xy_radius, z_radius); +} + } + //****************************************************************************** // @@ -234,10 +325,10 @@ void setup_ellipsoid(patch_system& ps, fp x_radius, fp y_radius, fp z_radius) { CCTK_VInfo(CCTK_THORNSTRING, - " h = ellipsoid: center=(%g,%g,%g)", + "setting h = ellipsoid: center=(%g,%g,%g)", double(x_center), double(y_center), double(z_center)); CCTK_VInfo(CCTK_THORNSTRING, - " radius=(%g,%g,%g)", + " radius=(%g,%g,%g)", double(x_radius), double(y_radius), double(z_radius)); for (int pn = 0 ; pn < ps.N_patches() ; ++pn) @@ -275,7 +366,7 @@ CCTK_VInfo(CCTK_THORNSTRING, then r = r_minus; else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" expected exactly one r>0 solution, got 0 or 2!\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" @@ -293,62 +384,3 @@ CCTK_VInfo(CCTK_THORNSTRING, } } } - -//****************************************************************************** - -// -// This function sets up the horizon of a Kerr black hole in 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 of my AHFinderDirect notes. -// -// Arguments: -// [xyz]_center = The position of the Kerr black hole. -// (mass,spin) = Describe the Kerr black hole. -// -// FIXME FIXME: at present [xyz}_center are ignored. -// -namespace { -void setup_Kerr_KerrSchild_horizon(patch_system& ps, - fp x_center, fp y_center, fp z_center, - fp mass, fp spin) -{ -CCTK_VInfo(CCTK_THORNSTRING, - " setting up Kerr/Kerr-Schild horizon"); -CCTK_VInfo(CCTK_THORNSTRING, - " mass=%g, spin=%g, posn=(%g,%g,%g)", - double(mass), double(spin), - double(x_center), double(y_center), double(z_center)); - -// horizon in Kerr-Schild is coordinate sphere $r = (1 + \sqrt{1-a^2}) m$ -const fp r = (1.0 + sqrt(1.0 - spin*spin)) * mass; - - for (int pn = 0 ; pn < ps.N_patches() ; ++pn) - { - patch& p = ps.ith_patch(pn); - - for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) - { - for (int isigma = p.min_isigma() ; - isigma <= p.max_isigma() ; - ++isigma) - { - const fp rho = p.rho_of_irho(irho); - const fp sigma = p.sigma_of_isigma(isigma); - - fp theta, phi; - p.theta_phi_of_rho_sigma(rho,sigma, theta,phi); - fp Kerr_x, Kerr_y, Kerr_z; - p.xyz_of_r_rho_sigma(r,rho,sigma, Kerr_x,Kerr_y,Kerr_z); - - const fp KS_x = Kerr_x - mass*spin*sin(theta)*sin(phi); - const fp KS_y = Kerr_y + mass*spin*sin(theta)*cos(phi); - const fp KS_z = Kerr_z; - p.ghosted_gridfn(ghosted_gfns::gfn__h,irho,isigma) - = jtutil::hypot3(KS_x, KS_y, KS_z); - } - } - } -} - } diff --git a/src/gr/gfn.hh b/src/gr/gfn.hh index 24c56fe..9f54b5c 100644 --- a/src/gr/gfn.hh +++ b/src/gr/gfn.hh @@ -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, - max_gfn = gfn__partial_H_wrt_partial_dd_h_22 // no comma + 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 7718796..bcaec78 100644 --- a/src/gr/horizon_Jacobian.cc +++ b/src/gr/horizon_Jacobian.cc @@ -3,9 +3,11 @@ // // <<>> // create_Jacobian - create a Jacobian matrix of the appropriate type -// horizon_Jacobian - top-level driver to compute the Jacobian +// +// horizon_Jacobian_SD - compute the Jacobian by symbolic differentiation /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian // +// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation #include #include @@ -53,12 +55,18 @@ void add_ghost_zone_Jacobian(const patch_system& ps, } //****************************************************************************** +//****************************************************************************** +//****************************************************************************** // // This function decodes the Jacobian_type parameter and creates the // appropriate type of Jacobian matrix object. // -Jacobian& create_Jacobian(const patch_system& ps, +// FIXME: the patch system shouldn't really have to be non-const, but +// the Jacobian constructors all require this to allow the +// linear solvers to directly update gridfns +// +Jacobian& create_Jacobian(patch_system& ps, const char Jacobian_type[]) { if (STRING_EQUAL(Jacobian_type, "dense matrix")) @@ -68,13 +76,15 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, Jacobian_type); /*NOTREACHED*/ } +//****************************************************************************** +//****************************************************************************** //****************************************************************************** // // This function computes the Jacobian matrix of the LHS function H(h) -// 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. +// 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. // // Inputs (angular gridfns, on ghosted grid): // h # shape of trial surface @@ -84,9 +94,10 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // Outputs: // The Jacobian matrix is stored in the Jacobian object Jac. // -void horizon_Jacobian(patch_system& ps, Jacobian& Jac) +void horizon_Jacobian_SD(patch_system& ps, + Jacobian& Jac) { -CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_SD"); ps.compute_synchronize_Jacobian(); Jac.zero(); @@ -189,7 +200,7 @@ Jac.zero(); 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, + (ps, xp, xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0, xm_irho, xm_isigma), II, xm_irho, xm_isigma, @@ -254,3 +265,91 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn } } } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function computes the Jacobian matrix of the LHS function H(h) +// by numerical perturbation of the H(h) function. The algorithm is as +// follows: +// +// evaluate H(h) +// array-copy H(h) to scratch array save_H +// for each point (y,JJ) +// { +// const fp save_h_y = h at y; +// h at y += perturbation_amplitude; +// evaluate H(h) +// 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) +// +void horizon_Jacobian_NP(patch_system& ps, + const struct cactus_grid_info& cgi, + const struct geometry_interpolator_info& ii, + Jacobian& Jac, + fp perturbation_amplitude) +{ +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian_NP"); + +// evaluate H(h) +jtutil::norm 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); + + for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn) + { + patch& yp = ps.ith_patch(ypn); + CCTK_VInfo(CCTK_THORNSTRING, " perturbing in %s patch", yp.name()); + + for (int y_irho = yp.min_irho() ; y_irho <= yp.max_irho() ; ++y_irho) + { + for (int y_isigma = yp.min_isigma() ; + y_isigma <= yp.max_isigma() ; + ++y_isigma) + { + const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma); + + const fp save_h_y = yp.ghosted_gridfn(ghosted_gfns::gfn__h, + y_irho,y_isigma); + yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma) + += perturbation_amplitude; + H_norms.reset(); + horizon_function(ps, cgi, ii, false, H_norms, false); + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) + { + patch& xp = ps.ith_patch(xpn); + + for (int x_irho = xp.min_irho() ; + x_irho <= xp.max_irho() ; + ++x_irho) + { + for (int x_isigma = xp.min_isigma() ; + x_isigma <= xp.max_isigma() ; + ++x_isigma) + { + const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma); + const fp new_xH = xp.gridfn(nominal_gfns::gfn__H, + x_irho,x_isigma); + Jac(II,JJ) = (new_xH - save_H[II]) / perturbation_amplitude; + } + } + } + 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; +} diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 82f5f31..9c6cb17 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -46,11 +46,15 @@ using jtutil::error_exit; // namespace { -void setup_xyz_posns(patch_system& ps); +void setup_xyz_posns(patch_system& ps, bool msg_flag); void interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& ii); -void compute_H(patch_system& ps, bool Jacobian_flag); + const struct geometry_interpolator_info& ii, + bool msg_flag); +void compute_H(patch_system& ps, + bool Jacobian_flag, + jtutil::norm& H_norms, + bool msg_flag); } //****************************************************************************** @@ -92,26 +96,30 @@ void compute_H(patch_system& ps, bool Jacobian_flag); // 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. // void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, const struct geometry_interpolator_info& ii, - bool Jacobian_flag) + bool Jacobian_flag, + jtutil::norm& H_norms, + bool msg_flag = true) { -CCTK_VInfo(CCTK_THORNSTRING, " horizon function"); +if (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); +setup_xyz_posns(ps, msg_flag); // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$ -interpolate_geometry(ps, cgi, ii); +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); +compute_H(ps, Jacobian_flag, H_norms, msg_flag); } //****************************************************************************** @@ -121,9 +129,11 @@ compute_H(ps, Jacobian_flag); // in the gridfns global_[xyz]. These will be used by interplate_geometry(). // namespace { -void setup_xyz_posns(patch_system& ps) +void setup_xyz_posns(patch_system& ps, bool msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " xyz positions and derivative coefficients"); +if (msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " xyz positions and derivative coefficients"); for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { @@ -188,10 +198,12 @@ CCTK_VInfo(CCTK_THORNSTRING, " xyz positions and derivative coefficients"); namespace { void interpolate_geometry(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii) + const struct geometry_interpolator_info& gii, + bool msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, - " interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid"); +if (msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid"); int status; @@ -313,8 +325,9 @@ if (first_time) first_time = false; // store derivative info in interpolator parameter table - CCTK_VInfo(CCTK_THORNSTRING, - " setting up derivative info for interpolator"); + if (msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " setting up derivative info for interpolator"); status = Util_TableSetIntArray(gii.param_table_handle, N_OUTPUT_ARRAYS, operand_indices, @@ -339,9 +352,10 @@ if (first_time) status); /*NOTREACHED*/ } -CCTK_VInfo(CCTK_THORNSTRING, - " calling interpolator (%d points)", - N_interp_points); +if (msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, + " calling interpolator (%d points)", + N_interp_points); status = CCTK_InterpLocalUniform(N_GRID_DIMS, gii.operator_handle, gii.param_table_handle, cgi.coord_origin, cgi.coord_delta, @@ -421,11 +435,16 @@ 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) +void compute_H(patch_system& ps, + bool Jacobian_flag, + jtutil::norm& H_norms, + bool msg_flag) { -CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); +if (msg_flag) + then CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { @@ -504,6 +523,9 @@ CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); #include "../gr.cg/horizon_function.c" } + // update running norms of H(h) function + H_norms.data(H); + if (Jacobian_flag) then { // partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn index 7def0d5..8e9745a 100644 --- a/src/gr/make.code.defn +++ b/src/gr/make.code.defn @@ -4,7 +4,8 @@ # Source files in this directory SRCS = driver.cc \ horizon_function.cc \ - horizon_Jacobian.cc + horizon_Jacobian.cc \ + Newton.cc # Subdirectories containing source files SUBDIRS = diff --git a/src/include/config.hh b/src/include/config.hh index dcbf2c9..6cf3693 100644 --- a/src/include/config.hh +++ b/src/include/config.hh @@ -8,7 +8,13 @@ typedef CCTK_REAL fp; +// Fortran 'integer' type +typedef CCTK_INT integer; + // FIXME: this assumes fp == C 'double' +// (CCTK_REAL_PRECISION_{4,8,16} are helpful, but not quite enough) +#undef FP_IS_FLOAT +#define FP_IS_DOUBLE #define FP_SCANF_FORMAT "%lf" // diff --git a/src/jtutil/README b/src/jtutil/README index 102e1b1..1b6ea5e 100644 --- a/src/jtutil/README +++ b/src/jtutil/README @@ -13,12 +13,15 @@ linear_map is a template class (templated on the floating-point type) representing a linear map between a contiguous interval of integers, and floating-point numbers. + cpm_map is a template class (templated on the floating-point type) representing a mapping from the integers to the integers, of the form i --> constant +/- i . (The name abbreviates "Constant Plus or Minus MAP".) + fuzzy + (declarations in "util.hh") is a template class (templated on the floating-point type) providing fuzzy arithmetic, to try to paper over some of the effects of floating-point rounding errors. For example, @@ -30,8 +33,20 @@ fuzzy and have the loop execute as one would naievly expect, even if rounding errors cause the final value of \code{x} to be 0.9999999999999999 or 1.000000000000001 or suchlike. + round + (declarations in "util.hh") is a template class (templated on the floating-point type) to do machine-independent rounding of floating point values +norm + (declarations in "util.hh") + is a template class (templated on the floating-point type) + to the compute 2-, rms-, and/or infinity-norms of a set of + numbers + +miscfp.cc + (declarations in "util.hh") + contains various misc floating-point routines + There are also a number of test drivers in the files test_*.cc . diff --git a/src/jtutil/miscfp.cc b/src/jtutil/miscfp.cc index 5ccda53..4263d2c 100644 --- a/src/jtutil/miscfp.cc +++ b/src/jtutil/miscfp.cc @@ -6,11 +6,19 @@ // 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 #include "stdc.h" #include "util.hh" +//****************************************************************************** +//****************************************************************************** //****************************************************************************** // @@ -113,3 +121,78 @@ if (! (fuzzy::GE(xx, xmin) && fuzzy::LE(xx, xmax)) ) return xx; } } // namespace jtutil:: + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function zeros the size-N array dst[]. +// +namespace jtutil + { +template + 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 + 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 + 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(int, float []); +template void array_zero(int, double[]); + +template void array_copy(int, const float [], float []); +template void array_copy(int, const double[], double[]); + +template void array_scale(int, float , float []); +template void array_scale(int, double, double[]); + } // namespace jtutil:: diff --git a/src/jtutil/util.hh b/src/jtutil/util.hh index 6873e62..89d6bde 100644 --- a/src/jtutil/util.hh +++ b/src/jtutil/util.hh @@ -1,13 +1,5 @@ // util.hh -- stuff for my C++ utility library // $Id$ -// -// jtutil::{how_many_in_range,is_even,is_odd} - misc integer manipulation -// jtutil::abs - absolute value template -// jtutil::pow* - raise floating point value to small-integer power -// jtutil::norm - compute 2-, rms-, and infinity-norms -// jtutil::fuzzy:: - fuzzy floating point comparisons and other operations -// jtutil::round:: - floating point rounding -// // // prerequisites: @@ -27,10 +19,14 @@ inline int is_even(int i) { return !(i & 0x1); } inline int is_odd (int i) { return (i & 0x1); } // -// absolute value template +// min/max/absolute value template // FIXME: puts this in std::, but older g++ versions // (which don't fully implement namespaces) get confused // +template + inline fp min(fp x, fp y) { return (x < y) ? x : y; } +template + inline fp max(fp x, fp y) { return (x > y) ? x : y; } template inline fp abs(fp x) { return (x > 0.0) ? x : -x; } @@ -65,6 +61,16 @@ double arctan_xy(double x, double y); // or error_exit() if no such value exists double modulo_reduce(double x, double xmod, double xmin, double xmax); +// +// simple BLAS-like array arithmetic routines +// +template + void array_zero(int N, fp dst[]); // dst[] = 0 +template + void array_copy(int N, const fp src[], fp dst[]); // dst[] = src[] +template + void array_scale(int N, fp alpha, fp dst[]); // dst[] *= alpha + // // more misc math stuff, valid only if has been #included // diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index dcc54c4..fa0d28d 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -207,8 +207,11 @@ void patch_system::create_patches(const struct patch_info patch_info_in[], fp delta_drho_dsigma) { CCTK_VInfo(CCTK_THORNSTRING, - "constructing %s patch system", + " constructing %s patch system", name_of_type(type())); +CCTK_VInfo(CCTK_THORNSTRING, + " at origin (%g,%g,%g)", + double(origin_x()), double(origin_y()), double(origin_z())); N_grid_points_ = 0; ghosted_N_grid_points_ = 0; for (int pn = 0 ; pn < N_patches() ; ++pn) @@ -216,7 +219,7 @@ ghosted_N_grid_points_ = 0; const struct patch_info& pi = patch_info_in[pn]; CCTK_VInfo(CCTK_THORNSTRING, - " constructing %s patch", + " constructing %s patch", pi.name); struct patch *p; @@ -322,13 +325,13 @@ const int ghosted_gfn_stride = ghosted_N_grid_points(); const int N_storage = gfn_stride * N_gridfns_in; const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in; -CCTK_VInfo(CCTK_THORNSTRING, "setting up gridfn storage"); +CCTK_VInfo(CCTK_THORNSTRING, " setting up gridfn storage"); CCTK_VInfo(CCTK_THORNSTRING, - " gfn=[%d,%d] ghosted_gfn=[%d,%d]", + " gfn=[%d,%d] ghosted_gfn=[%d,%d]", min_gfn_in, max_gfn_in, ghosted_min_gfn_in, ghosted_max_gfn_in); CCTK_VInfo(CCTK_THORNSTRING, - " N_grid_points()=%d ghosted_N_grid_points()=%d", + " N_grid_points()=%d ghosted_N_grid_points()=%d", N_grid_points(), ghosted_N_grid_points()); // storage arrays for all gridfns @@ -359,7 +362,7 @@ ghosted_gridfn_storage_ = new fp[ghosted_N_storage]; } CCTK_VInfo(CCTK_THORNSTRING, - " checking that storage is partitioned properly"); + " checking that storage is partitioned properly"); // check to make sure storage for distinct gridfns // forms a partition of the overall storage array @@ -498,7 +501,7 @@ void patch_system::interlink_full_sphere_patch_system (int N_overlap_points, int interp_handle, int interp_par_table_handle) { -CCTK_VInfo(CCTK_THORNSTRING, "interlinking full sphere patch system"); +CCTK_VInfo(CCTK_THORNSTRING, " interlinking full sphere patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -508,7 +511,7 @@ patch& my = ith_patch(patch_number_of_name("-y")); patch& mz = ith_patch(patch_number_of_name("-z")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(pz, mx, N_overlap_points); @@ -523,7 +526,7 @@ create_interpatch_ghost_zones(mz, mx, N_overlap_points); create_interpatch_ghost_zones(mz, my, N_overlap_points); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -574,7 +577,7 @@ void patch_system::interlink_plus_z_hemisphere_patch_system (int N_overlap_points, int interp_handle, int interp_par_table_handle) { -CCTK_VInfo(CCTK_THORNSTRING, "interlinking +z hemisphere patch system"); +CCTK_VInfo(CCTK_THORNSTRING, " interlinking +z hemisphere patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -583,7 +586,7 @@ patch& mx = ith_patch(patch_number_of_name("-x")); patch& my = ith_patch(patch_number_of_name("-y")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(pz, mx, N_overlap_points); @@ -598,7 +601,7 @@ mx.create_mirror_symmetry_ghost_zone(mx.min_rho_patch_edge()); my.create_mirror_symmetry_ghost_zone(my.min_rho_patch_edge()); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -637,7 +640,7 @@ void patch_system::interlink_plus_xy_quadrant_patch_system (int N_overlap_points, int interp_handle, int interp_par_table_handle) { -CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xy quadrant patch system"); +CCTK_VInfo(CCTK_THORNSTRING, " interlinking +xy quadrant patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -645,7 +648,7 @@ patch& py = ith_patch(patch_number_of_name("+y")); patch& mz = ith_patch(patch_number_of_name("-z")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(px, py, N_overlap_points); @@ -662,7 +665,7 @@ create_periodic_symmetry_ghost_zones(mz.max_rho_patch_edge(), true); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -692,7 +695,7 @@ void patch_system::interlink_plus_xz_quadrant_patch_system (int N_overlap_points, int interp_handle, int interp_par_table_handle) { -CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xz quadrant patch system"); +CCTK_VInfo(CCTK_THORNSTRING, " interlinking +xz quadrant patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); @@ -700,7 +703,7 @@ patch& py = ith_patch(patch_number_of_name("+y")); patch& my = ith_patch(patch_number_of_name("-y")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(pz, my, N_overlap_points); @@ -717,7 +720,7 @@ create_periodic_symmetry_ghost_zones(py.max_sigma_patch_edge(), false); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); @@ -747,14 +750,14 @@ void patch_system::interlink_plus_xyz_octant_patch_system (int N_overlap_points, int interp_handle, int interp_par_table_handle) { -CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xyz octant patch system"); +CCTK_VInfo(CCTK_THORNSTRING, " interlinking +xyz octant patch system"); patch& pz = ith_patch(patch_number_of_name("+z")); patch& px = ith_patch(patch_number_of_name("+x")); patch& py = ith_patch(patch_number_of_name("+y")); // create the ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); +CCTK_VInfo(CCTK_THORNSTRING, " creating ghost zones"); create_interpatch_ghost_zones(pz, px, N_overlap_points); create_interpatch_ghost_zones(pz, py, N_overlap_points); create_interpatch_ghost_zones(px, py, N_overlap_points); @@ -768,7 +771,7 @@ create_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(), true); // finish setting up the interpatch ghost zones -CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); +CCTK_VInfo(CCTK_THORNSTRING, " finishing interpatch setup"); finish_interpatch_setup(pz, px, N_overlap_points, interp_handle, interp_par_table_handle); -- cgit v1.2.3