Merge pull request #6284 from knelli2/strahlkorper_coords_txt_file

Add executable to write CCE worldtube coords to a file

src/Evolution/Systems/Cce/Callbacks/DumpBondiSachsOnWorldtube.hpp

@@ -24,13 +24,13 @@
 #include "IO/Observer/Actions/GetLockPointer.hpp"
 #include "IO/Observer/ObserverComponent.hpp"
 #include "IO/Observer/ReductionActions.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
 #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCoefficients.hpp"
 #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
 #include "Parallel/GlobalCache.hpp"
 #include "Parallel/Invoke.hpp"
 #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
 #include "ParallelAlgorithms/Interpolation/Protocols/PostInterpolationCallback.hpp"
-#include "ParallelAlgorithms/Interpolation/Targets/AngularOrdering.hpp"
 #include "ParallelAlgorithms/Interpolation/Targets/Sphere.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 #include "Utilities/ConstantExpressions.hpp"
@@ -133,7 +133,7 @@ struct DumpBondiSachsOnWorldtube
         Parallel::get<Tags::Sphere<InterpolationTargetTag>>(cache);
     const auto& filename_prefix = Parallel::get<Cce::Tags::FilePrefix>(cache);
 
-    if (sphere.angular_ordering != intrp::AngularOrdering::Cce) {
+    if (sphere.angular_ordering != ylm::AngularOrdering::Cce) {
       ERROR(
           "To use the DumpBondiSachsOnWorldtube post interpolation callback, "
           "the angular ordering of the Spheres must be Cce, not "

src/Executables/CMakeLists.txt

@@ -10,3 +10,4 @@ add_subdirectory(ExportCoordinates)
 add_subdirectory(ParallelInfo)
 add_subdirectory(ReduceCceWorldtube)
 add_subdirectory(TimeStepperSummary)
+add_subdirectory(WriteCceWorldtubeCoordsToFile)

src/Executables/WriteCceWorldtubeCoordsToFile/CMakeLists.txt

@@ -0,0 +1,24 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(EXECUTABLE WriteCceWorldtubeCoordsToFile)
+
+add_spectre_executable(
+  ${EXECUTABLE}
+  EXCLUDE_FROM_ALL
+  WriteCceWorldtubeCoordsToFile.cpp
+  )
+
+target_link_libraries(
+  ${EXECUTABLE}
+  PRIVATE
+  Boost::boost
+  Boost::program_options
+  Printf
+  SphericalHarmonics
+  SphericalHarmonicsIO
+  )
+
+if(BUILD_TESTING)
+  add_dependencies(test-executables ${EXECUTABLE})
+endif()

src/Executables/WriteCceWorldtubeCoordsToFile/WriteCceWorldtubeCoordsToFile.cpp

@@ -0,0 +1,79 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include <array>
+#include <boost/program_options.hpp>
+#include <boost/program_options/value_semantic.hpp>
+#include <cstddef>
+#include <exception>
+#include <string>
+
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.hpp"
+#include "Parallel/Printf/Printf.hpp"
+
+// Charm looks for this function but since we build without a main function or
+// main module we just have it be empty
+extern "C" void CkRegisterMainModule(void) {}
+
+/*
+ * This executable is used for writing the collocation points of a CCE
+ * Strahlkorper to a text file. This is useful for other NR codes so they can
+ * write data to certain points and we do the necessary transformations.
+ */
+int main(int argc, char** argv) {
+  boost::program_options::options_description desc(
+      "Program for writing the collocation points (coordinates) of a worldtube "
+      "sphere that SpECTRE CCE is able to read and interpret to a "
+      "file.\nDetails about the output file:\n"
+      " * There are no header or comment lines\n"
+      " * Each point is written to a new line of the output file\n"
+      " * The delimiter for the x, y, z components of each point is a space\n"
+      " * The points are written in scientific notation\n"
+      " * The sphere is centered on the origin (0, 0, 0)");
+  desc.add_options()("help,h", "show this help message")(
+      "radius,r", boost::program_options::value<double>()->required(),
+      "Radius of the worltube.")(
+      "lmax,L", boost::program_options::value<size_t>()->required(),
+      "The spherical harmonic L of the surface. The surface will have (L + 1) "
+      "* (2L + 1) total points")(
+      "output_file,o", boost::program_options::value<std::string>()->required(),
+      "Output filename for the points. No extension will be added.")(
+      "force,f", boost::program_options::bool_switch(),
+      "Overwrite the output file if it already exists");
+
+  boost::program_options::variables_map vars;
+
+  boost::program_options::store(
+      boost::program_options::command_line_parser(argc, argv)
+          .options(desc)
+          .run(),
+      vars);
+
+  if (vars.contains("help")) {
+    Parallel::printf("%s\n", desc);
+    return 0;
+  }
+
+  for (const auto& option : {"radius", "lmax", "output_file"}) {
+    if (not vars.contains(option)) {
+      Parallel::printf("Missing option: %s\n\n%s", option, desc);
+      return 1;
+    }
+  }
+
+  const double radius = vars["radius"].as<double>();
+  const size_t l_max = vars["lmax"].as<size_t>();
+  const std::array<double, 3> center{0.0, 0.0, 0.0};
+  const std::string output_file = vars["output_file"].as<std::string>();
+  const bool overwrite = vars["force"].as<bool>();
+
+  try {
+    ylm::write_strahlkorper_coords_to_text_file(
+        radius, l_max, center, output_file, ylm::AngularOrdering::Cce,
+        overwrite);
+  } catch (const std::exception& exception) {
+    Parallel::printf("%s\n", exception.what());
+    return 1;
+  }
+}

src/ParallelAlgorithms/Interpolation/Targets/AngularOrdering.cpp → src/NumericalAlgorithms/SphericalHarmonics/AngularOrdering.cpp

@@ -1,15 +1,15 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
-#include "ParallelAlgorithms/Interpolation/Targets/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
 
 #include <string>
 
 #include "Options/Options.hpp"
 #include "Options/ParseOptions.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
 
-namespace intrp {
+namespace ylm {
 std::ostream& operator<<(std::ostream& os, const AngularOrdering ordering) {
   switch (ordering) {
     case AngularOrdering::Strahlkorper:
@@ -20,17 +20,17 @@ std::ostream& operator<<(std::ostream& os, const AngularOrdering ordering) {
       ERROR("Unknown AngularOrdering type");
   }
 }
-}  // namespace intrp
+}  // namespace ylm
 
 template <>
-intrp::AngularOrdering
-Options::create_from_yaml<intrp::AngularOrdering>::create<void>(
+ylm::AngularOrdering
+Options::create_from_yaml<ylm::AngularOrdering>::create<void>(
     const Options::Option& options) {
   const auto ordering = options.parse_as<std::string>();
   if (ordering == "Strahlkorper") {
-    return intrp::AngularOrdering::Strahlkorper;
+    return ylm::AngularOrdering::Strahlkorper;
   } else if (ordering == "Cce") {
-    return intrp::AngularOrdering::Cce;
+    return ylm::AngularOrdering::Cce;
   }
   PARSE_ERROR(options.context(),
               "AngularOrdering must be 'Strahlkorper' or 'Cce'");

src/ParallelAlgorithms/Interpolation/Targets/AngularOrdering.hpp → src/NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp

@@ -13,7 +13,7 @@ struct create_from_yaml;
 }  // namespace Options
 /// \endcond
 
-namespace intrp {
+namespace ylm {
 /*!
  * \brief Label for the ordering of spherical harmonic points on a sphere
  *
@@ -23,17 +23,17 @@ namespace intrp {
  */
 enum class AngularOrdering { Strahlkorper, Cce };
 std::ostream& operator<<(std::ostream& os, AngularOrdering ordering);
-}  // namespace intrp
+}  // namespace ylm
 
 template <>
-struct Options::create_from_yaml<intrp::AngularOrdering> {
+struct Options::create_from_yaml<ylm::AngularOrdering> {
   template <typename Metavariables>
-  static intrp::AngularOrdering create(const Options::Option& options) {
+  static ylm::AngularOrdering create(const Options::Option& options) {
     return create<void>(options);
   }
 };
 
 template <>
-intrp::AngularOrdering
-Options::create_from_yaml<intrp::AngularOrdering>::create<void>(
+ylm::AngularOrdering
+Options::create_from_yaml<ylm::AngularOrdering>::create<void>(
     const Options::Option& options);

src/NumericalAlgorithms/SphericalHarmonics/CMakeLists.txt

@@ -8,6 +8,7 @@ add_spectre_library(${LIBRARY})
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
+  AngularOrdering.cpp
   ChangeCenterOfStrahlkorper.cpp
   RealSphericalHarmonics.cpp
   SpherepackIterator.cpp
@@ -23,6 +24,7 @@ spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
+  AngularOrdering.hpp
   ChangeCenterOfStrahlkorper.hpp
   RealSphericalHarmonics.hpp
   SpherepackIterator.hpp

src/NumericalAlgorithms/SphericalHarmonics/IO/CMakeLists.txt

@@ -10,6 +10,7 @@ spectre_target_sources(
   PRIVATE
   FillYlmLegendAndData.cpp
   ReadSurfaceYlm.cpp
+  StrahlkorperCoordsToTextFile.cpp
   )
 
 spectre_target_headers(
@@ -18,6 +19,7 @@ spectre_target_headers(
   HEADERS
   FillYlmLegendAndData.hpp
   ReadSurfaceYlm.hpp
+  StrahlkorperCoordsToTextFile.hpp
   )
 
 target_link_libraries(

src/NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.cpp

@@ -0,0 +1,95 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.hpp"
+
+#include <fstream>
+#include <iomanip>
+#include <limits>
+#include <ostream>
+#include <string>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Transpose.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/StrahlkorperFunctions.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/FileSystem.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+
+namespace Frame {
+struct Inertial;
+struct Distorted;
+struct Grid;
+}  // namespace Frame
+
+namespace ylm {
+template <typename Frame>
+void write_strahlkorper_coords_to_text_file(
+    const Strahlkorper<Frame>& strahlkorper,
+    const std::string& output_file_name, const AngularOrdering ordering,
+    const bool overwrite_file) {
+  if (not overwrite_file and
+      file_system::check_if_file_exists(output_file_name)) {
+    ERROR_NO_TRACE("The output file " << output_file_name
+                                      << " already exists.");
+  }
+
+  tnsr::I<DataVector, 3, Frame> cartesian_coords =
+      ylm::cartesian_coords(strahlkorper);
+
+  // Cce expects coordinates in a different order than a typical Strahlkorper
+  if (ordering == AngularOrdering::Cce) {
+    const auto physical_extents =
+        strahlkorper.ylm_spherepack().physical_extents();
+    auto transposed_coords =
+        tnsr::I<DataVector, 3, Frame>(get<0>(cartesian_coords).size());
+    for (size_t i = 0; i < 3; ++i) {
+      transpose(make_not_null(&transposed_coords.get(i)),
+                cartesian_coords.get(i), physical_extents[0],
+                physical_extents[1]);
+    }
+
+    cartesian_coords = std::move(transposed_coords);
+  }
+
+  std::ofstream output_file(output_file_name);
+  output_file << std::fixed
+              << std::setprecision(std::numeric_limits<double>::digits10 + 4)
+              << std::scientific;
+
+  const size_t num_points = get<0>(cartesian_coords).size();
+  for (size_t i = 0; i < num_points; i++) {
+    output_file << get<0>(cartesian_coords)[i] << " "
+                << get<1>(cartesian_coords)[i] << " "
+                << get<2>(cartesian_coords)[i] << std::endl;
+  }
+}
+
+void write_strahlkorper_coords_to_text_file(const double radius,
+                                            const size_t l_max,
+                                            const std::array<double, 3>& center,
+                                            const std::string& output_file_name,
+                                            const AngularOrdering ordering,
+                                            const bool overwrite_file) {
+  const Strahlkorper<Frame::Inertial> strahlkorper{l_max, radius, center};
+  write_strahlkorper_coords_to_text_file(strahlkorper, output_file_name,
+                                         ordering, overwrite_file);
+}
+
+#define FRAME(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                                               \
+  template void write_strahlkorper_coords_to_text_file(                    \
+      const Strahlkorper<FRAME(data)>& strahlkorper,                       \
+      const std::string& output_file_name, const AngularOrdering ordering, \
+      const bool overwrite_file);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE,
+                        (Frame::Grid, Frame::Distorted, Frame::Inertial))
+
+#undef INSTANTIATE
+#undef FRAME
+}  // namespace ylm

src/NumericalAlgorithms/SphericalHarmonics/IO/StrahlkorperCoordsToTextFile.hpp

@@ -0,0 +1,49 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <string>
+
+#include "NumericalAlgorithms/SphericalHarmonics/AngularOrdering.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+
+namespace ylm {
+/// @{
+/*!
+ * \brief Writes the collocation points of a `ylm::Strahlkorper` to an output
+ * text file.
+ *
+ * \details The ordering of the points can be either the typical
+ * `ylm::Spherepack` ordering or CCE ordering that works with libsharp. Also, an
+ * error will occur if the output file already exists, but the output file can
+ * be overwritten with the \p overwrite_file option.
+ *
+ * The second overload will construct a spherical `ylm::Strahlkorper` with the
+ * given \p radius, \p l_max, and \p center.
+ *
+ * The output file format will be as follows with no comment or header lines,
+ * a space as the delimiter, and decimals written in scientific notation:
+ *
+ * ```
+ * x0 y0 z0
+ * x1 y1 z1
+ * x2 y2 z2
+ * ...
+ * ```
+ */
+template <typename Frame>
+void write_strahlkorper_coords_to_text_file(
+    const Strahlkorper<Frame>& strahlkorper,
+    const std::string& output_file_name, AngularOrdering ordering,
+    bool overwrite_file = false);
+
+void write_strahlkorper_coords_to_text_file(double radius, size_t l_max,
+                                            const std::array<double, 3>& center,
+                                            const std::string& output_file_name,
+                                            AngularOrdering ordering,
+                                            bool overwrite_file = false);
+/// @}
+}  // namespace ylm