add support (not tested yet :)

for CactusEinstein/StaticConformal conformal metric


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

3 files changed, 297 insertions, 60 deletions

diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh
index 5404f46..f984f4f 100644
--- a/src/gr/gfns.hh
+++ b/src/gr/gfns.hh
@@ -15,22 +15,21 @@ enum	{
 // nominal gridfns
 enum	{
 	nominal_min_gfn = 1,
-	// n.b. no macros for next 3 gridfns in "cg.hh"
-	gfn__global_x = nominal_min_gfn,
-	gfn__global_y,
-	gfn__global_z,
+
+	//
+	// most of these gridfns have access macros in "cg.hh";
+	// the ones that don't are marked explicitly
+	//
+	gfn__global_x = nominal_min_gfn,	// no access macro
+	gfn__global_y,				// no access macro
+	gfn__global_z,				// no access macro
+
 	gfn__g_dd_11,
 	gfn__g_dd_12,
 	gfn__g_dd_13,
 	gfn__g_dd_22,
 	gfn__g_dd_23,
 	gfn__g_dd_33,
-	gfn__K_dd_11,
-	gfn__K_dd_12,
-	gfn__K_dd_13,
-	gfn__K_dd_22,
-	gfn__K_dd_23,
-	gfn__K_dd_33,
 	gfn__partial_d_g_dd_111,
 	gfn__partial_d_g_dd_112,
 	gfn__partial_d_g_dd_113,
@@ -49,6 +48,18 @@ enum	{
 	gfn__partial_d_g_dd_322,
 	gfn__partial_d_g_dd_323,
 	gfn__partial_d_g_dd_333,
+	gfn__K_dd_11,
+	gfn__K_dd_12,
+	gfn__K_dd_13,
+	gfn__K_dd_22,
+	gfn__K_dd_23,
+	gfn__K_dd_33,
+
+	gfn__psi,				// no access macro
+	gfn__partial_d_psi_1,			// no access macro
+	gfn__partial_d_psi_2,			// no access macro
+	gfn__partial_d_psi_3,			// no access macro
+
 	gfn__H,
 	gfn__partial_H_wrt_partial_d_h_1,
 	gfn__partial_H_wrt_partial_d_h_2,
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 53f88c4..beb1866 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -41,6 +41,12 @@ struct	cactus_grid_info
 	//	i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
 	CCTK_INT gridfn_dims[N_GRID_DIMS];
 
+	// this describes the semantics of the Cactus gij gridfns:
+	// true ==> the Cactus gij are conformal values as per
+	//          CactusEinstein/StaticConformal and the  psi  gridfn
+	// false ==> the Cactus gij are the physical metric
+	bool Cactus_conformal_metric;
+
 	// --> Cactus gridfn data for grid posn (0,0,0)
 	const fp* g_dd_11_data;
 	const fp* g_dd_12_data;
@@ -54,6 +60,8 @@ struct	cactus_grid_info
 	const fp* K_dd_22_data;
 	const fp* K_dd_23_data;
 	const fp* K_dd_33_data;
+
+	const fp* psi_data;		// NULL if !Cactus_conformal_metric
 	};
 
 //
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 3804be2..9b1e0a6 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -5,6 +5,7 @@
 // horizon_function - top-level driver
 /// setup_xyz_posns - setup global xyz posns of grid points
 /// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid
+/// convert_conformal_to_physical - convert conformal gij to physical
 /// compute_H - compute H(h) given earlier setup
 ///
 
@@ -53,6 +54,7 @@ bool interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
 			  bool print_msg_flag);
+void convert_conformal_to_physical(patch_system& ps, bool print_msg_flag);
 void Schwarzschild_EF_geometry(patch_system& ps,
 			       const struct cactus_grid_info& cgi,
 			       const struct geometry_info& gi,
@@ -135,6 +137,8 @@ switch	(gi.geometry_method)
 case geometry__interpolate_from_Cactus_grid:
 	if (! interpolate_geometry(ps, cgi, gi, print_msg_flag))
 	   then return false;				// *** ERROR RETURN ***
+	if (cgi.Cactus_conformal_metric)
+	   then convert_conformal_to_physical(ps, print_msg_flag);
 	break;
 case geometry__Schwarzschild_EF:
 	Schwarzschild_EF_geometry(ps, gi, print_msg_flag);
@@ -200,34 +204,45 @@ if (print_msg_flag)
 //******************************************************************************
 
 //
-// This function interpolates the Cactus gridfns $g_{ij}$ and $K_{ij}$
-// to determine $g_{ij}$, $K_{ij}$, and the spatial derivatives
-// $\partial_k g_{ij}$, on the trial horizon surface position given by h.
-//
-// Inputs (angular gridfns, on ghosted grid):
-// ... defined on ghosted grid
-// ... only values on nominal grid are actually used as input
-//	h				# shape of trial surface
+// This function interpolates the Cactus gridfns
+//	gxx...gzz
+//	kxx...kzz
+//	psi			# optional
+// to determine the nominal-grid angular gridfns
+//	g_dd_ij
+//	partial_d_g_dd_kij
+//	K_dd_ij
+//	psi			# optional
+//	partial_d_psi_k		# optional
+// at the nominal-grid trial horizon surface positions given by the
+// global_(x,y,z) angular gridfns.  The psi interpolation is done if and
+// only if the cgi.Cactus_conformal_metric flag is set.  Note that this
+// function ignores the physical-vs-conformal semantics of the gridfns;
+// it just interpolates and takes derivatives of the stored gridfn values.
 //
 // Inputs (angular gridfns, all on the nominal grid):
 //	global_[xyz]			# xyz positions of grid points
 //
 // Inputs (Cactus 3-D gridfns):
 //	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
+//					# (may be either physical or conformal)
 //	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
+//	psi				# optional conformal factor $\psi$
 //
 // Outputs (angular gridfns, all on the nominal grid):
-//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
-//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
-//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
+//	g_dd_{11,12,13,22,23,33}		# $\stored{g}_{ij}$
+//	K_dd_{11,12,13,22,23,33}		# $K_{ij}$
+//	partial_d_g_dd_[123]{11,12,13,22,23,33}	# $\partial_k \stored{g}_{ij}$
+//	psi					# (optional) $\psi$
+//	partial_d_psi_[123]			# (optional) $\partial_k \psi$
 //
-// On the first call this function also modifies the interpolator
-// parameter table.
+// This function may also modify the interpolator parameter table.
 //
 // Results:
 // This function returns true for a successful computation, or false
 // if the computation failed because the geometry interpolation would
 // need data outside the Cactus grid, or data from an excised region.
+//
 // FIXME: excision isn't implemented yet :(
 //
 namespace {
@@ -236,12 +251,32 @@ bool interpolate_geometry(patch_system& ps,
 			  const struct geometry_info& gi,
 			  bool print_msg_flag)
 {
+//
+// Implementation Notes:
+//
+// To handle the optional interpolation of psi, we set up all the data
+// type and pointer arrays to include psi, but with the psi entries at
+// the end.  We then choose the array sizes passed to the interpolator
+// to either include or exclude the psi entries as appropriate.
+//
+// We remember whether or not psi was interpolated on the previous call,
+// and only modify the interpolator parameter table if this changes (and
+// on our first call).
+//
+
+const bool psi_flag = cgi.Cactus_conformal_metric;
+
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "      interpolating g_ij and Kij from Cactus 3-D grid");
+		   "      interpolating %s from Cactus 3-D grid",
+		   psi_flag ? "{g_ij, K_ij, psi}" : "{g_ij, K_ij}");
 
 int status;
 
+
+//
+// ***** interpolation points *****
+//
 const int N_interp_points = ps.N_grid_points();
 const int interp_coords_type_code = CCTK_VARIABLE_REAL;
 const void* const interp_coords[N_GRID_DIMS]
@@ -251,8 +286,12 @@ const void* const interp_coords[N_GRID_DIMS]
     static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_z)),
     };
 
