diff --git a/src/Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.cpp b/src/Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.cpp
index 417d78e7a4f6..cdfd343d3a18 100644
--- a/src/Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.cpp
+++ b/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,
diff --git a/src/Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp b/src/Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp
index 818d691dbde1..dc308174aa3f 100644
--- a/src/Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp
+++ b/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;
 
diff --git a/src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.cpp b/src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.cpp
index 9abc02af2632..725e4d49b66d 100644
--- a/src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.cpp
+++ b/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);
diff --git a/src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp b/src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp
index e9202009c880..8c1ff57ebbce 100644
--- a/src/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp
+++ b/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
diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/SetInitialData.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/SetInitialData.hpp
index 64160e06de77..32a65d6e1545 100644
--- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/SetInitialData.hpp
+++ b/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]() {
diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/SetPiAndPhiFromConstraints.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/SetPiAndPhiFromConstraints.hpp
index 185a34db5e33..3962a54f2b41 100644
--- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/SetPiAndPhiFromConstraints.hpp
+++ b/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);
     }
   }
 };
diff --git a/src/ParallelAlgorithms/Actions/MutateApply.hpp b/src/ParallelAlgorithms/Actions/MutateApply.hpp
index d25fa19b448e..cbdc8c147475 100644
--- a/src/ParallelAlgorithms/Actions/MutateApply.hpp
+++ b/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>
diff --git a/support/Pipelines/Bbh/Ringdown.yaml b/support/Pipelines/Bbh/Ringdown.yaml
index e28129b36c7c..2b062228ace8 100644
--- a/support/Pipelines/Bbh/Ringdown.yaml
+++ b/support/Pipelines/Bbh/Ringdown.yaml
@@ -22,6 +22,7 @@ InitialData:
     Variables:
       SpacetimeMetric: SpacetimeMetric
       Pi: Pi
+      Phi: Phi
 
 DomainCreator:
   Sphere:
diff --git a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp
index 87ae2b775bbe..76e12b2ea498 100644
--- a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp
+++ b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp
@@ -19,6 +19,7 @@
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Tags.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
+#include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/System.hpp"
 #include "Framework/ActionTesting.hpp"
 #include "Framework/TestCreation.hpp"
@@ -50,6 +51,8 @@ struct MockElementArray {
   using metavariables = Metavariables;
   using chare_type = ActionTesting::MockArrayChare;
   using array_index = ElementId<3>;
+  using mutable_global_cache_tags =
+      tmpl::list<gh::Tags::SetPiAndPhiFromConstraints>;
   using phase_dependent_action_list = tmpl::list<
       Parallel::PhaseActions<
           Parallel::Phase::Initialization,
@@ -87,6 +90,8 @@ struct MockReadVolumeData {
             "CustomSpacetimeMetric");
       CHECK(get<importers::Tags::Selected<Tags::Pi<DataVector, 3>>>(
                 selected_fields) == "CustomPi");
+      CHECK(get<importers::Tags::Selected<Tags::Phi<DataVector, 3>>>(
+                selected_fields) == "CustomPhi");
       CHECK_FALSE(
           get<importers::Tags::Selected<
               gr::Tags::SpatialMetric<DataVector, 3>>>(selected_fields));
@@ -112,6 +117,8 @@ struct MockReadVolumeData {
               gr::Tags::SpacetimeMetric<DataVector, 3>>>(selected_fields));
       CHECK_FALSE(get<importers::Tags::Selected<Tags::Pi<DataVector, 3>>>(
           selected_fields));
+      CHECK_FALSE(get<importers::Tags::Selected<Tags::Phi<DataVector, 3>>>(
+          selected_fields));
     } else {
       REQUIRE(false);
     }
@@ -172,14 +179,14 @@ void test_set_initial_data(
   using element_array = MockElementArray<Metavariables>;
 
   ActionTesting::MockRuntimeSystem<Metavariables> runner{
-      initial_data.get_clone()};
+      {initial_data.get_clone()}, {true}};
 
   // Setup mock data file reader
   ActionTesting::emplace_nodegroup_component<reader_component>(
       make_not_null(&runner));
 
   // Setup element
