[ET Trac] #2944: CCE_Export Issues Ticket

14 May 2026


      #2944: CCE_Export Issues Ticket
Reporter: Zach Etienne
   Status: new
Milestone: 
  Version: 
     Type: bug
 Priority: major
Component:
@{5a73cd832d61371e861f3270}
The below review of CCE_Export was the product of a number of LLM worker agents, including Qwen3.6-27b, Claude Sonnet-4.7-XHigh, as well as a supervisor GPT-5.5-XHigh agent.
```
* Issue 1
 - Severity high
 - Description of the error
   The interpolation helpers ignore too many error/status paths. The ADMBase
   variable names used here are fixed by the thorn and inherited from ADMBase,
   so "bad user-supplied variable name" is not the main risk. However,
   CCTK_VarIndex, CCTK_InterpHandle, Util_TableCreate, the table setters,
   CCTK_CoordSystemHandle, CCTK_InterpGridArrays, and Util_TableDestroy can all
   report failure, and the output buffers are consumed immediately afterward.
   This is especially important for out-of-domain extraction radii and missing
   or misconfigured interpolation support.
 - Filename + line number(s)
   repos/CCE_Export/src/interpolate.cc:11, 33-60, 68, 85-112
   repos/einsteinanalysis/Multipole/src/interpolate.cc:13-19, 69-79
   repos/einsteinanalysis/ET_BHaHAHA/src/BHaHAHA_interpolate_metric_data.c:142-151
 - Original code snippet
   CCTK_INT variable_index = CCTK_VarIndex(name.c_str());
   const int operator_handle =
       CCTK_InterpHandle("Hermite polynomial interpolation");
   int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
   Util_TableSetFromString(param_table_handle, "order=3 ...");
   Util_TableSetIntArray(param_table_handle, 4, operand_indices,
                         "operand_indices");
   Util_TableSetIntArray(param_table_handle, 4, operation_codes,
                         "operation_codes");
   const int coord_system_handle = CCTK_CoordSystemHandle("cart3d");
   CCTK_InterpGridArrays(cctkGH, num_dims, operator_handle,
                         param_table_handle, coord_system_handle,
                         CCTK_MyProc(cctkGH) == 0 ? array_size : 0,
                         CCTK_VARIABLE_REAL, interp_coords, num_input_arrays,
                         input_array_indices, num_output_arrays,
                         output_array_types, output_arrays);
 - Suggested patch
   Add abort-level checks before any interpolated data are consumed:
CCTK_INT variable_index = CCTK_VarIndex(name.c_str());
   if (variable_index < 0)
     CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Could not find variable '%s'", name.c_str());
const int operator_handle =
       CCTK_InterpHandle("Hermite polynomial interpolation");
   if (operator_handle < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not find Hermite interpolator");
int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
   if (param_table_handle < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not create interpolation table");
int ierr = Util_TableSetFromString(param_table_handle, "order=3 ...");
   if (ierr < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not populate interpolation table");
ierr = Util_TableSetIntArray(param_table_handle, 4, operand_indices,
                                "operand_indices");
   if (ierr < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not set operand_indices");
ierr = Util_TableSetIntArray(param_table_handle, 4, operation_codes,
                                "operation_codes");
   if (ierr < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not set operation_codes");
const int coord_system_handle = CCTK_CoordSystemHandle("cart3d");
   if (coord_system_handle < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not find cart3d coordinate system");
ierr = CCTK_InterpGridArrays(...);
   if (ierr < 0)
     CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Interpolation failed for '%s' with status %d",
                name.c_str(), ierr);
ierr = Util_TableDestroy(param_table_handle);
   if (ierr < 0)
     CCTK_WARN(CCTK_WARN_ABORT, "Could not destroy interpolation table");
* Issue 2
 - Severity medium
 - Description of the error
   The interpolation call must remain collective, but the spherical-harmonic
   work after interpolation is only useful on rank 0. Cactus schedule option
   GLOBAL means once per grid hierarchy/component scheduling level, not only on
   MPI rank 0. The existing rank-0 guard around output confirms that all ranks
   enter the routine. Nonzero ranks request zero interpolation points and then
   perform decomposition work whose results are discarded.
 - Filename + line number(s)
   repos/CCE_Export/schedule.ccl:3-8
   repos/CCE_Export/src/cce_export.cc:135-183
   repos/flesh/doc/UsersGuide/Appendices.tex:1034-1085
 - Original code snippet
   Extract_Metric_Shift_Lapse_On_Sphere(...);
   Compute_Ylms(th, ph, re_ylms, im_ylms, lmax, array_size);
   for (int i = 0; i < 3; i++) {
     for (int j = i; j < 3; j++) {
       Decompose_Spherical_Harmonics(...);
     }
   }
   ...
   if (CCTK_MyProc(cctkGH) == 0) {
     Output_Decomposed_Metric_Data(...);
   }
 - Suggested patch
   Keep the collective interpolation outside the rank guard, then perform mode
   decomposition and output only on rank 0:
Extract_Metric_Shift_Lapse_On_Sphere(...);
if (CCTK_MyProc(cctkGH) == 0) {
     Compute_Ylms(th, ph, re_ylms, im_ylms, lmax, array_size);
for (int i = 0; i < 3; i++) {
       for (int j = i; j < 3; j++) {
         Decompose_Spherical_Harmonics(...);
       }
     }
for (int i = 0; i < 3; i++) {
       Decompose_Spherical_Harmonics(... beta ...);
       Decompose_Spherical_Harmonics(... dr_beta ...);
       Decompose_Spherical_Harmonics(... dt_beta ...);
     }
Decompose_Spherical_Harmonics(... alpha ...);
     Decompose_Spherical_Harmonics(... dr_alpha ...);
     Decompose_Spherical_Harmonics(... dt_alpha ...);
Output_Decomposed_Metric_Data(...);
   }
* Issue 3
 - Severity high
 - Description of the error
   CCTK_REAL precision is configurable in Cactus/ET. The writer stores data in
   vector<CCTK_REAL> but tells HDF5 that the memory buffer is H5T_NATIVE_DOUBLE.
   This is only safe for CCTK_REAL_PRECISION_8. CarpetIOHDF5 maps CCTK_REAL to
   H5T_NATIVE_FLOAT, H5T_NATIVE_DOUBLE, or H5T_NATIVE_LDOUBLE according to the
   configured precision, which is the relevant ET convention. The utility reader
   also needs a deliberate memory datatype if the on-disk datatype changes.
 - Filename + line number(s)
   repos/CCE_Export/src/h5_export.cc:65, 132
   repos/CCE_Export/src/util/cce_export_hdf5_to_ascii.c:50, 76-77
   repos/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh:22-27
   repos/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc:249-302
 - Original code snippet
   dataset = H5Dcreate(file, datasetname.c_str(), H5T_NATIVE_DOUBLE, dataspace,
                       cparms);
   ...
   HDF5_ERROR(H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memdataspace, filespace,
                       H5P_DEFAULT, data));
 - Suggested patch
   Use a CCTK_REAL-aware memory datatype for writes. If SpECTRE compatibility
   requires double-precision datasets on disk, keep the dataset datatype double
   and let HDF5 convert from the correct memory datatype:
#if defined(CCTK_REAL_PRECISION_16)
   #define CCE_HDF5_REAL_MEMORY H5T_NATIVE_LDOUBLE
   #elif defined(CCTK_REAL_PRECISION_8)
   #define CCE_HDF5_REAL_MEMORY H5T_NATIVE_DOUBLE
   #elif defined(CCTK_REAL_PRECISION_4)
   #define CCE_HDF5_REAL_MEMORY H5T_NATIVE_FLOAT
   #else
   #error "Unsupported CCTK_REAL precision"
   #endif
const hid_t file_real_type = H5T_NATIVE_DOUBLE;
HDF5_ERROR(dataset = H5Dcreate(file, datasetname.c_str(), file_real_type,
                                  dataspace, cparms));
   HDF5_ERROR(H5Dwrite(dataset, CCE_HDF5_REAL_MEMORY, memdataspace, filespace,
                       H5P_DEFAULT, data));
In the utility reader, allocate the requested output type deliberately and
   pass the matching HDF5 memory datatype, for example:
double *data = NULL;
   status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                    H5P_DEFAULT, data);
* Issue 4
 - Severity medium
 - Description of the error
   This thorn has C++14/C++17 dependencies without an explicit local contract.
   Full ET context does not make this disappear: Cactus documentation still says
   C++11 is required, and several simfactory optionlists still select C++11 or
   C++14, even though newer generic and machine optionlists commonly use C++17.
   The code uses std::make_unique, comma-separated using-declarations, and
   std::filesystem. The filesystem feature-test check is also unreliable because
   it tests __cpp_lib_filesystem before including a header that would define it.
 - Filename + line number(s)
   repos/CCE_Export/src/h5_export.cc:3-8, 13
   repos/CCE_Export/src/cce_export.cc:12
   repos/CCE_Export/src/cce_export.hh:15
   repos/CCE_Export/src/spherical_harmonic_decomposition.cc:5, 81-82
   repos/flesh/doc/UsersGuide/Notes.tex:47-49
   repos/simfactory2/mdb/optionlists/generic.cfg:30-32
   repos/simfactory2/mdb/optionlists/debian-cuda.cfg:44
   repos/simfactory2/mdb/optionlists/wheeler.cfg:19
 - Original code snippet
   #if defined __cpp_lib_filesystem && __cpp_lib_filesystem < 201703L
   #include <experimental/filesystem>
   namespace fs = std::experimental::filesystem;
   #else
   #include <filesystem>
   namespace fs = std::filesystem;
   #endif
   using std::string, std::ostringstream, std::map, std::ios, std::setprecision;
   static auto fr = std::make_unique<CCTK_REAL []>(array_size);
 - Suggested patch
   Either make C++17 an explicit requirement, or remove the C++17-only code. A
   minimal explicit-contract patch is:
#if __cplusplus < 201703L
   #error "CCE_Export requires C++17 or later"
   #endif
#include <filesystem>
   namespace fs = std::filesystem;
using std::string;
   using std::ostringstream;
   using std::map;
   using std::ios;
   using std::setprecision;
Also document the C++17 requirement in README/ThornGuide and ensure the
   build configuration used for this thorn selects a C++17-capable compiler
   mode. If C++11 compatibility is desired instead, replace std::filesystem
   with Cactus/local path handling and replace comma-separated using-declarations
   with individual declarations.
* Issue 5
 - Severity high
 - Description of the error
   For gravitational-wave extraction/CCE, r = 0 is not a valid extraction
   radius or worldtube. The default configuration activates one radius at
   radius[0] = 0.0, so a user who merely activates the thorn with defaults is
   requesting mathematically invalid data. The implementation then compounds
   this by computing radial derivatives through x^i/r partial_i, which divides
   by radius[rad_index]. Unlike Multipole, CCE_Export explicitly exports radial
   derivatives of worldtube data, so zero radius must be rejected, not treated
   as an allowed lower endpoint.
 - Filename + line number(s)
   repos/CCE_Export/param.ccl:9-17
   repos/CCE_Export/src/interpolate.cc:173-179, 197-202, 215-218
 - Original code snippet
   CCTK_INT nradii "How many radii to export the data on" STEERABLE=recover
   {
     0:100 :: "number of radii"
   } 1
CCTK_REAL radius[101] "Radii to export the metric data on" STEERABLE=recover
   {
     0.0:* :: "Extraction radius"
   } 0.0
dr_g.at(i).at(j).at(array_index) =
       (xs.at(array_index) / radius[rad_index]) *
           dx_g.at(i).at(j).at(array_index) +
       (ys.at(array_index) / radius[rad_index]) *
           dy_g.at(i).at(j).at(array_index) +
       (zs.at(array_index) / radius[rad_index]) *
           dz_g.at(i).at(j).at(array_index);
 - Suggested patch
   Make the default inactive or positive, and reject every active non-positive
   extraction radius before interpolation:
CCTK_INT nradii "How many radii to export the data on" STEERABLE=recover
   {
     0:100 :: "number of radii"
   } 0
CCTK_REAL radius[101] "Radii to export the metric data on" STEERABLE=recover
   {
     (0.0:* :: "Positive extraction radius"
   } 1.0
for (int ir = 0; ir < nradii; ++ir) {
     if (radius[ir] <= 0.0) {
       CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "CCE_Export::radius[%d] = %.17g is invalid: CCE/worldtube "
                  "extraction radii must be strictly positive",
                  ir, static_cast<double>(radius[ir]));
     }
   }
* Issue 6
 - Severity medium
 - Description of the error
   Output is documented and implemented as one file per extraction radius, but
   the filename uses a float argument and fixed two-decimal formatting. Distinct
   CCTK_REAL radii can therefore collide in the same output filename. The
   parameter documentation mentions two decimal places, but that does not
   protect the data from accidental collisions and is inconsistent with "a
   separate h5 file for each desired extraction radius."
 - Filename + line number(s)
   repos/CCE_Export/src/h5_export.hh:25-39
   repos/CCE_Export/src/h5_export.cc:154, 168-171
   repos/CCE_Export/src/cce_export.cc:178-182
   repos/CCE_Export/param.ccl:25-29
   repos/CCE_Export/doc/documentation.tex:49-50
 - Original code snippet
   void Output_Decomposed_Metric_Data(..., float rad, int lmax) {
     ...
     basename << base_file_name << "R" << setiosflags(ios::fixed)
              << setprecision(2) << rad << "." << extension;
     string output_name = (fs::path(my_out_dir) / basename.str()).string();
   }
 - Suggested patch
   Preserve CCTK_REAL precision through the output path, increase filename
   precision, and validate that active radii produce unique names:
void Output_Decomposed_Metric_Data(..., CCTK_REAL rad, int lmax) {
     ...
     basename << base_file_name << "R" << setiosflags(ios::fixed)
              << setprecision(12) << static_cast<double>(rad)
              << "." << extension;
     string output_name = (fs::path(my_out_dir) / basename.str()).string();
   }
std::setstd::string output_names;
   for (int ir = 0; ir < nradii; ++ir) {
     std::ostringstream name;
     name << base_file_name << "R" << std::fixed << std::setprecision(12)
          << static_cast<double>(radius[ir]) << "." << extension;
     if (!output_names.insert(name.str()).second) {
       CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Multiple CCE_Export radii map to output file '%s'",
                  name.str().c_str());
     }
   }
* Issue 7
 - Severity medium
 - Description of the error
   The documentation is wrong in two places. The metric components are
   interpolated directly from ADMBase; the code uses the ADM extrinsic-curvature
   relation to compute partial_t gamma_ij, not to compute gamma_ij itself. The
   spherical-harmonic coefficient integral implemented here is over dOmega using
   sin(theta) dtheta dphi and does not multiply by r^2.
 - Filename + line number(s)
   repos/CCE_Export/doc/documentation.tex:76-84, 93-100
   repos/CCE_Export/src/interpolate.cc:166-170, 221-311
   repos/CCE_Export/src/spherical_harmonic_decomposition.cc:84-96
 - Original code snippet
   The metric is computed from the extrinsic curvature, lapse, and shift as:
\begin{equation}
   K_{ij} = \frac{1}{2 \alpha} \left ( -\partial_0 \gamma_{ij}
     + \beta^k \partial_k \gamma_{ij}
     + \gamma_{ki} \partial_j \beta^k
     + \gamma_{kj} \partial_i \beta^k \right ) , .
   \end{equation}
C^{lm}(t, r) = \int {}_0 Y_{lm}^* u(t, r, \theta, \varphi) r^2 d \Omega
 - Suggested patch
   Update the ThornGuide to describe the implemented mathematics:
The metric components are interpolated directly from ADMBase. The metric
   time derivative is computed by inverting the ADM definition of extrinsic
   curvature:
\begin{equation}
   \partial_t \gamma_{ij} =
     -2\alpha K_{ij}
     + \beta^k \partial_k \gamma_{ij}
     + \gamma_{ki} \partial_j \beta^k
     + \gamma_{kj} \partial_i \beta^k .
   \end{equation}
The spherical-harmonic coefficients are
\begin{eqnarray}
     C^{lm}(t, r) =
       \int {}_0 Y_{lm}^* u(t, r, \theta, \varphi) d\Omega .
   \end{eqnarray}
* Issue 8
 - Severity medium
 - Description o
--
Ticket URL: https://bitbucket.org/einsteintoolkit/tickets/issues/2944/cce_export-issues-...

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

[ET Trac] #2944: CCE_Export Issues Ticket