Allow to read in numerical initial data for CurvedScalarWave

src/Evolution/Systems/CurvedScalarWave/Actions/CMakeLists.txt

@@ -0,0 +1,15 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+spectre_target_sources(
+  ${LIBRARY}
+  PRIVATE
+  SetInitialData.cpp
+  )
+
+spectre_target_headers(
+  ${LIBRARY}
+  INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
+  HEADERS
+  SetInitialData.hpp
+  )

src/Evolution/Systems/CurvedScalarWave/Actions/SetInitialData.cpp

@@ -0,0 +1,52 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/CurvedScalarWave/Actions/SetInitialData.hpp"
+
+#include <boost/functional/hash.hpp>
+#include <optional>
+#include <string>
+#include <utility>
+#include <variant>
+
+#include "Utilities/PrettyType.hpp"
+
+namespace CurvedScalarWave {
+
+NumericInitialData::NumericInitialData(
+    std::string file_glob, std::string subfile_name,
+    std::variant<double, importers::ObservationSelector> observation_value,
+    std::optional<double> observation_value_epsilon, bool enable_interpolation,
+    ScalarVars selected_variables
+    )
+    : importer_options_(
+          std::move(file_glob), std::move(subfile_name), observation_value,
+          observation_value_epsilon.value_or(1.0e-12), enable_interpolation),
+      selected_variables_(std::move(selected_variables)) {}
+
+NumericInitialData::NumericInitialData(CkMigrateMessage* msg)
+    : InitialData(msg) {}
+
+PUP::able::PUP_ID NumericInitialData::my_PUP_ID = 0;
+
+size_t NumericInitialData::volume_data_id() const {
+  size_t hash = 0;
+  boost::hash_combine(hash, pretty_type::get_name<NumericInitialData>());
+  boost::hash_combine(hash,
+                      get<importers::OptionTags::FileGlob>(importer_options_));
+  boost::hash_combine(hash,
+                      get<importers::OptionTags::Subgroup>(importer_options_));
+  return hash;
+}
+
+void NumericInitialData::pup(PUP::er& p) {
+  p | importer_options_;
+  p | selected_variables_;
+}
+
+bool operator==(const NumericInitialData& lhs, const NumericInitialData& rhs) {
+  return lhs.importer_options_ == rhs.importer_options_ and
+         lhs.selected_variables_ == rhs.selected_variables_;
+}
+
+}  // namespace CurvedScalarWave

src/Evolution/Systems/CurvedScalarWave/Actions/SetInitialData.hpp

