diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-18 17:42:03 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-18 17:42:03 +0000 |
commit | 07784c803efe3a065851d5e0e273cac82c12f0cb (patch) | |
tree | 60bff0605f54ecf0f9906df412947acbb6bf0d5b /src/gr/driver.cc | |
parent | 695c3f5b78149e18eeb0cc74f5ae971f5ff9cd32 (diff) |
move Jacobian.* to ../elliptic/
update other files to compute more of Jacobian coeffs
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@633 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src/gr/driver.cc')
-rw-r--r-- | src/gr/driver.cc | 170 |
1 files changed, 61 insertions, 109 deletions
diff --git a/src/gr/driver.cc b/src/gr/driver.cc index 3982be2..0c04dfe 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -3,7 +3,7 @@ // // <<<prototypes for functions local to this file>>> // AHFinderDirect_driver - top-level driver -/// setup_ellipsoid_initial_guess - set up ellipsoid in h gridfn +/// setup_Kerr_KerrSchild_initial_guess - set up Kerr/Kerr-Schild initial guess // #include <stdio.h> @@ -32,6 +32,7 @@ using jtutil::error_exit; #include "../util/patch_interp.hh" #include "../util/ghost_zone.hh" #include "../util/patch_system.hh" +#include "../util/Jacobian.hh" #include "gfn.hh" #include "AHFinderDirect.hh" @@ -43,10 +44,9 @@ using jtutil::error_exit; // namespace { -void setup_ellipsoid_initial_guess - (patch_system& ps, - fp global_center_x, fp global_center_y, fp global_center_z, - fp radius_x, fp radius_y, fp radius_z); +void setup_Kerr_KerrSchild_horizon(patch_system& ps, + fp mass, fp spin, + fp xposn, fp yposn, fp zposn); }; //****************************************************************************** @@ -113,42 +113,22 @@ cgi.coord_delta[2] = cctk_delta_space[2]; cgi.gridfn_dims[0] = cctk_lsh[0]; cgi.gridfn_dims[1] = cctk_lsh[1]; cgi.gridfn_dims[2] = cctk_lsh[2]; -cgi.g_dd_11_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gxx") - ); -cgi.g_dd_12_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gxy") - ); -cgi.g_dd_13_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gxz") - ); -cgi.g_dd_22_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gyy") - ); -cgi.g_dd_23_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gyz") - ); -cgi.g_dd_33_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::gzz") - ); -cgi.K_dd_11_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kxx") - ); -cgi.K_dd_12_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kxy") - ); -cgi.K_dd_13_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kxz") - ); -cgi.K_dd_22_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kyy") - ); -cgi.K_dd_23_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kyz") - ); -cgi.K_dd_33_data = static_cast<const fp *>( - CCTK_VarDataPtr(cctkGH, 0, "einstein::kzz") - ); +// n.b. The cgi.[gK]_dd_??_data are actually const fp * pointers, +// since we won't modify the 3-D gridfn data! But static_cast<...> +// won't change const modifiers, so we just cast to fp* and let +// the assignment take care of the const part... +cgi.g_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxx")); +cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxy")); +cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxz")); +cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyy")); +cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyz")); +cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gzz")); +cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxx")); +cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxy")); +cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxz")); +cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyy")); +cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyz")); +cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kzz")); // @@ -175,15 +155,14 @@ if (STRING_EQUAL(initial_guess_method, "read from file")) initial_guess__read_from_file__file_name, false); // no ghost zones } -else if (STRING_EQUAL(initial_guess_method, "ellipsoid")) +else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild")) then { - setup_ellipsoid_initial_guess(ps, - initial_guess__ellipsoid__center_global_x, - initial_guess__ellipsoid__center_global_y, - initial_guess__ellipsoid__center_global_z, - initial_guess__ellipsoid__radius_x, - initial_guess__ellipsoid__radius_y, - initial_guess__ellipsoid__radius_z); + setup_Kerr_KerrSchild_horizon(ps, + initial_guess__Kerr_KerrSchild__mass, + initial_guess__Kerr_KerrSchild__spin, + initial_guess__Kerr_KerrSchild__xposn, + initial_guess__Kerr_KerrSchild__yposn, + initial_guess__Kerr_KerrSchild__zposn); ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, true, ghosted_gfns::gfn__h, "h.dat", @@ -204,6 +183,12 @@ if (STRING_EQUAL(method, "horizon")) true, ghosted_gfns::gfn__h, "H.dat"); } +if (STRING_EQUAL(method, "horizon Jacobian")) + then { + Jacobian& Jac = create_Jacobian(ps, Jacobian_type); + horizon_function(ps, cgi, gii, true); + horizon_Jacobian(ps, Jac); + } else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, "unknown method=\"%s\"!", method); /*NOTREACHED*/ @@ -212,38 +197,30 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, //****************************************************************************** // -// This function sets up an ellipsoid in the gridfn h, using the -// formulas in "ellipsoid.maple" and the Maple-generated C code in -// "ellipsoid.c": -// -// ellipsoid has center (A,B,C), radius (a,b,c) -// angular coordinate system has center (U,V,W) -// -// direction cosines wrt angular coordinate center are (xcos,ycos,zcos) -// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) +// This function sets up the horizon of a Kerr black hole in Kerr-Schild +// coordinates, on the nominal grid, in the h gridfn. // -// then the equation of the ellipsoid is -// (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 -// ----------------- + ---------------- + ----------------- = 1 -// a^2 b^2 c^2 +// 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. // -// to solve this, we introduce intermediate variables -// AU = A - U -// BV = B - V -// CW = C - W +// Arguments: +// (mass,spin) = Describe the Kerr black hole. +// [xyz]_posn = The position of the Kerr black hole. // namespace { -void setup_ellipsoid_initial_guess - (patch_system& ps, - fp global_center_x, fp global_center_y, fp global_center_z, - fp radius_x, fp radius_y, fp radius_z) +void setup_Kerr_KerrSchild_horizon(patch_system& ps, + fp mass, fp spin, + fp xposn, fp yposn, fp zposn) { CCTK_VInfo(CCTK_THORNSTRING, - " h = ellipsoid: global_center=(%g,%g,%g)", - global_center_x, global_center_y, global_center_z); + " setting up Kerr/Kerr-Schild horizon"); CCTK_VInfo(CCTK_THORNSTRING, - " radius=(%g,%g,%g)", - radius_x, radius_y, radius_z); + " mass=%g, spin=%g, posn=(%g,%g,%g)", + double(mass), double(spin), + double(xposn), double(yposn), double(zposn)); + +// 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) { @@ -255,44 +232,19 @@ CCTK_VInfo(CCTK_THORNSTRING, isigma <= p.max_isigma() ; ++isigma) { - const fp rho = p.rho_of_irho(irho); + const fp rho = p.rho_of_irho(irho); const fp sigma = p.sigma_of_isigma(isigma); - fp xcos, ycos, zcos; - p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos); - - // set up variables used by Maple-generated code - const fp AU = global_center_x - ps.origin_x(); - const fp BV = global_center_y - ps.origin_y(); - const fp CW = global_center_z - ps.origin_z(); - const fp a = radius_x; - const fp b = radius_y; - const fp c = radius_z; - - // compute the solutions r_plus and r_minus - fp r_plus, r_minus; - #include "ellipsoid.c" - // exactly one of the solutions (call it r) should be positive - fp r; - if ((r_plus > 0.0) && (r_minus < 0.0)) - then r = r_plus; - else if ((r_plus < 0.0) && (r_minus > 0.0)) - then r = r_minus; - else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "\n" -" expected exactly one r>0 solution, 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" - , - p.name(), irho, isigma, - double(rho), double(sigma), - double(xcos), double(ycos), double(zcos), - double(r_plus), double(r_minus)); - /*NOTREACHED*/ + 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); - // r = horizon radius at this grid point - p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r; + const fp KS_x = Kerr_x - spin*sin(theta)*sin(phi); + const fp KS_y = Kerr_y + 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); } } } |