Merge pull request #6113 from knelli2/dc_time_dep_opts_2

Add several common time dependent options

src/Domain/Creators/TimeDependentOptions/CMakeLists.txt

@@ -5,17 +5,23 @@ spectre_target_sources(
   ${LIBRARY}
   PRIVATE
   BinaryCompactObject.cpp
+  ExpansionMap.cpp
   FromVolumeFile.cpp
+  RotationMap.cpp
   ShapeMap.cpp
   Sphere.cpp
+  TranslationMap.cpp
 )
 
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   BinaryCompactObject.hpp
+  ExpansionMap.hpp
   FromVolumeFile.hpp
+  RotationMap.hpp
   ShapeMap.hpp
   Sphere.hpp
+  TranslationMap.hpp
   )

src/Domain/Creators/TimeDependentOptions/ExpansionMap.cpp

@@ -0,0 +1,90 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp"
+
+#include <array>
+#include <string>
+#include <variant>
+
+#include "DataStructures/DataVector.hpp"
+#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp"
+#include "Options/Context.hpp"
+#include "Options/ParseError.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace domain::creators::time_dependent_options {
+ExpansionMapOptions::ExpansionMapOptions(
+    const std::variant<std::array<double, 3>, FromVolumeFile<names::Expansion>>&
+        expansion_values,
+    const std::variant<std::array<double, 3>, FromVolumeFile<names::Expansion>>&
+        expansion_outer_boundary_values,
+    std::optional<double> decay_timescale_outer_boundary_in,
+    std::optional<double> decay_timescale_in,
+    std::optional<double> asymptotic_velocity_outer_boundary_in,
+    const Options::Context& context)
+    : decay_timescale(std::move(decay_timescale_in)) {
+  const auto set_values =
+      [&, this](const gsl::not_null<std::array<DataVector, 3>*> to_set,
+                const auto& input_values, const bool is_outer_boundary) {
+        if (std::holds_alternative<std::array<double, 3>>(input_values)) {
+          asymptotic_velocity_outer_boundary =
+              asymptotic_velocity_outer_boundary_in;
+          if (decay_timescale.has_value() ==
+              asymptotic_velocity_outer_boundary.has_value()) {
+            PARSE_ERROR(
+                context,
+                "When specifying the ExpansionMap initial outer boundary "
+                "values directly, you must specify one of DecayTimescale or "
+                "AsymptoticVelocityOuterBoundary, but not both.");
+          }
+          auto& values = std::get<std::array<double, 3>>(input_values);
+          for (size_t i = 0; i < to_set->size(); i++) {
+            gsl::at(*to_set, i) = DataVector{1, gsl::at(values, i)};
+          }
+
+          if (is_outer_boundary) {
+            if (not decay_timescale_outer_boundary_in.has_value()) {
+              PARSE_ERROR(context,
+                          "When specifying the ExpansionMap initial outer "
+                          "boundary values directly, you must also specify a "
+                          "'DecayTimescaleOuterBoundary'.");
+            }
+
+            decay_timescale_outer_boundary =
+                decay_timescale_outer_boundary_in.value();
+          }
+
+        } else if (std::holds_alternative<FromVolumeFile<names::Expansion>>(
+                       input_values)) {
+          if (decay_timescale.has_value()) {
+            PARSE_ERROR(
+                context,
+                "When specifying the initial values from a volume file, "
+                "the decay timescale must be 'Auto'.");
+          }
+          auto& values_from_file =
+              std::get<FromVolumeFile<names::Expansion>>(input_values);
+          *to_set = is_outer_boundary
+                        ? values_from_file.expansion_values_outer_boundary
+                        : values_from_file.expansion_values;
+
+          if (is_outer_boundary) {
+            decay_timescale_outer_boundary =
+                decay_timescale_outer_boundary_in.value_or(
+                    values_from_file.decay_timescale_outer_boundary);
+            asymptotic_velocity_outer_boundary =
+                asymptotic_velocity_outer_boundary_in.value_or(
+                    values_from_file.velocity_outer_boundary);
+          }
+        }
+      };
+
+  // Expansion values
+  set_values(make_not_null(&initial_values), expansion_values, false);
+  // Outer boundary
+  set_values(make_not_null(&initial_values_outer_boundary),
+             expansion_outer_boundary_values, true);
+}
+}  // namespace domain::creators::time_dependent_options

src/Domain/Creators/TimeDependentOptions/ExpansionMap.hpp

