more code to compute the Jacobian matrix

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

3 files changed, 123 insertions, 52 deletions

diff --git a/src/gr/Jacobian.hh b/src/gr/Jacobian.hh
index 626d080..e49187a 100644
--- a/src/gr/Jacobian.hh
+++ b/src/gr/Jacobian.hh
@@ -17,6 +17,10 @@
 class	Jacobian
 	{
 public:
+	// basic meta-info
+	patch_system& my_patch_system();
+
+public:
 	// access a matrix element
 	virtual const fp& operator()(int II, int JJ) const;	// rvalue
 	virtual       fp& operator()(int II, int JJ) const;	// lvalue
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index 40fada6..39db258 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -43,7 +43,11 @@ using jtutil::error_exit;
 //
 
 namespace {
-	  };
+void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
+				    int x_irho, int x_isigma,
+				    int xm_irho, int xm_isigma,
+				    fp mol);
+	  }
 
 //******************************************************************************
 //******************************************************************************
@@ -67,7 +71,8 @@ void horizon_Jacobian(const patch_system& ps,
 {
 CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian");
 
-Jac.zero(II);
+ps.compute_synchronize_Jacobian();
+Jac.zero();
 
     for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
     {
@@ -92,9 +97,6 @@ Jac.zero(II);
 	//		Newton's method really wants the latter
 	//
 
-	// Jacobian row index for this point
-	const int II = ps.gpn_of_patch_irho_isigma(xp,x_irho,x_isigma);
-
 	// Jacobian coefficients for this point
 	const fp Jacobian_coeff_rho
 	   = xp.gridfn(gfn__partial_H_wrt_partial_d_h_1, x_irho,x_isigma);
@@ -112,50 +114,115 @@ Jac.zero(II);
 		 m_irho <= xp.molecule_max_m() ;
 		 ++m_irho)
 	    {
-	    const fp Jac_rho = Jacobian_coeff_rho * xp.partial_rho_coeff(m_rho);
+	    const int xm_irho = x_irho + m_irho;
+	    const fp Jac_rho     = Jacobian_coeff_rho
+				   * xp.partial_rho_coeff(m_rho);
 	    const fp Jac_rho_rho = Jacobian_coeff_rho_rho
 				   * xp.partial_rho_rho_coeff(m_rho);
-	    const fp Jac_sum = Jac_rho + Jac_rho_rho;
-	    if (xp.is_in_nominal_grid(x_irho+m_irho, x_isigma))
-	       then {
-		    const int JJ
-			= ps.gpn_of_patch_irho_isigma(xp, x_irho+m_irho,
-							  x_isigma);
-		    Jac.add_to_element(II,JJ, Jac_sum);
-		    }
-	       else {
-		    const patch_edge& xe = (m_irho < 0)
-					   ? xp.min_rho_patch_edge()
-					   : xp.max_rho_patch_edge()
-		    const ghost_zone& gz = xp.ghost_zone_on_edge(xe);
-		    const int x_iperp = xe.iperp_of_irho_isigma(x_irho,
-								x_isigma);
-		    const int x_ipar = xe.ipar_of_irho_isigma(x_irho, x_isigma);
-		    gz.compute_Jacobian(gfn__h, gfn__h);
-		    const patch&      yp = gz.Jacobian_y_patch();
-		    const patch_edge& ye = gz.Jacobian_y_edge();
-		    const int min_y_ipar_m = gz.Jacobian_min_y_ipar_m();
-		    const int max_y_ipar_m = gz.Jacobian_max_y_ipar_m();
-		    const int y_iperp = gz.Jacobian_y_iperp(x_iperp);
-		    const int y_ipar_posn
-			= gz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
-			for (int y_ipar_m = min_y_ipar_m ;
-			     y_ipar_m <= max_y_ipar_m ;
-			     ++y_ipar_m)
-			{
-			const int y_ipar = y_ipar_posn + y_ipar_m;
-			const int y_irho
-			   = ye.irho_of_iperp_ipar(y_iperp, y_ipar);
-			const int y_isigma
-			   = ye.isigma_of_iperp_ipar(y_iperp, y_ipar);
-			const int JJ
-			   = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma);
-			const fp gz_Jac = gz.Jacobian(x_iperp,x_ipar, y_ipar_m);
-			Jac.add_to_element(II,JJ, Jac_sum*gz_Jac);
-			}
-			}
-		}
+	    add_molecule_point_to_Jacobian(ps, xp,
+					   x_irho, x_isigma,
+					   xm_irho, x_isigma,
+					   Jac_rho + Jac_rho_rho);
+	    }
+
+	// partial_sigma, partial_sigma_sigma
+	    for (int m_isigma = xp.molecule_min_m() ;
+		 m_isigma <= xp.molecule_max_m() ;
+		 ++m_isigma)
+	    {
+	    const int xm_isigma = x_isigma + m_isigma;
+	    const fp Jac_sigma       = Jacobian_coeff_sigma
+				       * xp.partial_sigma_coeff(m_sigma);
+	    const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma
+				       * xp.partial_sigma_sigma_coeff(m_sigma);
+	    add_molecule_point_to_Jacobian(ps, xp,
+					   x_irho, x_isigma,
+					   x_irho, xm_isigma,
+					   Jac_sigma + Jac_sigma_sigma);
+	    }
+
+	// partial_rho_sigma
+	    for (int m_irho = xp.molecule_min_m() ;
+		 m_irho <= xp.molecule_max_m() ;
+		 ++m_irho)
+	    {
+	    for (int m_isigma = xp.molecule_min_m() ;
+		 m_isigma <= xp.molecule_max_m() ;
+		 ++m_isigma)
+	    {
+	    const int xm_irho   = x_irho   + m_irho;
+	    const int xm_isigma = x_isigma + m_isigma;
+	    const fp Jac_rho_sigma
+		= Jacobian_coeff_rho_sigma
+		  * xp.partial_rho_sigma_coeff(m_rho, m_sigma);
+	    add_molecule_point_to_Jacobian(ps, xp,
+					   x_irho, x_isigma,
+					   xm_irho, xm_isigma,
+					   Jac_rho_sigma);
+	    }
+	    }
+
+	}
 	}
