Merge pull request #6430 from knelli2/no_pi_phi_constraints

Allow choice of reading Pi and Phi from numeric ID

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

@@ -54,7 +54,8 @@ void initial_gh_variables_from_adm(
 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,
+    std::optional<double> observation_value_epsilon,
+    const bool enable_interpolation,
     std::variant<AdmVars, GhVars> selected_variables)
     : importer_options_(
           std::move(file_glob), std::move(subfile_name), observation_value,

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

@@ -9,11 +9,13 @@
 #include <variant>
 
 #include "DataStructures/DataBox/DataBox.hpp"
+#include "DataStructures/DataBox/Tag.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Tags.hpp"
 #include "Evolution/Initialization/InitialData.hpp"
+#include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
 #include "IO/Importers/Actions/ReadVolumeData.hpp"
 #include "IO/Importers/ElementDataReader.hpp"
@@ -96,21 +98,21 @@ class NumericInitialData : public evolution::initial_data::InitialData {
         "ADM variables: 'Lapse', 'Shift', 'SpatialMetric' and "
         "'ExtrinsicCurvature'. The initial GH variables will be computed "
         "from these numeric fields, as well as their numeric spatial "
-        "derivatives on the computational grid.";
+        "derivatives on the computational grid. The GH variable Pi will be set "
+        "to satisfy the gauge constraint using the evolution gauge. The GH "
+        "variable Phi will be set to satisfy the 3-index constraint.";
     using options = tags_list;
     using TaggedTuple::TaggedTuple;
   };
 
   // - Generalized harmonic variables
   using gh_vars = tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>,
-                             Tags::Pi<DataVector, 3>>;
+                             Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>;
   struct GhVars
       : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<VarName, gh_vars>> {
     static constexpr Options::String help =
-        "GH variables: 'SpacetimeMetric' and 'Pi'. These variables are "
-        "used to set the initial data directly; Phi is then set to the "
-        "numerical derivative of SpacetimeMetric, to enforce the 3-index "
-        "constraint.";
+        "GH variables: 'SpacetimeMetric', 'Pi', and 'Phi'. These variables are "
+        "used to set the initial data directly.";
     using options = tags_list;
     using TaggedTuple::TaggedTuple;
   };
@@ -232,8 +234,7 @@ class NumericInitialData : public evolution::initial_data::InitialData {
       *spacetime_metric = std::move(
           get<gr::Tags::SpacetimeMetric<DataVector, 3>>(*numeric_data));
       *pi = std::move(get<Tags::Pi<DataVector, 3>>(*numeric_data));
-      // Set Phi to the numerical spatial derivative of spacetime_metric
-      partial_derivative(phi, *spacetime_metric, mesh, inv_jacobian);
+      *phi = get<Tags::Phi<DataVector, 3>>(*numeric_data);
     } else if (std::holds_alternative<NumericInitialData::AdmVars>(
                    selected_variables_)) {
       // We have loaded ADM variables from the file. Convert to GH variables.
@@ -288,7 +289,7 @@ struct SetInitialData {
       db::DataBox<DbTagsList>& box,
       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
       Parallel::GlobalCache<Metavariables>& cache,
-      const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
+      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 =
@@ -297,21 +298,33 @@ struct SetInitialData {
     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,
+        [&box, &cache, &array_index,
+         &parallel_component](const auto* const initial_data) {
+          return apply(make_not_null(&box), *initial_data, cache, array_index,
                        parallel_component);
         });
   }
 
  private:
   // Numeric initial data
-  template <typename DbTagsList, typename Metavariables,
+  template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
             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*/) {
+      const ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
+    // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
+    // we are reading them in. Also we only need to mutate this tag once so do
+    // it on the first element.
+    if (is_zeroth_element(array_index) and
+        std::holds_alternative<NumericInitialData::GhVars>(
+            initial_data.selected_variables())) {
+      Parallel::mutate<Tags::SetPiAndPhiFromConstraints,
+                       gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
+          cache, false);
+    }
+
     // 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<
@@ -332,11 +345,12 @@ struct SetInitialData {
 
   // "AnalyticData"-type initial data
   template <typename DbTagsList, typename InitialData, typename Metavariables,
-            typename ParallelComponent>
+            typename ArrayIndex, 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 ArrayIndex& /*array_index*/,
       const ParallelComponent* const /*meta*/) {
     static constexpr size_t Dim = Metavariables::volume_dim;
 

src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.cpp

@@ -33,6 +33,11 @@
 #include "Utilities/GenerateInstantiations.hpp"
 
 namespace gh::gauges {
+void SetPiAndPhiFromConstraintsCacheMutator::apply(
+    const gsl::not_null<bool*> value, const bool new_value) {
+  *value = new_value;
+}
+
 template <size_t Dim>
 void SetPiAndPhiFromConstraints<Dim>::apply(
     const gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*> pi,
@@ -46,7 +51,11 @@ void SetPiAndPhiFromConstraints<Dim>::apply(
         functions_of_time,
     const tnsr::I<DataVector, Dim, Frame::ElementLogical>& logical_coordinates,
     const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
-    const gauges::GaugeCondition& gauge_condition) {
+    const gauges::GaugeCondition& gauge_condition,
+    const bool set_pi_and_phi_from_constraints) {
+  if (not set_pi_and_phi_from_constraints) {
+    return;
+  }
   const auto grid_coords = logical_to_grid_map(logical_coordinates);
   const auto inv_jac_logical_to_grid =
       logical_to_grid_map.inv_jacobian(logical_coordinates);

src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp

@@ -21,6 +21,7 @@
 #include "Evolution/Initialization/Tags.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Tags/GaugeCondition.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
+#include "Parallel/GlobalCache.hpp"
 #include "Utilities/TMPL.hpp"
 
 /// \cond
@@ -29,7 +30,29 @@ struct Time;
 }  // namespace Tags
 /// \endcond
 
-namespace gh::gauges {
+namespace gh {
+namespace Tags {
+/// DataBox tag for holding whether or not to set GH variables $\Pi$ and $\Phi$
+/// from constraints.
+struct SetPiAndPhiFromConstraints : db::SimpleTag {
+  using type = bool;
+
+  using option_tags = tmpl::list<>;
+  static constexpr bool pass_metavariables = false;
+
+  static bool create_from_options() { return true; }
+};
+}  // namespace Tags
+
+namespace gauges {
+/*!
+ * \brief GlobalCache mutator to set the value of the
+ * `gh::Tags::SetPiAndPhiFromConstraints` tag.
+ */
+struct SetPiAndPhiFromConstraintsCacheMutator {
+  static void apply(gsl::not_null<bool*> value, bool new_value);
+};
+
 /*!
  * \brief Set \f$\Pi_{ab}\f$ from the gauge source function (or 1-index
  * constraint) and \f$\Phi_{iab}\f$ from the 3-index constraint.
@@ -50,9 +73,14 @@ struct SetPiAndPhiFromConstraints {
                  domain::Tags::FunctionsOfTime,
                  domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
                  gr::Tags::SpacetimeMetric<DataVector, Dim>,
-                 gh::gauges::Tags::GaugeCondition>;
+                 gh::gauges::Tags::GaugeCondition,
+                 gh::Tags::SetPiAndPhiFromConstraints>;
 
+  using compute_tags = tmpl::list<
+      Parallel::Tags::FromGlobalCache<gh::Tags::SetPiAndPhiFromConstraints>>;
   using const_global_cache_tags = tmpl::list<gh::gauges::Tags::GaugeCondition>;
+  using mutable_global_cache_tags =
+      tmpl::list<gh::Tags::SetPiAndPhiFromConstraints>;
 
   static void apply(
       gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*> pi,
@@ -68,6 +96,8 @@ struct SetPiAndPhiFromConstraints {
       const tnsr::I<DataVector, Dim, Frame::ElementLogical>&
           logical_coordinates,
       const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
-      const gauges::GaugeCondition& gauge_condition);
+      const gauges::GaugeCondition& gauge_condition,
+      bool set_pi_and_phi_from_constraints = true);
 };
-}  // namespace gh::gauges
+}  // namespace gauges
+}  // namespace gh

src/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/SetInitialData.hpp

@@ -20,6 +20,7 @@
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
 #include "Evolution/NumericInitialData.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
+#include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Actions/NumericInitialData.hpp"
 #include "IO/Importers/Actions/ReadVolumeData.hpp"
@@ -186,7 +187,7 @@ struct SetInitialData {
       db::DataBox<DbTagsList>& box,
       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
       Parallel::GlobalCache<Metavariables>& cache,
-      const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
+      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 =
@@ -195,8 +196,9 @@ struct SetInitialData {
     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,
+        [&box, &cache, &array_index,
+         &parallel_component](const auto* const initial_data) {
+          return apply(make_not_null(&box), *initial_data, cache, array_index,
                        parallel_component);
         });
   }
@@ -205,13 +207,23 @@ struct SetInitialData {
   static constexpr size_t Dim = 3;
 
   // Numeric initial data
-  template <typename DbTagsList, typename Metavariables,
+  template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
             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*/) {
+      const ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
+    // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
+    // we are reading them in. Also we only need to mutate this tag once so do
+    // it on the first element.
+    if (is_zeroth_element(array_index) and
+        std::holds_alternative<gh::NumericInitialData::GhVars>(
+            initial_data.gh_numeric_id().selected_variables())) {
+      Parallel::mutate<gh::Tags::SetPiAndPhiFromConstraints,
+                       gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
+          cache, false);
+    }
     // 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<
@@ -232,11 +244,12 @@ struct SetInitialData {
 
   // "AnalyticData"-type initial data
   template <typename DbTagsList, typename InitialData, typename Metavariables,
-            typename ParallelComponent>
+            typename ArrayIndex, 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 ArrayIndex& /*array_index*/,
       const ParallelComponent* const /*meta*/) {
     // Get ADM + hydro variables from analytic data / solution
     const auto& [coords, mesh, inv_jacobian] = [&box]() {

src/Evolution/Systems/GrMhd/GhValenciaDivClean/SetPiAndPhiFromConstraints.hpp

@@ -45,27 +45,27 @@ namespace grmhd::GhValenciaDivClean {
 struct SetPiAndPhiFromConstraints {
  public:
   using return_tags =
-      tmpl::list<gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>>;
-  using argument_tags = tmpl::list<
-      ::Tags::Time, domain::Tags::Mesh<3>,
+      typename gh::gauges::SetPiAndPhiFromConstraints<3>::return_tags;
+
+  using argument_tags = tmpl::push_back<
+      typename gh::gauges::SetPiAndPhiFromConstraints<3>::argument_tags,
       evolution::dg::subcell::Tags::Mesh<3>,
-      domain::Tags::ElementMap<3, Frame::Grid>,
-      domain::CoordinateMaps::Tags::CoordinateMap<3, Frame::Grid,
-                                                  Frame::Inertial>,
-      domain::Tags::FunctionsOfTime,
-      domain::Tags::Coordinates<3, Frame::ElementLogical>,
       evolution::dg::subcell::Tags::Coordinates<3, Frame::ElementLogical>,
-      gr::Tags::SpacetimeMetric<DataVector, 3>,
-      gh::gauges::Tags::GaugeCondition,
       evolution::dg::subcell::Tags::ActiveGrid>;
 
-  using const_global_cache_tags = tmpl::list<gh::gauges::Tags::GaugeCondition>;
+  using compute_tags =
+      typename gh::gauges::SetPiAndPhiFromConstraints<3>::compute_tags;
+  using const_global_cache_tags =
+      typename gh::gauges::SetPiAndPhiFromConstraints<
+          3>::const_global_cache_tags;
+  using mutable_global_cache_tags =
+      typename gh::gauges::SetPiAndPhiFromConstraints<
+          3>::mutable_global_cache_tags;
 
   static void apply(
       const gsl::not_null<tnsr::aa<DataVector, 3, Frame::Inertial>*> pi,
       const gsl::not_null<tnsr::iaa<DataVector, 3, Frame::Inertial>*> phi,
       const double initial_time, const Mesh<3>& dg_mesh,
-      const Mesh<3>& subcell_mesh,
       const ElementMap<3, Frame::Grid>& logical_to_grid_map,
       const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, 3>&
           grid_to_inertial_map,
@@ -75,21 +75,22 @@ struct SetPiAndPhiFromConstraints {
           functions_of_time,
       const tnsr::I<DataVector, 3, Frame::ElementLogical>&
           dg_logical_coordinates,
-      const tnsr::I<DataVector, 3, Frame::ElementLogical>&
-          subcell_logical_coordinates,
       const tnsr::aa<DataVector, 3, Frame::Inertial>& spacetime_metric,
       const gh::gauges::GaugeCondition& gauge_condition,
+      const bool set_pi_and_phi_from_constraints, const Mesh<3>& subcell_mesh,
+      const tnsr::I<DataVector, 3, Frame::ElementLogical>&
+          subcell_logical_coordinates,
       const evolution::dg::subcell::ActiveGrid active_grid) {
     if (active_grid == evolution::dg::subcell::ActiveGrid::Dg) {
       gh::gauges::SetPiAndPhiFromConstraints<3>::apply(
           pi, phi, initial_time, dg_mesh, logical_to_grid_map,
           grid_to_inertial_map, functions_of_time, dg_logical_coordinates,
-          spacetime_metric, gauge_condition);
+          spacetime_metric, gauge_condition, set_pi_and_phi_from_constraints);
     } else {
       gh::gauges::SetPiAndPhiFromConstraints<3>::apply(
           pi, phi, initial_time, subcell_mesh, logical_to_grid_map,
           grid_to_inertial_map, functions_of_time, subcell_logical_coordinates,
-          spacetime_metric, gauge_condition);
+          spacetime_metric, gauge_condition, set_pi_and_phi_from_constraints);
     }
   }
 };

src/ParallelAlgorithms/Actions/MutateApply.hpp

@@ -9,6 +9,7 @@
 #include "DataStructures/DataBox/DataBox.hpp"
 #include "Parallel/AlgorithmExecution.hpp"
 #include "Utilities/TMPL.hpp"
+#include "Utilities/TypeTraits/CreateGetTypeAliasOrDefault.hpp"
 
 /// \cond
 namespace tuples {
@@ -23,6 +24,12 @@ class GlobalCache;
 /// \endcond
 
 namespace Actions {
+namespace detail {
+CREATE_GET_TYPE_ALIAS_OR_DEFAULT(simple_tags)
+CREATE_GET_TYPE_ALIAS_OR_DEFAULT(compute_tags)
+CREATE_GET_TYPE_ALIAS_OR_DEFAULT(const_global_cache_tags)
+CREATE_GET_TYPE_ALIAS_OR_DEFAULT(mutable_global_cache_tags)
+}  // namespace detail
 /*!
  * \ingroup ActionsGroup
  * \brief Apply the function `Mutator::apply` to the DataBox
@@ -40,6 +47,15 @@ namespace Actions {
  */
 template <typename Mutator>
 struct MutateApply {
+  using simple_tags =
+      detail::get_simple_tags_or_default_t<Mutator, tmpl::list<>>;
+  using compute_tags =
+      detail::get_compute_tags_or_default_t<Mutator, tmpl::list<>>;
+  using const_global_cache_tags =
+      detail::get_const_global_cache_tags_or_default_t<Mutator, tmpl::list<>>;
+  using mutable_global_cache_tags =
+      detail::get_mutable_global_cache_tags_or_default_t<Mutator, tmpl::list<>>;
+
   template <typename DataBox, typename... InboxTags, typename Metavariables,
             typename ArrayIndex, typename ActionList,
             typename ParallelComponent>

support/Pipelines/Bbh/Ringdown.yaml

@@ -22,6 +22,7 @@ InitialData:
     Variables:
       SpacetimeMetric: SpacetimeMetric
       Pi: Pi
+      Phi: Phi
 
 DomainCreator:
   Sphere: