Add ext curvature to PreprocessCceWorldtube

cmake/AddInputFileTests.cmake

@@ -277,6 +277,7 @@ configure_file(
 # These paths should be relative to the input file directory passed to
 # `add_input_file_tests`
 set(INPUT_FILE_WHITELIST
-    "PreprocessCceWorldtube/PreprocessCceWorldtube.yaml")
+    "PreprocessCceWorldtube/PreprocessCceWorldtube.yaml"
+    "PreprocessCceWorldtube/AdmPreprocessCceWorldtube.yaml")
 
-add_input_file_tests("${CMAKE_SOURCE_DIR}/tests/InputFiles/" ${INPUT_FILE_WHITELIST})
+add_input_file_tests("${CMAKE_SOURCE_DIR}/tests/InputFiles/" "${INPUT_FILE_WHITELIST}")

docs/Tutorials/CCE.md

@@ -181,13 +181,13 @@ complex modes in m-varies-fastest format. That is,
 Each dataset in the H5 file must also have an attribute
 named `Legend` which is an ASCII-encoded null-terminated variable-length string.
 
-##### Spherical harmonic nodes {#spherical_nodes}
+#### Spherical harmonic nodes {#spherical_nodes}
 
 When we refer to a "nodal" data format, we mean that the worldtube data are
-stored as complex values at specially chosen collocation points (a.k.a. grid
-points or nodes). This allows SpECTRE to perform integrals, derivatives, and
-interpolation exactly on the input data. These grid points are Gauss-Legendre in
-$cos(\theta)$ and equally spaced in $\phi$.
+stored as values at specially chosen collocation points (a.k.a. grid points or
+nodes). This allows SpECTRE to perform integrals, derivatives, and interpolation
+exactly on the input data. These grid points are Gauss-Legendre in $cos(\theta)$
+and equally spaced in $\phi$.
 
 Below is a routine for computing the spherical
 harmonic $\theta$ and $\phi$ values. These can be used to compute the Cartesian
@@ -223,6 +223,46 @@ named `Legend` which is an ASCII-encoded null-terminated variable-length string.
 \note Nodal data is likely the easiest to write out since no conversion to
 spherical harmonic coefficients is necessary.
 
