Convert M1Grey to use runtime initial data

src/Evolution/Executables/RadiationTransport/M1Grey/CMakeLists.txt

@@ -1,7 +1,7 @@
 # Distributed under the MIT License.
 # See LICENSE.txt for details.
 
-set(EXECUTABLE EvolveM1GreyConstantM1)
+set(EXECUTABLE EvolveM1Grey)
 
 add_spectre_executable(
   ${EXECUTABLE}

src/Evolution/Executables/RadiationTransport/M1Grey/EvolveM1Grey.hpp

@@ -70,9 +70,10 @@
 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
 #include "PointwiseFunctions/AnalyticData/Tags.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
-#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/ConstantM1.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/Factory.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
 #include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
 #include "Time/Actions/AdvanceTime.hpp"
 #include "Time/Actions/CleanHistory.hpp"
 #include "Time/Actions/RecordTimeStepperData.hpp"
@@ -109,13 +110,8 @@ class CProxy_GlobalCache;
 struct EvolutionMetavars {
   static constexpr size_t volume_dim = 3;
 
-  // To switch which initial data is evolved you only need to change the
-  // line `using initial_data = ...;` and include the header file for the
-  // solution.
-  using initial_data = RadiationTransport::M1Grey::Solutions::ConstantM1;
-  static_assert(
-      is_analytic_data_v<initial_data> xor is_analytic_solution_v<initial_data>,
-      "initial_data must be either an analytic_data or an analytic_solution");
+  using initial_data_list =
+      RadiationTransport::M1Grey::Solutions::all_solutions;
 
   // Set list of neutrino species to be used by M1 code
   using neutrino_species = tmpl::list<neutrinos::ElectronNeutrinos<1>>;
@@ -128,17 +124,12 @@ struct EvolutionMetavars {
       TimeStepperBase::local_time_stepping;
   static constexpr bool use_dg_element_collection = false;
 
-  using initial_data_tag =
-      tmpl::conditional_t<is_analytic_solution_v<initial_data>,
-                          Tags::AnalyticSolution<initial_data>,
-                          Tags::AnalyticData<initial_data>>;
   using analytic_variables_tags = typename system::variables_tag::tags_list;
   using limiter = Tags::Limiter<
       Limiters::Minmod<3, typename system::variables_tag::tags_list>>;
 
-  using analytic_compute =
-      evolution::Tags::AnalyticSolutionsCompute<volume_dim,
-                                                analytic_variables_tags, false>;
+  using analytic_compute = evolution::Tags::AnalyticSolutionsCompute<
+      volume_dim, analytic_variables_tags, false, initial_data_list>;
   using error_compute = Tags::ErrorsCompute<analytic_variables_tags>;
   using error_tags = db::wrap_tags_in<Tags::Error, analytic_variables_tags>;
   using observe_fields = tmpl::push_back<
@@ -165,7 +156,7 @@ struct EvolutionMetavars {
                            volume_dim, observe_fields, non_tensor_compute_tags>,
                        Events::time_events<system>>>>,
         tmpl::pair<ImexTimeStepper, TimeSteppers::imex_time_steppers>,
-
+        tmpl::pair<evolution::initial_data::InitialData, initial_data_list>,
         tmpl::pair<PhaseChange, PhaseControl::factory_creatable_classes>,
         tmpl::pair<RadiationTransport::M1Grey::BoundaryConditions::
                        BoundaryCondition<neutrino_species>,
@@ -192,7 +183,7 @@ struct EvolutionMetavars {
   static_assert(not local_time_stepping);
   using step_actions = tmpl::flatten<tmpl::list<
       Actions::MutateApply<
-          evolution::dg::BackgroundGrVars<system, EvolutionMetavars, false>>,
+          evolution::dg::BackgroundGrVars<system, EvolutionMetavars, true>>,
       evolution::dg::Actions::ComputeTimeDerivative<
           volume_dim, system, AllStepChoosers, local_time_stepping,
           use_dg_element_collection>,
@@ -227,7 +218,7 @@ struct EvolutionMetavars {
           evolution::dg::Initialization::Domain<volume_dim>,
           Initialization::TimeStepperHistory<EvolutionMetavars>>,
       Initialization::Actions::AddSimpleTags<
-          evolution::dg::BackgroundGrVars<system, EvolutionMetavars, false>>,
+          evolution::dg::BackgroundGrVars<system, EvolutionMetavars, true>>,
       Initialization::Actions::ConservativeSystem<system>,
       evolution::Initialization::Actions::SetVariables<
           domain::Tags::Coordinates<volume_dim, Frame::ElementLogical>>,
@@ -275,7 +266,8 @@ struct EvolutionMetavars {
                  observers::ObserverWriter<EvolutionMetavars>,
                  dg_element_array>;
 
-  using const_global_cache_tags = tmpl::list<initial_data_tag>;
+  using const_global_cache_tags =
+      tmpl::list<evolution::initial_data::Tags::InitialData>;
 
   static constexpr Options::String help{
       "Evolve the M1Grey system (without coupling to hydro).\n\n"};

src/Evolution/Systems/RadiationTransport/M1Grey/BoundaryConditions/DirichletAnalytic.hpp

@@ -22,8 +22,11 @@
 #include "Options/String.hpp"
 #include "PointwiseFunctions/AnalyticData/Tags.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/Factory.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 #include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
+#include "Utilities/CallWithDynamicType.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
 #include "Utilities/TMPL.hpp"
@@ -53,18 +56,42 @@ template <typename... NeutrinoSpecies>
 class DirichletAnalytic<tmpl::list<NeutrinoSpecies...>> final
     : public BoundaryCondition<tmpl::list<NeutrinoSpecies...>> {
  public:
-  using options = tmpl::list<>;
+  /// \brief What analytic solution/data to prescribe.
+  struct AnalyticPrescription {
+    static constexpr Options::String help =
+        "What analytic solution/data to prescribe.";
+    using type = std::unique_ptr<evolution::initial_data::InitialData>;
+  };
+
+  using options = tmpl::list<AnalyticPrescription>;
   static constexpr Options::String help{
       "DirichletAnalytic boundary conditions using either analytic solution or "
       "analytic data."};
 
   DirichletAnalytic() = default;
   DirichletAnalytic(DirichletAnalytic&&) = default;
   DirichletAnalytic& operator=(DirichletAnalytic&&) = default;
-  DirichletAnalytic(const DirichletAnalytic&) = default;
-  DirichletAnalytic& operator=(const DirichletAnalytic&) = default;
+
+  DirichletAnalytic(const DirichletAnalytic& rhs)
+      : BoundaryCondition<tmpl::list<NeutrinoSpecies...>>{dynamic_cast<
+            const BoundaryCondition<tmpl::list<NeutrinoSpecies...>>&>(rhs)},
+        analytic_prescription_(rhs.analytic_prescription_->get_clone()) {}
+
+  DirichletAnalytic& operator=(const DirichletAnalytic& rhs) {
+    if (&rhs == this) {
+      return *this;
+    }
+    analytic_prescription_ = rhs.analytic_prescription_->get_clone();
+    return *this;
+  }
+
   ~DirichletAnalytic() override = default;
 
+  explicit DirichletAnalytic(
+      std::unique_ptr<evolution::initial_data::InitialData>
+          analytic_prescription)
+      : analytic_prescription_(std::move(analytic_prescription)) {}
+
   explicit DirichletAnalytic(CkMigrateMessage* msg)
       : BoundaryCondition<tmpl::list<NeutrinoSpecies...>>(msg) {}
 
@@ -81,16 +108,15 @@ class DirichletAnalytic<tmpl::list<NeutrinoSpecies...>> final
 
   void pup(PUP::er& p) override {
     BoundaryCondition<tmpl::list<NeutrinoSpecies...>>::pup(p);
+    p | analytic_prescription_;
   }
 
   using dg_interior_evolved_variables_tags = tmpl::list<>;
   using dg_interior_temporary_tags =
       tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>>;
   using dg_interior_primitive_variables_tags = tmpl::list<>;
-  using dg_gridless_tags =
-      tmpl::list<::Tags::Time, ::Tags::AnalyticSolutionOrData>;
+  using dg_gridless_tags = tmpl::list<::Tags::Time>;
 
-  template <typename AnalyticSolutionOrData>
   std::optional<std::string> dg_ghost(
       const gsl::not_null<typename Tags::TildeE<
           Frame::Inertial, NeutrinoSpecies>::type*>... tilde_e,
@@ -111,35 +137,47 @@ class DirichletAnalytic<tmpl::list<NeutrinoSpecies...>> final
           tnsr::I<DataVector, 3, Frame::Inertial>>& /*face_mesh_velocity*/,
       const tnsr::i<DataVector, 3, Frame::Inertial>& /*normal_covector*/,
       const tnsr::I<DataVector, 3, Frame::Inertial>& /*normal_vector*/,
-      const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
-      const AnalyticSolutionOrData& analytic_solution_or_data) const {
-    auto boundary_values = [&analytic_solution_or_data, &coords, &time]() {
-      if constexpr (is_analytic_solution_v<AnalyticSolutionOrData>) {
-        return analytic_solution_or_data.variables(
-            coords, time,
-            tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
-                       Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
-                       hydro::Tags::LorentzFactor<DataVector>,
-                       hydro::Tags::SpatialVelocity<DataVector, 3>,
-                       gr::Tags::Lapse<DataVector>,
-                       gr::Tags::Shift<DataVector, 3>,
-                       gr::Tags::SpatialMetric<DataVector, 3>,
-                       gr::Tags::InverseSpatialMetric<DataVector, 3>>{});
-
-      } else {
-        (void)time;
-        return analytic_solution_or_data.variables(
-            coords,
-            tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
-                       Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
-                       hydro::Tags::LorentzFactor<DataVector>,
-                       hydro::Tags::SpatialVelocity<DataVector, 3>,
-                       gr::Tags::Lapse<DataVector>,
-                       gr::Tags::Shift<DataVector, 3>,
-                       gr::Tags::SpatialMetric<DataVector, 3>,
-                       gr::Tags::InverseSpatialMetric<DataVector, 3>>{});
-      }
-    }();
+      const tnsr::I<DataVector, 3, Frame::Inertial>& coords,
+      const double time) const {
+    auto boundary_values = call_with_dynamic_type<
+        tuples::TaggedTuple<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
+                            Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
+                            hydro::Tags::LorentzFactor<DataVector>,
+                            hydro::Tags::SpatialVelocity<DataVector, 3>,
+                            gr::Tags::Lapse<DataVector>,
+                            gr::Tags::Shift<DataVector, 3>,
+                            gr::Tags::SpatialMetric<DataVector, 3>,
+                            gr::Tags::InverseSpatialMetric<DataVector, 3>>,
+        RadiationTransport::M1Grey::Solutions::all_solutions>(
+        analytic_prescription_.get(),
+        [&coords, &time](const auto* const analytic_solution_or_data) {
+          if constexpr (is_analytic_solution_v<std::decay_t<
+                            decltype(*analytic_solution_or_data)>>) {
+            return analytic_solution_or_data->variables(
+                coords, time,
+                tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
+                           Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
+                           hydro::Tags::LorentzFactor<DataVector>,
+                           hydro::Tags::SpatialVelocity<DataVector, 3>,
+                           gr::Tags::Lapse<DataVector>,
+                           gr::Tags::Shift<DataVector, 3>,
+                           gr::Tags::SpatialMetric<DataVector, 3>,
+                           gr::Tags::InverseSpatialMetric<DataVector, 3>>{});
+
+          } else {
+            (void)time;
+            return analytic_solution_or_data->variables(
+                coords,
+                tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
+                           Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
+                           hydro::Tags::LorentzFactor<DataVector>,
+                           hydro::Tags::SpatialVelocity<DataVector, 3>,
+                           gr::Tags::Lapse<DataVector>,
+                           gr::Tags::Shift<DataVector, 3>,
+                           gr::Tags::SpatialMetric<DataVector, 3>,
+                           gr::Tags::InverseSpatialMetric<DataVector, 3>>{});
+          }
+        });
 
     *inv_spatial_metric =
         get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(boundary_values);
@@ -214,6 +252,9 @@ class DirichletAnalytic<tmpl::list<NeutrinoSpecies...>> final
         tilde_e, tilde_s, flux_tilde_e, flux_tilde_s, NeutrinoSpecies{})...);
     return {};
   }
+
+ private:
+  std::unique_ptr<evolution::initial_data::InitialData> analytic_prescription_;
 };
 
 /// \cond

src/Evolution/Systems/RadiationTransport/M1Grey/Initialize.hpp

@@ -13,6 +13,7 @@
 #include "Parallel/AlgorithmExecution.hpp"
 #include "Parallel/GlobalCache.hpp"
 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
@@ -21,8 +22,10 @@
 namespace Frame {
 struct Inertial;
 }  // namespace Frame
+namespace evolution::initial_data::Tags {
+struct InitialData;
+}  // namespace evolution::initial_data::Tags
 namespace Tags {
-struct AnalyticSolutionOrData;
 struct Time;
 }  // namespace Tags
 namespace domain {
@@ -57,41 +60,46 @@ struct InitializeM1Tags {
   static Parallel::iterable_action_return_t apply(
       db::DataBox<DbTagsList>& box,
       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
-      const Parallel::GlobalCache<Metavariables>& cache,
+      const Parallel::GlobalCache<Metavariables>& /*cache*/,
       const ArrayIndex& /*array_index*/, ActionList /*meta*/,
       const ParallelComponent* const /*meta*/) {
     using EvolvedVars = typename evolved_variables_tag::type;
     using HydroVars = typename hydro_variables_tag::type;
     using M1Vars = typename m1_variables_tag::type;
+    using derived_classes =
+        tmpl::at<typename Metavariables::factory_creation::factory_classes,
+                 evolution::initial_data::InitialData>;
+    call_with_dynamic_type<void, derived_classes>(
+        &db::get<evolution::initial_data::Tags::InitialData>(box),
+        [&box](const auto* const data_or_solution) {
+          static constexpr size_t dim = System::volume_dim;
+          const double initial_time = db::get<::Tags::Time>(box);
+          const size_t num_grid_points =
+              db::get<domain::Tags::Mesh<dim>>(box).number_of_grid_points();
+          const auto& inertial_coords =
+              db::get<domain::Tags::Coordinates<dim, Frame::Inertial>>(box);
 
-    static constexpr size_t dim = System::volume_dim;
-    const double initial_time = db::get<::Tags::Time>(box);
-    const size_t num_grid_points =
-        db::get<domain::Tags::Mesh<dim>>(box).number_of_grid_points();
-    const auto& inertial_coords =
-        db::get<domain::Tags::Coordinates<dim, Frame::Inertial>>(box);
+          db::mutate<evolved_variables_tag>(
+              [initial_time, &inertial_coords, &data_or_solution](
+                  const gsl::not_null<EvolvedVars*> evolved_vars) {
+                evolved_vars->assign_subset(
+                    evolution::Initialization::initial_data(
+                        *data_or_solution, inertial_coords, initial_time,
+                        typename evolved_variables_tag::tags_list{}));
+              },
+              make_not_null(&box));
 
-    db::mutate<evolved_variables_tag>(
-        [&cache, initial_time,
-         &inertial_coords](const gsl::not_null<EvolvedVars*> evolved_vars) {
-          evolved_vars->assign_subset(evolution::Initialization::initial_data(
-              Parallel::get<::Tags::AnalyticSolutionOrData>(cache),
-              inertial_coords, initial_time,
-              typename evolved_variables_tag::tags_list{}));
-        },
-        make_not_null(&box));
-
-    // Get hydro variables
-    HydroVars hydro_variables{num_grid_points};
-    hydro_variables.assign_subset(evolution::Initialization::initial_data(
-        Parallel::get<::Tags::AnalyticSolutionOrData>(cache), inertial_coords,
-        initial_time, typename hydro_variables_tag::tags_list{}));
-
-    M1Vars m1_variables{num_grid_points, -1.};
-    Initialization::mutate_assign<simple_tags>(make_not_null(&box),
-                                               std::move(hydro_variables),
-                                               std::move(m1_variables));
+          // Get hydro variables
+          HydroVars hydro_variables{num_grid_points};
+          hydro_variables.assign_subset(evolution::Initialization::initial_data(
+              *data_or_solution, inertial_coords, initial_time,
+              typename hydro_variables_tag::tags_list{}));
 
+          M1Vars m1_variables{num_grid_points, -1.};
+          Initialization::mutate_assign<simple_tags>(make_not_null(&box),
+                                                     std::move(hydro_variables),
+                                                     std::move(m1_variables));
+        });
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };

src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/CMakeLists.txt

@@ -16,6 +16,7 @@ spectre_target_headers(
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   ConstantM1.hpp
+  Factory.hpp
   Solutions.hpp
   )
 

src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/ConstantM1.cpp

@@ -24,12 +24,22 @@ ConstantM1::ConstantM1(const std::array<double, 3>& mean_velocity,
       mean_velocity_(std::move(mean_velocity)),  // NOLINT
       comoving_energy_density_(comoving_energy_density) {}
 
+std::unique_ptr<evolution::initial_data::InitialData> ConstantM1::get_clone()
+    const {
+  return std::make_unique<ConstantM1>(*this);
+}
+
+ConstantM1::ConstantM1(CkMigrateMessage* msg) : InitialData(msg) {}
+
 void ConstantM1::pup(PUP::er& p) {
+  InitialData::pup(p);
   p | mean_velocity_;
   p | comoving_energy_density_;
   p | background_spacetime_;
 }
 
+PUP::able::PUP_ID ConstantM1::my_PUP_ID = 0;  // NOLINT
+
 // Variables templated on neutrino species.
 template <typename NeutrinoSpecies>
 tuples::TaggedTuple<