-const int N_INPUT_ARRAYS = 12;
-const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
+
+//
+// ***** input arrays *****
+//
+
+const CCTK_INT input_array_type_codes[]
 	= {
 	  // g_ij
 	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
@@ -262,8 +301,10 @@ const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
 	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 						  CCTK_VARIABLE_REAL,
+	  // psi
+	  CCTK_VARIABLE_REAL,
 	  };
-const void* const input_arrays[N_INPUT_ARRAYS]
+const void* const input_arrays[]
 	= {
 	  static_cast<const void*>(cgi.g_dd_11_data),
 	  static_cast<const void*>(cgi.g_dd_12_data),
@@ -277,11 +318,21 @@ const void* const input_arrays[N_INPUT_ARRAYS]
 	  static_cast<const void*>(cgi.K_dd_22_data),
 	  static_cast<const void*>(cgi.K_dd_23_data),
 	  static_cast<const void*>(cgi.K_dd_33_data),
+	  static_cast<const void*>(cgi.psi_data),
 	  };
 
-const int N_OUTPUT_ARRAYS = 30;
+const int N_input_arrays_for_psi = 1;
+const int N_input_arrays_dim = sizeof(input_arrays) / sizeof(input_arrays[0]);
+const int N_input_arrays_use
+	= psi_flag ? N_input_arrays_dim
+		   : N_input_arrays_dim - N_input_arrays_for_psi;
 
