SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/GeneralizedHarmonic/Actions - SetInitialData.hpp Hit Total Coverage
Commit: 664546099c4dbf27a1b708fac45e39c82dd743d2 Lines: 9 50 18.0 %
Date: 2024-04-19 16:28:01
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <pup.h>
       8             : #include <string>
       9             : #include <variant>
      10             : 
      11             : #include "DataStructures/DataBox/DataBox.hpp"
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "Domain/Structure/ElementId.hpp"
      15             : #include "Domain/Tags.hpp"
      16             : #include "Evolution/Initialization/InitialData.hpp"
      17             : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
      18             : #include "IO/Importers/Actions/ReadVolumeData.hpp"
      19             : #include "IO/Importers/ElementDataReader.hpp"
      20             : #include "IO/Importers/Tags.hpp"
      21             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      22             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      23             : #include "Parallel/AlgorithmExecution.hpp"
      24             : #include "Parallel/GlobalCache.hpp"
      25             : #include "Parallel/Invoke.hpp"
      26             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      27             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      28             : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
      29             : #include "Utilities/CallWithDynamicType.hpp"
      30             : #include "Utilities/ErrorHandling/Error.hpp"
      31             : #include "Utilities/Gsl.hpp"
      32             : #include "Utilities/Serialization/CharmPupable.hpp"
      33             : #include "Utilities/TMPL.hpp"
      34             : #include "Utilities/TaggedTuple.hpp"
      35             : 
      36             : /// \cond
      37             : namespace Tags {
      38             : struct Time;
      39             : }  // namespace Tags
      40             : /// \endcond
      41             : 
      42           1 : namespace gh {
      43             : 
      44             : /*!
      45             :  * \brief Compute initial GH variables from ADM variables
      46             :  *
      47             :  * - The spacetime metric is assembled from the spatial metric, lapse, and
      48             :  *   shift. See `gr::spacetime_metric` for details.
      49             :  * - Phi is set to the numerical derivative of the spacetime metric. This
      50             :  *   ensures that the 3-index constraint is initially satisfied.
      51             :  * - Pi is computed by choosing the time derivatives of lapse and shift to be
      52             :  *   zero. The `gh::gauges::SetPiAndPhiFromConstraints` mutator exists to
      53             :  *   override Pi later in the algorithm (it should be combined with this
      54             :  *   function).
      55             :  */
      56             : template <size_t Dim>
      57           1 : void initial_gh_variables_from_adm(
      58             :     gsl::not_null<tnsr::aa<DataVector, Dim>*> spacetime_metric,
      59             :     gsl::not_null<tnsr::aa<DataVector, Dim>*> pi,
      60             :     gsl::not_null<tnsr::iaa<DataVector, Dim>*> phi,
      61             :     const tnsr::ii<DataVector, Dim>& spatial_metric,
      62             :     const Scalar<DataVector>& lapse, const tnsr::I<DataVector, Dim>& shift,
      63             :     const tnsr::ii<DataVector, Dim>& extrinsic_curvature, const Mesh<Dim>& mesh,
      64             :     const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
      65             :                           Frame::Inertial>& inv_jacobian);
      66             : 
      67             : /*!
      68             :  * \brief Numeric initial data loaded from volume data files
      69             :  *
      70             :  * This class can be factory-created in the input file to start an evolution
      71             :  * from numeric initial data. It selects the set of variables to load from
      72             :  * the volume data file (ADM or GH variables).
      73             :  */
      74           1 : class NumericInitialData : public evolution::initial_data::InitialData {
      75             :  public:
      76             :   /// Name of a variable in the volume data file
      77             :   template <typename Tag>
      78           1 :   struct VarName {
      79           0 :     using tag = Tag;
      80           0 :     static std::string name() { return db::tag_name<Tag>(); }
      81           0 :     using type = std::string;
      82           0 :     static constexpr Options::String help =
      83             :         "Name of the variable in the volume data file";
      84             :   };
      85             : 
      86             :   // These are the sets of variables that we support loading from volume data
      87             :   // files:
      88             :   // - ADM variables
      89           0 :   using adm_vars =
      90             :       tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
      91             :                  gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
      92             :                  gr::Tags::ExtrinsicCurvature<DataVector, 3>>;
      93           0 :   struct AdmVars : tuples::tagged_tuple_from_typelist<
      94             :                        db::wrap_tags_in<VarName, adm_vars>> {
      95           0 :     static constexpr Options::String help =
      96             :         "ADM variables: 'Lapse', 'Shift', 'SpatialMetric' and "
      97             :         "'ExtrinsicCurvature'. The initial GH variables will be computed "
      98             :         "from these numeric fields, as well as their numeric spatial "
      99             :         "derivatives on the computational grid.";
     100           0 :     using options = tags_list;
     101             :     using TaggedTuple::TaggedTuple;
     102             :   };
     103             : 
     104             :   // - Generalized harmonic variables
     105           0 :   using gh_vars = tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>,
     106             :                              Tags::Pi<DataVector, 3>>;
     107           0 :   struct GhVars
     108             :       : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<VarName, gh_vars>> {
     109           0 :     static constexpr Options::String help =
     110             :         "GH variables: 'SpacetimeMetric' and 'Pi'. These variables are "
     111             :         "used to set the initial data directly; Phi is then set to the "
     112             :         "numerical derivative of SpacetimeMetric, to enforce the 3-index "
     113             :         "constraint.";
     114           0 :     using options = tags_list;
     115             :     using TaggedTuple::TaggedTuple;
     116             :   };
     117             : 
     118             :   // Collect all variables that we support loading from volume data files.
     119             :   // Remember to `tmpl::remove_duplicates` when adding overlapping sets of
     120             :   // vars.
     121           0 :   using all_vars = tmpl::append<adm_vars, gh_vars>;
     122             : 
     123             :   // Input-file options
     124           0 :   struct Variables {
     125             :     // The user can supply any of these choices of variables in the input
     126             :     // file
     127           0 :     using type = std::variant<AdmVars, GhVars>;
     128           0 :     static constexpr Options::String help =
     129             :         "Set of initial data variables from which the generalized harmonic "
     130             :         "system variables are computed.";
     131             :   };
     132             : 
     133           0 :   using options =
     134             :       tmpl::list<importers::OptionTags::FileGlob,
     135             :                  importers::OptionTags::Subgroup,
     136             :                  importers::OptionTags::ObservationValue,
     137             :                  importers::OptionTags::EnableInterpolation, Variables>;
     138             : 
     139           0 :   static constexpr Options::String help =
     140             :       "Numeric initial data loaded from volume data files";
     141             : 
     142           0 :   NumericInitialData() = default;
     143           0 :   NumericInitialData(const NumericInitialData& rhs) = default;
     144           0 :   NumericInitialData& operator=(const NumericInitialData& rhs) = default;
     145           0 :   NumericInitialData(NumericInitialData&& /*rhs*/) = default;
     146           0 :   NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
     147           0 :   ~NumericInitialData() = default;
     148             : 
     149             :   /// \cond
     150             :   explicit NumericInitialData(CkMigrateMessage* msg);
     151             :   using PUP::able::register_constructor;
     152             :   WRAPPED_PUPable_decl_template(NumericInitialData);
     153             :   /// \endcond
     154             : 
     155           0 :   std::unique_ptr<evolution::initial_data::InitialData> get_clone()
     156             :       const override {
     157             :     return std::make_unique<NumericInitialData>(*this);
     158             :   }
     159             : 
     160           0 :   NumericInitialData(
     161             :       std::string file_glob, std::string subfile_name,
     162             :       std::variant<double, importers::ObservationSelector> observation_value,
     163             :       bool enable_interpolation,
     164             :       std::variant<AdmVars, GhVars> selected_variables);
     165             : 
     166           0 :   const importers::ImporterOptions& importer_options() const {
     167             :     return importer_options_;
     168             :   }
     169             : 
     170           0 :   const std::variant<AdmVars, GhVars>& selected_variables() const {
     171             :     return selected_variables_;
     172             :   }
     173             : 
     174             :   /*!
     175             :    * \brief Unique identifier for loading this volume data
     176             :    *
     177             :    * Involves a hash of the type name and the volume data file names.
     178             :    */
     179           1 :   size_t volume_data_id() const;
     180             : 
     181             :   /*!
     182             :    * \brief Selects which of the `fields` to import based on the choices in the
     183             :    * input-file options
     184             :    *
     185             :    * The `fields` are all datasets that are available to import, represented by
     186             :    * `importers::Tags::Selected<Tag>` tags. We select only those that we need by
     187             :    * setting their dataset name.
     188             :    */
     189             :   template <typename... AllTags>
     190           1 :   void select_for_import(
     191             :       const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
     192             :     // Select the subset of the available variables that we want to read from
     193             :     // the volume data file
     194             :     std::visit(
     195             :         [&fields](const auto& vars) {
     196             :           // This lambda is invoked with the set of vars selected in the input
     197             :           // file, which map to the tensor names that should be read from the H5
     198             :           // file
     199             :           using selected_vars = std::decay_t<decltype(vars)>;
     200             :           // Get the mapped tensor name from the input file and select it in the
     201             :           // set of all possible vars.
     202             :           tmpl::for_each<typename selected_vars::tags_list>(
     203             :               [&fields, &vars](const auto tag_v) {
     204             :                 using tag = typename std::decay_t<decltype(tag_v)>::type::tag;
     205             :                 get<importers::Tags::Selected<tag>>(*fields) =
     206             :                     get<VarName<tag>>(vars);
     207             :               });
     208             :         },
     209             :         selected_variables_);
     210             :   }
     211             : 
     212             :   /*!
     213             :    * \brief Set GH initial data given numeric data loaded from files
     214             :    *
     215             :    * The `numeric_data` contains the datasets selected above (and possibly
     216             :    * more). We either set the GH variables directly, or compute them from the
     217             :    * ADM variables.
     218             :    */
     219             :   template <typename... AllTags>
     220           1 :   void set_initial_data(
     221             :       const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
     222             :       const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
     223             :       const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
     224             :       const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
     225             :       const Mesh<3>& mesh,
     226             :       const InverseJacobian<DataVector, 3, Frame::ElementLogical,
     227             :                             Frame::Inertial>& inv_jacobian) const {
     228             :     if (std::holds_alternative<NumericInitialData::GhVars>(
     229             :             selected_variables_)) {
     230             :       // We have loaded the GH system variables from the file, so just move the
     231             :       // data for spacetime_metric and Pi into the DataBox directly, with no
     232             :       // conversion needed. Set Phi to the spatial derivative of the spacetime
     233             :       // metric to enforce the 3-index constraint.
     234             :       *spacetime_metric = std::move(
     235             :           get<gr::Tags::SpacetimeMetric<DataVector, 3>>(*numeric_data));
     236             :       *pi = std::move(get<Tags::Pi<DataVector, 3>>(*numeric_data));
     237             :       // Set Phi to the numerical spatial derivative of spacetime_metric
     238             :       partial_derivative(phi, *spacetime_metric, mesh, inv_jacobian);
     239             :     } else if (std::holds_alternative<NumericInitialData::AdmVars>(
     240             :                    selected_variables_)) {
     241             :       // We have loaded ADM variables from the file. Convert to GH variables.
     242             :       const auto& spatial_metric =
     243             :           get<gr::Tags::SpatialMetric<DataVector, 3>>(*numeric_data);
     244             :       const auto& lapse = get<gr::Tags::Lapse<DataVector>>(*numeric_data);
     245             :       const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(*numeric_data);
     246             :       const auto& extrinsic_curvature =
     247             :           get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(*numeric_data);
     248             : 
     249             :       initial_gh_variables_from_adm(spacetime_metric, pi, phi, spatial_metric,
     250             :                                     lapse, shift, extrinsic_curvature, mesh,
     251             :                                     inv_jacobian);
     252             :     } else {
     253             :       ERROR(
     254             :           "These initial data variables are not implemented yet. Please add "
     255             :           "an implementation to gh::NumericInitialData.");
     256             :     }
     257             :   }
     258             : 
     259           0 :   void pup(PUP::er& p) override;
     260             : 
     261           0 :   friend bool operator==(const NumericInitialData& lhs,
     262             :                          const NumericInitialData& rhs);
     263             : 
     264             :  private:
     265           0 :   importers::ImporterOptions importer_options_;
     266           0 :   std::variant<AdmVars, GhVars> selected_variables_{};
     267             : };
     268             : 
     269           0 : namespace Actions {
     270             : 
     271             : /*!
     272             :  * \brief Dispatch loading numeric initial data from files or set analytic
     273             :  * initial data.
     274             :  *
     275             :  * Place this action before
     276             :  * gh::Actions::ReceiveNumericInitialData in the action list.
     277             :  * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
     278             :  * invoked by this action.
     279             :  * Analytic initial data is set directly by this action and terminates the
     280             :  * phase.
     281             :  */
     282           1 : struct SetInitialData {
     283           0 :   using const_global_cache_tags =
     284             :       tmpl::list<evolution::initial_data::Tags::InitialData>;
     285             : 
     286             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
     287             :             typename ArrayIndex, typename ActionList,
     288             :             typename ParallelComponent>
     289           0 :   static Parallel::iterable_action_return_t apply(
     290             :       db::DataBox<DbTagsList>& box,
     291             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
     292             :       Parallel::GlobalCache<Metavariables>& cache,
     293             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
     294             :       const ParallelComponent* const parallel_component) {
     295             :     // Dispatch to the correct `apply` overload based on type of initial data
     296             :     using initial_data_classes =
     297             :         tmpl::at<typename Metavariables::factory_creation::factory_classes,
     298             :                  evolution::initial_data::InitialData>;
     299             :     return call_with_dynamic_type<Parallel::iterable_action_return_t,
     300             :                                   initial_data_classes>(
     301             :         &db::get<evolution::initial_data::Tags::InitialData>(box),
     302             :         [&box, &cache, &parallel_component](const auto* const initial_data) {
     303             :           return apply(make_not_null(&box), *initial_data, cache,
     304             :                        parallel_component);
     305             :         });
     306             :   }
     307             : 
     308             :  private:
     309             :   // Numeric initial data
     310             :   template <typename DbTagsList, typename Metavariables,
     311             :             typename ParallelComponent>
     312           0 :   static Parallel::iterable_action_return_t apply(
     313             :       const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
     314             :       const NumericInitialData& initial_data,
     315             :       Parallel::GlobalCache<Metavariables>& cache,
     316             :       const ParallelComponent* const /*meta*/) {
     317             :     // Select the subset of the available variables that we want to read from
     318             :     // the volume data file
     319             :     tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
     320             :         importers::Tags::Selected, NumericInitialData::all_vars>>
     321             :         selected_fields{};
     322             :     initial_data.select_for_import(make_not_null(&selected_fields));
     323             :     // Dispatch loading the variables from the volume data file
     324             :     // - Not using `ckLocalBranch` here to make sure the simple action
     325             :     //   invocation is asynchronous.
     326             :     auto& reader_component = Parallel::get_parallel_component<
     327             :         importers::ElementDataReader<Metavariables>>(cache);
     328             :     Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
     329             :         3, NumericInitialData::all_vars, ParallelComponent>>(
     330             :         reader_component, initial_data.importer_options(),
     331             :         initial_data.volume_data_id(), std::move(selected_fields));
     332             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     333             :   }
     334             : 
     335             :   // "AnalyticData"-type initial data
     336             :   template <typename DbTagsList, typename InitialData, typename Metavariables,
     337             :             typename ParallelComponent>
     338           0 :   static Parallel::iterable_action_return_t apply(
     339             :       const gsl::not_null<db::DataBox<DbTagsList>*> box,
     340             :       const InitialData& initial_data,
     341             :       Parallel::GlobalCache<Metavariables>& /*cache*/,
     342             :       const ParallelComponent* const /*meta*/) {
     343             :     static constexpr size_t Dim = Metavariables::volume_dim;
     344             : 
     345             :     // Get ADM variables from analytic data / solution
     346             :     const auto& x =
     347             :         db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
     348             :     const auto adm_vars = evolution::Initialization::initial_data(
     349             :         initial_data, x, db::get<::Tags::Time>(*box),
     350             :         tmpl::list<gr::Tags::SpatialMetric<DataVector, Dim>,
     351             :                    gr::Tags::Lapse<DataVector>,
     352             :                    gr::Tags::Shift<DataVector, Dim>,
     353             :                    gr::Tags::ExtrinsicCurvature<DataVector, Dim>>{});
     354             :     const auto& spatial_metric =
     355             :         get<gr::Tags::SpatialMetric<DataVector, Dim>>(adm_vars);
     356             :     const auto& lapse = get<gr::Tags::Lapse<DataVector>>(adm_vars);
     357             :     const auto& shift = get<gr::Tags::Shift<DataVector, Dim>>(adm_vars);
     358             :     const auto& extrinsic_curvature =
     359             :         get<gr::Tags::ExtrinsicCurvature<DataVector, Dim>>(adm_vars);
     360             : 
     361             :     // Compute GH vars from ADM vars
     362             :     const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
     363             :     const auto& inv_jacobian =
     364             :         db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     365             :                                               Frame::Inertial>>(*box);
     366             :     db::mutate<gr::Tags::SpacetimeMetric<DataVector, Dim>,
     367             :                Tags::Pi<DataVector, Dim>, Tags::Phi<DataVector, Dim>>(
     368             :         &gh::initial_gh_variables_from_adm<Dim>, box, spatial_metric, lapse,
     369             :         shift, extrinsic_curvature, mesh, inv_jacobian);
     370             : 
     371             :     // No need to import numeric initial data, so we terminate the phase by
     372             :     // pausing the algorithm on this element
     373             :     return {Parallel::AlgorithmExecution::Pause, std::nullopt};
     374             :   }
     375             : };
     376             : 
     377             : /*!
     378             :  * \brief Receive numeric initial data loaded by gh::Actions::SetInitialData.
     379             :  *
     380             :  * Place this action in the action list after
     381             :  * gh::Actions::SetInitialData to wait until the data
     382             :  * for this element has arrived, and then transform the data to GH variables and
     383             :  * store it in the DataBox to be used as initial data.
     384             :  *
     385             :  * This action modifies the following tags in the DataBox:
     386             :  * - gr::Tags::SpacetimeMetric<DataVector, 3>
     387             :  * - gh::Tags::Pi<DataVector, 3>
     388             :  * - gh::Tags::Phi<DataVector, 3>
     389             :  */
     390           1 : struct ReceiveNumericInitialData {
     391           0 :   static constexpr size_t Dim = 3;
     392           0 :   using inbox_tags =
     393             :       tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
     394             : 
     395             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
     396             :             typename ActionList, typename ParallelComponent>
     397           0 :   static Parallel::iterable_action_return_t apply(
     398             :       db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
     399             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
     400             :       const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
     401             :       const ParallelComponent* const /*meta*/) {
     402             :     auto& inbox =
     403             :         tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
     404             :             inboxes);
     405             :     const auto& initial_data = dynamic_cast<const NumericInitialData&>(
     406             :         db::get<evolution::initial_data::Tags::InitialData>(box));
     407             :     const auto& volume_data_id = initial_data.volume_data_id();
     408             :     if (inbox.find(volume_data_id) == inbox.end()) {
     409             :       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
     410             :     }
     411             :     auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
     412             : 
     413             :     const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
     414             :     const auto& inv_jacobian =
     415             :         db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     416             :                                               Frame::Inertial>>(box);
     417             : 
     418             :     db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
     419             :                Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>(
     420             :         [&initial_data, &numeric_data, &mesh, &inv_jacobian](
     421             :             const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
     422             :             const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
     423             :             const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi) {
     424             :           initial_data.set_initial_data(spacetime_metric, pi, phi,
     425             :                                         make_not_null(&numeric_data), mesh,
     426             :                                         inv_jacobian);
     427             :         },
     428             :         make_not_null(&box));
     429             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     430             :   }
     431             : };
     432             : 
     433             : }  // namespace Actions
     434             : }  // namespace gh

Generated by: LCOV version 1.14