+    }
+}
+
+//******************************************************************************
+
+//
+// This function adds all the elements for a single molecule point to
+// a Jacobian matrix, including any contributions from other patches if
+// the specified point lies in a ghost zone.
+//
+// Arguments:
+// ps = The patch system.
+// xp = The patch containing the center point of the molecule.
+// x_(irho,isigma) = The coordinates of the center point of the molecule.
+// xm_(irho,isigma) = The coordinates of the x+m point of the molecule.
+// mol = The molecule coefficient.
+//
+namespace {
+void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
+				    int x_irho, int x_isigma,
+				    int xm_irho, int xm_isigma,
+				    fp mol)
+{
+// Jacobian row
+const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);
+
+if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
+   then {
+	const int JJ = ps.gpn_of_patch_irho_isigma(xp, xm_irho, xm_isigma);
+	Jac.add_to_element(II,JJ, Jac_coeff);
+	}
+   else {
+	const ghost_zone& xmgz
+		= xp.ghost_zone_containing_point(xm_irho, xm_isigma);
+	const patch_edge& xme = xmgz.my_edge();
+	Const int xm_iperp = xme.iperp_of_irho_isigma(xm_irho, xm_isigma);
+	const int xm_ipar  = xme. ipar_of_irho_isigma(xm_irho, xm_isigma);
+
+	// on what other points ym does this molecule point xm depend
+	// via the patch_system::synchronize() operation?
+	const patch&      ymp = ps.synchronize_Jacobian_y_patch(xmgz);
+	const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz);
+	const int min_ym_ipar_m = ps.synchronize_Jacobian_min_y_ipar_m(xmgz);
+	const int max_ym_ipar_m = ps.synchronize_Jacobian_max_y_ipar_m(xmgz);
+	const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp);
+	const int ym_ipar_posn
+	   = ps.synchronize_Jacobian_y_ipar_posn(xmgz, xm_iperp, xm_ipar);
+	    for (int ym_ipar_m = min_ym_ipar_m ;
+		 ym_ipar_m <= max_ym_ipar_m ;
+		 ++ym_ipar_m)
+	    {
+	    const int ym_ipar = ym_ipar_posn + ym_ipar_m;
+	    const int ym_irho   = yme.  irho_of_iperp_ipar(ym_iperp,ym_ipar);
+	    const int ym_isigma = yme.isigma_of_iperp_ipar(ym_iperp,ym_ipar);
+	    const int JJ = ps.gpn_of_patch_irho_isigma(ymp, ym_irho, ym_isigma);
+	    const fp sync_Jac = ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar,
+							ym_ipar_m);
+	    Jac.add_to_element(II,JJ, mol*sync_Jac);
+	    }
 	}
-     }
 }
+	  }
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 745760b..400516b 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -50,7 +50,7 @@ void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_interpolator_info& ii);
 void compute_H(patch_system& ps);
-	  };
+	  }
 
 //******************************************************************************
 //******************************************************************************
@@ -100,7 +100,7 @@ void horizon_function(patch_system& ps,
 CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 
 // fill in values of ghosted gridfns in ghost zones
-ps.synchronize_ghost_zones(ghosted_gfns::gfn__h, ghosted_gfns::gfn__h);
+ps.synchronize(ghosted_gfns::gfn__h, ghosted_gfns::gfn__h);
 
 // set up xyz positions of grid points
 setup_xyz_posns(ps);
@@ -379,11 +379,11 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
 	assert(out_of_range_pt >= 0);
 	assert(out_of_range_pt < ps.N_grid_points());
 	const double global_x = ps.gridfn_data(nominal_gfns::gfn__global_x)
-				[out_of_range_pt];
+				   [out_of_range_pt];
 	const double global_y = ps.gridfn_data(nominal_gfns::gfn__global_y)
-				[out_of_range_pt];
+				   [out_of_range_pt];
 	const double global_z = ps.gridfn_data(nominal_gfns::gfn__global_z)
-				[out_of_range_pt];
+				   [out_of_range_pt];
 
 	assert(out_of_range_axis >= 0);
 	assert(out_of_range_axis < N_GRID_DIMS);