-const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
+
+//
+// ***** output arrays *****
+//
+
+const CCTK_INT output_array_type_codes[]
 	= {
  // g_ij             partial_x g_ij      partial_y g_ij      partial_z g_ij
  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
@@ -294,33 +345,36 @@ const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 		     CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 					 CCTK_VARIABLE_REAL,
+ // psi              partial_x psi       partial_y psi       partial_z psi
+ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 	  };
-const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
+
+const CCTK_INT operand_indices[]
 	= {
-	  0, 0, 0, 0,		// g_dd_11
-	  1, 1, 1, 1,		// g_dd_12
-	  2, 2, 2, 2,		// g_dd_13
-	  3, 3, 3, 3,		// g_dd_22
-	  4, 4, 4, 4,		// g_dd_23
-	  5, 5, 5, 5,		// g_dd_33
-	   6,  7,  8,		// K_dd_{11,12,13}
-	       9, 10,		// K_dd_{22,23}
-		  11,		// K_dd_33
+	  0, 0, 0, 0,		// g_dd_11, partial_[xyz] g_dd_11
+	  1, 1, 1, 1,		// g_dd_12, partial_[xyz] g_dd_12
+	  2, 2, 2, 2,		// g_dd_13, partial_[xyz] g_dd_13
+	  3, 3, 3, 3,		// g_dd_22, partial_[xyz] g_dd_22
+	  4, 4, 4, 4,		// g_dd_23, partial_[xyz] g_dd_23
+	  5, 5, 5, 5,		// g_dd_33, partial_[xyz] g_dd_33
+	  6, 7, 8, 9, 10, 11,	// K_dd_{11,12,13,22,23,33}
+	  12, 12, 12, 12,	// psi, partial_[xyz] psi
 	  };
 #define DERIV(x)	x
-const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
-	= {
-	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_11
-	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_12
-	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_13
-	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_22
-	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_23
-	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_33
-	  DERIV(0), DERIV(0), DERIV(0),			// K_dd_{11,12,13}
-		    DERIV(0), DERIV(0),			// K_dd_{22,23}
-			      DERIV(0),			// K_dd_{33}
-	  };
-void* const output_arrays[N_OUTPUT_ARRAYS]
+const CCTK_INT operation_codes[]
+  = {
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_11, partial_[xyz] g_dd_11
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_12, partial_[xyz] g_dd_12
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_13, partial_[xyz] g_dd_13
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_22, partial_[xyz] g_dd_22
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_23, partial_[xyz] g_dd_23
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_33, partial_[xyz] g_dd_33
+    DERIV(0), DERIV(0), DERIV(0), DERIV(0), DERIV(0), DERIV(0),
+					    // K_dd_{11,12,13,22,23,33}
+    DERIV(0), DERIV(1), DERIV(2), DERIV(3), // psi, partial_[xyz] psi
+    };
+
+void* const output_arrays[]
   = {
     static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_11)),
     static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)),
@@ -352,20 +406,48 @@ void* const output_arrays[N_OUTPUT_ARRAYS]
     static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_22)),
     static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_23)),
     static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_33)),
+    static_cast<void*>(ps.gridfn_data(gfns::gfn__psi)),
+    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_1)),
+    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_2)),
+    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_3)),
     };
 
