switch to molecule-at-a-time interface to access synchronize_Jacobian()

results


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

1 files changed, 18 insertions, 17 deletions

diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index d4d11f9..0dd7a31 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -269,17 +269,11 @@ CCTK_VInfo(CCTK_THORNSTRING,
 	++x_isigma)
 	{
 	//
-	// compute the Jacobian row for this grid point, i.e.
+	// compute the main Jacobian terms for this grid point, i.e.
 	//	  partial H(this point x, Jacobian row II)
 	//	---------------------------------------------
 	//	partial h(other points y, Jacobian column JJ)
 	//
-	// FIXME FIXME: we still have to take into account the
-	//		position-dependence of the coefficients,
-	//		cf the difference between J[3H(h)] and J[2H(h)];
-	//		here we're sort of computing the former, but
-	//		Newton's method really wants the latter
-	//
 
 	// Jacobian row index
 	const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);
@@ -397,15 +391,24 @@ 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);
 
+// FIXME: this won't change from one call to another
+//        ==> it would be more efficient to reuse the same buffer
+//            across multiple calls on this function
+int global_min_m, global_max_m;
+ps.synchronize_Jacobian_global_minmax_y_ipar_m(global_min_m, global_max_m);
+jtutil::array1d<fp> Jacobian_buffer(global_min_m, global_max_m);
+
 // 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);
+const patch&      ymp = ps.synchronize_Jacobian_y_patch(xmgz, xm_iperp,xm_ipar);
+const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz, xm_iperp,xm_ipar);
+const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp,xm_ipar);
+const int ym_ipar_posn
+	       = ps.synchronize_Jacobian_y_ipar_posn(xmgz, xm_iperp,xm_ipar);
+int min_ym_ipar_m, max_ym_ipar_m;
+ps.synchronize_Jacobian_minmax_y_ipar_m(xmgz, xm_iperp, xm_ipar,
+					min_ym_ipar_m, max_ym_ipar_m);
+ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar, Jacobian_buffer);
 
 // add the Jacobian contributions from the ym points
 	for (int ym_ipar_m = min_ym_ipar_m ;
@@ -416,9 +419,7 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn
 	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(x_II,JJ) += mol*sync_Jac;
+	Jac(x_II,JJ) += mol*Jacobian_buffer(ym_ipar_m);
 	}
 }
 	  }