Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <array>
       8             : #include <cmath>
       9             : #include <cstddef>
      10             : #include <iterator>
      11             : #include <type_traits>
      12             : 
      13             : #include "DataStructures/Tensor/EagerMath/Determinant.hpp"
      14             : #include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
      15             : #include "DataStructures/Tensor/Tensor.hpp"
      16             : #include "Options/Options.hpp"
      17             : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      18             : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
      19             : #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp"
      20             : #include "PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp"
      21             : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
      22             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      23             : #include "Utilities/Gsl.hpp"
      24             : #include "Utilities/Serialization/CharmPupable.hpp"
      25             : #include "Utilities/TMPL.hpp"
      26             : #include "Utilities/TaggedTuple.hpp"
      27             : 
      28             : /// \cond
      29             : namespace PUP {
      30             : class er;
      31             : }  // namespace PUP
      32             : /// \endcond
      33             : 
      34             : namespace grmhd::AnalyticData {
      35             : /*!
      36             :  * \brief Magnetized fluid disk orbiting a Kerr black hole in the Kerr-Schild
      37             :  * Cartesian coordinates, but in a toroidal mesh defined from a torus map
      38             :  * (see grmhd::AnalyticData::SphericalTorus).
      39             :  */
      40           1 : class PolarMagnetizedFmDisk
      41             :     : public virtual evolution::initial_data::InitialData,
      42             :       public MarkAsAnalyticData,
      43             :       public AnalyticDataBase {
      44             :  public:
      45           0 :   struct DiskParameters {
      46           0 :     using type = MagnetizedFmDisk;
      47           0 :     static constexpr Options::String help = "Parameters for the disk.";
      48             :   };
      49             : 
      50           0 :   struct TorusParameters {
      51           0 :     using type = grmhd::AnalyticData::SphericalTorus;
      52           0 :     static constexpr Options::String help =
      53             :         "Parameters for the evolution region.";
      54             :   };
      55             : 
      56           0 :   using options = tmpl::list<DiskParameters, TorusParameters>;
      57             : 
      58           0 :   static constexpr Options::String help =
      59             :       "Magnetized Fishbone-Moncrief disk in polar coordinates.";
      60             : 
      61           0 :   PolarMagnetizedFmDisk() = default;
      62           0 :   PolarMagnetizedFmDisk(const PolarMagnetizedFmDisk& /*rhs*/) = default;
      63           0 :   PolarMagnetizedFmDisk& operator=(const PolarMagnetizedFmDisk& /*rhs*/) =
      64             :       default;
      65           0 :   PolarMagnetizedFmDisk(PolarMagnetizedFmDisk&& /*rhs*/) = default;
      66           0 :   PolarMagnetizedFmDisk& operator=(PolarMagnetizedFmDisk&& /*rhs*/) = default;
      67           0 :   ~PolarMagnetizedFmDisk() override = default;
      68             : 
      69           0 :   PolarMagnetizedFmDisk(MagnetizedFmDisk fm_disk,
      70             :                         grmhd::AnalyticData::SphericalTorus torus_map);
      71             : 
      72           0 :   auto get_clone() const
      73             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
      74             : 
      75             :   /// \cond
      76             :   explicit PolarMagnetizedFmDisk(CkMigrateMessage* msg);
      77             :   using PUP::able::register_constructor;
      78             :   WRAPPED_PUPable_decl_template(PolarMagnetizedFmDisk);
      79             :   /// \endcond
      80             : 
      81             :   /// The grmhd variables.
      82             :   ///
      83             :   /// \note The functions are optimized for retrieving the hydro variables
      84             :   /// before the metric variables.
      85             :   template <typename DataType, typename... Tags>
      86           1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
      87             :                                          tmpl::list<Tags...> /*meta*/) const {
      88             :     // In this function, we label the coordinates this solution works
      89             :     // in with Frame::BlockLogical, and the coordinates the wrapped
      90             :     // solution uses Frame::Inertial.  This means the input and output
      91             :     // have to be converted to the correct label.
      92             : 
      93             :     const tnsr::I<DataType, 3> observation_coordinates(torus_map_(x));
      94             : 
      95             :     using dependencies = tmpl::map<
      96             :         tmpl::pair<gr::AnalyticSolution<3>::DerivShift<DataType>,
      97             :                    gr::Tags::Shift<DataType, 3, Frame::Inertial>>,
      98             :         tmpl::pair<gr::AnalyticSolution<3>::DerivSpatialMetric<DataType>,
      99             :                    gr::Tags::SpatialMetric<DataType, 3, Frame::Inertial>>>;
     100             :     using required_tags = tmpl::remove_duplicates<
     101             :         tmpl::remove<tmpl::list<Tags..., tmpl::at<dependencies, Tags>...>,
     102             :                      tmpl::no_such_type_>>;
     103             : 
     104             :     auto observation_data =
     105             :         fm_disk_.variables(observation_coordinates, required_tags{});
     106             : 
     107             :     const auto jacobian = torus_map_.jacobian(x);
     108             :     const auto inv_jacobian = torus_map_.inv_jacobian(x);
     109             : 
     110             :     const auto change_frame = [this, &inv_jacobian, &jacobian, &x](
     111             :                                   const auto& data, auto tag) {
     112             :       using Tag = decltype(tag);
     113             :       auto result =
     114             :           transform::to_different_frame(get<Tag>(data), jacobian, inv_jacobian);
     115             : 
     116             :       if constexpr (std::is_same_v<
     117             :                         Tag, gr::AnalyticSolution<3>::DerivShift<DataType>>) {
     118             :         const auto deriv_inv_jacobian =
     119             :             torus_map_.derivative_of_inv_jacobian(x);
     120             :         const auto& shift =
     121             :             get<gr::Tags::Shift<DataType, 3, Frame::Inertial>>(data);
     122             :         for (size_t i = 0; i < 3; ++i) {
     123             :           for (size_t j = 0; j < 3; ++j) {
     124             :             for (size_t k = 0; k < 3; ++k) {
     125             :               result.get(i, j) +=
     126             :                   deriv_inv_jacobian.get(j, k, i) * shift.get(k);
     127             :             }
     128             :           }
     129             :         }
     130             :       } else if constexpr (std::is_same_v<
     131             :                                Tag, gr::AnalyticSolution<3>::DerivSpatialMetric<
     132             :                                         DataType>>) {
     133             :         const auto hessian = torus_map_.hessian(x);
     134             :         const auto& spatial_metric =
     135             :             get<gr::Tags::SpatialMetric<DataType, 3, Frame::Inertial>>(data);
     136             :         for (size_t i = 0; i < 3; ++i) {
     137             :           for (size_t j = 0; j < 3; ++j) {
     138             :             for (size_t k = j; k < 3; ++k) {
     139             :               for (size_t l = 0; l < 3; ++l) {
     140             :                 for (size_t m = 0; m < 3; ++m) {
     141             :                   result.get(i, j, k) +=
     142             :                       (hessian.get(l, j, i) * jacobian.get(m, k) +
     143             :                        hessian.get(l, k, i) * jacobian.get(m, j)) *
     144             :                       spatial_metric.get(l, m);
     145             :                 }
     146             :               }
     147             :             }
     148             :           }
     149             :         }
     150             :       } else if constexpr (std::is_same_v<
     151             :                                Tag, gr::Tags::SqrtDetSpatialMetric<DataType>>) {
     152             :         get(result) *= abs(get(determinant(jacobian)));
     153             :       }
     154             : 
     155             :       typename Tag::type result_with_replaced_frame{};
     156             :       std::copy(std::move_iterator(result.begin()),
     157             :                 std::move_iterator(result.end()),
     158             :                 result_with_replaced_frame.begin());
     159             :       return result_with_replaced_frame;
     160             :     };
     161             : 
     162             :     return {change_frame(observation_data, Tags{})...};
     163             :   }
     164             : 
     165           0 :   using equation_of_state_type = MagnetizedFmDisk::equation_of_state_type;
     166           0 :   const equation_of_state_type& equation_of_state() const {
     167             :     return fm_disk_.equation_of_state();
     168             :   }
     169             : 
     170             :   // NOLINTNEXTLINE(google-runtime-references)
     171           0 :   void pup(PUP::er& p) override;
     172             : 
     173             :  private:
     174           0 :   friend bool operator==(const PolarMagnetizedFmDisk& lhs,
     175             :                          const PolarMagnetizedFmDisk& rhs);
     176             : 
     177           0 :   MagnetizedFmDisk fm_disk_;
     178           0 :   grmhd::AnalyticData::SphericalTorus torus_map_;
     179             : };
     180             : 
     181           0 : bool operator!=(const PolarMagnetizedFmDisk& lhs,
     182             :                 const PolarMagnetizedFmDisk& rhs);
     183             : }  // namespace grmhd::AnalyticData