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/gr/driver.cc | 192 ++++++++++++++++++++++++++++++++----------------------- 1 file changed, 112 insertions(+), 80 deletions(-) (limited to 'src/gr/driver.cc') 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); - } - } - } -} - } -- cgit v1.2.3