@@ -0,0 +1,108 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <optional>
+#include <string>
+#include <variant>
+
+#include "DataStructures/DataVector.hpp"
+#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp"
+#include "Options/Auto.hpp"
+#include "Options/Context.hpp"
+#include "Options/String.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace domain::creators::time_dependent_options {
+/*!
+ * \brief Class to be used as an option for initializing expansion map
+ * coefficients.
+ */
+struct ExpansionMapOptions {
+  using type = Options::Auto<ExpansionMapOptions, Options::AutoLabel::None>;
+  static std::string name() { return "ExpansionMap"; }
+  static constexpr Options::String help = {
+      "Options for a time-dependent expansion of the coordinates. Specify "
+      "'None' to not use this map."};
+
+  struct InitialValues {
+    using type =
+        std::variant<std::array<double, 3>, FromVolumeFile<names::Expansion>>;
+    static constexpr Options::String help = {
+        "Initial values for the expansion map, its velocity and "
+        "acceleration."};
+  };
+
+  struct InitialValuesOuterBoundary {
+    using type =
+        std::variant<std::array<double, 3>, FromVolumeFile<names::Expansion>>;
+    static constexpr Options::String help = {
+        "Initial values for the expansion map, its velocity and "
+        "acceleration at the outer boundary. Unless you are starting from a "
+        "checkpoint or continuing an evolution, this option should likely be "
+        "[1.0, 0.0, 0.0] at the start of an evolution"};
+  };
+
+  struct DecayTimescaleOuterBoundary {
+    using type = Options::Auto<double>;
+    static constexpr Options::String help = {
+        "A timescale for how fast the outer boundary expansion approaches its "
+        "asymptotic value. Can optionally specify 'Auto' when reading the "
+        "initial values 'FromVolumeFile' to use the decay timescale from the "
+        "function of time in the volume file. Cannot specify 'Auto' when "
+        "initial values are specified directly."};
+  };
+
+  struct DecayTimescale {
+    using type = Options::Auto<double>;
+    static constexpr Options::String help = {
+        "If specified, a SettleToConstant function of time will be used for "
+        "the expansion map and this number will determine the timescale that "
+        "the expansion approaches its asymptotic value. If 'Auto' is "
+        "specified, a PiecewisePolynomial function of time will be used for "
+        "the expansion map. Note that if you are reading the initial values "
+        "from a volume file, you must specify 'Auto' for this option."};
+  };
+
+  struct AsymptoticVelocityOuterBoundary {
+    using type = Options::Auto<double>;
+    static constexpr Options::String help = {
+        "There are two choices for this option. If a value is specified, a "
+        "FixedSpeedCubic function of time will be used for the expansion map "
+        "at the outer boundary and this number will determine its velocity. If "
+        "'Auto' is specified, the behavior will depend on what is chosen for "
+        "'InitialValuesOuterBoundary'. If values are specified for "
+        "'InitialValuesOuterBoundary', then 'Auto' here means a "
+        "SettleToConstant function of time will be used for the expansion map "
+        "at the outer boundary. If 'FromVolumeFile' is specified for "
+        "'InitialValuesOuterBoundary', then a FixedSpeedCubic function of time "
+        "will be used and the velocity from the function of "
+        "time in the volume file will be used."};
+  };
+
+  using options = tmpl::list<InitialValues, InitialValuesOuterBoundary,
+                             DecayTimescaleOuterBoundary, DecayTimescale,
+                             AsymptoticVelocityOuterBoundary>;
+
+  ExpansionMapOptions() = default;
+  ExpansionMapOptions(
+      const std::variant<std::array<double, 3>,
+                         FromVolumeFile<names::Expansion>>& expansion_values,
+      const std::variant<std::array<double, 3>,
+                         FromVolumeFile<names::Expansion>>&
+          expansion_outer_boundary_values,
+      std::optional<double> decay_timescale_outer_boundary_in,
+      std::optional<double> decay_timescale_in,
+      std::optional<double> asymptotic_velocity_outer_boundary_in,
+      const Options::Context& context = {});
+
+  std::array<DataVector, 3> initial_values{};
+  std::array<DataVector, 3> initial_values_outer_boundary{};
+  double decay_timescale_outer_boundary{};
+  std::optional<double> decay_timescale{};
+  std::optional<double> asymptotic_velocity_outer_boundary{};
+};
+}  // namespace domain::creators::time_dependent_options

src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp

@@ -12,6 +12,7 @@
 #include <vector>
 
 #include "DataStructures/DataVector.hpp"
+#include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp"
 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
 #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp"
 #include "Domain/FunctionsOfTime/QuaternionHelpers.hpp"