-bool first_time = true;
-if (first_time)
-   then {
-	first_time = false;
+const int N_output_arrays_for_psi = 4;
+const int N_output_arrays_dim = sizeof(output_arrays)
+				/ sizeof(output_arrays[0]);
+const int N_output_arrays_use
+	= psi_flag ? N_output_arrays_dim
+		   : N_output_arrays_dim - N_output_arrays_for_psi;
+
 
+//
+// ***** parameter table *****
+//
+
+// this flag is true if and only if the parameter table already has the
+//	operand_indices
+//	operand_codes
+// entries for  psi_flag == par_table_psi_flag .
+static bool par_table_setup = false;
+
+// if  par_table_setup,
+//    this flag is the value of  psi_flag  for the parameter table entries
+// otherwise this flag is ignored
+static bool par_table_psi_flag = false;
+
+if (par_table_setup && (psi_flag == par_table_psi_flag))
+   then {
+	// parameter table is already set to just what we need
+	// ==> no-op here
+	}
+   else {
 	// store derivative info in interpolator parameter table
 	if (print_msg_flag)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "         setting up interpolator derivative info");
 
 	status = Util_TableSetIntArray(gi.param_table_handle,
-				       N_OUTPUT_ARRAYS, operand_indices,
+				       N_output_arrays_use, operand_indices,
 				       "operand_indices");
 	if (status < 0)
 	   then error_exit(ERROR_EXIT,
@@ -376,7 +458,7 @@ if (first_time)
 			   status);				/*NOTREACHED*/
 
 	status = Util_TableSetIntArray(gi.param_table_handle,
-				       N_OUTPUT_ARRAYS, operation_codes,
+				       N_output_arrays_use, operation_codes,
 				       "operation_codes");
 	if (status < 0)
 	   then error_exit(ERROR_EXIT,
@@ -385,8 +467,15 @@ if (first_time)
 "        Util_TableSetIntArray() status=%d\n"
 ,
 			   status);				/*NOTREACHED*/
+
+	par_table_setup = true;
+	par_table_psi_flag = psi_flag;
 	}
 
+
+//
+// ***** the actual interpolation *****
+//
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "         calling interpolator (%d points)",
@@ -397,13 +486,18 @@ status = CCTK_InterpLocalUniform(N_GRID_DIMS,
 				 N_interp_points,
 				    interp_coords_type_code,
 				    interp_coords,
-				 N_INPUT_ARRAYS,
+				 N_input_arrays_use,
 				    cgi.gridfn_dims,
 				    input_array_type_codes,
 				    input_arrays,
-				 N_OUTPUT_ARRAYS,
+				 N_output_arrays_use,
 				    output_array_type_codes,
 				    output_arrays);
+
+
+//
+// ***** handle any interpolation errors *****
+//
 if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
    then {
 	// look in interpolator output table entries
@@ -461,6 +555,7 @@ if (status < 0)
 ,
 		   status);					/*NOTREACHED*/
 
+
 return true;						// *** NORMAL RETURN ***
 }
 	  }
@@ -468,6 +563,129 @@ return true;						// *** NORMAL RETURN ***
 //******************************************************************************
 
 //
