From 0cca329914063e0a34b472c77e81385d0a2b5560 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 3 Oct 2002 20:56:46 +0000 Subject: 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 --- src/gr/gfns.hh | 31 +++-- src/gr/gr.hh | 8 ++ src/gr/horizon_function.cc | 318 ++++++++++++++++++++++++++++++++++++++------- 3 files changed, 297 insertions(+), 60 deletions(-) (limited to 'src/gr') 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(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(cgi.g_dd_11_data), static_cast(cgi.g_dd_12_data), @@ -277,11 +318,21 @@ const void* const input_arrays[N_INPUT_ARRAYS] static_cast(cgi.K_dd_22_data), static_cast(cgi.K_dd_23_data), static_cast(cgi.K_dd_33_data), + static_cast(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(ps.gridfn_data(gfns::gfn__g_dd_11)), static_cast(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)), @@ -352,20 +406,48 @@ void* const output_arrays[N_OUTPUT_ARRAYS] static_cast(ps.gridfn_data(gfns::gfn__K_dd_22)), static_cast(ps.gridfn_data(gfns::gfn__K_dd_23)), static_cast(ps.gridfn_data(gfns::gfn__K_dd_33)), + static_cast(ps.gridfn_data(gfns::gfn__psi)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_psi_1)), + static_cast(ps.gridfn_data(gfns::gfn__partial_d_psi_2)), + static_cast(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,12 +555,136 @@ if (status < 0) , status); /*NOTREACHED*/ + 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 -- cgit v1.2.3