@@ -0,0 +1,321 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+#include <optional>
+#include <pup.h>
+#include <string>
+#include <variant>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
+#include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/Structure/ElementId.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/Initialization/InitialData.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
+#include "IO/Importers/Actions/ReadVolumeData.hpp"
+#include "IO/Importers/ElementDataReader.hpp"
+#include "IO/Importers/Tags.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "Options/String.hpp"
+#include "Parallel/AlgorithmExecution.hpp"
+#include "Parallel/GlobalCache.hpp"
+#include "Parallel/Invoke.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
+#include "Utilities/CallWithDynamicType.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/Serialization/CharmPupable.hpp"
+#include "Utilities/SetNumberOfGridPoints.hpp"
+#include "Utilities/TMPL.hpp"
+#include "Utilities/TaggedTuple.hpp"
+
+/// \cond
+namespace Tags {
+struct Time;
+}  // namespace Tags
+/// \endcond
+
+namespace CurvedScalarWave {
+/*!
+ * \brief Numeric initial data loaded from volume data files
+ */
+class NumericInitialData : public evolution::initial_data::InitialData {
+ public:
+  template <typename Tag>
+  struct VarName {
+    using tag = Tag;
+    static std::string name() { return db::tag_name<Tag>(); }
+    using type = std::string;
+    static constexpr Options::String help =
+        "Name of the variable in the volume data file";
+  };
+
+  // These are the scalar variables that we support loading from volume
+  // data files
+  using all_vars =
+      tmpl::list<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi>;
+  using optional_primitive_vars = tmpl::list<>;
+
+  struct ScalarVars : tuples::tagged_tuple_from_typelist<
+                          db::wrap_tags_in<VarName, all_vars>> {
+    static constexpr Options::String help =
+        "Scalar variables: 'Psi', 'Pi' and 'Phi'.";
+    using options = tags_list;
+    using TaggedTuple::TaggedTuple;
+  };
+
+  // Input-file options
+  struct Variables {
+    using type = ScalarVars;
+    static constexpr Options::String help =
+        "Set of initial data variables from which the Valencia evolution "
+        "variables are computed.";
+  };
+
+  using options =
+      tmpl::push_back<importers::ImporterOptions::tags_list, Variables>;
+
+  static constexpr Options::String help =
+      "Numeric initial data loaded from volume data files";
+
+  NumericInitialData() = default;
+  NumericInitialData(const NumericInitialData& rhs) = default;
+  NumericInitialData& operator=(const NumericInitialData& rhs) = default;
+  NumericInitialData(NumericInitialData&& /*rhs*/) = default;
+  NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
+  ~NumericInitialData() = default;
+
+  /// \cond
+  explicit NumericInitialData(CkMigrateMessage* msg);
+  using PUP::able::register_constructor;
+  WRAPPED_PUPable_decl_template(NumericInitialData);
+  /// \endcond
+
+  std::unique_ptr<evolution::initial_data::InitialData> get_clone()
+      const override {
+    return std::make_unique<NumericInitialData>(*this);
+  }
+
+  NumericInitialData(
+      std::string file_glob, std::string subfile_name,
+      std::variant<double, importers::ObservationSelector> observation_value,
+      std::optional<double> observation_value_epsilon,
+      bool enable_interpolation, ScalarVars selected_variables);
+
+  const importers::ImporterOptions& importer_options() const {
+    return importer_options_;
+  }
+
+  const ScalarVars& selected_variables() const { return selected_variables_; }
+
+  size_t volume_data_id() const;
+
+  template <typename... AllTags>
+  void select_for_import(
+      const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
+    // Select the subset of the available variables that we want to read from
+    // the volume data file
+    using selected_vars = std::decay_t<decltype(selected_variables_)>;
+    tmpl::for_each<typename selected_vars::tags_list>(
+        [&fields, this](const auto tag_v) {
+          using tag = typename std::decay_t<decltype(tag_v)>::type::tag;
+          get<importers::Tags::Selected<tag>>(*fields) =
+              get<VarName<tag>>(selected_variables_);
+        });
+  }
+
+  template <typename... AllTags>
+  void set_initial_data(
+      const gsl::not_null<Scalar<DataVector>*> psi_scalar,
+      const gsl::not_null<Scalar<DataVector>*> pi_scalar,
+      const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar,
+      const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
+      const Mesh<3>& mesh,
+      const InverseJacobian<DataVector, 3, Frame::ElementLogical,
+                            Frame::Inertial>& inv_jacobian) const {
+    *psi_scalar = std::move(get<CurvedScalarWave::Tags::Psi>(*numeric_data));
+    *pi_scalar = std::move(get<CurvedScalarWave::Tags::Pi>(*numeric_data));
+    // Set Phi to the numerical spatial derivative of the scalar
+    partial_derivative(phi_scalar, *psi_scalar, mesh, inv_jacobian);
+  }
+
+  void pup(PUP::er& p) override;
+
+  friend bool operator==(const NumericInitialData& lhs,
+                         const NumericInitialData& rhs);
+
+ private:
+  importers::ImporterOptions importer_options_{};
+  ScalarVars selected_variables_{};
+};
+
+namespace Actions {
+
+/*!
+ * \brief Dispatch loading numeric initial data from files.
+ *
+ * Place this action before
+ * ScalarTensor::Actions::SetNumericInitialData in the action list.
+ * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
+ * invoked by this action.
+ */
+struct SetInitialData {
+  using const_global_cache_tags =
+      tmpl::list<evolution::initial_data::Tags::InitialData>;
+
+  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
+            typename ArrayIndex, typename ActionList,
+            typename ParallelComponent>
+  static Parallel::iterable_action_return_t apply(
+      db::DataBox<DbTagsList>& box,
+      const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
+      Parallel::GlobalCache<Metavariables>& cache,
+      const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
+      const ParallelComponent* const parallel_component) {
+    // Dispatch to the correct `apply` overload based on type of initial data
+    using initial_data_classes =
+        tmpl::at<typename Metavariables::factory_creation::factory_classes,
+                 evolution::initial_data::InitialData>;
+    return call_with_dynamic_type<Parallel::iterable_action_return_t,
+                                  initial_data_classes>(
+        &db::get<evolution::initial_data::Tags::InitialData>(box),
+        [&box, &cache, &parallel_component](const auto* const initial_data) {
+          return apply(make_not_null(&box), *initial_data, cache,
+                       parallel_component);
+        });
+  }
+
+ private:
+  static constexpr size_t Dim = 3;
+
+  // Numeric initial data
+  template <typename DbTagsList, typename Metavariables,
+            typename ParallelComponent>
+  static Parallel::iterable_action_return_t apply(
+      const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
+      const NumericInitialData& initial_data,
+      Parallel::GlobalCache<Metavariables>& cache,
+      const ParallelComponent* const /*meta*/) {
+    // Select the subset of the available variables that we want to read from
+    // the volume data file
+    tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
+        importers::Tags::Selected, NumericInitialData::all_vars>>
+        selected_fields{};
+    initial_data.select_for_import(make_not_null(&selected_fields));
+    // Dispatch loading the variables from the volume data file
+    // - Not using `ckLocalBranch` here to make sure the simple action
+    //   invocation is asynchronous.
+    auto& reader_component = Parallel::get_parallel_component<
+        importers::ElementDataReader<Metavariables>>(cache);
+    Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
+        3, NumericInitialData::all_vars, ParallelComponent>>(
+        reader_component, initial_data.importer_options(),
+        initial_data.volume_data_id(), std::move(selected_fields));
+    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+  }
+
+  // "AnalyticData"-type initial data
+  template <typename DbTagsList, typename InitialData, typename Metavariables,
+            typename ParallelComponent>
+  static Parallel::iterable_action_return_t apply(
+      const gsl::not_null<db::DataBox<DbTagsList>*> box,
+      const InitialData& initial_data,
+      Parallel::GlobalCache<Metavariables>& /*cache*/,
+      const ParallelComponent* const /*meta*/) {
+    // Get scalar variables from analytic data / solution
+    const auto& [coords, mesh, inv_jacobian] = [&box]() {
+      return std::forward_as_tuple(
+          db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box),
+          db::get<domain::Tags::Mesh<Dim>>(*box),
+          db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
+                                                Frame::Inertial>>(*box));
+    }();
+    auto vars = evolution::Initialization::initial_data(
+        initial_data, coords, db::get<::Tags::Time>(*box),
+        tmpl::append<tmpl::list<>,
+                     // Don't use the scalar gradient
+                     tmpl::list<CurvedScalarWave::Tags::Psi,
+                                CurvedScalarWave::Tags::Pi>>{});
+
+    // Move scalar variables and compute gradient
+    db::mutate<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
+               CurvedScalarWave::Tags::Phi<3>>(
+        [&vars](const gsl::not_null<Scalar<DataVector>*> psi_scalar,
+                const gsl::not_null<Scalar<DataVector>*> pi_scalar,
+                const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar,
+                const Mesh<3>& local_mesh,
+                const InverseJacobian<DataVector, 3_st, Frame::ElementLogical,
+                                      Frame::Inertial>& local_inv_jacobian) {
+          *psi_scalar = std::move(get<CurvedScalarWave::Tags::Psi>(vars));
+          *pi_scalar = std::move(get<CurvedScalarWave::Tags::Pi>(vars));
+          // Set Phi to the numerical spatial derivative of the scalar
+          partial_derivative(phi_scalar, *psi_scalar, local_mesh,
+                             local_inv_jacobian);
+        },
+        box, mesh, inv_jacobian);
+
+    // No need to import numeric initial data, so we terminate the phase by
+    // pausing the algorithm on this element
+    return {Parallel::AlgorithmExecution::Pause, std::nullopt};
+  }
+};
+
+/*!
+ * \brief Receive numeric initial data loaded by
+ * CurvedScalarWave::Actions::SetInitialData.
+ */
+struct ReceiveNumericInitialData {
+  static constexpr size_t Dim = 3;
+  using inbox_tags =
+      tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
+
+  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
+            typename ActionList, typename ParallelComponent>
+  static Parallel::iterable_action_return_t apply(
+      db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
+      const Parallel::GlobalCache<Metavariables>& /*cache*/,
+      const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
+      const ParallelComponent* const /*meta*/) {
+    auto& inbox =
+        tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
+            inboxes);
+    const auto& initial_data = dynamic_cast<const NumericInitialData&>(
+        db::get<evolution::initial_data::Tags::InitialData>(box));
+    const size_t volume_data_id = initial_data.volume_data_id();
+    if (inbox.find(volume_data_id) == inbox.end()) {
+      return {Parallel::AlgorithmExecution::Retry, std::nullopt};
+    }
+    auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
+
+    const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
+    const auto& inv_jacobian =
+        db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
+                                              Frame::Inertial>>(box);
+
+    db::mutate<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
+               CurvedScalarWave::Tags::Phi<3>>(
+        [&initial_data, &numeric_data, &mesh, &inv_jacobian](
+            const gsl::not_null<Scalar<DataVector>*> psi_scalar,
+            const gsl::not_null<Scalar<DataVector>*> pi_scalar,
+            const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar) {
+          initial_data.set_initial_data(psi_scalar, pi_scalar, phi_scalar,
+                                        make_not_null(&numeric_data), mesh,
+                                        inv_jacobian);
+        },
+        make_not_null(&box));
+
+    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+  }
+};
+
+}  // namespace Actions
+
+}  // namespace CurvedScalarWave