Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines tags related to domain quantities
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : #include <cstddef>
      11             : #include <memory>
      12             : #include <string>
      13             : #include <unordered_map>
      14             : #include <unordered_set>
      15             : #include <vector>
      16             : 
      17             : #include "DataStructures/DataBox/Subitems.hpp"
      18             : #include "DataStructures/DataBox/Tag.hpp"
      19             : #include "DataStructures/DataBox/TagName.hpp"
      20             : #include "DataStructures/Tensor/EagerMath/Determinant.hpp"
      21             : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
      22             : #include "DataStructures/Tensor/TypeAliases.hpp"
      23             : #include "DataStructures/Variables.hpp"
      24             : #include "Domain/Structure/Direction.hpp"
      25             : #include "Domain/Structure/Element.hpp"
      26             : #include "Utilities/GetOutput.hpp"
      27             : #include "Utilities/Gsl.hpp"
      28             : #include "Utilities/NoSuchType.hpp"
      29             : #include "Utilities/Requires.hpp"
      30             : #include "Utilities/TMPL.hpp"
      31             : #include "Utilities/TypeTraits.hpp"
      32             : #include "Utilities/TypeTraits/IsA.hpp"
      33             : 
      34             : /// \cond
      35             : class DataVector;
      36             : template <size_t VolumeDim>
      37             : class Domain;
      38             : template <size_t VolumeDim, typename Frame>
      39             : class ElementMap;
      40             : template <size_t VolumeDim>
      41             : class Mesh;
      42             : /// \endcond
      43             : 
      44             : namespace domain {
      45             : /// \ingroup ComputationalDomainGroup
      46             : /// \brief %Tags for the domain.
      47             : namespace Tags {
      48             : /// \ingroup DataBoxTagsGroup
      49             : /// \ingroup ComputationalDomainGroup
      50             : /// The ::Element associated with the DataBox
      51             : template <size_t VolumeDim>
      52           1 : struct Element : db::SimpleTag {
      53           0 :   using type = ::Element<VolumeDim>;
      54             : };
      55             : 
      56             : /// \ingroup DataBoxTagsGroup
      57             : /// \ingroup ComputationalDomainGroup
      58             : /// \brief The computational grid of the Element in the DataBox
      59             : /// \details The corresponding interface tag uses Mesh::slice_through to compute
      60             : /// the mesh on the face of the element.
      61             : template <size_t VolumeDim>
      62           1 : struct Mesh : db::SimpleTag {
      63           0 :   using type = ::Mesh<VolumeDim>;
      64             : };
      65             : 
      66             : /// \ingroup DataBoxTagsGroup
      67             : /// \ingroup ComputationalDomainGroup
      68             : /// The coordinate map from the ElementLogical frame to the TargetFrame
      69             : template <size_t VolumeDim, typename TargetFrame = Frame::Inertial>
      70           1 : struct ElementMap : db::SimpleTag {
      71           0 :   static constexpr size_t dim = VolumeDim;
      72           0 :   using target_frame = TargetFrame;
      73           0 :   using source_frame = Frame::ElementLogical;
      74             : 
      75           0 :   static std::string name() {
      76             :     return "ElementMap(" + get_output(TargetFrame{}) + ")";
      77             :   }
      78           0 :   using type = ::ElementMap<VolumeDim, TargetFrame>;
      79             : };
      80             : 
      81             : /// \ingroup DataBoxTagsGroup
      82             : /// \ingroup ComputationalDomainGroup
      83             : /// The coordinates in a given frame.
      84             : template <size_t Dim, typename Frame>
      85           1 : struct Coordinates : db::SimpleTag {
      86           0 :   static std::string name() { return get_output(Frame{}) + "Coordinates"; }
      87           0 :   using type = tnsr::I<DataVector, Dim, Frame>;
      88             : };
      89             : 
      90             : /// \ingroup DataBoxTagsGroup
      91             : /// \ingroup ComputationalDomainGroup
      92             : /// The coordinates in the target frame of `MapTag`. The `SourceCoordsTag`'s
      93             : /// frame must be the source frame of `MapTag`
      94             : template <class MapTag, class SourceCoordsTag,
      95             :           template <size_t, class> class CoordinatesTag = Coordinates>
      96           1 : struct MappedCoordinates
      97             :     : CoordinatesTag<MapTag::dim, typename MapTag::target_frame>,
      98             :       db::ComputeTag {
      99           0 :   using base = CoordinatesTag<MapTag::dim, typename MapTag::target_frame>;
     100           0 :   using return_type = typename base::type;
     101           0 :   using argument_tags = tmpl::list<MapTag, SourceCoordsTag>;
     102           0 :   static constexpr auto function(
     103             :       const gsl::not_null<return_type*> target_coords,
     104             :       const typename MapTag::type& element_map,
     105             :       const tnsr::I<DataVector, MapTag::dim, typename MapTag::source_frame>&
     106             :           source_coords) {
     107             :     if (source_coords.get(0).size() == 0) {
     108             :       // When subcell is not supported on an element due to topology, the
     109             :       // subcell mesh is
     110             :       // {0, Basis::Uninitialized, Quadrature::Uninitialized}. It is never
     111             :       // used, but the compute items still need to run
     112             :       return;
     113             :     }
     114             :     *target_coords = element_map(source_coords);
     115             :   }
     116             : };
     117             : 
     118             : /// \ingroup DataBoxTagsGroup
     119             : /// \ingroup ComputationalDomainGroup
     120             : /// \brief The inverse Jacobian from the source frame to the target frame.
     121             : ///
     122             : /// Specifically, \f$\partial x^{\bar{i}} / \partial x^i\f$, where \f$\bar{i}\f$
     123             : /// denotes the source frame and \f$i\f$ denotes the target frame.
     124             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
     125           1 : struct InverseJacobian : db::SimpleTag {
     126           0 :   static std::string name() {
     127             :     return "InverseJacobian(" + get_output(SourceFrame{}) + "," +
     128             :            get_output(TargetFrame{}) + ")";
     129             :   }
     130           0 :   using type = ::InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>;
     131             : };
     132             : 
     133             : /// \ingroup DataBoxTagsGroup
     134             : /// \ingroup ComputationalDomainGroup
     135             : /// Computes the inverse Jacobian of the map held by `MapTag` at the coordinates
     136             : /// held by `SourceCoordsTag`. The coordinates must be in the source frame of
     137             : /// the map.
     138             : template <typename MapTag, typename SourceCoordsTag>
     139           1 : struct InverseJacobianCompute
     140             :     : InverseJacobian<MapTag::dim, typename MapTag::source_frame,
     141             :                       typename MapTag::target_frame>,
     142             :       db::ComputeTag {
     143           0 :   using base = InverseJacobian<MapTag::dim, typename MapTag::source_frame,
     144             :                                typename MapTag::target_frame>;
     145           0 :   using return_type = typename base::type;
     146           0 :   using argument_tags = tmpl::list<MapTag, SourceCoordsTag>;
     147           0 :   static constexpr auto function(
     148             :       const gsl::not_null<return_type*> inv_jacobian,
     149             :       const typename MapTag::type& element_map,
     150             :       const tnsr::I<DataVector, MapTag::dim, typename MapTag::source_frame>&
     151             :           source_coords) {
     152             :     *inv_jacobian = element_map.inv_jacobian(source_coords);
     153             :   }
     154             : };
     155             : 
     156             : /// \ingroup DataBoxTagsGroup
     157             : /// \ingroup ComputationalDomainGroup
     158             : /// \brief The Jacobian from the source frame to the target frame.
     159             : ///
     160             : /// Specifically, \f$\partial x^{i} / \partial \xi^{\bar{i}}\f$, where
     161             : /// \f$\xi^\bar{i}\f$ denotes the source frame and \f$x^i\f$ denotes the target
     162             : /// frame.
     163             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
     164           1 : struct Jacobian : db::SimpleTag {
     165           0 :   static std::string name() {
     166             :     return "Jacobian(" + get_output(SourceFrame{}) + "," +
     167             :            get_output(TargetFrame{}) + ")";
     168             :   }
     169           0 :   using type = ::Jacobian<DataVector, Dim, SourceFrame, TargetFrame>;
     170             : };
     171             : 
     172             : /// \ingroup DataBoxTagsGroup
     173             : /// \ingroup ComputationalDomainGroup
     174             : /// Computes the Jacobian of the map from the `InverseJacobian<Dim, SourceFrame,
     175             : /// TargetFrame>` tag.
     176             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
     177           1 : struct JacobianCompute : Jacobian<Dim, SourceFrame, TargetFrame>,
     178             :                          db::ComputeTag {
     179           0 :   using base = Jacobian<Dim, SourceFrame, TargetFrame>;
     180           0 :   using return_type = typename base::type;
     181           0 :   using argument_tags =
     182             :       tmpl::list<InverseJacobian<Dim, SourceFrame, TargetFrame>>;
     183           0 :   static constexpr auto function(
     184             :       const gsl::not_null<return_type*> jacobian,
     185             :       const ::InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>&
     186             :           inv_jac) {
     187             :     *jacobian = determinant_and_inverse(inv_jac).second;
     188             :   }
     189             : };
     190             : 
     191             : /// \ingroup DataBoxTagsGroup
     192             : /// \ingroup ComputationalDomainGroup
     193             : /// \brief The determinant of the inverse Jacobian from the source frame to the
     194             : /// target frame.
     195             : template <typename SourceFrame, typename TargetFrame>
     196           1 : struct DetInvJacobian : db::SimpleTag {
     197           0 :   using type = Scalar<DataVector>;
     198           0 :   static std::string name() {
     199             :     return "DetInvJacobian(" + get_output(SourceFrame{}) + "," +
     200             :            get_output(TargetFrame{}) + ")";
     201             :   }
     202             : };
     203             : 
     204             : /// \ingroup DataBoxTagsGroup
     205             : /// \ingroup ComputationalDomainGroup
     206             : /// Computes the determinant of the inverse Jacobian.
     207             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
     208           1 : struct DetInvJacobianCompute : db::ComputeTag,
     209             :                                DetInvJacobian<SourceFrame, TargetFrame> {
     210           0 :   using base = DetInvJacobian<SourceFrame, TargetFrame>;
     211           0 :   using return_type = typename base::type;
     212           0 :   using argument_tags =
     213             :       tmpl::list<InverseJacobian<Dim, SourceFrame, TargetFrame>>;
     214           0 :   static void function(const gsl::not_null<return_type*> det_inv_jac,
     215             :                        const ::InverseJacobian<DataVector, Dim, SourceFrame,
     216             :                                                TargetFrame>& inv_jac) {
     217             :     determinant(det_inv_jac, inv_jac);
     218             :   }
     219             : };
     220             : 
     221             : /// \ingroup DataBoxTagsGroup
     222             : /// \ingroup ComputationalDomainGroup
     223             : /// \brief The determinant of the Jacobian from the source frame to the target
     224             : /// frame.
     225             : template <typename SourceFrame, typename TargetFrame>
     226           1 : struct DetJacobian : db::SimpleTag {
     227           0 :   using type = Scalar<DataVector>;
     228           0 :   static std::string name() {
     229             :     return "DetJacobian(" + get_output(SourceFrame{}) + "," +
     230             :            get_output(TargetFrame{}) + ")";
     231             :   }
     232             : };
     233             : 
     234             : /// \ingroup DataBoxTagsGroup
     235             : /// \ingroup ComputationalDomainGroup
     236             : /// \brief The inverse Jacobian times the determinant of the Jacobian.
     237             : ///
     238             : /// This quantity is divergence-free analytically. See
     239             : /// `::dg::metric_identity_det_jac_times_inv_jac` for more information.
     240             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
     241           1 : struct DetTimesInvJacobian : db::SimpleTag {
     242           0 :   using type = ::InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>;
     243           0 :   static std::string name() {
     244             :     return "DetTimesInvJacobian(" + get_output(SourceFrame{}) + "," +
     245             :            get_output(TargetFrame{}) + ")";
     246             :   }
     247             : };
     248             : 
     249             : /// @{
     250             : /// \ingroup DataBoxTagsGroup
     251             : /// \ingroup ComputationalDomainGroup
     252             : /// The set of directions to neighboring Elements
     253             : template <size_t VolumeDim>
     254           1 : struct InternalDirections : db::SimpleTag {
     255           0 :   static constexpr size_t volume_dim = VolumeDim;
     256           0 :   using type = std::unordered_set<::Direction<VolumeDim>>;
     257             : };
     258             : 
     259             : template <size_t VolumeDim>
     260           0 : struct InternalDirectionsCompute : InternalDirections<VolumeDim>,
     261             :                                    db::ComputeTag {
     262           0 :   static constexpr size_t volume_dim = VolumeDim;
     263           0 :   using base = InternalDirections<VolumeDim>;
     264           0 :   using return_type = std::unordered_set<::Direction<VolumeDim>>;
     265           0 :   using argument_tags = tmpl::list<Element<VolumeDim>>;
     266           0 :   static void function(const gsl::not_null<return_type*> directions,
     267             :                        const ::Element<VolumeDim>& element) {
     268             :     for (const auto& direction_neighbors : element.neighbors()) {
     269             :       directions->insert(direction_neighbors.first);
     270             :     }
     271             :   }
     272             : };
     273             : /// @}
     274             : 
     275             : /// @{
     276             : /// \ingroup DataBoxTagsGroup
     277             : /// \ingroup ComputationalDomainGroup
     278             : /// The set of directions which correspond to external boundaries.
     279             : /// Used for representing data on the interior side of the external boundary
     280             : /// faces.
     281             : template <size_t VolumeDim>
     282           1 : struct BoundaryDirectionsInterior : db::SimpleTag {
     283           0 :   static constexpr size_t volume_dim = VolumeDim;
     284           0 :   using type = std::unordered_set<::Direction<VolumeDim>>;
     285             : };
     286             : 
     287             : template <size_t VolumeDim>
     288           0 : struct BoundaryDirectionsInteriorCompute
     289             :     : BoundaryDirectionsInterior<VolumeDim>,
     290             :       db::ComputeTag {
     291           0 :   static constexpr size_t volume_dim = VolumeDim;
     292           0 :   using base = BoundaryDirectionsInterior<VolumeDim>;
     293           0 :   using return_type = std::unordered_set<::Direction<VolumeDim>>;
     294           0 :   using argument_tags = tmpl::list<Element<VolumeDim>>;
     295           0 :   static void function(const gsl::not_null<return_type*> directions,
     296             :                        const ::Element<VolumeDim>& element) {
     297             :     *directions = element.external_boundaries();
     298             :   }
     299             : };
     300             : /// @}
     301             : 
     302             : /// @{
     303             : /// \ingroup DataBoxTagsGroup
     304             : /// \ingroup ComputationalDomainGroup
     305             : /// The set of directions which correspond to external boundaries. To be used
     306             : /// to represent data which exists on the exterior side of the external boundary
     307             : /// faces.
     308             : template <size_t VolumeDim>
     309           1 : struct BoundaryDirectionsExterior : db::SimpleTag {
     310           0 :   static constexpr size_t volume_dim = VolumeDim;
     311           0 :   using type = std::unordered_set<::Direction<VolumeDim>>;
     312             : };
     313             : 
     314             : template <size_t VolumeDim>
     315           0 : struct BoundaryDirectionsExteriorCompute
     316             :     : BoundaryDirectionsExterior<VolumeDim>,
     317             :       db::ComputeTag {
     318           0 :   static constexpr size_t volume_dim = VolumeDim;
     319           0 :   using base = BoundaryDirectionsExterior<VolumeDim>;
     320           0 :   using return_type = std::unordered_set<::Direction<VolumeDim>>;
     321           0 :   using argument_tags = tmpl::list<Element<VolumeDim>>;
     322           0 :   static constexpr auto function(const gsl::not_null<return_type*> directions,
     323             :                                  const ::Element<VolumeDim>& element) {
     324             :     *directions = element.external_boundaries();
     325             :   }
     326             : };
     327             : /// @}
     328             : 
     329             : /// \ingroup DataBoxTagsGroup
     330             : /// \ingroup ComputationalDomainGroup
     331             : /// \brief Tag which is either a SimpleTag for quantities on an
     332             : /// interface, base tag to a compute item which acts on tags on an interface, or
     333             : /// base tag to a compute item which slices a tag from the volume to an
     334             : /// interface.
     335             : ///
     336             : /// The contained object will be a map from ::Direction to the item type of
     337             : /// `Tag`, with the set of directions being those produced by `DirectionsTag`.
     338             : ///
     339             : /// If a SimpleTag is desired on the interface, then this tag can be added to
     340             : /// the DataBox directly. If a ComputeTag which acts on Tags on the interface is
     341             : /// desired, then the tag should be added using `InterfaceCompute`. If a
     342             : /// ComputeTag which slices a TensorTag or a VariablesTag in the volume to an
     343             : /// Interface is desired, then it should be added using `Slice`. In all cases,
     344             : /// the tag can then be retrieved using `Tags::Interface<DirectionsTag, Tag>`.
     345             : ///
     346             : /// \tparam DirectionsTag the item of directions
     347             : /// \tparam Tag the tag labeling the item
     348             : ///
     349             : /// \see InterfaceCompute, Slice
     350             : template <typename DirectionsTag, typename Tag>
     351           1 : struct Interface : db::SimpleTag {
     352             :   static_assert(db::is_simple_tag_v<DirectionsTag>);
     353             :   static_assert(db::is_simple_tag_v<Tag>);
     354           0 :   static std::string name() {
     355             :     return "Interface<" + db::tag_name<DirectionsTag>() + ", " +
     356             :            db::tag_name<Tag>() + ">";
     357             :   };
     358           0 :   using tag = Tag;
     359           0 :   using type = std::unordered_map<::Direction<DirectionsTag::volume_dim>,
     360             :                                   typename Tag::type>;
     361             : };
     362             : 
     363             : /// \ingroup DataBoxTagsGroup
     364             : /// \ingroup ComputationalDomainGroup
     365             : /// ::Direction to an interface
     366             : template <size_t VolumeDim>
     367           1 : struct Direction : db::SimpleTag {
     368           0 :   using type = ::Direction<VolumeDim>;
     369             : };
     370             : 
     371             : }  // namespace Tags
     372             : }  // namespace domain
     373             : 
     374             : namespace db {
     375             : namespace detail {
     376             : template <typename DirectionsTag, typename VariablesTag>
     377             : struct InterfaceSubitemsImpl {
     378             :   using type = tmpl::transform<
     379             :       typename VariablesTag::type::tags_list,
     380             :       tmpl::bind<domain::Tags::Interface, tmpl::pin<DirectionsTag>, tmpl::_1>>;
     381             : 
     382             :   using tag = domain::Tags::Interface<DirectionsTag, VariablesTag>;
     383             : 
     384             :   template <typename Subtag>
     385             :   static void create_item(
     386             :       const gsl::not_null<typename tag::type*> parent_value,
     387             :       const gsl::not_null<typename Subtag::type*> sub_value) {
     388             :     sub_value->clear();
     389             :     for (auto& direction_vars : *parent_value) {
     390             :       const auto& direction = direction_vars.first;
     391             :       auto& parent_vars = get<typename Subtag::tag>(direction_vars.second);
     392             :       auto& sub_var = (*sub_value)[direction];
     393             :       for (auto vars_it = parent_vars.begin(), sub_var_it = sub_var.begin();
     394             :            vars_it != parent_vars.end(); ++vars_it, ++sub_var_it) {
     395             :         sub_var_it->set_data_ref(&*vars_it);
     396             :       }
     397             :     }
     398             :   }
     399             : 
     400             :   // The `return_type` can be anything for Subitems because the DataBox figures
     401             :   // out the correct return type, we just use the `return_type` type alias to
     402             :   // signal to the DataBox we want mutating behavior.
     403             :   using return_type = NoSuchType;
     404             : 
     405             :   template <typename Subtag>
     406             :   static void create_compute_item(
     407             :       const gsl::not_null<typename Subtag::type*> sub_value,
     408             :       const typename tag::type& parent_value) {
     409             :     for (const auto& direction_vars : parent_value) {
     410             :       const auto& direction = direction_vars.first;
     411             :       const auto& parent_vars =
     412             :           get<typename Subtag::tag>(direction_vars.second);
     413             :       auto& sub_var = (*sub_value)[direction];
     414             :       auto sub_var_it = sub_var.begin();
     415             :       for (auto vars_it = parent_vars.begin();
     416             :            vars_it != parent_vars.end(); ++vars_it, ++sub_var_it) {
     417             :         // clang-tidy: do not use const_cast
     418             :         // The DataBox will only give out a const reference to the
     419             :         // result of a compute item.  Here, that is a reference to a
     420             :         // const map to Tensors of DataVectors.  There is no (publicly
     421             :         // visible) indirection there, so having the map const will
     422             :         // allow only allow const access to the contained DataVectors,
     423             :         // so no modification through the pointer cast here is
     424             :         // possible.
     425             :         sub_var_it->set_data_ref(const_cast<DataVector*>(&*vars_it));  // NOLINT
     426             :       }
     427             :     }
     428             :   }
     429             : };
     430             : }  // namespace detail
     431             : 
     432             : template <typename DirectionsTag, typename VariablesTag>
     433           0 : struct Subitems<domain::Tags::Interface<DirectionsTag, VariablesTag>,
     434             :                 Requires<tt::is_a_v<Variables, typename VariablesTag::type>>>
     435             :     : detail::InterfaceSubitemsImpl<DirectionsTag, VariablesTag> {};
     436             : 
     437             : }  // namespace db