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/Tensor/EagerMath/Determinant.hpp"
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Domain/ElementMap.hpp"
      13             : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
      14             : #include "Utilities/SetNumberOfGridPoints.hpp"
      15             : 
      16             : /// \cond
      17             : namespace Tags {
      18             : struct Time;
      19             : }  // namespace Tags
      20             : /// \endcond
      21             : 
      22           1 : namespace evolution::dg::subcell::fd::Tags {
      23             : /// \brief The inverse Jacobian from the element logical frame to the grid frame
      24             : /// at the cell centers.
      25             : ///
      26             : /// Specifically, \f$\partial x^{\bar{i}} / \partial x^i\f$, where \f$\bar{i}\f$
      27             : /// denotes the element logical frame and \f$i\f$ denotes the grid frame.
      28             : template <size_t Dim>
      29           1 : struct InverseJacobianLogicalToGrid : db::SimpleTag {
      30           0 :   static std::string name() { return "InverseJacobian(Logical,Grid)"; }
      31           0 :   using type =
      32             :       ::InverseJacobian<DataVector, Dim, Frame::ElementLogical, Frame::Grid>;
      33             : };
      34             : 
      35             : /// \brief The determinant of the inverse Jacobian from the element logical
      36             : /// frame to the grid frame at the cell centers.
      37           1 : struct DetInverseJacobianLogicalToGrid : db::SimpleTag {
      38           0 :   static std::string name() { return "Det(InverseJacobian(Logical,Grid))"; }
      39           0 :   using type = Scalar<DataVector>;
      40             : };
      41             : 
      42             : /// \brief The determinant of the inverse Jacobian from the element logical
      43             : /// frame to the inertial frame at the cell centers.
      44           1 : struct DetInverseJacobianLogicalToInertial : db::SimpleTag {
      45           0 :   static std::string name() { return "Det(InverseJacobian(Logical,Inertial))"; }
      46           0 :   using type = Scalar<DataVector>;
      47             : };
      48             : 
      49             : /// \brief The inverse Jacobian from the element logical frame to the inertial
      50             : /// frame at the cell centers.
      51             : ///
      52             : /// Specifically, \f$\partial x^{\bar{i}} / \partial x^i\f$, where \f$\bar{i}\f$
      53             : /// denotes the element logical frame and \f$i\f$ denotes the inertial frame.
      54             : template <size_t Dim>
      55           1 : struct InverseJacobianLogicalToInertial : db::SimpleTag {
      56           0 :   static std::string name() { return "InverseJacobian(Logical,Inertial)"; }
      57           0 :   using type = ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
      58             :                                  Frame::Inertial>;
      59             : };
      60             : 
      61             : /// \brief The inverse Hessian from the element logical frame to the grid frame
      62             : /// at the cell centers.
      63             : ///
      64             : /// Specifically, \f$\partial^2 x^{\bar{i}} / {\partial x^i \partial x^j}\f$,
      65             : /// where \f$\bar{i}\f$ denotes the element logical frame and \f$i\f$ denotes
      66             : /// the grid frame.
      67             : template <size_t Dim>
      68           1 : struct InverseHessianLogicalToGrid : db::SimpleTag {
      69           0 :   static std::string name() { return "InverseHessian(Logical,Grid)"; }
      70           0 :   using type =
      71             :       ::InverseHessian<DataVector, Dim, Frame::ElementLogical, Frame::Grid>;
      72             : };
      73             : 
      74             : /// \brief The inverse Hessian from the element logical frame to the inertial
      75             : /// frame at the cell centers.
      76             : ///
      77             : /// Specifically, \f$\partial^2 x^{\bar{i}} / {\partial x^i \partial x^j}\f$,
      78             : /// where \f$\bar{i}\f$ denotes the element logical frame and \f$i\f$ denotes
      79             : /// the inertial frame.
      80             : template <size_t Dim>
      81           1 : struct InverseHessianLogicalToInertial : db::SimpleTag {
      82           0 :   static std::string name() { return "InverseHessian(Logical,Inertial)"; }
      83           0 :   using type =
      84             :       ::InverseHessian<DataVector, Dim, Frame::ElementLogical, Frame::Inertial>;
      85             : };
      86             : 
      87             : /// Compute item for the inverse jacobian matrix from logical to
      88             : /// grid coordinates
      89             : template <typename MapTagLogicalToGrid, size_t Dim>
      90           1 : struct InverseJacobianLogicalToGridCompute : InverseJacobianLogicalToGrid<Dim>,
      91             :                                              db::ComputeTag {
      92           0 :   static constexpr size_t dim = Dim;
      93           0 :   using base = InverseJacobianLogicalToGrid<Dim>;
      94           0 :   using return_type = typename base::type;
      95           0 :   using argument_tags = tmpl::list<
      96             :       MapTagLogicalToGrid,
      97             :       evolution::dg::subcell::Tags::Coordinates<dim, Frame::ElementLogical>>;
      98           0 :   static void function(
      99             :       const gsl::not_null<return_type*> inverse_jacobian_logical_to_grid,
     100             :       const ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
     101             :       const tnsr::I<DataVector, dim, Frame::ElementLogical>& logical_coords) {
     102             :     *inverse_jacobian_logical_to_grid =
     103             :         logical_to_grid_map.inv_jacobian(logical_coords);
     104             :   }
     105             : };
     106             : 
     107             : /// Compute item for the determinant of the inverse jacobian matrix
     108             : /// from logical to grid coordinates
     109             : template <size_t Dim>
     110           1 : struct DetInverseJacobianLogicalToGridCompute
     111             :     : DetInverseJacobianLogicalToGrid,
     112             :       db::ComputeTag {
     113           0 :   static constexpr size_t dim = Dim;
     114           0 :   using base = DetInverseJacobianLogicalToGrid;
     115           0 :   using return_type = typename base::type;
     116           0 :   using argument_tags = tmpl::list<InverseJacobianLogicalToGrid<Dim>>;
     117           0 :   static void function(
     118             :       const gsl::not_null<return_type*> det_inverse_jacobian_logical_to_grid,
     119             :       const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     120             :                               Frame::Grid>& inverse_jacobian_logical_to_grid) {
     121             :     *det_inverse_jacobian_logical_to_grid =
     122             :         determinant(inverse_jacobian_logical_to_grid);
     123             :   }
     124             : };
     125             : 
     126             : /// Compute item for the inverse jacobian matrix from logical to
     127             : /// inertial coordinates
     128             : template <typename MapTagGridToInertial, size_t Dim>
     129           1 : struct InverseJacobianLogicalToInertialCompute
     130             :     : InverseJacobianLogicalToInertial<Dim>,
     131             :       db::ComputeTag {
     132           0 :   static constexpr size_t dim = Dim;
     133           0 :   using base = InverseJacobianLogicalToInertial<Dim>;
     134           0 :   using return_type = typename base::type;
     135           0 :   using argument_tags =
     136             :       tmpl::list<MapTagGridToInertial,
     137             :                  evolution::dg::subcell::Tags::Coordinates<dim, Frame::Grid>,
     138             :                  ::Tags::Time, ::domain::Tags::FunctionsOfTime,
     139             :                  InverseJacobianLogicalToGrid<Dim>>;
     140           0 :   static void function(
     141             :       const gsl::not_null<return_type*> inverse_jacobian_logical_to_inertial,
     142             :       const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>&
     143             :           grid_to_inertial_map,
     144             :       const tnsr::I<DataVector, dim, Frame::Grid>& grid_coords,
     145             :       const double time,
     146             :       const std::unordered_map<
     147             :           std::string,
     148             :           std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
     149             :           functions_of_time,
     150             :       const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     151             :                               Frame::Grid>& inverse_jacobian_logical_to_grid) {
     152             :     if (grid_to_inertial_map.is_identity()) {
     153             :       // Optimization for time-independent maps; we just point to the
     154             :       // logical-to-grid inverse jacobian.
     155             :       const size_t num_pts = inverse_jacobian_logical_to_grid[0].size();
     156             :       for (size_t storage_index = 0;
     157             :            storage_index < inverse_jacobian_logical_to_grid.size();
     158             :            ++storage_index) {
     159             :         make_const_view(
     160             :             make_not_null(&std::as_const(
     161             :                 (*inverse_jacobian_logical_to_inertial)[storage_index])),
     162             :             inverse_jacobian_logical_to_grid[storage_index], 0, num_pts);
     163             :       }
     164             :     } else {
     165             :       set_number_of_grid_points(inverse_jacobian_logical_to_inertial,
     166             :                                 get<0>(grid_coords).size());
     167             :       // Get grid to inertial inverse jacobian
     168             :       const auto& inverse_jacobian_grid_to_inertial =
     169             :           grid_to_inertial_map.inv_jacobian(grid_coords, time,
     170             :                                             functions_of_time);
     171             :       for (size_t i = 0; i < Dim; i++) {
     172             :         for (size_t j = 0; j < Dim; j++) {
     173             :           auto& inv_jacobian_component =
     174             :               inverse_jacobian_logical_to_inertial->get(i, j);
     175             :           inv_jacobian_component = 0.;
     176             :           for (size_t k = 0; k < Dim; k++) {
     177             :             inv_jacobian_component +=
     178             :                 inverse_jacobian_logical_to_grid.get(i, k) *
     179             :                 inverse_jacobian_grid_to_inertial.get(k, j);
     180             :           }
     181             :         }
     182             :       }
     183             :     }
     184             :   }
     185             : };
     186             : 
     187             : /// Compute item for the determinant of the inverse jacobian matrix
     188             : /// from logical to inertial coordinates
     189             : template <typename MapTagGridToInertial, size_t Dim>
     190           1 : struct DetInverseJacobianLogicalToInertialCompute
     191             :     : DetInverseJacobianLogicalToInertial,
     192             :       db::ComputeTag {
     193           0 :   static constexpr size_t dim = Dim;
     194           0 :   using base = DetInverseJacobianLogicalToInertial;
     195           0 :   using return_type = typename base::type;
     196           0 :   using argument_tags =
     197             :       tmpl::list<MapTagGridToInertial, InverseJacobianLogicalToInertial<Dim>,
     198             :                  DetInverseJacobianLogicalToGrid>;
     199           0 :   static void function(
     200             :       const gsl::not_null<return_type*>
     201             :           det_inverse_jacobian_logical_to_inertial,
     202             :       const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>&
     203             :           grid_to_inertial_map,
     204             :       const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     205             :                               Frame::Inertial>&
     206             :           inverse_jacobian_logical_to_inertial,
     207             :       const Scalar<DataVector>& det_inverse_jacobian_logical_to_grid) {
     208             :     if (grid_to_inertial_map.is_identity()) {
     209             :       const size_t num_pts = get(det_inverse_jacobian_logical_to_grid).size();
     210             :       make_const_view(make_not_null(&std::as_const(
     211             :                           get(*det_inverse_jacobian_logical_to_inertial))),
     212             :                       get(det_inverse_jacobian_logical_to_grid), 0, num_pts);
     213             :     } else {
     214             :       *det_inverse_jacobian_logical_to_inertial =
     215             :           determinant(inverse_jacobian_logical_to_inertial);
     216             :     }
     217             :   }
     218             : };
     219             : 
     220             : #ifdef SPECTRE_AUTODIFF
     221             : /// Compute item for the inverse hessian matrix from logical to
     222             : /// grid coordinates
     223             : template <typename MapTagLogicalToGrid, size_t Dim>
     224           1 : struct InverseHessianLogicalToGridCompute : InverseHessianLogicalToGrid<Dim>,
     225             :                                             db::ComputeTag {
     226           0 :   static constexpr size_t dim = Dim;
     227           0 :   using base = InverseHessianLogicalToGrid<Dim>;
     228           0 :   using return_type = typename base::type;
     229           0 :   using argument_tags = tmpl::list<
     230             :       MapTagLogicalToGrid,
     231             :       evolution::dg::subcell::Tags::Coordinates<dim, Frame::ElementLogical>>;
     232           0 :   static void function(
     233             :       const gsl::not_null<return_type*> inverse_hessian_logical_to_grid,
     234             :       const ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
     235             :       const tnsr::I<DataVector, dim, Frame::ElementLogical>& logical_coords) {
     236             :     *inverse_hessian_logical_to_grid =
     237             :         logical_to_grid_map.inv_hessian(logical_coords);
     238             :   }
     239             : };
     240             : 
     241             : /// Compute item for the inverse hessian matrix from logical to
     242             : /// inertial coordinates
     243             : template <typename MapTagGridToInertial, size_t Dim>
     244           1 : struct InverseHessianLogicalToInertialCompute
     245             :     : InverseHessianLogicalToInertial<Dim>,
     246             :       db::ComputeTag {
     247           0 :   static constexpr size_t dim = Dim;
     248           0 :   using base = InverseHessianLogicalToInertial<Dim>;
     249           0 :   using return_type = typename base::type;
     250           0 :   using argument_tags =
     251             :       tmpl::list<MapTagGridToInertial,
     252             :                  evolution::dg::subcell::Tags::Coordinates<dim, Frame::Grid>,
     253             :                  ::Tags::Time, ::domain::Tags::FunctionsOfTime,
     254             :                  InverseJacobianLogicalToGrid<Dim>,
     255             :                  InverseHessianLogicalToGrid<Dim>>;
     256           0 :   static void function(
     257             :       const gsl::not_null<return_type*> inverse_hessian_logical_to_inertial,
     258             :       const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>&
     259             :           grid_to_inertial_map,
     260             :       const tnsr::I<DataVector, dim, Frame::Grid>& grid_coords,
     261             :       const double time,
     262             :       const std::unordered_map<
     263             :           std::string,
     264             :           std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
     265             :           functions_of_time,
     266             :       const ::InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     267             :                               Frame::Grid>& inverse_jacobian_logical_to_grid,
     268             :       const ::InverseHessian<DataVector, Dim, Frame::ElementLogical,
     269             :                              Frame::Grid>& inverse_hessian_logical_to_grid) {
     270             :     if (grid_to_inertial_map.is_identity()) {
     271             :       // Optimization for time-independent maps; we just point to the
     272             :       // logical-to-grid inverse hessian.
     273             :       const size_t num_pts = inverse_hessian_logical_to_grid[0].size();
     274             :       for (size_t storage_index = 0;
     275             :            storage_index < inverse_hessian_logical_to_grid.size();
     276             :            ++storage_index) {
     277             :         make_const_view(
     278             :             make_not_null(&std::as_const(
     279             :                 (*inverse_hessian_logical_to_inertial)[storage_index])),
     280             :             inverse_hessian_logical_to_grid[storage_index], 0, num_pts);
     281             :       }
     282             :     } else {
     283             :       set_number_of_grid_points(inverse_hessian_logical_to_inertial,
     284             :                                 get<0>(grid_coords).size());
     285             :       // Get grid to inertial inverse jacobian
     286             :       const auto& inverse_jacobian_grid_to_inertial =
     287             :           grid_to_inertial_map.inv_jacobian(grid_coords, time,
     288             :                                             functions_of_time);
     289             :       // Get grid to inertial inverse hessian
     290             :       const auto& inverse_hessian_grid_to_inertial =
     291             :           grid_to_inertial_map.inv_hessian(grid_coords,
     292             :                                            inverse_jacobian_grid_to_inertial,
     293             :                                            time, functions_of_time);
     294             :       for (size_t i = 0; i < Dim; i++) {
     295             :         for (size_t j = 0; j < Dim; j++) {
     296             :           for (size_t k = j; k < Dim; k++) {
     297             :             auto& inv_hessian_component =
     298             :                 inverse_hessian_logical_to_inertial->get(i, j, k);
     299             :             inv_hessian_component = 0.;
     300             :             for (size_t n = 0; n < Dim; n++) {
     301             :               inv_hessian_component +=
     302             :                   inverse_hessian_grid_to_inertial.get(n, j, k) *
     303             :                   inverse_jacobian_logical_to_grid.get(i, n);
     304             :               for (size_t m = 0; m < Dim; m++) {
     305             :                 inv_hessian_component +=
     306             :                     inverse_jacobian_grid_to_inertial.get(n, k) *
     307             :                     inverse_hessian_logical_to_grid.get(i, m, n) *
     308             :                     inverse_jacobian_grid_to_inertial.get(m, j);
     309             :               }
     310             :             }
     311             :           }
     312             :         }
     313             :       }
     314             :     }
     315             :   }
     316             : };
     317             : #endif  // SPECTRE_AUTODIFF
     318             : }  // namespace evolution::dg::subcell::fd::Tags