+// This function converts the g_dd_ij and partial_d_g_dd_kij gridfns
+// from the Cactus conformal semantics to the physical semantics.  As
+// documented in the CactusEinstein/ConformalState thorn guide, the
+// Cactus conformal semantics are
+//	physical_gij = psi^4 * stored_gij
+// From this it's trivial to derive the transformation for the derivatives,
+//	partial_k physical_gij = 4 psi^3 (partial_k psi) stored_gij
+//				 + psi^4 partial_k stored_gij
+//
+// Inputs (angular gridfns, all on the nominal grid):
+//   psi					# $\psi$
+//   partial_d_psi_{1,2,3}			# $\partial_k \psi$
+//   g_dd_{11,12,13,22,23,33}			# $\stored{g}_{ij}$
+//   partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k \stored{g}_{ij}$
+//
+//
+// Outputs (angular gridfns, all on the nominal grid):
+//   g_dd_{11,12,13,22,23,33}			# $g_{ij}$
+//   partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
+//
+namespace {
+void convert_conformal_to_physical(patch_system& ps, bool print_msg_flag)
+{
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "      converting Cactus conformal gij --> physical gij");
+
+    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 psi = p.gridfn(gfns::gfn__psi, irho,isigma);
+	const fp psi3 = jtutil::pow3(psi);
+	const fp psi4 = jtutil::pow4(psi);
+
+	const fp partial_d_psi_1
+		= p.gridfn(gfns::gfn__partial_d_psi_1, irho,isigma);
+	const fp partial_d_psi_2
+		= p.gridfn(gfns::gfn__partial_d_psi_2, irho,isigma);
+	const fp partial_d_psi_3
+		= p.gridfn(gfns::gfn__partial_d_psi_3, irho,isigma);
+
+	const fp stored_g_dd_11 = p.gridfn(gfns::gfn__g_dd_11, irho, isigma);
+	const fp stored_g_dd_12 = p.gridfn(gfns::gfn__g_dd_12, irho, isigma);
+	const fp stored_g_dd_13 = p.gridfn(gfns::gfn__g_dd_13, irho, isigma);
+	const fp stored_g_dd_22 = p.gridfn(gfns::gfn__g_dd_22, irho, isigma);
+	const fp stored_g_dd_23 = p.gridfn(gfns::gfn__g_dd_23, irho, isigma);
+	const fp stored_g_dd_33 = p.gridfn(gfns::gfn__g_dd_33, irho, isigma);
+
+	p.gridfn(gfns::gfn__g_dd_11, irho, isigma) *= psi4;
+	p.gridfn(gfns::gfn__g_dd_12, irho, isigma) *= psi4;
+	p.gridfn(gfns::gfn__g_dd_13, irho, isigma) *= psi4;
+	p.gridfn(gfns::gfn__g_dd_22, irho, isigma) *= psi4;
+	p.gridfn(gfns::gfn__g_dd_23, irho, isigma) *= psi4;
+	p.gridfn(gfns::gfn__g_dd_33, irho, isigma) *= psi4;
+
+	p.gridfn(gfns::gfn__partial_d_g_dd_111, irho,isigma)
+		= 4.0*psi3*partial_d_psi_1*stored_g_dd_11
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_111, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_112, irho,isigma)
+		= 4.0*psi3*partial_d_psi_1*stored_g_dd_12
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_112, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_113, irho,isigma)
+		= 4.0*psi3*partial_d_psi_1*stored_g_dd_13
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_113, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_122, irho,isigma)
+		= 4.0*psi3*partial_d_psi_1*stored_g_dd_22
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_122, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_123, irho,isigma)
+		= 4.0*psi3*partial_d_psi_1*stored_g_dd_23
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_123, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_133, irho,isigma)
+		= 4.0*psi3*partial_d_psi_1*stored_g_dd_33
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_133, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_211, irho,isigma)
+		= 4.0*psi3*partial_d_psi_2*stored_g_dd_11
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_211, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_212, irho,isigma)
+		= 4.0*psi3*partial_d_psi_2*stored_g_dd_12
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_212, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_213, irho,isigma)
+		= 4.0*psi3*partial_d_psi_2*stored_g_dd_13
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_213, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_222, irho,isigma)
+		= 4.0*psi3*partial_d_psi_2*stored_g_dd_22
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_222, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_223, irho,isigma)
+		= 4.0*psi3*partial_d_psi_2*stored_g_dd_23
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_223, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_233, irho,isigma)
+		= 4.0*psi3*partial_d_psi_2*stored_g_dd_33
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_233, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_311, irho,isigma)
+		= 4.0*psi3*partial_d_psi_3*stored_g_dd_11
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_311, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_312, irho,isigma)
+		= 4.0*psi3*partial_d_psi_3*stored_g_dd_12
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_312, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_313, irho,isigma)
+		= 4.0*psi3*partial_d_psi_3*stored_g_dd_13
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_313, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_322, irho,isigma)
+		= 4.0*psi3*partial_d_psi_3*stored_g_dd_22
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_322, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_323, irho,isigma)
+		= 4.0*psi3*partial_d_psi_3*stored_g_dd_23
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_323, irho,isigma);
+	p.gridfn(gfns::gfn__partial_d_g_dd_333, irho,isigma)
+		= 4.0*psi3*partial_d_psi_3*stored_g_dd_33
+		  + psi4*p.gridfn(gfns::gfn__partial_d_g_dd_333, irho,isigma);
+	}
+	}
+    }
+}
+	  }
+
+//******************************************************************************
+
+//
 // This function computes H(h), and optionally its Jacobian coefficients,
 // (from which the Jacobian matrix may be computed later).  This function
 // uses a mixture of algebraic operations and (rho,sigma) finite differencing.