diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-18 19:36:42 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2002-07-18 19:36:42 +0000 |
commit | f6165749357e1a46fd09ac058c6f98e463ef2bdf (patch) | |
tree | 887d9073641d78ca2f7e0afaac7d335c9ab6497b /src/gr/driver.cc | |
parent | 3ffb8c841a8b727e98615c538d6d89295c734094 (diff) |
restore "ellipsoid" initial guess
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@645 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat (limited to 'src/gr/driver.cc')
-rw-r--r-- | src/gr/driver.cc | 127 |
1 files changed, 115 insertions, 12 deletions
diff --git a/src/gr/driver.cc b/src/gr/driver.cc index 5031e41..0783f43 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -3,7 +3,8 @@ // // <<<prototypes for functions local to this file>>> // AHFinderDirect_driver - top-level driver -/// setup_Kerr_KerrSchild_initial_guess - set up Kerr/Kerr-Schild initial guess +/// setup_ellipsoid - setup up a coordiante ellipsoid in h +/// setup_Kerr_KerrSchild_horizon - set up Kerr/Kerr-Schild horizon in h // #include <stdio.h> @@ -45,9 +46,12 @@ using jtutil::error_exit; // namespace { +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 mass, fp spin, - fp xposn, fp yposn, fp zposn); + fp x_center, fp y_center, fp z_center, + fp mass, fp spin); }; //****************************************************************************** @@ -156,14 +160,22 @@ 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")) + 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-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, - initial_guess__Kerr_KerrSchild__xposn, - initial_guess__Kerr_KerrSchild__yposn, - initial_guess__Kerr_KerrSchild__zposn); + initial_guess__Kerr_KerrSchild__spin); ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h, true, ghosted_gfns::gfn__h, "h.dat", @@ -184,7 +196,7 @@ if (STRING_EQUAL(method, "horizon")) true, ghosted_gfns::gfn__h, "H.dat"); } -if (STRING_EQUAL(method, "horizon Jacobian")) +else if (STRING_EQUAL(method, "horizon Jacobian")) then { Jacobian& Jac = create_Jacobian(ps, Jacobian_type); horizon_function(ps, cgi, gii, true); @@ -198,6 +210,95 @@ 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) +// +// 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 +// +// to solve this, we introduce intermediate variables +// AU = A - U +// BV = B - V +// CW = C - W +// +namespace { +void setup_ellipsoid(patch_system& ps, + fp x_center, fp y_center, fp z_center, + fp x_radius, fp y_radius, fp z_radius) +{ +CCTK_VInfo(CCTK_THORNSTRING, + " h = ellipsoid: center=(%g,%g,%g)", + double(x_center), double(y_center), double(z_center)); +CCTK_VInfo(CCTK_THORNSTRING, + " radius=(%g,%g,%g)", + double(x_radius), double(y_radius), double(z_radius)); + + 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 xcos, ycos, zcos; + p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos); + + // set up variables used by Maple-generated code + const fp AU = x_center - ps.origin_x(); + const fp BV = y_center - ps.origin_y(); + const fp CW = z_center - ps.origin_z(); + const fp a = x_radius; + const fp b = y_radius; + const fp c = z_radius; + + // 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*/ + + // r = horizon radius at this grid point + p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r; + } + } + } +} + } + +//****************************************************************************** + +// // This function sets up the horizon of a Kerr black hole in Kerr-Schild // coordinates, on the nominal grid, in the h gridfn. // @@ -205,20 +306,22 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // 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. -// [xyz]_posn = The position of the Kerr black hole. +// +// FIXME FIXME: at present [xyz}_center are ignored. // namespace { void setup_Kerr_KerrSchild_horizon(patch_system& ps, - fp mass, fp spin, - fp xposn, fp yposn, fp zposn) + 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(xposn), double(yposn), double(zposn)); + 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; |