restore "ellipsoid" initial guess

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@645 f88db872-0e4f-0410-b76b-b9085cfa78c5

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;