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             : 
       9             : #include "DataStructures/DataBox/Tag.hpp"
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
      12             : #include "DataStructures/Tensor/IndexType.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Domain/FunctionsOfTime/Tags.hpp"
      16             : #include "Domain/Tags.hpp"
      17             : #include "Evolution/DgSubcell/ActiveGrid.hpp"
      18             : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
      19             : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
      20             : #include "ParallelAlgorithms/Events/Tags.hpp"
      21             : #include "Utilities/ErrorHandling/Assert.hpp"
      22             : #include "Utilities/GetOutput.hpp"
      23             : #include "Utilities/Gsl.hpp"
      24             : #include "Utilities/TMPL.hpp"
      25             : 
      26             : /// \cond
      27             : namespace Tags {
      28             : struct Time;
      29             : }  // namespace Tags
      30             : /// \endcond
      31             : 
      32             : namespace evolution::dg::subcell::Tags {
      33             : /*!
      34             :  * \brief "Computes" the active coordinates by setting the `DataVector`s to
      35             :  * point into the coordinates of either the DG or subcell grid.
      36             :  */
      37             : template <size_t Dim, typename Fr>
      38           1 : struct ObserverCoordinatesCompute
      39             :     : db::ComputeTag,
      40             :       ::Events::Tags::ObserverCoordinates<Dim, Fr> {
      41           0 :   using base = ::Events::Tags::ObserverCoordinates<Dim, Fr>;
      42           0 :   using return_type = typename base::type;
      43           0 :   using argument_tags = tmpl::list<ActiveGrid, Coordinates<Dim, Fr>,
      44             :                                    ::domain::Tags::Coordinates<Dim, Fr>>;
      45           0 :   static void function(const gsl::not_null<return_type*> active_coords,
      46             :                        const subcell::ActiveGrid active_grid,
      47             :                        const tnsr::I<DataVector, Dim, Fr>& subcell_coords,
      48             :                        const tnsr::I<DataVector, Dim, Fr>& dg_coords) {
      49             :     const auto set_to_refs =
      50             :         [&active_coords](const tnsr::I<DataVector, Dim, Fr>& coords) {
      51             :           for (size_t i = 0; i < Dim; ++i) {
      52             :             active_coords->get(i).set_data_ref(
      53             :                 make_not_null(&const_cast<DataVector&>(coords.get(i))));
      54             :           }
      55             :         };
      56             :     if (active_grid == subcell::ActiveGrid::Dg) {
      57             :       set_to_refs(dg_coords);
      58             :     } else {
      59             :       ASSERT(active_grid == subcell::ActiveGrid::Subcell,
      60             :              "ActiveGrid should be subcell if it isn't DG. Maybe an extra enum "
      61             :              "entry was added?");
      62             :       set_to_refs(subcell_coords);
      63             :     }
      64             :   }
      65             : };
      66             : 
      67             : /*!
      68             :  * \brief Computes the active inverse Jacobian.
      69             :  */
      70             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
      71           1 : struct ObserverInverseJacobianCompute
      72             :     : ::Events::Tags::ObserverInverseJacobian<Dim, SourceFrame, TargetFrame>,
      73             :       db::ComputeTag {
      74           0 :   using base =
      75             :       ::Events::Tags::ObserverInverseJacobian<Dim, SourceFrame, TargetFrame>;
      76           0 :   using return_type = typename base::type;
      77           0 :   using argument_tags =
      78             :       tmpl::list<::Events::Tags::ObserverCoordinates<Dim, SourceFrame>,
      79             :                  ::Tags::Time, domain::Tags::FunctionsOfTime,
      80             :                  domain::Tags::ElementMap<Dim, Frame::Grid>,
      81             :                  domain::CoordinateMaps::Tags::CoordinateMap<Dim, Frame::Grid,
      82             :                                                              Frame::Inertial>>;
      83           0 :   static void function(
      84             :       const gsl::not_null<return_type*> observer_inverse_jacobian,
      85             :       const tnsr::I<DataVector, Dim, SourceFrame>& source_coords,
      86             :       [[maybe_unused]] const double time,
      87             :       const std::unordered_map<
      88             :           std::string,
      89             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
      90             :           functions_of_time,
      91             :       const ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
      92             :       const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>&
      93             :           grid_to_inertial_map) {
      94             :     static_assert(std::is_same_v<SourceFrame, Frame::ElementLogical> or
      95             :                   std::is_same_v<SourceFrame, Frame::Grid>);
      96             :     static_assert(std::is_same_v<TargetFrame, Frame::Inertial> or
      97             :                   std::is_same_v<TargetFrame, Frame::Grid>);
      98             :     static_assert(not std::is_same_v<SourceFrame, TargetFrame>);
      99             :     if constexpr (std::is_same_v<SourceFrame, Frame::ElementLogical>) {
     100             :       if constexpr (std::is_same_v<TargetFrame, Frame::Inertial>) {
     101             :         if (grid_to_inertial_map.is_identity()) {
     102             :           auto logical_to_grid_inv_jac =
     103             :               logical_to_grid_map.inv_jacobian(source_coords);
     104             :           for (size_t i = 0; i < observer_inverse_jacobian->size(); ++i) {
     105             :             (*observer_inverse_jacobian)[i] =
     106             :                 std::move(logical_to_grid_inv_jac[i]);
     107             :           }
     108             :         } else {
     109             :           const auto logical_to_grid_inv_jac =
     110             :               logical_to_grid_map.inv_jacobian(source_coords);
     111             :           const auto grid_to_inertial_inv_jac =
     112             :               grid_to_inertial_map.inv_jacobian(
     113             :                   logical_to_grid_map(source_coords), time, functions_of_time);
     114             :           for (size_t logical_i = 0; logical_i < Dim; ++logical_i) {
     115             :             for (size_t inertial_i = 0; inertial_i < Dim; ++inertial_i) {
     116             :               observer_inverse_jacobian->get(logical_i, inertial_i) =
     117             :                   logical_to_grid_inv_jac.get(logical_i, 0) *
     118             :                   grid_to_inertial_inv_jac.get(0, inertial_i);
     119             :               for (size_t grid_i = 1; grid_i < Dim; ++grid_i) {
     120             :                 observer_inverse_jacobian->get(logical_i, inertial_i) +=
     121             :                     logical_to_grid_inv_jac.get(logical_i, grid_i) *
     122             :                     grid_to_inertial_inv_jac.get(grid_i, inertial_i);
     123             :               }
     124             :             }
     125             :           }
     126             :         }
     127             :       } else {
     128             :         *observer_inverse_jacobian =
     129             :             logical_to_grid_map.inv_jacobian(source_coords);
     130             :       }
     131             :     } else {
     132             :       *observer_inverse_jacobian = grid_to_inertial_map.inv_jacobian(
     133             :           source_coords, time, functions_of_time);
     134             :     }
     135             :   }
     136             : };
     137             : 
     138             : /*!
     139             :  * \brief Computes the active Jacobian and determinant of the inverse Jacobian.
     140             :  */
     141             : template <size_t Dim, typename SourceFrame, typename TargetFrame>
     142           1 : struct ObserverJacobianAndDetInvJacobianCompute
     143             :     : ::Tags::Variables<tmpl::list<
     144             :           ::Events::Tags::ObserverDetInvJacobian<SourceFrame, TargetFrame>,
     145             :           ::Events::Tags::ObserverJacobian<Dim, SourceFrame, TargetFrame>>>,
     146             :       db::ComputeTag {
     147           0 :   using base = ::Tags::Variables<tmpl::list<
     148             :       ::Events::Tags::ObserverDetInvJacobian<SourceFrame, TargetFrame>,
     149             :       ::Events::Tags::ObserverJacobian<Dim, SourceFrame, TargetFrame>>>;
     150           0 :   using return_type = typename base::type;
     151           0 :   using argument_tags = tmpl::list<
     152             :       ::Events::Tags::ObserverInverseJacobian<Dim, SourceFrame, TargetFrame>>;
     153           0 :   static void function(const gsl::not_null<return_type*> result,
     154             :                        const InverseJacobian<DataVector, Dim, SourceFrame,
     155             :                                              TargetFrame>& inv_jacobian) {
     156             :     determinant_and_inverse(result, inv_jacobian);
     157             :   }
     158             : };
     159             : }  // namespace evolution::dg::subcell::Tags