SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/GrMhd/GhValenciaDivClean/Actions - SetInitialData.hpp Hit Total Coverage
Commit: 3c2e9d3ed337bca2146eee9de07432e292a38c3a Lines: 4 39 10.3 %
Date: 2024-06-11 22:56:19
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 <string>
       8             : #include <variant>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "Domain/Structure/ElementId.hpp"
      15             : #include "Evolution/DgSubcell/ActiveGrid.hpp"
      16             : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
      17             : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
      18             : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
      19             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      20             : #include "Evolution/NumericInitialData.hpp"
      21             : #include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
      22             : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
      23             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Actions/NumericInitialData.hpp"
      24             : #include "IO/Importers/Actions/ReadVolumeData.hpp"
      25             : #include "IO/Importers/ElementDataReader.hpp"
      26             : #include "IO/Importers/Tags.hpp"
      27             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      28             : #include "Options/String.hpp"
      29             : #include "Parallel/AlgorithmExecution.hpp"
      30             : #include "Parallel/GlobalCache.hpp"
      31             : #include "Parallel/Invoke.hpp"
      32             : #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
      33             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      34             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      35             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      36             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      37             : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
      38             : #include "Utilities/ErrorHandling/Error.hpp"
      39             : #include "Utilities/Gsl.hpp"
      40             : #include "Utilities/TMPL.hpp"
      41             : #include "Utilities/TaggedTuple.hpp"
      42             : 
      43             : /// \cond
      44             : namespace Tags {
      45             : struct Time;
      46             : }  // namespace Tags
      47             : /// \endcond
      48             : 
      49           1 : namespace grmhd::GhValenciaDivClean {
      50             : 
      51             : /*!
      52             :  * \brief Numeric initial data loaded from volume data files
      53             :  *
      54             :  * This class can be factory-created in the input file to start an evolution
      55             :  * from numeric initial data. It selects the set of GH variables to load from
      56             :  * the volume data file (ADM or GH variables), the set of hydro variables, and
      57             :  * allows to set constant values for some of the hydro variables. This class
      58             :  * mostly combines the `gh::NumericInitialData` and
      59             :  * `grmhd::ValenciaDivClean::NumericInitialData` classes.
      60             :  */
      61           1 : class NumericInitialData : public evolution::initial_data::InitialData,
      62             :                            public evolution::NumericInitialData {
      63             :  private:
      64           0 :   using GhNumericId = gh::NumericInitialData;
      65           0 :   using HydroNumericId = grmhd::ValenciaDivClean::NumericInitialData;
      66             : 
      67             :  public:
      68           0 :   using all_vars =
      69             :       tmpl::append<GhNumericId::all_vars, HydroNumericId::all_vars>;
      70             : 
      71           0 :   struct GhVariables : GhNumericId::Variables {};
      72           0 :   struct HydroVariables : HydroNumericId::Variables {};
      73             : 
      74           0 :   using options =
      75             :       tmpl::list<importers::OptionTags::FileGlob,
      76             :                  importers::OptionTags::Subgroup,
      77             :                  importers::OptionTags::ObservationValue,
      78             :                  importers::OptionTags::EnableInterpolation, GhVariables,
      79             :                  HydroVariables, HydroNumericId::DensityCutoff>;
      80             : 
      81           0 :   static constexpr Options::String help =
      82             :       "Numeric initial data loaded from volume data files";
      83             : 
      84           0 :   NumericInitialData() = default;
      85           0 :   NumericInitialData(const NumericInitialData& rhs) = default;
      86           0 :   NumericInitialData& operator=(const NumericInitialData& rhs) = default;
      87           0 :   NumericInitialData(NumericInitialData&& /*rhs*/) = default;
      88           0 :   NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
      89           0 :   ~NumericInitialData() = default;
      90             : 
      91             :   /// \cond
      92             :   explicit NumericInitialData(CkMigrateMessage* msg);
      93             :   using PUP::able::register_constructor;
      94             :   WRAPPED_PUPable_decl_template(NumericInitialData);
      95             :   /// \endcond
      96             : 
      97           0 :   std::unique_ptr<evolution::initial_data::InitialData> get_clone()
      98             :       const override {
      99             :     return std::make_unique<NumericInitialData>(*this);
     100             :   }
     101             : 
     102           0 :   NumericInitialData(
     103             :       std::string file_glob, std::string subfile_name,
     104             :       std::variant<double, importers::ObservationSelector> observation_value,
     105             :       bool enable_interpolation,
     106             :       typename GhNumericId::Variables::type gh_selected_variables,
     107             :       typename HydroNumericId::Variables::type hydro_selected_variables,
     108             :       double density_cutoff);
     109             : 
     110           0 :   const importers::ImporterOptions& importer_options() const {
     111             :     return gh_numeric_id_.importer_options();
     112             :   }
     113             : 
     114           0 :   const GhNumericId& gh_numeric_id() const { return gh_numeric_id_; }
     115             : 
     116           0 :   const HydroNumericId& hydro_numeric_id() const { return hydro_numeric_id_; }
     117             : 
     118           0 :   size_t volume_data_id() const;
     119             : 
     120             :   template <typename... AllTags>
     121           0 :   void select_for_import(
     122             :       const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
     123             :     gh_numeric_id_.select_for_import(fields);
     124             :     hydro_numeric_id_.select_for_import(fields);
     125             :   }
     126             : 
     127             :   template <typename... AllTags, size_t ThermodynamicDim>
     128           0 :   void set_initial_data(
     129             :       const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
     130             :       const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
     131             :       const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
     132             :       const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
     133             :       const gsl::not_null<Scalar<DataVector>*> electron_fraction,
     134             :       const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     135             :       const gsl::not_null<tnsr::I<DataVector, 3>*> spatial_velocity,
     136             :       const gsl::not_null<tnsr::I<DataVector, 3>*> magnetic_field,
     137             :       const gsl::not_null<Scalar<DataVector>*> div_cleaning_field,
     138             :       const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
     139             :       const gsl::not_null<Scalar<DataVector>*> pressure,
     140             :       const gsl::not_null<Scalar<DataVector>*> temperature,
     141             :       const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
     142             :       const Mesh<3>& mesh,
     143             :       const InverseJacobian<DataVector, 3, Frame::ElementLogical,
     144             :                             Frame::Inertial>& inv_jacobian,
     145             :       const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
     146             :           equation_of_state) const {
     147             :     gh_numeric_id_.set_initial_data(spacetime_metric, pi, phi, numeric_data,
     148             :                                     mesh, inv_jacobian);
     149             :     const auto spatial_metric = gr::spatial_metric(*spacetime_metric);
     150             :     const auto inv_spatial_metric =
     151             :         determinant_and_inverse(spatial_metric).second;
     152             :     hydro_numeric_id_.set_initial_data(
     153             :         rest_mass_density, electron_fraction, specific_internal_energy,
     154             :         spatial_velocity, magnetic_field, div_cleaning_field, lorentz_factor,
     155             :         pressure, temperature, numeric_data, inv_spatial_metric,
     156             :         equation_of_state);
     157             :   }
     158             : 
     159           0 :   void pup(PUP::er& p) override;
     160             : 
     161           0 :   friend bool operator==(const NumericInitialData& lhs,
     162             :                          const NumericInitialData& rhs);
     163             : 
     164             :  private:
     165           0 :   GhNumericId gh_numeric_id_{};
     166           0 :   HydroNumericId hydro_numeric_id_{};
     167             : };
     168             : 
     169           0 : namespace Actions {
     170             : 
     171             : /*!
     172             :  * \brief Dispatch loading numeric initial data from files.
     173             :  *
     174             :  * Place this action before
     175             :  * grmhd::GhValenciaDivClean::Actions::SetNumericInitialData in the action list.
     176             :  * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
     177             :  * invoked by this action.
     178             :  */
     179           1 : struct SetInitialData {
     180           0 :   using const_global_cache_tags =
     181             :       tmpl::list<evolution::initial_data::Tags::InitialData>;
     182             : 
     183             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
     184             :             typename ArrayIndex, typename ActionList,
     185             :             typename ParallelComponent>
     186           0 :   static Parallel::iterable_action_return_t apply(
     187             :       db::DataBox<DbTagsList>& box,
     188             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
     189             :       Parallel::GlobalCache<Metavariables>& cache,
     190             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
     191             :       const ParallelComponent* const parallel_component) {
     192             :     // Dispatch to the correct `apply` overload based on type of initial data
     193             :     using initial_data_classes =
     194             :         tmpl::at<typename Metavariables::factory_creation::factory_classes,
     195             :                  evolution::initial_data::InitialData>;
     196             :     return call_with_dynamic_type<Parallel::iterable_action_return_t,
     197             :                                   initial_data_classes>(
     198             :         &db::get<evolution::initial_data::Tags::InitialData>(box),
     199             :         [&box, &cache, &parallel_component](const auto* const initial_data) {
     200             :           return apply(make_not_null(&box), *initial_data, cache,
     201             :                        parallel_component);
     202             :         });
     203             :   }
     204             : 
     205             :  private:
     206           0 :   static constexpr size_t Dim = 3;
     207             : 
     208             :   // Numeric initial data
     209             :   template <typename DbTagsList, typename Metavariables,
     210             :             typename ParallelComponent>
     211           0 :   static Parallel::iterable_action_return_t apply(
     212             :       const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
     213             :       const NumericInitialData& initial_data,
     214             :       Parallel::GlobalCache<Metavariables>& cache,
     215             :       const ParallelComponent* const /*meta*/) {
     216             :     // Select the subset of the available variables that we want to read from
     217             :     // the volume data file
     218             :     tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
     219             :         importers::Tags::Selected, NumericInitialData::all_vars>>
     220             :         selected_fields{};
     221             :     initial_data.select_for_import(make_not_null(&selected_fields));
     222             :     // Dispatch loading the variables from the volume data file
     223             :     // - Not using `ckLocalBranch` here to make sure the simple action
     224             :     //   invocation is asynchronous.
     225             :     auto& reader_component = Parallel::get_parallel_component<
     226             :         importers::ElementDataReader<Metavariables>>(cache);
     227             :     Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
     228             :         3, NumericInitialData::all_vars, ParallelComponent>>(
     229             :         reader_component, initial_data.importer_options(),
     230             :         initial_data.volume_data_id(), std::move(selected_fields));
     231             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     232             :   }
     233             : 
     234             :   // "AnalyticData"-type initial data
     235             :   template <typename DbTagsList, typename InitialData, typename Metavariables,
     236             :             typename ParallelComponent>
     237           0 :   static Parallel::iterable_action_return_t apply(
     238             :       const gsl::not_null<db::DataBox<DbTagsList>*> box,
     239             :       const InitialData& initial_data,
     240             :       Parallel::GlobalCache<Metavariables>& /*cache*/,
     241             :       const ParallelComponent* const /*meta*/) {
     242             :     // Get ADM + hydro variables from analytic data / solution
     243             :     const auto& [coords, mesh, inv_jacobian] = [&box]() {
     244             :       if constexpr (db::tag_is_retrievable_v<
     245             :                         evolution::dg::subcell::Tags::ActiveGrid,
     246             :                         db::DataBox<DbTagsList>>) {
     247             :         const bool on_subcell =
     248             :             db::get<evolution::dg::subcell::Tags::ActiveGrid>(*box) ==
     249             :             evolution::dg::subcell::ActiveGrid::Subcell;
     250             :         if (on_subcell) {
     251             :           return std::forward_as_tuple(
     252             :               db::get<evolution::dg::subcell::Tags::Coordinates<
     253             :                   Dim, Frame::Inertial>>(*box),
     254             :               db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box),
     255             :               db::get<evolution::dg::subcell::fd::Tags::
     256             :                           InverseJacobianLogicalToInertial<Dim>>(*box));
     257             :         }
     258             :       }
     259             :       return std::forward_as_tuple(
     260             :           db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box),
     261             :           db::get<domain::Tags::Mesh<Dim>>(*box),
     262             :           db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     263             :                                                 Frame::Inertial>>(*box));
     264             :     }();
     265             :     auto vars = evolution::Initialization::initial_data(
     266             :         initial_data, coords, db::get<::Tags::Time>(*box),
     267             :         tmpl::append<tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
     268             :                                 gr::Tags::Lapse<DataVector>,
     269             :                                 gr::Tags::Shift<DataVector, 3>,
     270             :                                 gr::Tags::ExtrinsicCurvature<DataVector, 3>>,
     271             :                      hydro::grmhd_tags<DataVector>>{});
     272             :     const auto& spatial_metric =
     273             :         get<gr::Tags::SpatialMetric<DataVector, 3>>(vars);
     274             :     const auto& lapse = get<gr::Tags::Lapse<DataVector>>(vars);
     275             :     const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(vars);
     276             :     const auto& extrinsic_curvature =
     277             :         get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(vars);
     278             : 
     279             :     // Compute GH vars from ADM vars
     280             :     db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
     281             :                gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>>(
     282             :         &gh::initial_gh_variables_from_adm<3>, box, spatial_metric, lapse,
     283             :         shift, extrinsic_curvature, mesh, inv_jacobian);
     284             : 
     285             :     // Move hydro vars directly into the DataBox
     286             :     tmpl::for_each<hydro::grmhd_tags<DataVector>>(
     287             :         [&box, &vars](const auto tag_v) {
     288             :           using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     289             :           Initialization::mutate_assign<tmpl::list<tag>>(
     290             :               box, std::move(get<tag>(vars)));
     291             :         });
     292             : 
     293             :     // No need to import numeric initial data, so we terminate the phase by
     294             :     // pausing the algorithm on this element
     295             :     return {Parallel::AlgorithmExecution::Pause, std::nullopt};
     296             :   }
     297             : };
     298             : 
     299             : /*!
     300             :  * \brief Receive numeric initial data loaded by
     301             :  * grmhd::GhValenciaDivClean::Actions::SetInitialData.
     302             :  *
     303             :  * Place this action in the action list after
     304             :  * grmhd::GhValenciaDivClean::Actions::SetInitialData to wait until the
     305             :  * data for this element has arrived, and then compute the GH variables and the
     306             :  * remaining primitive variables and store them in the DataBox to be used as
     307             :  * initial data.
     308             :  *
     309             :  * This action modifies the GH system tags (spacetime metric, pi, phi) and the
     310             :  * tags listed in `hydro::grmhd_tags` in the DataBox (i.e., the hydro
     311             :  * primitives). It does not modify conservative variables, so it relies on a
     312             :  * primitive-to-conservative update in the action list before the evolution can
     313             :  * start.
     314             :  *
     315             :  * \requires This action requires an equation of state, which is retrieved from
     316             :  * the DataBox as `hydro::Tags::GrmhdEquationOfState`.
     317             :  */
     318           1 : struct ReceiveNumericInitialData {
     319           0 :   static constexpr size_t Dim = 3;
     320           0 :   using inbox_tags =
     321             :       tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
     322             : 
     323             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
     324             :             typename ActionList, typename ParallelComponent>
     325           0 :   static Parallel::iterable_action_return_t apply(
     326             :       db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
     327             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
     328             :       const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
     329             :       const ParallelComponent* const /*meta*/) {
     330             :     auto& inbox =
     331             :         tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
     332             :             inboxes);
     333             :     const auto& initial_data = dynamic_cast<const NumericInitialData&>(
     334             :         db::get<evolution::initial_data::Tags::InitialData>(box));
     335             :     const size_t volume_data_id = initial_data.volume_data_id();
     336             :     if (inbox.find(volume_data_id) == inbox.end()) {
     337             :       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
     338             :     }
     339             :     auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
     340             : 
     341             :     const auto& [mesh, inv_jacobian] = [&box]() {
     342             :       if constexpr (db::tag_is_retrievable_v<
     343             :                         evolution::dg::subcell::Tags::ActiveGrid,
     344             :                         db::DataBox<DbTagsList>>) {
     345             :         const bool on_subcell =
     346             :             db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) ==
     347             :             evolution::dg::subcell::ActiveGrid::Subcell;
     348             :         if (on_subcell) {
     349             :           return std::forward_as_tuple(
     350             :               db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box),
     351             :               db::get<evolution::dg::subcell::fd::Tags::
     352             :                           InverseJacobianLogicalToInertial<Dim>>(box));
     353             :         }
     354             :       }
     355             :       return std::forward_as_tuple(
     356             :           db::get<domain::Tags::Mesh<Dim>>(box),
     357             :           db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     358             :                                                 Frame::Inertial>>(box));
     359             :     }();
     360             :     const auto& equation_of_state =
     361             :         db::get<hydro::Tags::GrmhdEquationOfState>(box);
     362             : 
     363             :     db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
     364             :                gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>,
     365             :                hydro::Tags::RestMassDensity<DataVector>,
     366             :                hydro::Tags::ElectronFraction<DataVector>,
     367             :                hydro::Tags::SpecificInternalEnergy<DataVector>,
     368             :                hydro::Tags::SpatialVelocity<DataVector, 3>,
     369             :                hydro::Tags::MagneticField<DataVector, 3>,
     370             :                hydro::Tags::DivergenceCleaningField<DataVector>,
     371             :                hydro::Tags::LorentzFactor<DataVector>,
     372             :                hydro::Tags::Pressure<DataVector>,
     373             :                hydro::Tags::Temperature<DataVector>>(
     374             :         [&initial_data, &numeric_data, &equation_of_state](
     375             :             const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
     376             :             const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
     377             :             const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
     378             :             const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
     379             :             const gsl::not_null<Scalar<DataVector>*> electron_fraction,
     380             :             const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     381             :             const gsl::not_null<tnsr::I<DataVector, 3>*> spatial_velocity,
     382             :             const gsl::not_null<tnsr::I<DataVector, 3>*> magnetic_field,
     383             :             const gsl::not_null<Scalar<DataVector>*> div_cleaning_field,
     384             :             const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
     385             :             const gsl::not_null<Scalar<DataVector>*> pressure,
     386             :             const gsl::not_null<Scalar<DataVector>*> temperature,
     387             :             const auto& local_mesh, const auto& local_inv_jacobian) {
     388             :           initial_data.set_initial_data(
     389             :               spacetime_metric, pi, phi, rest_mass_density, electron_fraction,
     390             :               specific_internal_energy, spatial_velocity, magnetic_field,
     391             :               div_cleaning_field, lorentz_factor, pressure, temperature,
     392             :               make_not_null(&numeric_data), local_mesh, local_inv_jacobian,
     393             :               equation_of_state);
     394             :         },
     395             :         make_not_null(&box), mesh, inv_jacobian);
     396             : 
     397             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     398             :   }
     399             : };
     400             : 
     401             : }  // namespace Actions
     402             : 
     403             : }  // namespace grmhd::GhValenciaDivClean

Generated by: LCOV version 1.14