-  const ElementId<3> element_id{0, {{{1, 0}, {1, 0}, {1, 0}}}};
+  const ElementId<3> element_id{0};
   const Mesh<3> mesh{8, Spectral::Basis::Legendre,
                      Spectral::Quadrature::GaussLobatto};
   const auto map =
@@ -208,6 +215,8 @@ void test_set_initial_data(
 
   ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing);
 
+  const auto& cache = ActionTesting::cache<element_array>(runner, element_id);
+
   // SetInitialData
   ActionTesting::next_action<element_array>(make_not_null(&runner), element_id);
 
@@ -216,6 +225,10 @@ void test_set_initial_data(
     const auto& numeric_id =
         dynamic_cast<const NumericInitialData&>(initial_data);
 
+    CHECK(Parallel::get<gh::Tags::SetPiAndPhiFromConstraints>(cache) ==
+          std::holds_alternative<NumericInitialData::AdmVars>(
+              numeric_id.selected_variables()));
+
     REQUIRE_FALSE(ActionTesting::next_action_if_ready<element_array>(
         make_not_null(&runner), element_id));
 
@@ -235,6 +248,8 @@ void test_set_initial_data(
           get<gr::Tags::SpacetimeMetric<DataVector, 3>>(kerr_gh_vars);
       get<Tags::Pi<DataVector, 3>>(inbox) =
           get<Tags::Pi<DataVector, 3>>(kerr_gh_vars);
+      get<Tags::Phi<DataVector, 3>>(inbox) =
+          get<Tags::Phi<DataVector, 3>>(kerr_gh_vars);
     } else if (std::holds_alternative<NumericInitialData::AdmVars>(
                    selected_vars)) {
       const auto kerr_adm_vars = kerr.variables(
@@ -264,6 +279,8 @@ void test_set_initial_data(
     // ReceiveNumericInitialData
     ActionTesting::next_action<element_array>(make_not_null(&runner),
                                               element_id);
+  } else {
+    CHECK(Parallel::get<gh::Tags::SetPiAndPhiFromConstraints>(cache));
   }
 
   // Check result. These variables are not particularly precise because we are
@@ -286,13 +303,13 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.Gh.NumericInitialData",
                   "[Unit][Evolution]") {
   register_factory_classes_with_charm<Metavariables>();
   test_set_initial_data(
-      NumericInitialData{
-          "TestInitialData.h5",
-          "VolumeData",
-          0.,
-          {1.0e-9},
-          false,
-          NumericInitialData::GhVars{"CustomSpacetimeMetric", "CustomPi"}},
+      NumericInitialData{"TestInitialData.h5",
+                         "VolumeData",
+                         0.,
+                         {1.0e-9},
+                         false,
+                         NumericInitialData::GhVars{"CustomSpacetimeMetric",
+                                                    "CustomPi", "CustomPhi"}},
       "NumericInitialData:\n"
       "  FileGlob: TestInitialData.h5\n"
       "  Subgroup: VolumeData\n"
@@ -301,7 +318,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.Gh.NumericInitialData",
       "  ElementsAreIdentical: False\n"
       "  Variables:\n"
       "    SpacetimeMetric: CustomSpacetimeMetric\n"
-      "    Pi: CustomPi\n",
+      "    Pi: CustomPi\n"
+      "    Phi: CustomPhi\n",
       true);
   test_set_initial_data(
       NumericInitialData{
diff --git a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_SetPiAndPhiFromConstraints.cpp b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_SetPiAndPhiFromConstraints.cpp
index a4ab966cfb7e..90b1f44967e2 100644
--- a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_SetPiAndPhiFromConstraints.cpp
+++ b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_SetPiAndPhiFromConstraints.cpp
@@ -25,6 +25,7 @@
 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
 #include "Domain/Tags.hpp"
 #include "Evolution/Initialization/Tags.hpp"
+#include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/Constraints.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/DampedHarmonic.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Dispatch.hpp"
@@ -84,6 +85,44 @@ void test(const gsl::not_null<std::mt19937*> generator) {
         .get(i + 1, 0) *= 0.01;
   }
 
+  // Testing SetPiAndPhiFromConstraints = False. Note that we put the
+  // SetPiAndPhiFromConstraints tag in the box here because we don't have a
+  // cache.
+  {
+    auto box = db::create<db::AddSimpleTags<
+        ::Tags::Time, ::Tags::Variables<evolved_vars_tags>,
+        domain::Tags::Mesh<Dim>, domain::Tags::ElementMap<Dim, Frame::Grid>,
+        domain::CoordinateMaps::Tags::CoordinateMap<Dim, Frame::Grid,
+                                                    Frame::Inertial>,
+        domain::Tags::FunctionsOfTimeInitialize,
+        domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
+        gh::gauges::Tags::GaugeCondition,
+        gh::Tags::SetPiAndPhiFromConstraints>>(
+        0., evolved_vars, mesh,
+        ElementMap<Dim, Frame::Grid>{
+            ElementId<Dim>{0},
+            domain::make_coordinate_map_base<Frame::BlockLogical, Frame::Grid>(
+                domain::CoordinateMaps::Identity<Dim>{})},
+        domain::make_coordinate_map_base<Frame::Grid, Frame::Inertial>(
+            domain::CoordinateMaps::Identity<Dim>{}),
+        std::unordered_map<
+            std::string,
+            std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>{},
+        logical_coordinates(mesh),
+        std::unique_ptr<gh::gauges::GaugeCondition>(
+            std::make_unique<gh::gauges::DampedHarmonic>(
+                100., std::array{1.2, 1.5, 1.7}, std::array{2, 4, 6})),
+        false);
+    db::mutate_apply<gh::gauges::SetPiAndPhiFromConstraints<Dim>>(
+        make_not_null(&box));
+
+    // Should be exact since we didn't compute anything
+    CHECK(get<gh::Tags::Pi<DataVector, Dim>>(evolved_vars) ==
+          db::get<gh::Tags::Pi<DataVector, Dim>>(box));
+    CHECK(get<gh::Tags::Phi<DataVector, Dim>>(evolved_vars) ==
+          db::get<gh::Tags::Phi<DataVector, Dim>>(box));
+  }
+
   auto box = db::create<db::AddSimpleTags<
       ::Tags::Time, ::Tags::Variables<evolved_vars_tags>,
       domain::Tags::Mesh<Dim>, domain::Tags::ElementMap<Dim, Frame::Grid>,
@@ -91,7 +130,7 @@ void test(const gsl::not_null<std::mt19937*> generator) {
                                                   Frame::Inertial>,
       domain::Tags::FunctionsOfTimeInitialize,
       domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
-      gh::gauges::Tags::GaugeCondition>>(
+      gh::gauges::Tags::GaugeCondition, gh::Tags::SetPiAndPhiFromConstraints>>(
       0., evolved_vars, mesh,
       ElementMap<Dim, Frame::Grid>{
           ElementId<Dim>{0},
@@ -105,7 +144,8 @@ void test(const gsl::not_null<std::mt19937*> generator) {
       logical_coordinates(mesh),
       std::unique_ptr<gh::gauges::GaugeCondition>(
           std::make_unique<gh::gauges::DampedHarmonic>(
-              100., std::array{1.2, 1.5, 1.7}, std::array{2, 4, 6})));
+              100., std::array{1.2, 1.5, 1.7}, std::array{2, 4, 6})),
+      true);
   db::mutate_apply<gh::gauges::SetPiAndPhiFromConstraints<Dim>>(
       make_not_null(&box));
 
diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/Test_SetInitialData.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/Test_SetInitialData.cpp
index 084599ed3aec..6ec8ebbbe839 100644
--- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/Test_SetInitialData.cpp
+++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/Test_SetInitialData.cpp
@@ -60,6 +60,8 @@ struct MockElementArray {
   using metavariables = Metavariables;
   using chare_type = ActionTesting::MockArrayChare;
   using array_index = ElementId<3>;
+  using mutable_global_cache_tags =
+      tmpl::list<gh::Tags::SetPiAndPhiFromConstraints>;
   using phase_dependent_action_list = tmpl::list<
       Parallel::PhaseActions<
           Parallel::Phase::Initialization,
@@ -157,7 +159,7 @@ void test_set_initial_data(
   using element_array = MockElementArray<Metavariables>;
 
   ActionTesting::MockRuntimeSystem<Metavariables> runner{
-      initial_data.get_clone()};
+      {initial_data.get_clone()}, {true}};
 
   // We get test data from a TOV solution
   TovStar tov_star{
@@ -211,6 +213,8 @@ void test_set_initial_data(
 
   ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing);
 
+  const auto& cache = ActionTesting::cache<element_array>(runner, element_id);
+
   // SetInitialData
   ActionTesting::next_action<element_array>(make_not_null(&runner), element_id);
 
@@ -219,6 +223,10 @@ void test_set_initial_data(
     const auto& numeric_id =
         dynamic_cast<const NumericInitialData&>(initial_data);
 
+    CHECK(Parallel::get<gh::Tags::SetPiAndPhiFromConstraints>(cache) ==
+          std::holds_alternative<gh::NumericInitialData::AdmVars>(
+              numeric_id.gh_numeric_id().selected_variables()));
+
     REQUIRE_FALSE(ActionTesting::next_action_if_ready<element_array>(
         make_not_null(&runner), element_id));
 
@@ -240,6 +248,8 @@ void test_set_initial_data(
           get<gr::Tags::SpacetimeMetric<DataVector, 3>>(tov_vars);
       get<gh::Tags::Pi<DataVector, 3>>(inbox) =
           get<gh::Tags::Pi<DataVector, 3>>(tov_vars);
+      get<gh::Tags::Phi<DataVector, 3>>(inbox) =
+          get<gh::Tags::Phi<DataVector, 3>>(tov_vars);
     } else if (std::holds_alternative<gh::NumericInitialData::AdmVars>(
                    gh_selected_vars)) {
       const auto tov_adm_vars = tov_star.variables(
@@ -302,6 +312,8 @@ void test_set_initial_data(
     // ReceiveNumericInitialData
     ActionTesting::next_action<element_array>(make_not_null(&runner),
                                               element_id);
+  } else {
+    CHECK(Parallel::get<gh::Tags::SetPiAndPhiFromConstraints>(cache));
   }
 
   // Check result. The GH variables are not particularly precise because we
@@ -337,7 +349,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.GhValenciaDivClean.SetInitialData",
             0.,
             {1.0e-9},
             false,
-            gh::NumericInitialData::GhVars{"CustomSpacetimeMetric", "CustomPi"},
+            gh::NumericInitialData::GhVars{"CustomSpacetimeMetric", "CustomPi",
+                                           "CustomPhi"},
             {"CustomRho", "CustomUi", "CustomYe", "CustomB"},
             1.e-14},
         "NumericInitialData:\n"
@@ -349,6 +362,7 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.GhValenciaDivClean.SetInitialData",
         "  GhVariables:\n"
         "    SpacetimeMetric: CustomSpacetimeMetric\n"
         "    Pi: CustomPi\n"
+        "    Phi: CustomPhi\n"
         "  HydroVariables:\n"
         "    RestMassDensity: CustomRho\n"
         "    LowerSpatialFourVelocity: CustomUi\n"
diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Test_SetPiAndPhiFromConstraints.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Test_SetPiAndPhiFromConstraints.cpp
index 1654c6045c8b..da241813f577 100644
--- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Test_SetPiAndPhiFromConstraints.cpp
+++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Test_SetPiAndPhiFromConstraints.cpp
@@ -29,6 +29,7 @@
 #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
 #include "Evolution/Initialization/Tags.hpp"
+#include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/DampedHarmonic.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Dispatch.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Gauges.hpp"
@@ -48,8 +49,9 @@
 #include "Utilities/TMPL.hpp"
 
 SPECTRE_TEST_CASE(
- "Unit.Evolution.Systems.GrMhd.GhValenciaDivClean.SetPiAndPhiFromConstraints",
- "[Unit][Evolution][Actions]") {
+    "Unit.Evolution.Systems.GrMhd.GhValenciaDivClean."
+    "SetPiAndPhiFromConstraints",
+    "[Unit][Evolution][Actions]") {
   MAKE_GENERATOR(generator);
   std::uniform_real_distribution<> metric_dist(0.1, 1.);
   std::uniform_real_distribution<> deriv_dist(-1.e-5, 1.e-5);
@@ -86,6 +88,50 @@ SPECTRE_TEST_CASE(
 
   const auto initial_dg_vars = make_vars(dg_mesh);
 
+  // Testing SetPiAndPhiFromConstraints = False. Note that we put the
+  // SetPiAndPhiFromConstraints tag in the box here because we don't have a
+  // cache.
+  {
+    auto box = db::create<
+        db::AddSimpleTags<::Tags::Time, ::Tags::Variables<evolved_vars_tags>,
+                          domain::Tags::Mesh<3>,
+                          domain::Tags::ElementMap<3, Frame::Grid>,
+                          domain::CoordinateMaps::Tags::CoordinateMap<
+                              3, Frame::Grid, Frame::Inertial>,
+                          domain::Tags::FunctionsOfTimeInitialize,
+                          domain::Tags::Coordinates<3, Frame::ElementLogical>,
+                          gh::gauges::Tags::GaugeCondition,
+                          gh::Tags::SetPiAndPhiFromConstraints,
+                          evolution::dg::subcell::Tags::ActiveGrid>,
+        db::AddComputeTags<
+            evolution::dg::subcell::Tags::MeshCompute<3>,
+            evolution::dg::subcell::Tags::LogicalCoordinatesCompute<3>>>(
+        0., initial_dg_vars, dg_mesh,
+        ElementMap<3, Frame::Grid>{
+            ElementId<3>{0},
+            domain::make_coordinate_map_base<Frame::BlockLogical, Frame::Grid>(
+                domain::CoordinateMaps::Identity<3>{})},
+        domain::make_coordinate_map_base<Frame::Grid, Frame::Inertial>(
+            domain::CoordinateMaps::Identity<3>{}),
+        std::unordered_map<
+            std::string,
+            std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>{},
+        logical_coordinates(dg_mesh),
+        std::unique_ptr<gh::gauges::GaugeCondition>(
+            std::make_unique<gh::gauges::DampedHarmonic>(
+                100., std::array{1.2, 1.5, 1.7}, std::array{2, 4, 6})),
+        false, evolution::dg::subcell::ActiveGrid::Dg);
+
+    db::mutate_apply<grmhd::GhValenciaDivClean::SetPiAndPhiFromConstraints>(
+        make_not_null(&box));
+
+    // Should be exact since we didn't compute anything
+    CHECK(get<gh::Tags::Pi<DataVector, 3>>(initial_dg_vars) ==
+          db::get<gh::Tags::Pi<DataVector, 3>>(box));
+    CHECK(get<gh::Tags::Phi<DataVector, 3>>(initial_dg_vars) ==
+          db::get<gh::Tags::Phi<DataVector, 3>>(box));
+  }
+
   auto box = db::create<
       db::AddSimpleTags<::Tags::Time, ::Tags::Variables<evolved_vars_tags>,
                         domain::Tags::Mesh<3>,
@@ -95,6 +141,7 @@ SPECTRE_TEST_CASE(
                         domain::Tags::FunctionsOfTimeInitialize,
                         domain::Tags::Coordinates<3, Frame::ElementLogical>,
                         gh::gauges::Tags::GaugeCondition,
+                        gh::Tags::SetPiAndPhiFromConstraints,
                         evolution::dg::subcell::Tags::ActiveGrid>,
       db::AddComputeTags<
           evolution::dg::subcell::Tags::MeshCompute<3>,
@@ -113,7 +160,7 @@ SPECTRE_TEST_CASE(
       std::unique_ptr<gh::gauges::GaugeCondition>(
           std::make_unique<gh::gauges::DampedHarmonic>(
               100., std::array{1.2, 1.5, 1.7}, std::array{2, 4, 6})),
-      evolution::dg::subcell::ActiveGrid::Dg);
+      true, evolution::dg::subcell::ActiveGrid::Dg);
   const auto initial_subcell_vars =
       make_vars(db::get<evolution::dg::subcell::Tags::Mesh<3>>(box));
 
@@ -134,7 +181,7 @@ SPECTRE_TEST_CASE(
             box),
         db::get<domain::Tags::FunctionsOfTime>(box), logical_coordinates(mesh),
         get<gr::Tags::SpacetimeMetric<DataVector, 3>>(initial_vars),
-        db::get<gh::gauges::Tags::GaugeCondition>(box));
+        db::get<gh::gauges::Tags::GaugeCondition>(box), true);
 
     const auto& pi = db::get<gh::Tags::Pi<DataVector, 3>>(box);
     CHECK(pi == expected_pi);