@@ -110,6 +111,20 @@ FromVolumeFile<names::Expansion>::FromVolumeFile(
   expansion_values = exp_function_of_time->func_and_2_derivs(time);
   expansion_values_outer_boundary =
       exp_outer_boundary_function_of_time->func_and_2_derivs(time);
+
+  const auto* fixed_speed_cubic =
+      dynamic_cast<domain::FunctionsOfTime::FixedSpeedCubic*>(
+          exp_outer_boundary_function_of_time.get());
+
+  if (fixed_speed_cubic == nullptr) {
+    PARSE_ERROR(
+        context,
+        "When reading the Expansion map parameters from a volume file, the "
+        "ExpansionOuterBoundary function of time must be a FixedSpeedCubic.");
+  }
+
+  velocity_outer_boundary = fixed_speed_cubic->velocity();
+  decay_timescale_outer_boundary = fixed_speed_cubic->decay_timescale();
 }
 
 FromVolumeFile<names::Rotation>::FromVolumeFile(

src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp

@@ -94,6 +94,8 @@ struct FromVolumeFile<names::Expansion> : public detail::FromVolumeFileBase {
 
   std::array<DataVector, 3> expansion_values{};
   std::array<DataVector, 3> expansion_values_outer_boundary{};
+  double velocity_outer_boundary{};
+  double decay_timescale_outer_boundary{};
 };
 
 template <>

src/Domain/Creators/TimeDependentOptions/RotationMap.cpp

@@ -0,0 +1,105 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp"
+
+#include <array>
+#include <optional>
+#include <string>
+#include <utility>
+#include <variant>
+
+#include "DataStructures/DataVector.hpp"
+#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp"
+#include "Options/Context.hpp"
+#include "Options/ParseError.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/MakeArray.hpp"
+
+namespace domain::creators::time_dependent_options {
+template <size_t NumDerivs>
+RotationMapOptions<NumDerivs>::RotationMapOptions(
+    std::variant<std::vector<std::array<double, 4>>,
+                 FromVolumeFile<names::Rotation>>
+        initial_quaternions,
+    std::optional<std::vector<std::array<double, 3>>> initial_angles,
+    std::optional<double> decay_timescale_in, const Options::Context& context)
+    : decay_timescale(decay_timescale_in) {
+  quaternions = make_array<NumDerivs + 1, DataVector>(DataVector{4, 0.0});
+  angles = make_array<NumDerivs + 1, DataVector>(DataVector{3, 0.0});
+
+  if (std::holds_alternative<std::vector<std::array<double, 4>>>(
+          initial_quaternions)) {
+    auto& values =
+        std::get<std::vector<std::array<double, 4>>>(initial_quaternions);
+    if (values.empty() or values.size() > quaternions.size()) {
+      PARSE_ERROR(
+          context,
+          "Must specify at least the value of the quaternion, and optionally "
+          "up to "
+              << NumDerivs << " time derivatives.");
+    }
+    for (size_t i = 0; i < values.size(); i++) {
+      gsl::at(quaternions, i) =
+          DataVector{values[i][0], values[i][1], values[i][2], values[i][3]};
+    }
+
+    if (initial_angles.has_value()) {
+      auto& angle_values = initial_angles.value();
+      if (angle_values.empty() or angle_values.size() > angles.size()) {
+        PARSE_ERROR(
+            context,
+            "When specifying the angle, you must specify at least the value, "
+            "and optionally up to "
+                << NumDerivs << " time derivatives.");
+      }
+      for (size_t i = 0; i < angle_values.size(); i++) {
+        gsl::at(angles, i) = DataVector{angle_values[i][0], angle_values[i][1],
+                                        angle_values[i][2]};
+      }
+    }
+  } else if (std::holds_alternative<FromVolumeFile<names::Rotation>>(
+                 initial_quaternions)) {
+    if (decay_timescale.has_value()) {
+      PARSE_ERROR(context,
+                  "When specifying the initial quaternions from a volume file, "
+                  "the decay timescale must be 'Auto'.");
+    }
+    auto& values_from_file =
+        std::get<FromVolumeFile<names::Rotation>>(initial_quaternions);
+
+    for (size_t i = 0; i < values_from_file.quaternions.size(); i++) {
+      gsl::at(quaternions, i) = gsl::at(values_from_file.quaternions, i);
+      gsl::at(angles, i) = gsl::at(values_from_file.angle_values, i);
+    }
+
+    if (initial_angles.has_value()) {
+      // Reset angle func so derivs that weren't specified are zero
+      angles = make_array<NumDerivs + 1, DataVector>(DataVector{3, 0.0});
+      auto& angle_values = initial_angles.value();
+      if (angle_values.empty() or angle_values.size() > angles.size()) {
+        PARSE_ERROR(
+            context,
+            "When specifying the angle, you must specify at least the value, "
+            "and optionally up to "
+                << NumDerivs << " time derivatives.");
+      }
+      for (size_t i = 0; i < angle_values.size(); i++) {
+        gsl::at(angles, i) = DataVector{angle_values[i][0], angle_values[i][1],
+                                        angle_values[i][2]};
+      }
+    }
+  }
+}
+
+#define NUMDERIVS(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data) \
+  template struct RotationMapOptions<NUMDERIVS(data)>;
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (2, 3))
+
+#undef INSTANTIATE
+#undef NUMDERIVS
+}  // namespace domain::creators::time_dependent_options