+#### ADM Cartesian metric and derivatives {#adm_cartesian_metric_and_derivatives}
+
+For worldtube data stored in an H5 file in the "ADM metric nodal" format, there
+must be the following datasets with these exact names (including the `.dat`
+suffix):
+
+- `gxx.dat`, `gxy.dat`, `gxz.dat`, `gyy.dat`, `gyz.dat`, `gzz.dat`
+- `Dxgxx.dat`, `Dxgxy.dat`, `Dxgxz.dat`, `Dxgyy.dat`, `Dxgyz.dat`, `Dxgzz.dat`
+- `Dygxx.dat`, `Dygxy.dat`, `Dygxz.dat`, `Dygyy.dat`, `Dygyz.dat`, `Dygzz.dat`
+- `Dzgxx.dat`, `Dzgxy.dat`, `Dzgxz.dat`, `Dzgyy.dat`, `Dzgyz.dat`, `Dzgzz.dat`
+- `Shiftx.dat`, `Shifty.dat`, `Shiftz.dat`
+- `DxShiftx.dat`, `DxShifty.dat`, `DxShiftz.dat`
+- `DyShiftx.dat`, `DyShifty.dat`, `DyShiftz.dat`
+- `DzShiftx.dat`, `DzShifty.dat`, `DzShiftz.dat`
+- `Lapse.dat`, `DxLapse.dat`, `DyLapse.dat`, `DzLapse.dat`
+- `Kxx.dat`, `Kxy.dat`, `Kxz.dat`, `Kyy.dat`, `Kyz.dat`, `Kzz.dat`
+- `AuxiliaryShiftx.dat`, `AuxiliaryShifty.dat`, `AuxiliaryShiftz.dat`
+
+Here `g` represents the spacetime metric, but we only require the spatial
+components (e.g. `gxx.dat`, `gxy.dat`, etc...) so in practice, those are the
+tensor components of the spatial metric. The temporal components of the
+spacetime metric are stored separately in the lapse and shift. Each of the
+spatial metric, lapse, and shift must also have their cartesian derivatives. `K`
+is the extrinsic curvature and `AuxiliaryShift` is the auxiliary shift vector
+used in the Gamma-driver condition. We will compute the time derivative of the
+spatial metric using Eq. (2.134) of \cite BaumgarteShapiro, the time derivative
+of the lapse using the `1+log` slicing condition from Eq. (4.87) of
+\cite BaumgarteShapiro, and the time derivative of the shift using the
+Gamma-driver condition from Eq. (4.89) of \cite BaumgarteShapiro.
+
+\warning If your worldtube data is in the ADM metric nodal format but you have
+not used `1+log` slicing and the Gamma-driver condition, your time derivatives
+will be **wrong**. If you'd like us to support other commonly used gauge
+conditions, please open an issue on our
+[GitHub](https://github.com/sxs-collaboration/spectre) or email the core
+developers at [[email protected]](mailto:[email protected]).
+
+The layout of each of these datasets must be
+[spherical harmonic nodes](#spherical_nodes).
+
 #### Cartesian metric and derivatives {#cartesian_metric_and_derivatives}
 
 For worldtube data stored in an H5 file in either the "metric nodal" or "metric
@@ -251,7 +291,8 @@ each of these datasets must be in either
 
 In the "bondi nodal" format, you must have the same Bondi variables as the
 [required format](#required_h5_worldtube_data_format), but each variable layout
-must be the [spherical harmonic nodal layout](#spherical_nodes).
+must be the [spherical harmonic nodal layout](#spherical_nodes) with complex
+values interleaved as `Re`, `Im`, `Re`, `Im`, ...
 
 If you already have data in the
 [required "bondi modal" format](#required_h5_worldtube_data_format), then
@@ -297,6 +338,9 @@ Here are some notes about the different options in the YAML input file:
   *ascending* m.
 - `BufferDepth` is an advanced option that lets you load more data into RAM at
   once so there are fewer filesystem accesses.
+- If using the ADM metric nodal input data format, specify the
+  `InputDataFormat:` like so:
+\snippet AdmPreprocessCceWorldtube.yaml adm_input_data_format_example
 
 ### What Worldtube data "should" look like {#worldtube_data_looks}
 
@@ -310,7 +354,8 @@ The 2,2 modes are oscillatory and capture the orbits of the two objects. The
 real part of the 2,0 mode contains the gravitational memory of the system. Then
 for this system, all the other modes are subdominant.
 
-If you are using the [cartesian metric](#cartesian_metric_and_derivatives)
+If you are using the [cartesian metric](#cartesian_metric_and_derivatives) or
+[adm cartesian metric](#adm_cartesian_metric_and_derivatives)
 worldtube format, here is a plot of the imaginary part of the 2,2 mode of the
 lapse and its radial and time derivative during inspiral.
 

src/Executables/PreprocessCceWorldtube/PreprocessCceWorldtube.cpp

@@ -56,6 +56,8 @@ using modal_bondi_input_tags = Cce::Tags::worldtube_boundary_tags_for_writing<
 using nodal_bondi_input_tags =
     Cce::Tags::worldtube_boundary_tags_for_writing<Cce::Tags::BoundaryValue>;
 
+using AdmOptions = Cce::MetricWorldtubeH5BufferUpdater<DataVector>::AdmOptions;
+
 // from a data-varies-fastest set of buffers provided by
 // `MetricWorldtubeH5BufferUpdater` extract the set of coefficients for a
 // particular time given by `buffer_time_offset` into the `time_span` size of
@@ -326,12 +328,13 @@ void bondi_nodal_to_bondi_modal(
   }
 }
 
-void metric_nodal_to_bondi_modal(
-    const std::string& input_file, const std::string& output_file,
-    const size_t input_buffer_depth,
-    const std::optional<double>& extraction_radius) {
+void metric_nodal_to_bondi_modal(const std::string& input_file,
+                                 const std::string& output_file,
+                                 const size_t input_buffer_depth,
+                                 const std::optional<double>& extraction_radius,
+                                 const std::optional<AdmOptions>& adm_options) {
   Cce::MetricWorldtubeH5BufferUpdater<DataVector> buffer_updater{
-      input_file, extraction_radius, false};
+      input_file, extraction_radius, false, adm_options};
   const size_t l_max = buffer_updater.get_l_max();
 
   const size_t number_of_angular_points =
@@ -377,7 +380,50 @@ void metric_nodal_to_bondi_modal(
   }
 }
 
-enum class InputDataFormat { MetricNodal, MetricModal, BondiNodal, BondiModal };
+// If we use the AdmOptions class directly in an option tag, the input file
+// would look like:
+//
+// InputDataFormat:
+//   AdvectiveLapse: True
+//   AdvectiveShift: True
+//   AuxiliaryShiftFactor: False
+//
+// However, this doesn't have the name `AdmMetricNodal`. This is consistent with
+// how our options work, but is confusing (especially for users who are
+// unfamiliar with spectre) since the actual input data format isn't shown. This
+// class now changes the input file to look like:
+//
+// InputDataFormat:
+//   AdmMetricNodal:
+//     AdvectiveLapse: True
+//     AdvectiveShift: True
+//     AuxiliaryShiftFactor: False
+//
+// which is much clearer.
+struct AdmMetricNodalOptions {
+  struct AdmMetricNodal {
+    using type = ::AdmOptions;
+    static constexpr Options::String help = ::AdmOptions::help;
+  };
+
+  using options = tmpl::list<AdmMetricNodal>;
+  static constexpr Options::String help = ::AdmOptions::help;
+
+  AdmMetricNodalOptions() = default;
+  // NOLINTNEXTLINE
+  AdmMetricNodalOptions(::AdmOptions adm_metric_nodal_in)
+      : adm_metric_nodal(adm_metric_nodal_in) {}
+
+  ::AdmOptions adm_metric_nodal;
+};
+
+enum class InputDataFormat {
+  AdmMetricNodal,
+  MetricNodal,
+  MetricModal,
+  BondiNodal,
+  BondiModal
+};
 
 std::ostream& operator<<(std::ostream& os,
                          const InputDataFormat input_data_format) {
@@ -407,10 +453,11 @@ struct InputH5Files {
 };
 
 struct InputDataFormat {
-  using type = ::InputDataFormat;
+  using type = std::variant<::InputDataFormat, AdmMetricNodalOptions>;
   static constexpr Options::String help =
-      "The type of data stored in the 'InputH5Files'. Can be  'MetricNodal', "
-      "'MetricModal', 'BondiNodal', or 'BondiModal'.";
+      "The type of data stored in the 'InputH5Files'. Can be 'AdmMetricNodal' "
+      "with additional options, 'MetricNodal', 'MetricModal', 'BondiNodal', or "
+      "'BondiModal'.";
 };
 
 struct OutputH5File {
@@ -496,11 +543,23 @@ struct InputH5Files : db::SimpleTag {
 };
 
 struct InputDataFormat : db::SimpleTag {
-  using type = ::InputDataFormat;
+ private:
+  using option_type = std::variant<::InputDataFormat, AdmMetricNodalOptions>;
+
+ public:
+  using type = std::variant<::InputDataFormat, AdmOptions>;
   using option_tags = tmpl::list<OptionTags::InputDataFormat>;
   static constexpr bool pass_metavariables = false;
-  static type create_from_options(type input_data_format) {
-    return input_data_format;
+  static type create_from_options(option_type input_data_format) {
+    type result{};
+    if (std::holds_alternative<::InputDataFormat>(input_data_format)) {
+      result = std::get<::InputDataFormat>(input_data_format);
+    } else {
+      result =
+          std::get<AdmMetricNodalOptions>(input_data_format).adm_metric_nodal;
+    }
+
+    return result;
   }
 };
 
@@ -582,8 +641,8 @@ struct Options::create_from_yaml<InputDataFormat> {
       return InputDataFormat::BondiModal;
     }
     PARSE_ERROR(options.context(),
-                "InputDataFormat must be 'MetricNodal', 'MetricModal', "
-                "'BondiNodal', or 'BondiModal'");
+                "InputDataFormat must be 'AdmMetricNodal', 'MetricNodal', "
+                "'MetricModal', 'BondiNodal', or 'BondiModal'");
   }
 };
 
@@ -638,8 +697,12 @@ int main(int argc, char** argv) {
     const TagsTuple inputs =
         Parallel::create_from_options<void>(options, tags{});
 
-    const InputDataFormat input_data_format =
+    const auto& input_data_format =
         tuples::get<ReduceCceTags::InputDataFormat>(inputs);
+    const InputDataFormat input_data_format_enum =
+        std::holds_alternative<InputDataFormat>(input_data_format)
+            ? std::get<InputDataFormat>(input_data_format)
+            : InputDataFormat::AdmMetricNodal;
     const std::vector<std::string>& input_files =
         tuples::get<ReduceCceTags::InputH5Files>(inputs);
     const std::string& output_h5_file =
@@ -651,7 +714,7 @@ int main(int argc, char** argv) {
       // If the input format is BondiModal, then we don't actually have to do
       // any transformations, only combining H5 files. So the temporary file
       // name is just the output file
-      if (input_data_format == InputDataFormat::BondiModal) {
+      if (input_data_format_enum == InputDataFormat::BondiModal) {
         temporary_combined_h5_file = output_h5_file;
       } else {
         // Otherwise we have to do a transformation so a temporary H5 file is
@@ -668,7 +731,7 @@ int main(int argc, char** argv) {
       // Now combine the h5 files into a single file
       h5::combine_h5_dat(input_files, temporary_combined_h5_file.value(),
                          Verbosity::Quiet);
-    } else if (input_data_format == InputDataFormat::BondiModal) {
+    } else if (input_data_format_enum == InputDataFormat::BondiModal) {
       // Error here if the input data format is BondiModal since there's nothing
       // to do
       ERROR_NO_TRACE(
@@ -678,7 +741,7 @@ int main(int argc, char** argv) {
     }
 
     if (tuples::get<ReduceCceTags::FixSpecNormalization>(inputs)) {
-      if (input_data_format != InputDataFormat::MetricModal) {
+      if (input_data_format_enum != InputDataFormat::MetricModal) {
         ERROR_NO_TRACE(
             "The option FixSpecNormalization can only be 'true' when the input "
             "data format is MetricModal. Otherwise, it must be 'false'");
@@ -702,20 +765,22 @@ int main(int argc, char** argv) {
       }
     };
 
-    switch (input_data_format) {
-      case InputDataFormat::BondiModal:
+    switch (input_data_format_enum) {
+      case InputDataFormat::BondiModal: {
         // Nothing to do here because this is the desired output format and the
         // H5 files were combined above
         return 0;
-      case InputDataFormat::BondiNodal:
+      }
+      case InputDataFormat::BondiNodal: {
         bondi_nodal_to_bondi_modal(
             input_worldtube_filename(), output_h5_file,
             tuples::get<ReduceCceTags::BufferDepth>(inputs),
             tuples::get<ReduceCceTags::ExtractionRadius>(inputs));
 
         clean_temporary_file();
         return 0;
-      case InputDataFormat::MetricModal:
+      }
+      case InputDataFormat::MetricModal: {
         perform_cce_worldtube_reduction(
             input_worldtube_filename(), output_h5_file,
             tuples::get<ReduceCceTags::BufferDepth>(inputs),
@@ -726,16 +791,25 @@ int main(int argc, char** argv) {
 
         clean_temporary_file();
         return 0;
-      case InputDataFormat::MetricNodal:
+      }
+      case InputDataFormat::AdmMetricNodal: {
+        [[fallthrough]];
+      }
+      case InputDataFormat::MetricNodal: {
+        const std::optional<AdmOptions> adm_options =
+            std::holds_alternative<AdmOptions>(input_data_format)
+                ? std::optional{std::get<AdmOptions>(input_data_format)}
+                : std::nullopt;
         metric_nodal_to_bondi_modal(
             input_worldtube_filename(), output_h5_file,
             tuples::get<ReduceCceTags::BufferDepth>(inputs),
-            tuples::get<ReduceCceTags::ExtractionRadius>(inputs));
+            tuples::get<ReduceCceTags::ExtractionRadius>(inputs), adm_options);
 
         clean_temporary_file();
         return 0;
+      }
       default:
-        ERROR("Unknown input data format " << input_data_format);
+        ERROR("Unknown input data format " << input_data_format_enum);
     }
   } catch (const std::exception& exception) {
     Parallel::printf("%s\n", exception.what());

tests/InputFiles/PreprocessCceWorldtube/AdmPreprocessCceWorldtube.yaml

@@ -0,0 +1,17 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+InputH5File: InputFilename.h5
+OutputH5File: ReducedWorldtubeR0292.h5
+# [adm_input_data_format_example]
+InputDataFormat:
+  AdmMetricNodal:
+    AdvectiveLapse: True
+    AdvectiveShift: True
+    AuxiliaryShiftFactor: 0.75 # Eq. 4.89 B&S
+# [adm_input_data_format_example]
+ExtractionRadius: 292
+FixSpecNormalization: False
+DescendingM: False
+BufferDepth: Auto
+LMaxFactor: 3